text:
using Distributions
using Klara
using GAMCSampler

CURRENTDIR, CURRENTFILE = splitdir(@__FILE__)
ROOTDIR = splitdir(splitdir(splitdir(CURRENTDIR)[1])[1])[1]
SRCDIR = joinpath(ROOTDIR, "src")
DATADIR = joinpath(ROOTDIR, "data")
OUTDIR = joinpath(ROOTDIR, "output", "two_planets")

# SRCDIR = "../../src"
# DATADIR = "../../data"
# OUTDIR = "../../output/two_planets"

SUBOUTDIR = "SMMALA"

include(joinpath(SRCDIR, "rv_model.jl"))
include(joinpath(SRCDIR, "utils_ex.jl"))

using RvModelKeplerian

nchains = 10
nmcmc = 110000
nburnin = 10000

dataset = readdlm(joinpath(DATADIR, "two_planets.csv"), ',', header=false); # read observational data
obs_times = dataset[:, 1]
obs_rv = dataset[:, 2]
sigma_obs = dataset[:, 3]
set_times(obs_times); # set data to use for model evaluation
set_obs( obs_rv);
set_sigma_obs(sigma_obs);

param_true = make_param_true_ex2()
param_perturb_scale = make_param_perturb_scale(param_true)

p = BasicContMuvParameter(:p, logtarget=plogtarget, diffopts=DiffOptions(mode=:forward, order=2))

model = likelihood_model(p, false)

sampler = SMMALA(0.02, H -> softabs(H, 1000.))

mcrange = BasicMCRange(nsteps=nmcmc, burnin=nburnin)

mctuner = AcceptanceRateMCTuner(0.7, score=x -> logistic_rate_score(x, 3.), verbose=false)

outopts = Dict{Symbol, Any}(:monitor=>[:value], :diagnostics=>[:accept])

times = Array{Float64}(nchains)
stepsizes = Array{Float64}(nchains)
i = 1

while i <= nchains
  param_init = param_true.+0.01*param_perturb_scale.*randn(length(param_true))
  v0 = Dict(:p=>param_init)

  job = BasicMCJob(model, sampler, mcrange, v0, tuner=mctuner, outopts=outopts)

  tic()
  run(job)
  runtime = toc()

  chain = output(job)
  ratio = acceptance(chain)

  if 0.45 < ratio < 0.85
    writedlm(joinpath(OUTDIR, SUBOUTDIR, "chain"*lpad(string(i), 2, 0)*".csv"), chain.value, ',')
    writedlm(joinpath(OUTDIR, SUBOUTDIR, "diagnostics"*lpad(string(i), 2, 0)*".csv"), vec(chain.diagnosticvalues), ',')

    times[i] = runtime
    stepsizes[i] = job.sstate.tune.step

    println("Iteration ", i, " of ", nchains, " completed with acceptance ratio ", ratio)
    i += 1
  end
end

writedlm(joinpath(OUTDIR, SUBOUTDIR, "times.csv"), times, ',')
writedlm(joinpath(OUTDIR, SUBOUTDIR, "stepsizes.csv"), stepsizes, ',')


text:
@testset "bisect_interval" begin
    a, b = Arf(0), Arf(2)
    @test isequal(ArbExtras.bisect_interval(a, b), ((a, Arf(1)), (Arf(1), b)))
    a, b = Arf(0), Arf(Inf)
    @test isequal(ArbExtras.bisect_interval(a, b), ((a, Arf(Inf)), (Arf(Inf), b)))
    a, b, = Arf(0, prec = 80), Arf(2, prec = 90)
    @test precision(ArbExtras.bisect_interval(a, b)[1][2]) == 90

    a, b = Arf(1), Arf(2)
    @test ArbExtras.bisect_interval(a, b, log_midpoint = true)[1][2] ≈ exp(log(2) / 2)
    a, b = Arf(-2), Arf(-1)
    @test ArbExtras.bisect_interval(a, b, log_midpoint = true)[1][2] ≈ -exp(log(2) / 2)
    a, b = Arf(0), Arf(2)
    @test isequal(
        ArbExtras.bisect_interval(a, b, log_midpoint = true),
        ((a, Arf(1)), (Arf(1), b)),
    )
    a, b = Arf(-2), Arf(0)
    @test isequal(
        ArbExtras.bisect_interval(a, b, log_midpoint = true),
        ((a, Arf(-1)), (Arf(-1), b)),
    )
    a, b = Arf(-2), Arf(2)
    @test isequal(
        ArbExtras.bisect_interval(a, b, log_midpoint = true),
        ((a, Arf(0)), (Arf(0), b)),
    )

    a, b = Arb(0), Arb(2)
    @test isequal(ArbExtras.bisect_interval(a, b), ((a, Arb(1)), (Arb(1), b)))
    a, b = Arb(0), Arb(Inf)
    @test isequal(ArbExtras.bisect_interval(a, b), ((a, Arb(Inf)), (Arb(Inf), b)))
    a, b, = Arb(0, prec = 80), Arb(2, prec = 90)
    @test precision(ArbExtras.bisect_interval(a, b)[1][2]) == 90
end

@testset "bisect_interval_recursive" begin
    a, b = Arf(0), Arf(8)
    @test ArbExtras.bisect_interval_recursive(a, b, 0) == [(0, 8)]
    @test ArbExtras.bisect_interval_recursive(a, b, 1) == [(0, 4), (4, 8)]
    @test ArbExtras.bisect_interval_recursive(a, b, 2) == [(0, 2), (2, 4), (4, 6), (6, 8)]
    @test ArbExtras.bisect_interval_recursive(a, b, 3) ==
          [(0, 1), (1, 2), (2, 3), (3, 4), (4, 5), (5, 6), (6, 7), (7, 8)]

    a, b = Arb(0), Arb(8)
    @test ArbExtras.bisect_interval_recursive(a, b, 0) == [(0, 8)]
    @test ArbExtras.bisect_interval_recursive(a, b, 1) == [(0, 4), (4, 8)]
    @test ArbExtras.bisect_interval_recursive(a, b, 2) == [(0, 2), (2, 4), (4, 6), (6, 8)]
    @test ArbExtras.bisect_interval_recursive(a, b, 3) ==
          [(0, 1), (1, 2), (2, 3), (3, 4), (4, 5), (5, 6), (6, 7), (7, 8)]
end

@testset "check_tolerance" begin
    @test ArbExtras.check_tolerance(Arb((1, 1.1)), atol = 0.2, rtol = 0.2)
    @test ArbExtras.check_tolerance(Arb((1, 1.1)), atol = 0.2, rtol = 0.05)
    @test ArbExtras.check_tolerance(Arb((1, 1.1)), atol = 0.05, rtol = 0.2)
    @test !ArbExtras.check_tolerance(Arb((1, 1.1)), atol = 0.05, rtol = 0.05)

    @test ArbExtras.check_tolerance(Arb((3, 3.1)), atol = 0.2)
    @test !ArbExtras.check_tolerance(Arb((3, 3.1)), atol = 0.05)
    @test ArbExtras.check_tolerance(Arb(1), atol = 0)

    @test ArbExtras.check_tolerance(Arb((10, 10.1)), rtol = 0.02)
    @test !ArbExtras.check_tolerance(Arb((10, 10.1)), rtol = 0.005)
    @test ArbExtras.check_tolerance(Arb(1), rtol = 0)
    @test ArbExtras.check_tolerance(Arb(0), rtol = 0)
    @test !ArbExtras.check_tolerance(Arb(π) - Arb(π), rtol = 1)

    @test !ArbExtras.check_tolerance(Arb(NaN))
    @test !ArbExtras.check_tolerance(Arb(Inf))
    @test ArbExtras.check_tolerance(Arb((-1000, 1000)))
end

@testset "check_interval" begin
    @test ArbExtras.check_interval(Bool, Arf(1), Arf(2))
    @test ArbExtras.check_interval(Bool, Arf(-5), Arf(5))
    @test ArbExtras.check_interval(Bool, Arf(1), Arf(1))
    @test !ArbExtras.check_interval(Bool, Arf(2), Arf(1))
    @test !ArbExtras.check_interval(Bool, Arf(1), Arf(Inf))
    @test !ArbExtras.check_interval(Bool, Arf(-Inf), Arf(2))
    @test !ArbExtras.check_interval(Bool, Arf(-Inf), Arf(Inf))
    @test !ArbExtras.check_interval(Bool, Arf(Inf), Arf(2))
    @test !ArbExtras.check_interval(Bool, Arf(1), Arf(-Inf))
    @test !ArbExtras.check_interval(Bool, Arf(1), Arf(NaN))
    @test !ArbExtras.check_interval(Bool, Arf(NaN), Arf(2))

    @test isnothing(ArbExtras.check_interval(Arf(1), Arf(2)))
    @test isnothing(ArbExtras.check_interval(Arf(-5), Arf(5)))
    @test isnothing(ArbExtras.check_interval(Arf(1), Arf(1)))
    @test_throws ArgumentError ArbExtras.check_interval(Arf(2), Arf(1))
    @test_throws ArgumentError ArbExtras.check_interval(Arf(1), Arf(Inf))
    @test_throws ArgumentError ArbExtras.check_interval(Arf(-Inf), Arf(2))
    @test_throws ArgumentError ArbExtras.check_interval(Arf(-Inf), Arf(Inf))
    @test_throws ArgumentError ArbExtras.check_interval(Arf(Inf), Arf(2))
    @test_throws ArgumentError ArbExtras.check_interval(Arf(1), Arf(-Inf))
    @test_throws ArgumentError ArbExtras.check_interval(Arf(1), Arf(NaN))
    @test_throws ArgumentError ArbExtras.check_interval(Arf(NaN), Arf(2))
end

@testset "format_interval" begin
    @test ArbExtras.format_interval(Arf(1), Arf(2)) == "[1.0, 2.0]"
    @test ArbExtras.format_interval(Arf(1), Arf(Inf)) == "[1.0, Inf]"
    @test ArbExtras.format_interval(Arf(-Inf), Arf(2)) == "[-Inf, 2.0]"
    @test ArbExtras.format_interval(Arf(-Inf), Arf(Inf)) == "[-Inf, Inf]"

    @test ArbExtras.format_interval(Arf(3.1415), ubound(Arb(π))) == "[3.1415, 3.14159]"
    @test ArbExtras.format_interval(Arf(3.14159), ubound(Arb(π))) == "[3.14159 +/- 2.66e-6]"
    @test ArbExtras.format_interval(lbound(Arb(π)), Arf(Inf)) == "[3.14159, Inf]"
    @test ArbExtras.format_interval(Arf(0), ubound(Arb(π))) == "[0.0, 3.14159]"

    @test ArbExtras.format_interval(Arf(1.23456789), Arf(1.234567891), digits = 10) ==
          "[1.23456789, 1.234567891]"
end

@testset "taylor_remainder" begin
    @test let x = Arb((-1, 1))
        Arblib.overlaps(
            ArbExtras.taylor_remainder(cos(ArbSeries((x, 1), degree = 4)), x),
            Arb(1 // factorial(4)),
        )
    end

    for f in (cos, exp, atan)
        for x in (Arb((1, 2)), Arb(NaN), Arb(5), Arb((-10, 10)))
            for degree = 1:10
                p = f(ArbSeries((x, 1); degree))
                @test isequal(
                    ArbExtras.taylor_remainder(p, x),
                    p[degree] * (x - midpoint(x))^degree,
                )
            end
        end
    end
end

ext: jl


text:
function plot_logistic_coef(df::DataFrame, pars::Vector{Symbol}, mname="";
    fig="", title="")

    name = deepcopy(mname)
    if length(name) > 9
      name = name[1:9]
    end

    s = Vector{NamedTuple}(undef, length(pars))
    mindx = 1
    dftmp = DataFrame()
    for par in pars
        dftmp[!, par] = logistic.(df[:, par])
    end

    m, l, u = estimparam(dftmp)
    d = Dict{Symbol, NamedTuple}()
    for (indx, par) in enumerate(pars)
        d[Symbol(par)] = (mean=m[indx], lower=l[indx], upper=u[indx])
    end
    s =   (; d...)

    plot_model_coef(NamedTuple[s], pars; mnames=[name])
end

export
    plot_logistic_coef


text:
# Autogenerated wrapper script for ElasticFusionLib_jll for x86_64-linux-gnu-cxx11
export libefusion

using OpenBLAS32_jll
using SuiteSparse_jll
using GLEW_jll
using GLU_jll
using Zlib_jll
using JpegTurbo_jll
using Pangolin_jll
JLLWrappers.@generate_wrapper_header("ElasticFusionLib")
JLLWrappers.@declare_library_product(libefusion, "libefusion.so")
function __init__()
    JLLWrappers.@generate_init_header(OpenBLAS32_jll, SuiteSparse_jll, GLEW_jll, GLU_jll, Zlib_jll, JpegTurbo_jll, Pangolin_jll)
    JLLWrappers.@init_library_product(
        libefusion,
        "lib/libefusion.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# implementation of the GPUCompiler interfaces for generating eBPF code

## target

export BPFCompilerTarget

Base.@kwdef struct BPFCompilerTarget <: AbstractCompilerTarget
    prog_section::String="prog" # section for kernel to be placed in
    license::String="" # license for kernel and source code
    function_pointers::UnitRange{Int}=1:1000 # set of valid function "pointers"
end

llvm_triple(::BPFCompilerTarget) = "bpf-bpf-bpf"
llvm_datalayout(::BPFCompilerTarget) = "e-m:e-p:64:64-i64:64-n32:64-S128"

function llvm_machine(target::BPFCompilerTarget)
    triple = llvm_triple(target)
    t = Target(;triple=triple)

    cpu = ""
    feat = ""
    tm = TargetMachine(t, triple, cpu, feat)
    asm_verbosity!(tm, true)

    return tm
end

## job

runtime_slug(job::CompilerJob{BPFCompilerTarget}) = "bpf"

const bpf_intrinsics = () # TODO
isintrinsic(::CompilerJob{BPFCompilerTarget}, fn::String) = in(fn, bpf_intrinsics)

function finish_module!(job::CompilerJob{BPFCompilerTarget}, mod::LLVM.Module)
    for func in LLVM.functions(mod)
        if LLVM.name(func) == "gpu_signal_exception"
            throw(KernelError(job, "eBPF does not support exceptions"))
        end
        # Set entry section for loaders like libbpf
        LLVM.section!(func, job.target.prog_section)
    end

    # Set license
    license = job.target.license
    if license != ""
        ctx = LLVM.context(mod)
        i8 = LLVM.Int8Type(ctx)
        glob = GlobalVariable(mod, LLVM.ArrayType(i8, length(license)+1), "_license")
        linkage!(glob, LLVM.API.LLVMExternalLinkage)
        constant!(glob, true)
        section!(glob, "license")
        initializer!(glob, ConstantArray(i8, vcat(Vector{UInt8}(license),0x0)))
    end

    # Set all map definitions as external linkage
    for gv in filter(x->(section(x)=="maps")||(section(x)==".maps"), collect(LLVM.globals(mod)))
        linkage!(gv, LLVM.API.LLVMExternalLinkage)
    end
end

valid_function_pointer(job::CompilerJob{BPFCompilerTarget}, ptr::Ptr{Cvoid}) =
    reinterpret(UInt, ptr) in job.target.function_pointers


text:
# Isomap
# ------
# A Global Geometric Framework for Nonlinear Dimensionality Reduction,
# J. B. Tenenbaum, V. de Silva and J. C. Langford, Science 290 (5500): 2319-2323, 22 December 2000

#### Isomap type
struct Isomap{T <: Real} <: AbstractDimensionalityReduction
    k::Int
    model::KernelPCA
    component::AbstractVector{Int}

    Isomap{T}(k::Int, model::KernelPCA) where T = new(k, model)
    Isomap{T}(k::Int, model::KernelPCA, cc::AbstractVector{Int}) where T = new(k, model, cc)
end

## properties
outdim(R::Isomap) = outdim(R.model)
eigvals(R::Isomap) = principalvars(R.model)
neighbors(R::Isomap) = R.k
vertices(R::Isomap) = R.component

## show
summary(io::IO, R::Isomap) = print(io, "Isomap(outdim = $(outdim(R)), neighbors = $(neighbors(R)))")

## interface functions
function fit(::Type{Isomap}, X::AbstractMatrix{T}; k::Int=12, maxoutdim::Int=2) where {T<:Real}
    # Construct NN graph
    D, E = find_nn(X, k)
    G, C = largest_component(SimpleWeightedGraph(adjmat(D,E)))

    # Compute shortest path for every point
    n = length(C)
    DD = zeros(T, n, n)
    for i in 1:n
        dj = dijkstra_shortest_paths(G, i)
        DD[i,:] = dj.dists
    end

    broadcast!(x->-x*x/2, DD, DD)
    broadcast!((x,y)->(x+y)/2, DD, DD, DD') # remove roundoff error
    M = fit(KernelPCA, DD, kernel=nothing, maxoutdim=maxoutdim)

    return Isomap{T}(k, M, C)
end

transform(R::Isomap) = transform(R.model)


text:
my_I(T, N) = Matrix{T}(I, N, N)

Zygote._pullback(::AContext, ::typeof(my_I), args...) = my_I(args...), nograd_pullback

function lgssm_components(k::Separable, x::SpaceTimeGrid, storage)

    # Compute spatial covariance, and temporal GaussMarkovModel.
    r, t = x.xl, x.xr
    kr, kt = k.l, k.r
    Kr = kernelmatrix(kr, r)
    As_t, as_t, Qs_t, emission_proj_t, x0_t = lgssm_components(kt, t, storage)

    # Compute components of complete LGSSM.
    Nr = length(r)
    ident = my_I(eltype(storage), Nr)
    As = map(A -> kron(ident, A), As_t)
    as = map(a -> repeat(a, Nr), as_t)
    Qs = map(Q -> kron(Kr + ident_eps(1e-12), Q), Qs_t)
    emission_proj = _build_st_proj(emission_proj_t, Nr, ident)
    x0 = Gaussian(repeat(x0_t.m, Nr), kron(Kr, x0_t.P))
    return As, as, Qs, emission_proj, x0
end

function _build_st_proj((Hs, hs)::Tuple{AbstractVector, AbstractVector}, Nr::Integer, ident)
    return (map(H -> kron(ident, H), Hs), map(h -> fill(h, Nr), hs))
end

function build_Σs(x::RectilinearGrid, Σ::Diagonal{<:Real})
    return Diagonal.(collect.(eachcol(reshape(Σ.diag, :, length(x.xr)))))
end

function Zygote._pullback(
    ::AContext, ::typeof(build_Σs), x::RectilinearGrid, Σ::Diagonal{<:Real},
)
    build_Σs_pullback(Δ) = nothing, nothing, (diag=vcat(getfield.(Δ, :diag)...), )
    return build_Σs(x, Σ), build_Σs_pullback
end

function build_Σs(x::RegularInTime, Σ::Diagonal{<:Real})
    vs = restructure(Σ.diag, length.(x.vs))
    return Diagonal.(collect.(vs))
end

function restructure(y::AbstractVector{<:Real}, lengths::AbstractVector{<:Integer})
    idxs_start = cumsum(vcat(0, lengths)) .+ 1
    idxs_end = idxs_start[1:end-1] .+ lengths .- 1
    return map(n -> y[idxs_start[n]:idxs_end[n]], eachindex(lengths))
end

function restructure(y::AbstractVector{T}, lengths::AbstractVector{<:Integer}) where {T}
    idxs_start = cumsum(vcat(0, lengths)) .+ 1
    idxs_end = idxs_start[1:end-1] .+ lengths .- 1
    return map(eachindex(lengths)) do n
        y_missing = Vector{T}(undef, lengths[n])
        y_missing .= y[idxs_start[n]:idxs_end[n]]
        return y_missing
    end
end

function Zygote._pullback(
    ::AContext, ::typeof(restructure), y::Vector, lengths::AbstractVector{<:Integer},
)
    restructure_pullback(Δ::Vector) = nothing, vcat(Δ...), nothing
    return restructure(y, lengths), restructure_pullback
end

# Implement specific to Fills for AD's sake.
function restructure(y::Fill{<:Real}, lengths::AbstractVector{<:Integer})
    return map(l -> Fill(y.value, l), Zygote.dropgrad(lengths))
end

function restructure(y::AbstractVector, emissions::StructArray)
    return restructure(y, Zygote.dropgrad(map(dim_out, emissions)))
end

destructure(y::AbstractVector{<:AbstractVector{<:Real}}) = vcat(y...)


text:
## Dense matrix interface

dense_matrix_rows(s::CDenseMatrix) = ccall((:dense_matrix_rows, libsymengine), Int, (Ptr{Cvoid}, ), s.ptr)
dense_matrix_cols(s::CDenseMatrix) = ccall((:dense_matrix_cols, libsymengine), Int, (Ptr{Cvoid}, ), s.ptr)
function dense_matrix_rows_cols(mat::CDenseMatrix, r::UInt, c::UInt)
    ccall((:dense_matrix_rows_cols, libsymengine), Int, (Ptr{Cvoid}, UInt, UInt), s.ptr, r, c)
end

function dense_matrix_get_basic(s::CDenseMatrix, r::Int, c::Int)
    result = Basic()
    ccall((:dense_matrix_get_basic, libsymengine), Nothing, (Ref{Basic}, Ptr{Cvoid}, UInt, UInt), result, s.ptr, UInt(r), UInt(c))
    result
end

function dense_matrix_set_basic(s::CDenseMatrix, val, r::Int, c::Int)
    value = Basic(val)
    ccall((:dense_matrix_set_basic, libsymengine), Nothing, (Ptr{Cvoid}, UInt, UInt, Ref{Basic}), s.ptr, UInt(r), UInt(c), value)
    value
end

## Basic operations det, inv, transpose
function dense_matrix_det(s::CDenseMatrix)
    result = Basic()
    ccall((:dense_matrix_det, libsymengine), Nothing, (Ref{Basic}, Ptr{Cvoid}), result, s.ptr)
    result
end
function dense_matrix_inv(s::CDenseMatrix)
    result = CDenseMatrix()
    ccall((:dense_matrix_inv, libsymengine), Nothing, (Ptr{Cvoid}, Ptr{Cvoid}), result.ptr, s.ptr)
    result
end

function dense_matrix_transpose(s::CDenseMatrix)
    result = CDenseMatrix()
    ccall((:dense_matrix_transpose, libsymengine), Nothing, (Ptr{Cvoid}, Ptr{Cvoid}), result.ptr, s.ptr)
    result
end

function dense_matrix_submatrix(mat::CDenseMatrix, r1::Int, c1::Int, r2::Int, c2::Int, r::Int, c::Int)
    s = CDenseMatrix()
    ccall((:dense_matrix_inv, libsymengine), Nothing, (Ptr{Cvoid}, Ptr{Cvoid}, UInt, UInt, UInt, UInt, UInt, UInt),
          s.ptr, mat.ptr, UInt(r1), UInt(c1), UInt(r2), UInt(c2), UInt(r), UInt(c))
    s
end

## some matrix arithmetic methods
## These are unncecessary, as +,-,*,^ are inhertied by the AbstractArray Interface. Those return Array{Basic,2}. These
## return CDenseMatrix objects.
function dense_matrix_add_matrix(a::CDenseMatrix, b::CDenseMatrix)
    result = CDenseMatrix()
    ccall((:dense_matrix_add_matrix, libsymengine), Nothing, (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}), result.ptr, a.ptr, b.ptr)
    result
end

function dense_matrix_mul_matrix(a::CDenseMatrix, b::CDenseMatrix)
    result = CDenseMatrix()
    ccall((:dense_matrix_mul_matrix, libsymengine), Nothing, (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}), result.ptr, a.ptr, b.ptr)
    result
end

function dense_matrix_add_scalar(a::CDenseMatrix, b::Basic)
    result = CDenseMatrix()
    ccall((:dense_matrix_add_scalar, libsymengine), Nothing, (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}), result.ptr, a.ptr, b.ptr)
    result
end

function dense_matrix_mul_scalar(a::CDenseMatrix, b::Basic)
    result = CDenseMatrix()
    ccall((:dense_matrix_mul_scalar, libsymengine), Nothing, (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}), result.ptr, a.ptr, b.ptr)
    result
end

## Factorizations ##################################################
function dense_matrix_LU(mat::CDenseMatrix)
    ## need square?
    L = CDenseMatrix()
    U = CDenseMatrix()
    ccall((:dense_matrix_LU, libsymengine), Nothing, (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}), L.ptr, U.ptr,  mat.ptr)
    (L, U)
end

## LDL decomposition
function dense_matrix_LDL(mat::CDenseMatrix)
    L = CDenseMatrix()
    D = CDenseMatrix()
    ccall((:dense_matrix_LDL, libsymengine), Nothing, (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}), L.ptr, D.ptr,  mat.ptr)
    (L, D)
end

##  Fraction free LU factorization
function dense_matrix_FFLU(mat::CDenseMatrix)
    LU = CDenseMatrix()
    ccall((:dense_matrix_FFLU, libsymengine), Nothing, (Ptr{Cvoid}, Ptr{Cvoid}), LU.ptr, mat.ptr)
    LU
end

## Fraction free LDU factorization
function dense_matrix_FFLDU(mat::CDenseMatrix)
    L = CDenseMatrix()
    D = CDenseMatrix()
    U = CDenseMatrix()

    ccall((:dense_matrix_FFLDU, libsymengine), Nothing, (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}), L.ptr, D.ptr, U.ptr, mat.ptr)
    (L, D, U)
end

function dense_matrix_LU_solve(A::CDenseMatrix, b::CDenseMatrix)
    x = CDenseMatrix()
    ccall((:dense_matrix_LU_solve, libsymengine), Nothing, (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}), x.ptr, A.ptr, b.ptr)
    x
end

function dense_matrix_zeros(::Type{CDenseMatrix}, r::Int, c::Int)
    result = CDenseMatrix()
    ccall((:dense_matrix_zeros, libsymengine), Nothing, (Ptr{Cvoid}, UInt, UInt), result.ptr, UInt(r), UInt(c))
    result
end

function dense_matrix_ones(r::Int, c::Int)
    result = CDenseMatrix()
    ccall((:dense_matrix_ones, libsymengine), Nothing, (Ptr{Cvoid}, UInt, UInt), result.ptr, UInt(r), UInt(c))
    result
end

## dense_matrix_diag XXX don't have CVecBasic constructor

function dense_matrix_eye(N::Int, M::Int, k::Int=0)
    s = CDenseMatrix()
    ccall((:dense_matrix_eye, libsymengine), Nothing, (Ptr{Cvoid}, UInt, UInt, Int), s.ptr, UInt(N), UInt(M), k)
    s
end

dense_matrix_eq(a::CDenseMatrix, b::CDenseMatrix) = ccall((:dense_matrix_eq, libsymengine), Int, (Ptr{Cvoid},Ptr{Cvoid}), a.ptr, b.ptr)

## Plug into Julia's interface ##################################################

import Base: ==
==(a::CDenseMatrix, b::CDenseMatrix)  = dense_matrix_eq(a, b) == 1

## Abstract Array Interface
Base.size(s::CDenseMatrix) = (dense_matrix_rows(s), dense_matrix_cols(s))
Base.getindex(s::CDenseMatrix, r::Int, c::Int) = dense_matrix_get_basic(s, r-1, c-1)
Base.setindex!(s::CDenseMatrix, val, r::Int, c::Int) = dense_matrix_set_basic(s, val, r-1, c-1)

## special matrices
Base.zeros(::Type{CDenseMatrix}, r::Int, c::Int) = dense_matrix_zeros(s, r-1, c-1)
Base.ones(::Type{CDenseMatrix}, r::Int, c::Int) = dense_matrix_ones(r-1, c-1)

## basic functions
LinearAlgebra.det(s::CDenseMatrix) = dense_matrix_det(s)
Base.inv(s::CDenseMatrix) = dense_matrix_inv(s)
Base.transpose(s::CDenseMatrix) = dense_matrix_transpose(s)

LinearAlgebra.factorize(M::CDenseMatrix) = factorize(convert(Matrix, M))

"""
LU decomposition for CDenseMatrix, dense matrices of symbolic values

Also: lufact(a, val{:false}) for non-pivoting lu factorization
"""
function LinearAlgebra.lu(a::CDenseMatrix)
    l, u = dense_matrix_LU(a)
    convert(Matrix, l), convert(Matrix, u), Matrix{Basic}(LinearAlgebra.I, size(l)[1], size(l)[1])
end

# solve using LU_solve
import Base: \
\(A::CDenseMatrix, b::CDenseMatrix) = dense_matrix_LU_solve(A, b)
\(A::CDenseMatrix, b::Vector) = A \ convert(CDenseMatrix, [convert(Basic,u) for u in b])


text:
#=
We can determine how "out of order" an array A is by counting the number of inversions it has. Two elements A[i] and A[j] form an inversion if A[i] > A[j] but i < j. That is, a smaller element appears after a larger element.

Given an array, count the number of inversions it has. Do this faster than O(N^2) time.

You may assume each element in the array is distinct.

For example, a sorted list has zero inversions. The array [2, 4, 1, 3, 5] has three inversions: (2, 1), (4, 1), and (4, 3). The array [5, 4, 3, 2, 1] has ten inversions: every distinct pair forms an inversion.
=#

function inversions(A::Array{Int, 1})
    B = copy(A)
    number_of_inversions = 0
    while length(B) != 1
        if findmin(B)[2] != 1
            number_of_inversions = number_of_inversions + sum(B[1] .> B[2:end])
        end
        B = B[2:end]
    end
    return number_of_inversions
end

max_stars_repo_name: DominiqueCaron/daily-coding-problem

max_issues_repo_name: DominiqueCaron/daily-coding-problem

max_forks_repo_name: DominiqueCaron/daily-coding-problem


text:
module GeoRegions

## Modules Used
using Dates
using DelimitedFiles
using GeometryBasics
using Logging
using PolygonOps
using PrettyTables

import Base: show, read

## Exporting the following functions:
export
        GeoRegion, RectRegion, PolyRegion, RegionGrid, RectGrid, PolyGrid,
        resetGeoRegions, templateGeoRegions, listGeoRegions, readGeoRegions,
        isGeoRegion, addGeoRegions, removeGeoRegion, coordGeoRegion, isinGeoRegion,
        tableGeoRegions, tableRectRegions, tablePolyRegions,
        Point2

## Abstract types
"""
    GeoRegion

Abstract supertype for geographical regions, with the following subtypes:

    RectRegion{ST<:AbstractString, FT<:Real} <: GeoRegion
    PolyRegion{ST<:AbstractString, FT<:Real} <: GeoRegion

Both `RectRegion` and `PolyRegion` types contain the following fields:
* `regID` - A `String` Type, the identifier for the GeoRegion
* `parID` - A `String` Type, the identifier for the parent GeoRegion
* `name` - A `String` Type, the full name of the GeoRegion
* `N` - A `Float` Type, the north boundary of the GeoRegion
* `S` - A `Float` Type, the south boundary of the GeoRegion
* `E` - A `Float` Type, the east boundary of the GeoRegion
* `W` - A `Float` Type, the est boundary of the GeoRegion
* `is180` - A `Bool` Type, is `W` < 0
* `is360` - A `Bool` Type, is `E` > 180

A `PolyRegion` type will also contain the following field:
* `shape` - A vector of `Point2` Types, defining a non-rectilinear shape of the GeoRegion
"""
abstract type GeoRegion end

struct RectRegion{ST<:AbstractString, FT<:Real} <: GeoRegion
    regID :: ST
    parID :: ST
    name  :: ST
    N     :: FT
    S     :: FT
    E     :: FT
    W     :: FT
    is180 :: Bool
    is360 :: Bool
end

struct PolyRegion{ST<:AbstractString, FT<:Real} <: GeoRegion
    regID :: ST
    parID :: ST
    name  :: ST
    N     :: FT
    S     :: FT
    E     :: FT
    W     :: FT
    shape :: Vector{Point2{FT}}
    is180 :: Bool
    is360 :: Bool
end

"""
    RegionGrid

Abstract supertype for geographical regions, with the following subtypes:

    RectGrid{FT<:Real} <: RegionGrid
    PolyGrid{FT<:Real} <: RegionGrid

Both `RectGrid` and `PolyGrid` types contain the following fields:
* `igrid` - A vector of `Int`s defining the gridpoint indices of the [N,S,E,W] points respectively
* `ilon` - A vector of `Int`s defining indices of the parent longitude vector describing the region
* `ilat` - A vector of `Int`s defining indices of the parent latitude vector describing the region
* `glon` - A vector of `Float`s defining the latitude vector describing the region
* `glat` - A vector of `Float`s defining the latitude vector describing the region

A `PolyGrid` type will also contain the following field:
* `mask` - An array of 0s and 1s defining a non-rectlinear shape within a rectilinear grid where data is valid (only available in PolyGrid types)
"""
abstract type RegionGrid end

struct RectGrid{FT<:Real} <: RegionGrid
    igrid :: Vector{Int}
    ilon  :: Vector{Int}
    ilat  :: Vector{Int}
    glon  :: Vector{FT}
    glat  :: Vector{FT}
end

struct PolyGrid{FT<:Real} <: RegionGrid
    igrid :: Vector{Int}
    ilon  :: Vector{Int}
    ilat  :: Vector{Int}
    glon  :: Vector{FT}
    glat  :: Vector{FT}
    mask  :: Array{FT,2}
end

modulelog() = "$(now()) - GeoRegions.jl"

function __init__()
    jfol = joinpath(DEPOT_PATH[1],"files","GeoRegions"); mkpath(jfol);
    flist   = ["rectlist.txt","polylist.txt","giorgi.txt","srex.txt","ar6.txt"]

    for fname in flist
        if !isfile(joinpath(jfol,fname))
            copygeoregions(fname)
            @info "$(modulelog()) - $(fname) does not exist in $(jfol), copying ..."
        end
    end
end

## Including other files in the module
include("Read.jl")
include("Create.jl")
include("Query.jl")
include("IsIn.jl")
include("IsInGeoRegion.jl")
include("Extract.jl")
include("Show.jl")
include("Tables.jl")

end # module


text:
# Relativistic CSFs

#
# type CSF
# ------------------------------------------------------------------------------------------
"""
    $(TYPEDEF)

Represents a configuration state function.

[`CSF`](@ref) is iterable, returning a tuple with the elements
`(orbital::Orbital, n_electrons, orbital_coupling, csf_coupling)` on each
iteration.
"""
struct CSF
    angularsym :: AngularSymmetry
    orbitals :: Vector{RelativisticOrbital{Int}}
    occupations :: Vector{Int}
    orbcouplings :: Vector{AngularMomentum}
    csfcouplings :: Vector{AngularMomentum}

    function CSF(
            orbitals::Vector{<:RelativisticOrbital},
            electrons::Vector{Int},
            couplings_orbitals::Vector{AngularMomentum},
            couplings_csf::Vector{AngularMomentum},
            parity,
        )
        let norbitals = length(orbitals)
            @assert length(electrons) == norbitals
            @assert length(couplings_orbitals) == norbitals
            @assert length(couplings_csf) == norbitals
        end
        @assert issorted(orbitals)
        @assert all(x -> x!=0, electrons)
        new(
            AngularSymmetry(last(couplings_csf), parity),
            orbitals, electrons, couplings_orbitals, couplings_csf,
        )
    end
end

import Base: ==
function ==(csf1::CSF, csf2::CSF)
    (csf1.angularsym == csf2.angularsym)     &&
    (csf1.orbitals == csf2.orbitals)         &&
    (csf1.occupations == csf2.occupations)   &&
    (csf1.orbcouplings == csf2.orbcouplings) &&
    (csf1.csfcouplings == csf2.csfcouplings)
end

function Base.iterate(csf::CSF, s=1)
    s > length(csf) && return nothing
    item = (csf.orbitals[s], csf.occupations[s], csf.orbcouplings[s], csf.csfcouplings[s])
    return item, s + 1
end

Base.length(csf::CSF) = length(csf.orbitals)

parity(csf::CSF) = parity(csf.angularsym)
angularmomentum(csf::CSF) = angularmomentum(csf.angularsym)
nelectrons(csf :: CSF) = sum(csf.occupations)

"""
    nexcitations(csf_from::CSF, csf_to::CSF)

Calculates the excitations between relativistic orbitals. I.e. difference of occupation
number between `nl-` to `nl` counts as an excitation.
"""
function nexcitations(csf_from::CSF, csf_to::CSF)
    @assert nelectrons(csf_from) == nelectrons(csf_to)
    let from_orbs = [(i, orb.n, kappa(orb)) for (i, orb) in enumerate(csf_from.orbitals)],
        to_orbs = [(i, orb.n, kappa(orb)) for (i, orb) in enumerate(csf_to.orbitals)]
        @assert length(from_orbs) == length(unique(from_orbs))
        @assert length(to_orbs) == length(unique(to_orbs))
    end

    to_unchecked = fill(true, length(csf_to.orbitals))

    nexcitations = 0
    nremoved = 0 # for sanity checking only

    # Let's count all the electrons that have been excited into some of
    # the occupied orbitals of csf_from.
    for (i, orb) in enumerate(csf_from.orbitals)
        toids = findall(o -> o.n == orb.n && kappa(o) == kappa(orb), csf_to.orbitals)
        to_nelectrons = if length(toids) == 0
            0
        elseif length(toids) == 1
            # Mark this csf_to orbital checked.
            toid = first(toids)
            @assert to_unchecked[toid]
            to_unchecked[toid] = false
            csf_to.occupations[toid]
        else
            error("Too many matches ($(length(toids))) among csf_to.orbitals.")
        end

        diff = to_nelectrons - csf_from.occupations[i]
        # If diff < 0, electrons have been excited from this orbital
        (diff < 0) && (nremoved += abs(diff)) # for sanity checking
        nexcitations += (diff > 0) ? diff : 0
    end

    # Now let's loop over all the remaning csf_to orbitals
    for (i, orb) in enumerate(csf_to.orbitals)
        to_unchecked[i] || continue
        nexcitations += csf_to.occupations[i]
    end

    @assert nexcitations == nremoved # sanity check
    return nexcitations
end

function Base.show(io::IO, csf::CSF)
    write(io, '"')
    print(io, csf)
    write(io, '"')
end

function Base.print(io::IO, csf::CSF)
    first = true
    for (orb, nelec, orbcoupling, csfcoupling) in csf
        !first && write(io, ' ')
        print(io, orb)
        write(io, '(', string(nelec), '|', string(orbcoupling), '|', string(csfcoupling), ')')
        first = false
    end
    write(io, " | ", string(csf.angularsym))
    nothing
end

function kappa2rso(kappa :: Integer)
    lstr = specname(κ2ℓ(kappa))
    (kappa < 0) ? lstr : lstr*"-"
end


text:
# IJulia stopped working after switching to external Conda
# Start

using Conda
using PyCall
import Pkg

env = "py38"

path = Conda.bin_dir(dirname(PyCall.pyprogramname))
if isnothing(findfirst(path, ENV["PATH"])) && Sys.iswindows()
    ENV["PATH"] = path * ";" * ENV["PATH"]
    ENV["PATH"] = Conda.ROOTENV * ";" * ENV["PATH"]
end

#using IJulia
#jupyterlab(dir=pwd())

# Use external Conda
if false # does not work
    ENV["JUPYTER"] = joinpath(Conda.ROOTENV, "envs", env, "Scripts", "jupyter.exe")
    #ENV["JUPYTER"] = joinpath(Conda.ROOTENV, "envs", env, "Scripts")
    Pkg.build("IJulia")
end

# Use internal Conda
if false # does not work
    ENV["JUPYTER"] = ""
    Pkg.build("IJulia")
end


text:
"""
JohnsonSu(ξ,λ,γ,δ)

The *Johnson distribution* with location `ξ` and scale `λ` has probability density function

```math
f(x; \\xi,\\lambda,\\gamma,\\delta) = \\frac{\\delta}{\\sqrt(2*\\pi)} + \\frac{1}{sqrt(x^2+1)} * exp(-\\frac{1}{2}\\lambda+\\delta*log(x+sqrt(x^2 +1))^-2))
```

```
julia
```
"""
struct JohnsonSu{T<:Real} <: ContinuousUnivariateDistribution
    ξ::T
    λ::T
    γ::T
    δ::T
    JohnsonSu{T}(ξ::T, λ::T, γ::T, δ::T) where {T<:Real} = new{T}(ξ,λ,γ,δ)
end

function JohnsonSu(ξ::T, λ::T, γ::T, δ::T) where {T <: Real}
    @check_args(JohnsonSu, -Inf < ξ < Inf && λ > 0 && -Inf < γ < Inf && δ > 0)
    return JohnsonSu{T}(ξ,λ,γ,δ)
end
#using Distributions
#using Distributions: Normal
#### Outer constructors
#JohnsonSu{T<:Real}(ξ::T, λ::T, γ::T, δ::T) = JohnsonSu{T}(ξ, λ, γ, δ)
JohnsonSu(ξ::Real, λ::Real, γ::Real, δ::Real) = JohnsonSu(promote(ξ, λ, γ, δ)...)
JohnsonSu(ξ::Integer, λ::Integer, γ::Integer, δ::Integer) = JohnsonSu(Float64(ξ), Float64(λ), Float64(γ), Float64(δ))
JohnsonSu() = JohnsonSu(1.0, 1.0, 1.0, 1.0)

@distr_support JohnsonSu -Inf Inf

#### Conversions
function convert(::Type{JohnsonSu{T}}, ξ::Real, λ::Real, γ::Real, δ::Real) where T<:Real
 JohnsonSu(T(ξ), T(λ), T(γ), T(δ))
end

#convert{T <: Real, S <: Float64}(::Type{JohnsonSu{T}}, d::JohnsonSu{S}) = JohnsonSu( T(d.ξ), T(d.λ),  T(d.γ), T(d.δ))

#### Parameters
peakness(d::JohnsonSu) = d.ξ
shape(d::JohnsonSu) = d.λ
location(d::JohnsonSu) = d.γ
scale(d::JohnsonSu) = d.δ

params(d::JohnsonSu) = (d.ξ, d.λ, d.γ, d.δ)
@inline partype(d::JohnsonSu{T}) where {T<:Real} = T

#### Statistics
mean(d::JohnsonSu) = d.ξ-d.λ*exp.(d.δ^2/2)*sinh(d.γ/d.δ)
median(d::JohnsonSu) = d.ξ+d.λ*sinh(-(d.γ/d.δ))
var(d::JohnsonSu) = d.λ^2/2*(exp.(d.δ^-2)-1)*(exp.(d.δ^-2)*cosh(2d.γ/d.δ)+1)
mode(d::JohnsonSu) = d.ξ-d.λ*exp.(d.δ^2/2)*sinh(d.γ/d.δ) # copied from mean function.
skewness(d::JohnsonSu) = -(exp.(d.δ^-2)^(1/2)*(exp.(d.δ^-2)-1)^2*(exp.(d.δ^-2)*(exp.(d.δ^-2)+2)*sinh(3*(d.λ/d.δ))+3*sinh((d.λ/d.δ)))/sqrt(2)*((exp.(d.δ^-2)-1)*(exp.(d.δ^-2)*cosh(2*(d.λ/d.δ))+1))^3/2)
kurtosis(d::JohnsonSu) = (exp.(d.δ^-2)^2*((exp.(d.δ^-2)^4+2*exp.(d.δ^-2)^3+3*exp.(d.δ^-2)^2-3)*cosh(4*(d.λ/d.δ)+(4*(exp.(d.δ^-2)+2))*cosh(2*(d.λ/d.δ))))+3*(2*exp.(d.δ^-2)+1)) / (2(exp.(d.δ^-2)*cosh(2*(d.λ/d.δ))+1)^2)

#### Evaluation
@_delegate_statsfuns JohnsonSu ξ, λ, γ, δ

mgf(d::JohnsonSu, t::Real) = exp(t * d.ξ + d.λ^2/2 * t^2)
cf(d::JohnsonSu, t::Real) = exp(im * t * d.ξ - d.λ^2/2 * t^2)

#### Sampling
rand(d::JohnsonSu) = rand(GLOBAL_RNG, d)
rand(rng::AbstractRNG, d::JohnsonSu) = d.ξ + d.λ * randn(rng)

#### Fitting
struct JohnsonSuStats <: SufficientStats
    s::Float64    # (weighted) sum of x
    m::Float64    # (weighted) mean of x
    s2::Float64   # (weighted) sum of (x - μ)^2
    s3::Float64    # total sample weight
end

function suffstats(::Type{JohnsonSu}, x::AbstractArray{T}) where T<:Real
    n = length(x)

    # compute m
    s = x[1]
    for i = 2:n
        @inbounds s += x[i]
    end

    m = s / n

    # compute s2
    s2 = abs2(x[1] - m)
    for i = 2:n
        @inbounds s2 += abs2(x[i] - m)
    end

    # compute s3
    s3 = kurtosis(x)

    JohnsonSuStats(m, s, s2, s3)
end

# fit_mle based on sufficient statistics
fit_mle(::Type{JohnsonSu}, ss::JohnsonSuStats) = JohnsonSu(ss.m, ss.s, ss.s2, ss.s3)

@quantile_newton(JohnsonSu)

pdf(d::JohnsonSu, x::Real) = (d.δ/d.λ*sqrt.(2π)*1/(sqrt.(1+(x-d.ξ/d.λ).^2 +1)))exp.(-1/2*(d.γ+d.δ*asinh.(x-d.ξ/d.λ).^2))
#cdf(d::JohnsonSu, x::Real) = cdf(Normal(0,1),x)*(d.γ + d.δ*sinh.((x-d.ξ)/d.λ))

function fit_mle(::Type{JohnsonSu}, x::AbstractArray{T}; xi::Float64=NaN, lambda::Float64=NaN, gamma::Float64=NaN, delta::Float64=NaN) where T<:Real
    if isnan(xi)
        if isnan(lambda)
            fit_mle(JohnsonSu, suffstats(JohnsonSu, x))
        else
            g = NormalKnownSigma(lambda)
            fit_mle(g, suffstats(g, x))
        end
    else
        if isnan(lambda)
            g = NormalKnownMu(gamma)
            fit_mle(g, suffstats(g, x))
        else
            JohnsonSu(xi, lambda, gamma, delta)
        end
    end
end

function fit_mle(::Type{JohnsonSu}, x::AbstractArray{T}, w::AbstractArray{Float64}; xi::Float64=NaN, lambda::Float64=NaN, gamma::Float64=NaN, delta::Float64=NaN) where T<:Real
    if isnan(mu)
        if isnan(sigma)
            fit_mle(JohnsonSu, suffstats(JohnsonSu, x, w))
        else
            g = NormalKnownSigma(sigma)
            fit_mle(g, suffstats(g, x, w))
        end
    else
        if isnan(sigma)
            g = NormalKnownMu(mu)
            fit_mle(g, suffstats(g, x, w))
        else
            JohnsonSu(mu, sigma)
        end
    end
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

# Collection of sources required to build MPFRBuilder
sources = [
    "https://www.mpfr.org/mpfr-current/mpfr-4.0.1.tar.xz" =>
    "67874a60826303ee2fb6affc6dc0ddd3e749e9bfcb4c8655e3953d0458a6e16e",
]
name = "MPFR"
version = v"4.0.1"

# Bash recipe for building across all platforms
script = raw"""
cd mpfr-4.0.1
UNAME=`uname`
./configure --prefix=$prefix --host=$target --enable-shared --disable-static --with-gmp=$prefix
make -j
make install

"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libmpfr", :libmpfr)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    "https://github.com/JuliaMath/GMPBuilder/releases/download/v6.1.2-2/build_GMP.v6.1.2.jl"
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)


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

using JuliaFEM, Test

#=

@testset "test continuum 3d linear elasticity with surface load" begin
    nodes = Dict{Int64, Node}(
    1 => [0.0, 0.0, 0.0],
    2 => [1.0, 0.0, 0.0],
    3 => [1.0, 1.0, 0.0],
    4 => [0.0, 1.0, 0.0],
    5 => [0.0, 0.0, 1.0],
    6 => [1.0, 0.0, 1.0],
    7 => [1.0, 1.0, 1.0],
    8 => [0.0, 1.0, 1.0])

    element = Element(Hex8, [1, 2, 3, 4, 5, 6, 7, 8])
    update!([element], "geometry", nodes)
    update!([element], "youngs modulus", 200e3)
    update!([element], "poissons ratio", 1/3)

    plastic_parameters = Dict{Any, Any}("type" => JuliaFEM.ideal_plasticity!,
                                        "yield_surface" => Val{:von_mises},
                                        "params" => Dict("yield_stress" => 400.0))
    to_integ_points = Dict()
    map(x-> to_integ_points[x] = plastic_parameters, get_connectivity(element))
    update!(element, "plasticity", to_integ_points)

    elasticity_problem = Problem(Elasticity, "solve continuum block", 3)
    elasticity_problem.properties.finite_strain = false
    elasticity_problem.properties.geometric_stiffness = false
    push!(elasticity_problem.properties.store_fields, :plastic_strain)
    push!(elasticity_problem, element)

    bc = Element(Quad4, [1,4,8,5])
    update!([bc], "geometry", nodes)
    bc["displacement 1"] = 0.0
    bc["displacement 2"] = 0.0
    bc["displacement 3"] = 0.0
    boundary_problem = Problem(Dirichlet, "symmetry boundary conditions", 3, "displacement")
    push!(boundary_problem, bc)

    disp = Element(Quad4, [2,3,7,6])
    update!([disp], "geometry", nodes)
    disp["displacement 1"] = 0.002
    boundary_motion = Problem(Dirichlet, "displacement bc", 3, "displacement")
    push!(boundary_motion, disp)

    solver = NonlinearSolver("solve block problem")
    solver.time = 1.0
    push!(solver, elasticity_problem)
    push!(solver, boundary_problem)
    push!(solver, boundary_motion)
    solver()

    disp = element("displacement", [1.0, 1.0, 1.0], 1.0)
    @info("displacement at tip: $disp")
    u_expected = 2.0 * [-1/3, -1/3, 1.0]
     @test isapprox(disp, u_expected)
end

=#

# function solve_rod_model_elasticity(eltype)
#     fn = @__DIR__() * "/testdata/rod_short.med"
#     mesh = aster_read_mesh(fn, eltype)
#     element_sets = join(keys(mesh.element_sets), ", ")
#     @info("element sets: $element_sets")
#     p1 = Problem(Elasticity, "rod", 3)
#     p2 = Problem(Elasticity, "trac", 3)
#     p3 = Problem(Dirichlet, "fixed", 3, "displacement")
#     p4 = Problem(Dirichlet, "fixed", 3, "displacement")
#     p5 = Problem(Dirichlet, "fixed", 3, "displacement")
#     p1.elements = create_elements(mesh, "ROD")
#     p2.elements = create_elements(mesh, "FACE2")
#     p3.elements = create_elements(mesh, "FACE1")
#     p4.elements = create_elements(mesh, "FACE3")
#     p5.elements = create_elements(mesh, "FACE5")
#     update!(p1, "youngs modulus", 96.0)
#     update!(p1, "poissons ratio", 1/3)
#     update!(p2, "displacement traction force 1", 96.0)
#     update!(p3, "displacement 1", 0.0)
#     update!(p4, "displacement 2", 0.0)
#     update!(p5, "displacement 3", 0.0)
#     solver = LinearSolver(p1, p2, p3, p4, p5)
#     solver()
#     u_max = maximum(p1.assembly.u)
#     @info("$eltype, u_max = $u_max")
#     return u_max
# end
# @testset "compare 3d rod to CA solution" begin
#     @test isapprox(solve_rod_model_elasticity("Tet4"), 0.2)
#     @test isapprox(solve_rod_model_elasticity("Tet10"), 0.2)
#     @test isapprox(solve_rod_model_elasticity("Hex8"), 0.2)
#     @test isapprox(solve_rod_model_elasticity("Hex20"), 0.2)
#     @test isapprox(solve_rod_model_elasticity("Hex27"), 0.2)
# end


text:
struct ProportionalControlSolver <: Solver
    k::Float64
end

ProportionalControlSolver() = ProportionalControlSolver(1.0)

struct ProportionalControlPlanner <: Policy
    k::Float64
end

function POMDPs.solve(pc::ProportionalControlSolver, pomdp)
    return ProportionalControlPlanner(pc.k)
end

function BasicPOMCP.convert_estimator(pc::ProportionalControlSolver, ::Any, pomdp)
    return solve(pc, pomdp)
end

function POMDPs.action(pc::ProportionalControlPlanner, s::CovidState)
    I = infected(s)
    inf_prop = I / (s.S + I + s.R)
    return CovidAction(min(pc.k*inf_prop, 1.0))
end

function POMDPs.action(pc::ProportionalControlPlanner, b::ParticleCollection{<:CovidState})
    s̄ = Statistics.mean(b.particles)
    return POMDPs.action(pc, s̄)
end

function POMDPs.action(pc::ProportionalControlPlanner, b)
    b = ParticleCollection([rand(b) for _ in 1:100]) # TODO: change default particle size
    return POMDPs.action(pc, b)
end

function BasicPOMCP.estimate_value(
    estimator::ProportionalControlPlanner,
    pomdp::POMDP{S},
    s::S,
    node,
    depth::Int) where {S}

    disc = 1.0
    r_total = 0.0
    step = 0
    γ = discount(pomdp)

    while !isterminal(pomdp, s) && step < depth

        a = action(estimator, s)

        sp,o,r = gen(pomdp, s, a)

        r_total += disc*r

        s = sp

        disc *= γ
        step += 1
    end

    return r_total
end


text:
using ImageMorphology
using ImageCore
using ImageCore.PaddedViews
using Test
using OffsetArrays
using ImageMetadata
using Suppressor

@test isempty(detect_ambiguities(ImageMorphology))

using Documenter
Base.VERSION >= v"1.8" && doctest(ImageMorphology; manual=false)

include("testutils.jl")
@testset "ImageMorphology" begin
    include("structuring_element.jl")
    include("utils.jl")
    include("extreme_filter.jl")
    include("convexhull.jl")
    include("connected.jl")
    include("ops/dilate.jl")
    include("ops/erode.jl")
    include("ops/opening.jl")
    include("ops/closing.jl")
    include("ops/tophat.jl")
    include("ops/bothat.jl")
    include("ops/mgradient.jl")
    include("ops/mlaplacian.jl")
    include("isboundary.jl")
    include("thinning.jl")
    include("imfill.jl")
    include("maxtree.jl")
    include("feature_transform.jl")
    include("clearborder.jl")
    @info "Beginning deprecation tests, warnings are expected"
    include("deprecations.jl")
end

# multithreaded.jl runs on CI


text:
module OutlierDetection
    using MLJModelInterface
    using NearestNeighbors
    using MacroTools
    using Distances
    using Requires:@require
    const MMI = MLJModelInterface
    const NN = NearestNeighbors
    const DI = Distances
    const CLASS_NORMAL = 1
    const CLASS_OUTLIER = -1
    const CLASS_UNKNOWN = 0

    # base
    export Detector,
           UnsupervisedDetector,
           SupervisedDetector,
           Model,
           Scores,
           Data,
           Labels,
           Fit,
           Result,
           fit,
           score,
           detect

    # models
    export DNN,
           DNNModel,
           KNN,
           KNNModel,
           LOF,
           LOFModel,
           COF,
           COFModel,
           ABOD,
           ABODModel,
           AE,
           AEModel,
           DeepSAD,
           DeepSADModel,
           ESAD,
           ESADModel

    export PyABOD,
           PyCBLOF,
           PyCOF,
           PyCOPOD,
           PyHBOS,
           PyIForest,
           PyKNN,
           PyLMDD,
           PyLODA,
           PyLOF,
           PyLOCI,
           PyMCD,
           PyOCSVM,
           PyPCA,
           PyROD,
           PySOD,
           PySOS,
           PyModel

    # evaluation
    export classify,
           normalize,
           unify,
           combine,
           Evaluator,
           Score,
           Class

    # basic types
    include("base.jl")

    # utilities
    include("utils/neighbors.jl")
    include("utils/neural.jl")

    # detectors
    include("models/abod.jl")
    include("models/ae.jl")
    include("models/cof.jl")
    include("models/deepsad.jl")
    include("models/dnn.jl")
    include("models/esad.jl")
    include("models/knn.jl")
    include("models/lof.jl")
    include("pymodels/utils.jl")
    include("pymodels/detectors.jl")

    # evaluation
    include("evaluate/utils.jl")
    include("evaluate/evaluate.jl")

    # examples
    include("examples/examples.jl")

    function __init__()
       @require MLJ="add582a8-e3ab-11e8-2d5e-e98b27df1bc7" begin
           include("extension/mlj.jl")
           include("extension/mlj_extra.jl")
       end
    end
end


text:
##
# sparse sampling solver environments

abstract type AbstractVSSOptions end

# accessors
maxdepth(opt::AbstractVSSOptions) = opt.maxdepth
width(opt::AbstractVSSOptions) = opt.width
action_width(opt::AbstractVSSOptions) = opt.action_width

@with_kw struct VOSSOptions{V<:VOOActionGenerator} <: AbstractVSSOptions
    maxdepth::Int
    width::Int
    action_width::Int
    action_width_decay::Float64
    last_action_null::Bool
    voo::V
end

@with_kw struct VOWSSOptions{V<:VOOActionGenerator} <: AbstractVSSOptions
    maxdepth::Int
    width::Int
    action_width::Int
    action_width_decay::Float64
    last_action_null::Bool
    voo::V
end

struct VOOSparseSamplingSolver
    opt::AbstractVSSOptions
    rng::AbstractRNG
end

VOOSparseSamplingSolver(opt::AbstractVSSOptions; rng=Random.GLOBAL_RNG) = VOOSparseSamplingSolver(opt, rng)
VOOSparseSamplingSolver(;rng=Random.GLOBAL_RNG, kwargs...) = VOOSparseSamplingSolver(VOWSSOptions(kwargs...), rng)

function POMDPs.solve(s::VOOSparseSamplingSolver, m::Union{MDP,POMDP})
    if s.opt isa VOWSSOptions
        @assert !(Nothing <: statetype(m)) "POWSS does not support problems where the state can be nothing (this can be fixed by adding Some in appropriate places)."
    end
    VSSPlanner(m, s.opt, s.rng)
end

struct VSSPlanner{M<:Union{MDP,POMDP},O<:AbstractVSSOptions,R<:AbstractRNG}
    m::M
    opt::O
    rng::R
end


text:
#=
testODE:
- Julia version:
- Author: ymocquar
- Date: 2019-11-13
=#
using HOODESolver
using DifferentialEquations
using LinearAlgebra
using Plots
using Random

function twoscales_solve( par_u0::PrepareU0, order, t, nb)

    parGen = PrepareTwoScalePureAB(nb, t, order, par_u0)

    return twoscales_pure_ab(parGen, only_end=true)

end

function getmindif(tab::Vector{Vector{BigFloat}})
    nmmin = Inf
    ret = (1, 2)
    for i=1:size(tab,1)
        for j=(i+1):size(tab,1)
            nm = norm(tab[i]-tab[j], Inf)
            if nm < nmmin
                nmmin = nm
                ret = i, j
            end
        end
    end
    return ret, min(nmmin,1.1)
end

function fctMain(n_tau)

 #   u0 =[big"0.12345678", big"0.1209182736", big"0.1290582671", big"0.1239681094" ]
    seed=12996
    Random.seed!(seed)
    u0=rand(BigFloat,4)
    println("seed = $seed")
 #   tab_eps = zeros(BigFloat,7)
 #   tab_eps= [big"0.5", big"0.2", big"0.1",big"0.05", big"0.02", big"0.01",big"0.005", big"0.002", big"0.001",big"0.0005", big"0.0002", big"0.0001"]
    tab_eps= [big"0.1",big"0.01", big"0.001",big"0.0001", big"0.00001", big"0.000001", big"0.0000001", big"0.00000001", big"0.000000001"]
    # epsilon=big"0.1"
    # for i=1:7
    #     tab_eps[i] = epsilon
    #     epsilon /= 100
    # end
    nbmaxtest=11
    t_max = big"1.0"
    A=[0 0 1 0; 0 0 0 0;-1 0 0 0; 0 0 0 0]
    tabtabsol = Vector{Vector{Vector{BigFloat}}}()
    for order=2:14
        ordprep=order+2
        y = ones(Float64, nbmaxtest, size(tab_eps,1) )
        y_big = ones(BigFloat, nbmaxtest, size(tab_eps,1) )
        x=zeros(Float64,nbmaxtest)
        ind=1
        for epsilon in tab_eps
    #        fct = u -> [ 1.11u[2]^2+0.8247u[4]^3+0.12647u[3]^4,
    #       0.4356789141u[1]-0.87u[3]*u[4],
    #       1.9898985/(1+1.1237u[4]^3), 0.97u[1]*u[3]+0.8111u[2]^2 ]
            fct = (u,p,t) ->  [
    u[3]^3 + sin(t^2),
    u[4] + u[1]^2 + 1/(1+exp(t))-0.35,
    u[1]*u[2]*u[4]*(1+t),
    u[2] - u[2]^2 + 1/(1+u[4]^2) + u[1]^2 + cos(exp(t))
]
            parphi = PreparePhi(epsilon, n_tau, A, fct)
            println("prepareU0 eps=$epsilon n_tau=$n_tau")
            nb=10
            @time par_u0 = PrepareU0(parphi, ordprep, u0)
            @time pargen = PrepareTwoScalesPureAB(nb*2^nbmaxtest, t_max, order, par_u0)
            @time solref = twoscales_pure_ab(pargen, only_end=true)
            if size(tabtabsol,1) < ind
                push!(tabtabsol,Vector{Vector{BigFloat}}())
            end
            tabsol = tabtabsol[ind]
            push!(tabsol, solref)
            res_gen = Array{ Array{BigFloat,1}, 1}(undef, nbmaxtest)
            indref = 1
            eps_v = convert(Float32,epsilon)
            println("epsilon = $eps_v solref=$solref")
            indc =1
            labels=Array{String,2}(undef, 1, ind)
            while indc <= nbmaxtest
                @time pargen = PrepareTwoScalesPureAB(nb, t_max, order, par_u0)
                @time sol= twoscales_pure_ab(pargen, only_end=true)
                push!(tabsol, sol)
                res_gen[indc] = sol
                diff=solref-sol
                println("tabsol=$tabsol")
                (a, b), nm = getmindif(tabsol)
                if a != indref
                    println("New solref !!!! a=$a, b=$b nm=$nm")
                    indref = a
                    solref = tabsol[a]
                    for i=1:indc
                        nm2 = min( norm(res_gen[i] - solref, Inf), 1.1)
                        y[i, ind] = nm2 == 0 ? nm : nm2
                        y_big[i, ind] = nm2 == 0 ? nm : nm2
                    end
                else
                    diff=solref-sol
                    y[indc,ind] = min(norm(diff,Inf), 1.1)
                    y_big[indc,ind] = min(norm(diff,Inf), 1.1)
                end
                println("solref=$solref")
                println("nb=$nb sol=$sol")
                x[indc] = 1.0/nb
                println("epsilon=$epsilon result=$y")
                println("epsilon=$epsilon reslog2=$(log2.(y))")
                println("epsilon=$epsilon result=$y_big")
                println("epsilon=$epsilon reslog2=$(log2.(y_big))")
                nb *= 2
                indc += 1
            end
            if ind == size(tab_eps,1)
                for i=1:ind
                    labels[1,i] = " epsilon,order=$(convert(Float32,tab_eps[i])),$order "
                end
                gr()
                p=Plots.plot(
            x,
            view(y,:,1:ind),
            xlabel="delta t",
            xaxis=:log,
            ylabel="error",
            yaxis=:log,
            legend=:bottomright,
            label=labels,
            marker=2
        )

                prec_v = precision(BigFloat)
                Plots.savefig(p,"out/r20_$(prec_v)_$(eps_v)_$(order)_$(ordprep)_$(n_tau)_epsilon_fct.pdf")
            end
            ind+= 1
        end
    end
end

# testODESolverEps()

# for i=3:9
#     fctMain(2^i)
# end
# setprecision(512)
fctMain(32)


text:
module NanchenData

    using Statistics
    using ProjAssistant
    @gen_sub_proj

    include("data.jl")

    function __init__()
        _populate_bundle()
        @create_proj_dirs
    end

end


text:
@testset "FluxGradient" begin
    ir = 24
    Tair,pressure,VPD,LE = tha48[ir, Cols(:Tair,:pressure,:VPD,:LE)]
    Gs24 = surface_conductance(FluxGradient(), Tair,pressure,VPD,LE)
    @test keys(Gs24) == (:Gs_ms, :Gs_mol)
    # compare with R
    @test all(isapprox.(values(Gs24), (0.00919248, 0.3751482), rtol=1e-3))
    #
    df = copy(tha48)
    surface_conductance!(df, FluxGradient())
    @test propertynames(df)[(end-1):end] == SA[:Gs_ms, :Gs_mol]
    @test df[ir, Cols(:Gs_ms, :Gs_mol)] == Gs24
end

@testset "PenmanMonteith" begin
    ir = 24
    Tair,pressure,VPD,LE,Rn,Ga_h,G = tha48[ir, Cols(:Tair,:pressure,:VPD,:LE, :Rn,:Ga_h,:G)]
    Gs24 = surface_conductance(InversePenmanMonteith(), Tair,pressure,VPD,LE,Rn,Ga_h;G)
    @test keys(Gs24) == (:Gs_ms, :Gs_mol)
    # compare with R
    @test all(isapprox.(values(Gs24), (0.006846274, 0.2793988), rtol=1e-3))
    #
    df = copy(tha48)
    surface_conductance!(df, InversePenmanMonteith(); G=df.G)
    @test propertynames(df)[(end-1):end] == SA[:Gs_ms, :Gs_mol]
    @test df[ir, Cols(:Gs_ms, :Gs_mol)] == Gs24
end


text:
# Model
include_model("particle")
include_constraints("contact_no_slip")
model = no_slip_model(model)

# Horizon
T = 51

# Time step
tf = 5.0
h = tf / (T-1)

# Bounds
_uu = Inf * ones(model.m)
_uu[model.idx_u] .= 0.0
_ul = zeros(model.m)
ul, uu = control_bounds(model, T, _ul, _uu)

# Initial and final states
q1 = [0.0, 0.0, 1.0]
v1 = [1.0, 0.0, 0.0]
v2 = v1 - G_func(model,q1) * h
q2 = q1 + 0.5 * h * (v1 + v2)

x1 = [q1; q2]

xl, xu = state_bounds(model, T, x1 = x1)

# Objective
obj = PenaltyObjective(1.0e5, model.m)

# Constraints
con_contact = contact_no_slip_constraints(model, T)

# Problem
prob = trajectory_optimization_problem(model,
               obj,
               T,
               h = h,
               xl = xl,
               xu = xu,
               ul = ul,
               uu = uu,
               con = con_contact
               )

# Trajectory initialization
x0 = [0.01 * rand(model.n) for t = 1:T] #linear_interpolation(x1, x1, T) # linear interpolation on state
u0 = [0.001 * rand(model.m) for t = 1:T-1] # random controls

# Pack trajectories into vector
z0 = pack(x0, u0, prob)

#NOTE: may need to run examples multiple times to get good trajectories
# Solve nominal problem
@time z̄ , info = solve(prob, copy(z0), tol = 1.0e-5, c_tol = 1.0e-5)

check_slack(z̄, prob)
x̄, ū = unpack(z̄, prob)

t = 25
q3 = view(x̄[t], model.nq .+ (1:model.nq))
q2 = view(x̄[t], 1:model.nq)
λ = view(ū[t], model.idx_λ)
s = view(ū[t], model.idx_s)
λ_stack = λ[1] * ones(2)
s[1] - (λ_stack' * _P_func(model, q3) * (q3 - q2) / h)[1]

include(joinpath(pwd(), "models/visualize.jl"))
vis = Visualizer()
open(vis)
visualize!(vis, model, state_to_configuration(x̄), Δt = h)


text:
using Glimpse
const Gl = Glimpse
import Glimpse: System, Component, ComponentData, Spatial, TimingData, valid_entities,
                RGBAf0
using LinearAlgebra

# define rotation component and system
# add new component and system to the diorama, alongside with the already present rendering systems/components
push_system(dio, Gl.Rotator());
# add some sphere entities and the new extra component which will automatically engage with the new system
geom = Gl.PolygonGeometry(Sphere(Point3f0(0.0), 1.0f0))
progtag = Gl.ProgramTag{Gl.PeelingProgram}()
for i in 1:2
    Entity(dio, Gl.UniformColor(rand(RGBAf0)),
           Gl.Spatial(50 * (1 .- rand(Point3f0)), rand(Vec3f0)), Gl.Material(), Gl.Shape(),
           Gl.Rotation(0.2f0, 50 * (1 .- rand(Point3f0)), normalize(1 .- rand(Vec3f0))),
           Gl.Dynamic(), progtag, geom)
end


text:
using CommonCrawl
using URIParser

cache_dir = joinpath(tempdir(), "cc")
!isdir(cache_dir) && mkdir(cache_dir)
cc = CrawlCorpus(cache_dir, true)

try
    segs = segments(cc)
    @assert !isempty(segs)
    for seg in segs[1:5]
        println("\t$(seg)")
    end
    println("\t... total $(length(segs))")

    println("fetching 10 archives...")
    arcs = archives(cc, 10)
    @assert length(arcs) == 10
    for arc in arcs
        println("\t$(arc)")
    end

    println("reading archive $(arcs[1]) for 10 entries")
    f = open(cc, arcs[1])
    text_entries = read_entries(cc, f, "text/", 10)
    @assert length(text_entries) == 10
    for entry in text_entries
        println("\t$(entry.mime) : $(length(entry.data)) : $(entry.uri)")
    end
    close(f)
finally
    clear_cache(cc)
    rm(cache_dir)
end


text:
using GLM, QuantileRegression

# Load data
url = "http://vincentarelbundock.github.io/Rdatasets/csv/quantreg/engel.csv"
# TODO: Make automatic URL downloads work
dat = readtable("engel.csv")

# Fit least absolute deviation model (quantile =.5)
res = qreg(:(foodexp~income), dat, .5)
describe(res)

# Fit quantile regression for a bunch of different quantiles
dat_plot = [summary(qreg(:(foodexp ~ income), dat, i/10))[2,:] for i in 1:9]
dat_plot = reduce(vcat, dat_plot)

# Fit OLS model to compare
res_lm = lm(:(foodexp~income), dat)
ols = coeftable(res_lm)[2,1]
ols = rep(ols, 9)

# Plot results
using Winston

x = [i/10 for i=1:9]
y = dat_plot["Estimate"]
p = FramedPlot()
add(p, Curve(x, dat_plot["Estimate"]))
add(p, Curve(x, dat_plot["2.5%"], "type", "dash"))
add(p, Curve(x, dat_plot["97.5%"], "type", "dash"))
add(p, Curve(x, ols, "color", "red"))
setattr(p, "title", "Quantile regression: Food Expenditure ~ Income\n Red bar: OLS coefficient")
setattr(p, "xlabel", "Quantiles")
setattr(p, "ylabel", "Coefficient on Income")
file(p, "/Users/USERNAME/qreg_example_plot.png")

# > summary(res)
# 2x5 DataFrame:
#         Estimate Std.Error t value     2.5%    97.5%
# [1,]     81.4823   14.6345 5.56783  52.7987  110.166
# [2,]    0.560181 0.0131756 42.5164 0.534356 0.586005
# insheet using engel.csv
# qreg foodexp income, vce(iid, kernel(epan2))

# . qreg foodexp income, vce(iid, kernel(epan2))
# Iteration  1:  WLS sum of weighted deviations =  18051.196
#
# Iteration  1: sum of abs. weighted deviations =  18035.127
# Iteration  2: sum of abs. weighted deviations =  17815.249
# Iteration  3: sum of abs. weighted deviations =  17567.085
# Iteration  4: sum of abs. weighted deviations =  17559.932
#
# Median regression                                    Number of obs =       235
#   Raw sum of deviations 46278.06 (about 582.54126)
#   Min sum of deviations 17559.93                     Pseudo R2     =    0.6206
#
# ------------------------------------------------------------------------------
#      foodexp |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
# -------------+----------------------------------------------------------------
#       income |   .5601805   .0131763    42.51   0.000     .5342206    .5861403
#        _cons |   81.48233   14.63518     5.57   0.000     52.64815    110.3165
# ------------------------------------------------------------------------------


text:
using BinDeps

@BinDeps.setup

rsvg = library_dependency("rsvg", aliases = ["librsvg","librsvg-2.so.2"])

@linux_only begin
    provides(AptGet, "librsvg2-2",rsvg)
end

@BinDeps.install Dict(:rsvg => :librsvg)


text:
using Test

using Distributed
addprocs(1)

using Logging: NullLogger, with_logger
silence(f) = with_logger(f, NullLogger())

using Outsource: Connector, outsource

@testset "Connector" begin
    con = Connector()
    rcon = reverse(con)

    @test isopen(con) && isopen(rcon)

    # input
    x = rand()
    res = @async put!(con, x)
    @test take!(rcon) == x
    @test fetch(res) === con

    # output
    y = rand()
    res = @async put!(rcon, y)
    @test take!(con) == y
    @test fetch(res) === rcon

    @test isnothing(close(con))
    @test !isopen(con) && !isopen(rcon)
end

f() = nothing
@testset "outsource" begin
    silence() do
        con = outsource(_ -> nothing, 1)
        @test isa(con, Connector)
    end
    @test_throws Exception outsource(f, 2)
    # TODO: test async remote error handling
end


text:
@testset "void types" begin
    @test scitype(nothing) == Nothing
    @test scitype(missing) == Missing
    @test scitype([nothing, nothing]) == AbstractVector{Nothing}
end

@testset "scitype" begin
    X = [1, 2, 3]
    @test scitype(X) == AbstractVector{Unknown}

    @test scitype(missing) == Missing
    @test scitype((5, 2)) == Tuple{Unknown,Unknown}
    anyv = Any[5]
    @test scitype(anyv[1]) == Unknown

    X = [missing, 1, 2, 3]
    @test scitype(X) == AbstractVector{Union{Missing, Unknown}}

    Xnm = X[2:end]
    @test scitype(Xnm) == AbstractVector{Union{Missing, Unknown}}
    @test scitype(Xnm, tight=true) == AbstractVector{Unknown}

    Xm = Any[missing, missing]
    @test scitype(Xm) == AbstractVector{Missing}

    @test scitype([missing, missing]) == AbstractVector{Missing}
end

@testset "scitype2" begin
    ScientificTypesBase.Scitype(::Type{<:Integer}, ::MockMLJ) = Count
    X = [1, 2, 3]
    @test scitype(X) == AbstractVector{Count}
    Xm = [missing, 1, 2, 3]
    @test scitype(Xm) == AbstractVector{Union{Missing,Count}}
    Xnm = Xm[2:end]
    @test scitype(Xnm) == AbstractVector{Union{Missing,Count}}
    @test scitype(Xnm; tight=true) == AbstractVector{Count}

    @test elscitype(X) == Count
    @test elscitype(Xm) == Union{Missing,Count}
    @test elscitype(Xnm) == Union{Missing,Count}
    @test elscitype(Xnm, tight=true) == Count
end

@testset "temporal types" begin
    @test ScientificDate <: ScientificTimeType
    @test ScientificDateTime <: ScientificTimeType
    @test ScientificTime <: ScientificTimeType
end

@testset "Empty array" begin
    set_convention(MockMLJ())
    ScientificTypesBase.Scitype(::Type{<:Integer}, ::MockMLJ) = Count
    ScientificTypesBase.Scitype(::Type{Missing}, ::MockMLJ) = Missing
    @test scitype(Int[]) == AbstractVector{Count}
    @test scitype(Any[]) == AbstractVector{Unknown}
    @test scitype(Missing[]) == AbstractVector{Missing}
    @test scitype(Vector{Union{Int,Missing}}()) == AbstractVector{Union{Missing,Count}}
end


text:
using Documenter, MegTools

makedocs(;
    modules=[MegTools],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
        "DocStrings" => "functions.md",
        "Simple Pipeline" => "pipelines/testingpreproc.md",

    ],
    repo="https://github.com/AsimHDar/MegTools.jl/blob/{commit}{path}#L{line}",
    sitename="MegTools.jl",
    authors="Asim H. Dar",
    assets=String[],
)

deploydocs(;
    repo="github.com/AsimHDar/MegTools.jl",
)

max_stars_repo_name: ElectronicTeaCup/MegTools

max_issues_repo_name: ElectronicTeaCup/MegTools

max_forks_repo_name: ElectronicTeaCup/MegTools


text:
using Revise, Parameters, Logging, Test
using CUDA
includet("C:\\GitHub\\GitHub\\NuclearMedEval\\src\\utils\\CUDAAtomicUtils.jl")
#includet("C:\\GitHub\\GitHub\\NuclearMedEval\\src\\kernelEvolutions.jl")
includet("C:\\GitHub\\GitHub\\NuclearMedEval\\src\\structs\\BasicStructs.jl")
includet("C:\\GitHub\\GitHub\\NuclearMedEval\\src\\utils\\CUDAGpuUtils.jl")
includet("C:\\GitHub\\GitHub\\NuclearMedEval\\src\\utils\\IterationUtils.jl")
includet("C:\\GitHub\\GitHub\\NuclearMedEval\\src\\utils\\ReductionUtils.jl")
includet("C:\\GitHub\\GitHub\\NuclearMedEval\\src\\utils\\MemoryUtils.jl")
#includet("C:\\GitHub\\GitHub\\NuclearMedEval\\src\\distanceMetrics\\MeansMahalinobis.jl")
includet("C:/GitHub/GitHub/NuclearMedEval/src/distanceMetrics/Housdorff/verB/PrepareArrtoBool.jl")
includet("C:/GitHub/GitHub/NuclearMedEval/src/distanceMetrics/Housdorff/verB/MetaDataUtils.jl")
using ..PrepareArrtoBool, ..CUDAGpuUtils, ..PrepareArrtoBool, ..CUDAAtomicUtils, ..MetaDataUtils
using CUDA

localQuesValues = CUDA.zeros(Float32,14)
threads=(32,5)
blocks =3

mainArrDims= (67,78,90)
dataBdim = (17,7,12)

metaDataDims= (cld(mainArrDims[1],dataBdim[1] ),cld(mainArrDims[2],dataBdim[2]),cld(mainArrDims[3],dataBdim[3]))
metaData = MetaDataUtils.allocateMetadata(mainArrDims,dataBdim);
size(metaData)


text:
# -*- encoding: utf-8 -*-
#
# The MIT License (MIT)
#
# Copyright © 2021 Matteo Foglieni and Riccardo Gervasoni
#

function check_is_positive(string::String="")
	var = filter(x -> !isspace(x) && x≠"\"", string)
     (var == "") && (return true)

	!isnothing(tryparse(Float64, var)) || (return false)

     return parse(Float64, var)≥0.0 ? true : false
end

check_is_positive(number::Number) = check_is_positive(string(number))

"""
     check_is_positive(string::String="") :: Bool
     check_is_positive(number::Number) = check_is_positive(string(number))

Checks if the input `string` is a number that can be parsed as
a positive Float64.
"""
check_is_positive

function string2positive(string::String, uint::Bool=false)
     if check_is_positive(string)==false
          throw(ArgumentError(
               "invalid number; it must be parsed as a positive Int64"
          ))
     end
     var = filter(x -> !isspace(x) && x≠"\"", string)

     !(var=="") || (return 0)

	return parse(Float64, var)
end

string2positive(number::Number, uint::Bool=false) = string2positive(string(number), uint)

"""
     string2positive(string::String) :: Float64
     string2positive(number::Number, uint::Bool=false) =
          string2positive(string(number), uint)

Checks if the input `string` is a number that can be parsed as
a positive Float64 with [`check_is_positive`](@ref), and return it as
a `Float64`.
"""
string2positive

##########################################################################################92

function check_is_uint64(string::String="")
	var = filter(x -> !isspace(x) && x≠"\"", string)
     (var == "") && (return true)

	!isnothing(tryparse(Float64, var)) || (return false)
     number = parse(Float64, var)
     (number - floor(number) ≈ 0.0 ) || (return false)

     return convert(Int64, number)>0 ? true : false
end

check_is_uint64(number::Number) = check_is_uint64(string(number))

"""
     check_is_uint64(string::String="") :: Bool
     check_is_uint64(number::Number) = check_is_uint64(string(number))

Checks if the input `string` is a number that can be parsed as
a positive Int64.
"""
check_is_uint64

function string2int64(string::String, uint::Bool=false)
     if check_is_uint64(string)==false
          throw(ArgumentError(
               "invalid number; it must be parsed as a positive Int64"
          ))
     end
     var = filter(x -> !isspace(x) && x≠"\"", string)

     !(var=="") || (return uint==false ? Int64(0) : UInt64(0) )

	return uint==false ?
          convert(Int64,parse(Float64, var)) :
          convert(UInt64,parse(Float64, var))
end

string2int64(number::Number, uint::Bool=false) = string2int64(string(number), uint)

"""
     string2int64(string::String, uint::Bool=false) :: Union{Int64, UInt64}
     string2int64(number::Number, uint::Bool=false) =
          string2int64(string(number), uint)

Checks if the input `string` is a number that can be parsed as
a positive Int64 with [`check_is_uint64`](@ref), and return it as
a `Int64` if `uint==false`or as a `UInt64` if `uint==true`.
"""
string2int64

##########################################################################################92

function check_is_even_uint64(string::String="")
	var = filter(x -> !isspace(x) && x≠"\"", string)
     (var == "") && (return true)

	!isnothing(tryparse(Float64, var)) || (return false)
     number = parse(Float64, var)
     (number - floor(number) ≈ 0.0 ) || (return false)
     even = convert(Int64, number)

     return (even>=0 && iseven(even)) ? true : false
end

check_is_even_uint64(number::Number) = check_is_even_uint64(string(number))

"""
     check_is_even_uint64(string::String="") :: Bool
     check_is_even_uint64(number::Number) = check_is_even_uint64(string(number))

Checks if the input `string` is a number that can be parsed as
an even positive Int64.
"""
check_is_even_uint64

function string2evenint64(string::String, uint::Bool=false)
     if check_is_even_uint64(string)==false
          throw(ArgumentError(
               "invalid number; it must be parsed as an even positive Int64"
          ))
     end
     var = filter(x -> !isspace(x) && x≠"\"", string)

     !(var=="") || (return uint==false ? Int64(0) : UInt64(0))

	return uint==false ?
          convert(Int64,parse(Float64, var)) :
          convert(UInt64,parse(Float64, var))
end

string2evenint64(number::Number, uint::Bool=false) = string2evenint64(string(number), uint)

"""
     string2evenint64(string::String, uint::Bool=false) :: Union{Int64, UInt64}
     string2evenint64(number::Number, uint::Bool=false) = string2evenint64(string(number), uint)

Checks if the input `string` is a number that can be parsed as
an even positive Int64 with [`check_is_even_uint64`](@ref), and return it as
a `Int64` if `uint==false`or as a `UInt64` if `uint==true`.
"""
string2evenint64

##########################################################################################92

function check_is_square(string::String="")
	var = filter(x -> !isspace(x) && x≠"\"", string)
     (var == "") && (return true)

	!isnothing(tryparse(Float64, var)) || (return false)
     square = parse(Float64,var)
     (square >= 0) || (return false)

     return √square - floor(√square) ≈ 0. ? true : false
end

check_is_square(number::Number) = check_is_square(string(number))

"""
     check_is_square(string::String="") :: Bool
     check_is_square(number::Number) = check_is_square(string(number))

Checks if the input `string` is a number that can be parsed as
a squared positive Int64.
"""
check_is_square

function string2rootint64(string::String)
     if check_is_square(string)==false
          throw(ArgumentError(
               "invalid number; it must be parsed as a squared positive Int64"
          ))
     end

     var = filter(x -> !isspace(x) && x≠"\"", string)
     !(var=="") || (return 0)
     return convert(Int64, √parse(Float64,var))
end

string2rootint64(number::Number) = string2rootint64(string(number))

"""
     string2rootint64(string::String) :: Int64
     string2rootint64(number::Number) = string2rootint64(string(number))

Checks if the input `string` is a number that can be parsed as
a squared positive Int64 with [`check_is_square`](@ref), and return
the square root as a `Int64`.
"""
string2rootint64

##########################################################################################92

"""
     check_is_color(string::String="") :: Bool

Checks if the input `string` is a color written in RGB components
as "<R, G, B>".
"""
function check_is_color(string::String="")
	color = filter(x -> !isspace(x) && x≠"\"", string)
     (color == "") && (return true)

	(color[begin] == '<' && color[end] == '>') || (return false)

	color = color[begin+1:end-1]
	color = split(color, ",")
	(length(color)==3) || (return false)

	for c in color
		if isnothing(tryparse(Float64, c)) && c ∉ keys(SYM_NUM)
               return false
          end
     end

	return true
end

"""
     string2color(string::String="") :: RGB{Float32}

Checks if the input `string` is a color written in RGB components
as "<R, G, B>" with [`check_is_color`](@ref), and return it.
"""
function string2color(string::String)
     if check_is_color(string)==false
          throw(ArgumentError(
               "invalid color sintax; must be: <R, G, B>\n"*
               "Example: --background_color=<1,2,3>"
          ))
     end

     var = filter(x -> !isspace(x) && x≠"\"", string)

     !(var=="") || (return RGB{Float32}(0,0,0) )

	rgb = Vector{String}(split(var[begin+1:end-1], ","))

     R = isnothing(tryparse(Float64, rgb[1])) ? SYM_NUM[rgb[1]] : parse(Float64, rgb[1])
     G = isnothing(tryparse(Float64, rgb[2])) ? SYM_NUM[rgb[2]] : parse(Float64, rgb[2])
     B = isnothing(tryparse(Float64, rgb[3])) ? SYM_NUM[rgb[3]] : parse(Float64, rgb[3])

     return RGB{Float32}(R,G,B)
end

##########################################################################################92

"""
     check_is_vector(string::String="") :: Bool

Checks if the input `string` is a vector written in X,Y,Z components
as "[X, Y, Z]".
"""
function check_is_vector(string::String="")
	vector = filter(x -> !isspace(x), string)
     (vector == "") && (return true)

	(vector[begin] == '[' && vector[end] == ']') || (return false)

	vector = vector[begin+1:end-1]
	vector = split(vector, ",")
	(length(vector)==3) || (return false)

	for c in vector
		if isnothing(tryparse(Float64, c)) && c ∉ keys(SYM_NUM)
               return false
          end
	end

	return true
end

"""
     string2vector(string::String="") :: Union{Vec, Nothing}

Checks if the input `string` is  a vector written in X,Y,Z components
as "[X, Y, Z]" with [`check_is_vector`](@ref), and return `Vec(X,Y,Z)`.

See also: [`Vec`](@ref)
"""
function string2vector(string::String)
     if check_is_vector(string)==false
          throw(ArgumentError(
               "invalid vector sintax; must be: [1,2,3]\n"*
               "Example: --camera_position=[1,2,3]"
          ))
     end

     var = filter(x -> !isspace(x) && x≠"\"", string)

     !(var=="") || ( return Vec(0.0, 0, 0) )

     vector = filter(x -> !isspace(x) && x≠"\"", var)[begin+1:end-1]
	vec = split(vector, ",")
     x = isnothing(tryparse(Float64, vec[1])) ? SYM_NUM[vec[1]] : parse(Float64, vec[1])
     y = isnothing(tryparse(Float64, vec[2])) ? SYM_NUM[vec[2]] : parse(Float64, vec[2])
     z = isnothing(tryparse(Float64, vec[3])) ? SYM_NUM[vec[3]] : parse(Float64, vec[3])

	return Vec(x, y, z)
end

##########################################################################################92

"""
     check_is_declare_float(string::String="")

Checks if the input `string` is a declaration of one (or more) floats
in the form "NAME:VALUE" with [`check_is_declare_float`](@ref).
Examples:
```bash
    --declare_float=name:1.0
    --declare_float=name1:1.0,name2:2.0
    --declare_float=" name1 : 1.0 , name2: 2.0"
```
"""
function check_is_declare_float(string::String="")
	string_without_spaces = filter(x -> !isspace(x), string)
     (string_without_spaces == "") && (return true)

	vec_nameval = split.(split(string_without_spaces, ","), ":" )
	for declare_float ∈ vec_nameval
		if !(length(declare_float)==2 && !isnothing(tryparse(Float64, declare_float[2])))
			return false
		end
	end

	return true
end

"""
     declare_float2dict(string::String) :: Union{Dict{String, Float64}, Nothing}

Checks if the input `string` is a declaration of one (or more) floats
in the form "NAME:VALUE" with [`check_is_declare_float`](@ref).
Return a `Dict{String, Float64}` that associates each NAME (as keys) with
its `Float64` value, or nothing if `string==""`.
"""
function declare_float2dict(string::String)
     if check_is_declare_float(string)==false
          throw(ArgumentError(
               "invalid declare_float usage. Correct usage: \n"*
               "\t--declare_float=name:1.0\n"*
               "\t--declare_float=name1:1.0,name2:2.0\n"*
               """\t--declare_float=" name1 : 1.0 , name2: 2.0\n"""
          ))
     end

     string_without_spaces = filter(x -> !isspace(x), string)

     !(string_without_spaces=="") || (return nothing)

     vec_nameval = split.(split(string_without_spaces, ","), ":" )
     declare_float = Dict{String, Float64}([v[1]=>parse(Float64, v[2]) for v in vec_nameval]...)
     return declare_float
end

##########################################################################################92

"""
    check_is_one_of(string::String, vec::Vector{String}) :: Bool

Checks if the input `string` is inside one of the strings
contained in `vec`.
"""
function check_is_one_of(string::String, vec::Vector{String})
	var = filter(x -> !isspace(x) && x≠"\"", string)
	return var ∈ vec ? true : false
end

"""
    string2stringoneof(string::String, vec::Vector{String}) :: String

Checks if the input `string` is inside one of the strings
contained in `vec` with [`check_is_one_of`](@ref), and return it.
"""
function string2stringoneof(string::String, vec::Vector{String})
     if check_is_one_of(string, vec)==false
          throw(ArgumentError("$(string) is not in $(vec)"))
     end

     return filter(x -> !isspace(x) && x≠"\"", string)
end

##########################################################################################92

"""
     check_is_iterable(string::String, type::Union{Type, Nothing} = nothing) :: Bool
     check_is_iterable(object::T, type::Union{Type, Nothing} = nothing) where T<:Any
          = check_is_iterable(string(object), type)

Checks if the input `string` can be parsed in a iterable object,
returning `true` if it is, otherwise `false`.
If specified an input `type`, check also if all the elements contained
in `object` are of a type `T` such that `T <: type`.
"""
function check_is_iterable(string::String, type::Union{Type, Nothing} = nothing)
     object = eval(Meta.parse(string))
     #applicable(length, object) || (return false)
     applicable(iterate, object) || (return false)

     next = iterate(object)
     while !isnothing(next)
          (val, state) = next

          if !isnothing(type)
               (typeof(val) <: type) ||  (return false)
          end

          next = iterate(object, state)
     end

     return true
end

check_is_iterable(object::T, type::Union{Type, Nothing} = nothing) where T<:Any = check_is_iterable(string(object), type)

"""
     string2iterable(string::String="") :: Vector{String}

Checks if the input `string` is a vector of variable names written
as "[namevar1, namevar2, ...]" with [`check_is_vec_variables`](@ref),
and return a `Vector{String} = ["namevar1", "namevar2", ...]`.
"""
function string2iterable(string::String, type::Union{Type, Nothing} = nothing)
     if check_is_iterable(string)==false
          throw(ArgumentError(
               "invalid input: it is not an iterable object"
          ))
     end

     if check_is_iterable(string, Float64)==false && check_is_iterable(string, Int64)==false
          throw(ArgumentError(
               "invalid input iterable: its elements cannot be converted to Float64 "
          ))
     end

	return eval(Meta.parse(string))
end

##########################################################################################92

"""
     check_is_vec_variables(string::String="") :: Bool

Checks if the input `string` is a vector of variable names written
as "[namevar1, namevar2, ...]".
"""
function check_is_vec_variables(string::String="")
	vector = filter(x -> !isspace(x), string)
     (vector == "") && (return true)

	(vector[begin] == '[' && vector[end] == ']') || (return false)

     ':' ∉ vector || (return false)

	vector = vector[begin+1:end-1]
	vector = split(vector, ",")
	(length(vector)≥1) || (return false)

	return true
end

"""
     string2vector(string::String="") :: Union{Vector{String}, Nothing}

Checks if the input `string` is a vector of variable names written
as "[namevar1, namevar2, ...]" with [`check_is_vec_variables`](@ref),
and return a `Vector{String} = ["namevar1", "namevar2", ...]`.
"""
function string2vec_variables(string::String)
     if check_is_vec_variables(string)==false
          throw(ArgumentError(
               "invalid vector of variables syntax; must be: [namevar1, namevar2, ...]\n"*
               "Example: --vec_variables=\"[ namevar1 , namevar2 , ...]\" "
          ))
     end

     vector = filter(x -> !isspace(x) && x≠"\"", string)
     !(vector=="") || (return nothing)

     vector = vector[begin+1:end-1]
     vector = Vector{String}(split(vector, ","))

	return vector
end

##########################################################################################92

"""
     check_is_function(string::String="") :: Bool

Checks if the input `string` is a function name defined in `Raytracing`.
"""
function check_is_function(string::String="")
	name = filter(x -> !isspace(x), string)
     (name == "") && (return false)

     isdefined(Raytracing, Symbol(name)) || (return false)
     isa(eval(Symbol(name)),  Function) || (return false)
     try
          func::Function = eval(Symbol(name))
     catch
          return false
     end

	return true
end

"""
     string2function(string::String="") :: Function

Checks if the input `string` is a function name defined in `Raytracing`
with [`check_is_function`](@ref), and return it.
"""
function string2function(string::String)
     if check_is_function(string)==false
          throw(ArgumentError("$(string) is not a function defined in Raytracing module"))
     end

	return  eval(Symbol(string))
end

##########################################################################################92

"""
     from_CLI_to_vecstring(string::String) :: Vector{String}

Parse a `string` as it would be from the Command Line, and return the
`Vector{String}` that contains all the commands parsed.
"""
function from_CLI_to_vecstring(string::String)
     cmd = lstrip(rstrip(string))
     if lstrip(cmd) == ""
          return Vector{String}()
     end

     arglist = Vector{String}()
     parsestack = Vector{Char}()
     insideliteral = false

     for i in 1:length(cmd)
          #global insideliteral
          islast = i >= length(cmd)
          if cmd[i] == ' ' && (insideliteral==true)
               # Whitespace within literal is kept
               push!(parsestack, cmd[i])
          elseif cmd[i] == ' '
               # Whitespace delimits arguments
               push!(arglist, join(parsestack))
               empty!(parsestack)
          elseif !islast && cmd[i]  == '\\' && cmd[i+1] == '\"'
               # Escaped double quote
               push!(parsestack, cmd[i+1])
               i+=1
          elseif cmd[i] == '\"' && insideliteral==false
               # Begin literal
               insideliteral = true
          elseif cmd[i]=='\"' && insideliteral==true
               # End literal
               insideliteral = false
          else
               push!(parsestack, cmd[i])
          end
     end
     push!(arglist, join(parsestack))
     empty!(parsestack)

     println(arglist)
     if arglist[1] == "./Raytracer.jl"
          return arglist[begin+1:end]
     else
          throw(ArgumentError("first argument must be './Raytracer.jl', not $(vec[1])"))
     end
end


text:
using Test
using LinearAlgebra
using SparseArrays
using SimplePartitions
using SimplePolynomials # Polynomials
# using Pkg
# Pkg.add(PackageSpec(name = "SimpleTools", rev = "master"))
# Pkg.resolve()
# using SimpleTools
using SimpleGraphs

@testset "Core" begin
    G = Path(10)
    name(G, "Example")
    @test name(G) == "Example"
    @test NV(G) == 10
    @test NE(G) == 9
    @test G[2] == [1, 3]
    @test deg(G, 2) == 2
    @test sum(deg(G)) == 2NE(G)

    G = IntGraph(2)
    add!(G, 1, 2)
    add!(G, 1, 3)
    add!(G, 2, 3)
    @test NV(G) == NE(G)
    A = adjacency(G)
    H = SimpleGraph(A)
    @test G == H
    H = SimpleGraph(convert(BitMatrix, A))
    @test G == H

    @test get_edge(G, 1, 2) == get_edge(G, 2, 1)

    G = StringGraph()
    add!(G, "alpha", "beta")
    @test G["alpha", "beta"]

    G = SimpleGraph{Complex{Float64}}()
    H = SimpleGraph{Complex{Float64}}()
    add!(G, im, -im)
    add!(H, -im, im)
    @test G == H
end

@testset "Ops" begin
    G = IntGraph(3)
    ee = [(1, 2), (2, 3), (1, 3)]
    add_edges!(G, ee)
    delete!(G, 1, 2)
    @test NE(G) == 2
    G = Cycle(5)
    delete!(G, 5)
    @test G == Path(4)

    G = line_graph(Cycle(5))
    @test deg(G) == 2 * ones(5)

    H = G'
    @test G == H'

    G = Cycle(4)
    H = cartesian(G, G)
    @test NV(H) == NV(G)^2

    H = disjoint_union(Cycle(4), Cycle(5))
    @test NV(H) == 9
    @test NE(H) == 9
    @test H == Cycle(4) + Cycle(5)

    @test 2H == H + H

    G = join(IntGraph(5), IntGraph(4))
    @test NE(G) == 20

    H = IntGraph(5) ∨ IntGraph(4)
    @test G == H

    G = Cycle(9)
    H = Cycle(10)
    contract!(H, 1, 10)
    @test G == H

    G = Path(9)
    H = trim(G, 1)
    @test NV(H) == 0

    G = Complete(5)
    H = subdivide(G)
    @test NV(H) == NV(G) + NE(G)
    @test NE(H) == 2 * NE(G)

end

@testset "Constructors" begin
    G = Complete(5)
    @test NE(G) == 10
    G = Complete(4, 5)
    @test NE(G) == 20
    G = Complete([5, 5, 5])
    @test NE(G') == 30

    G = RandomTree(10)
    @test NE(G) == 9
    G = Grid(3, 3)
    @test NV(G) == 9
    G = Wheel(10)
    @test NV(G) == 10
    G = Cube(4)
    @test NV(G) == 16
    G = BuckyBall()
    @test NE(G) == 90
    G = Petersen()
    @test NE(G) == 15
    G = Paley(13)
    @test NE(G) == NE(G')
    G = RandomRegular(10, 3)
    @test NE(G) == 15
    G = Knight(5, 5)
    @test NV(G) == 25
    @test NE(HoffmanSingleton()) == 175
    p1 = char_poly(Hoffman())
    p2 = char_poly(Cube(4))
    @test p1 == p2
end

@testset "Platonics" begin
    G = Icosahedron()
    H = Dodecahedron()
    @test NE(G) == NE(H)
    G = Octahedron()
    H = Complete([2, 2, 2])
    @test NE(G) == NE(H)
    @test Tetrahedron() == Complete(4)
end

@testset "Unit distance" begin
    G = Spindle()
    @test is_unit_distance(G)
    G = Golomb()
    @test is_unit_distance(G)
end

@testset "Connectivity" begin
    G = Path(10)
    @test diam(G) == 9
    @test eccentricity(G, 2) == 8
    @test num_components(G) == 1
    @test spanning_forest(G) == G
    @test is_acyclic(G)
    @test is_cut_edge(G, 3, 4)
    @test radius(G) == 5
    @test graph_center(G) == Set([5, 6])

    G = Complete(5, 5)'
    @test num_components(G) == 2
    H = spanning_forest(G)
    @test num_components(H) == 2
    @test NE(H) == 8

    G = Cycle(5) + Cycle(8) + Cycle(10)
    A = max_component(G)
    @test length(A) == 10
end

@testset "Matrices" begin
    G = Petersen()
    M = incidence(G)
    L = laplace(G)
    @test L == M * M'
    A = adjacency(G)
    v = A * ones(10)
    @test v == 3 * ones(Int, 10)
end

@testset "Twins" begin
    G = Complete(5, 3)
    @test twins(G) == bipartition(G)
end

@testset "Coloring" begin
    G = RandomTree(10)
    d = two_color(G)
    @test Set(values(d)) == Set([1, 2])
    f = greedy_color(G)
    @test length(keys(f)) == NV(G)
end

@testset "Euler" begin
    G = Cube(4)
    tour = euler(G)
    @test length(tour) == NE(G) + 1
    @test tour[1] == tour[end]
end

@testset "Hamiltonian" begin
    G = Cube(4)
    tour = hamiltonian_cycle(G)
    @test length(tour) == NV(G)
end

@testset "Girth" begin
    G = Cube(4)
    @test girth(G) == 4
    G = RandomTree(10)
    @test girth(G) == 0
end

@testset "Bisect" begin
    G = RandomTree(10)
    A, B = bisect(G)
    @test G.V == union(A, B)
    @test length(intersect(A, B)) == 0
    @test length(cross_edges(G, A, B)) == 1
end

@testset "Transitive" begin
    G = RandomTree(10)
    D = transitive_orientation(G)
    @test G == simplify(D)
    @test num_trans_orientations(G) == 2
end

@testset "Prufer" begin
    for k = 1:10
        G = RandomTree(10)
        code = prufer_code(G)
        H = prufer_restore(code)
        @test G == H
    end
end

@testset "Polynomials" begin
    G = Complete(3, 4)
    p = indep_poly(G)
    @test coeffs(p) == [1; 7; 9; 5; 1]

    p = matching_poly(G)
    @test coeffs(p) == [0, -24, 0, 36, 0, -12, 0, 1]

    p = interlace(G)
    @test coeffs(p) == [0, 2, 3, 4, 3, 1]

end

# DIRECTED GRAPH STUFF

@testset "Basic Directed" begin
    G = StringDigraph()
    add!(G, "alpha", "bravo")
    add!(G, "bravo", "charlie")
    @test NV(G) == 3
    @test NE(G) == 2
    @test sum(out_deg(G)) == sum(in_deg(G))

    H = simplify(relabel(G))
    @test H == Path(3)
end

@testset "Digraph Constructors" begin
    G = DirectedPath(10)
    add!(G, 10, 1)
    @test G == DirectedCycle(10)

    G = DirectedComplete(10, false)
    @test NE(G) == 10 * 9
    G = DirectedComplete(10, true)
    @test NE(G) == 10 * 10

    G = RandomTournament(10)
    @test NE(G) == 45
end

@testset "Directed Distance" begin
    G = DirectedCycle(10)
    @test diam(G) == 9
end

@testset "Directed Euler" begin
    G = TorusDigraph(4, 4)
    P = euler(G)
    @test length(P) == NE(G) + 1
end

@testset "Directed Hamiltonian Cycle" begin
    G = TorusDigraph(4, 4)
    P = hamiltonian_cycle(G)
    @test length(P) == NV(G)
end

@testset "Strong Connectivity" begin
    G = TorusDigraph(5, 5)
    @test is_strongly_connected(G)
end

@testset "Directed Matrices" begin
    G = RandomTournament(10)
    A = adjacency(G)
    B = A + A'
    @test B == adjacency(Complete(10))

    G = RandomTournament(10)
    M = incidence(G)
    @test M * M' == laplace(Complete(10))
end

@testset "Hypergraphs" begin
    H = IntHypergraph(2)

    add!(H, [1, 2, 3])
    add!(H, [3, 4])

    @test eltype(H) == Int

    @test sort(vlist(H)) == [1, 2, 3, 4]
    @test length(elist(H)) == 2

    @test NV(H) == 4
    @test NE(H) == 2
    @test has(H, 3)
    @test has(H, [2, 3, 1])

    G = SimpleGraph(H)
    @test is_connected(G)

    @test deg(H, 3) == 2

    add!(H, [4, 5])
    @test delete!(H, 4)
    @test NV(H) == 4
    @test NE(H) == 1

    add!(H, 4, 5, 6)
    @test is_uniform(H)

    G = Cycle(8)
    H = SimpleHypergraph(G)
    @test is_uniform(H)
    @test G.V == H.V
    @test NE(H) == NE(G)

    A = rand(6, 10) .> 0.5
    H = SimpleHypergraph(A)
    B = Float64.(A)
    K = SimpleHypergraph(B)
    @test H == K

    H = CompleteHypergraph(6, 3)
    @test NE(H) == 20

end

@testset "Rotation Systems" begin
    G = Grid(3, 3)
    @test euler_char(G) == 2
    G = Cycle(5)
    set_rot(G)
    @test euler_char(G) == 2
    G = Dodecahedron()
    @test euler_char(G) == 2
    H = dual(G)
    @test char_poly(Icosahedron()) == char_poly(H)
    @test euler_char(H) == 2
end


text:
# Test function taken from https://discourse.julialang.org/t/how-do-you-unfold-a-nested-julia-array/2243/3
# Helps flatten nested arrays
function flattenall(a)
    while any(x -> typeof(x) <: AbstractArray, a)
        a = collect(Iterators.flatten(a))
    end
    return a
end

mutable struct WeightChange{L<:Learner, A} <: IntrinsicReward
    learner::L
    previous_weights::A
    previous_reward::Float64
end

function WeightChange(learner::Learner)
    current_ws = flattenall(get_weights(learner))
    WeightChange(learner, deepcopy(current_ws), 0.0)
end

function update_reward!(wc::WeightChange, agent)
    cw = flattenall(get_weights(agent.demon_learner))
    curiosity_reward = sum(abs.(cw .- wc.previous_weights))
    wc.previous_weights .= deepcopy(cw)
    wc.previous_reward = curiosity_reward
    return curiosity_reward
end

function get_reward(wc::WeightChange, agent)
    return wc.previous_reward
end

mutable struct NoReward <: IntrinsicReward end

function update_reward!(::NoReward, agent)
    return 0
end

function get_reward(NR::NoReward, agent)
    return 0.0
end

max_stars_repo_name: MatthewMcLeod/curiosity

max_issues_repo_name: MatthewMcLeod/curiosity

max_forks_repo_name: MatthewMcLeod/curiosity


text:
module MyExample

include("extra_file.jl")

export sum_values

end


text:
# Заполняет словарь d индексами индивидуальных значений
function indsdict!(d::Dict{T}, cdata::Tuple) where T
    @inbounds for (i, element) in enumerate(zip(cdata...))
        ind = ht_keyindex(d, element)
        if ind > 0
            push!(d.vals[ind], i)
        else
            v = Vector{Int}(undef, 1)
            v[1] = i
            d[element] = v
        end
    end
    d
end

function sortbyvec!(a, vec)
    sort!(a, by = x -> findfirst(y -> x == y, vec))
end

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

isnanormissing(x::Number) = isnan(x)
isnanormissing(x::AbstractFloat) = isnan(x)
isnanormissing(x::Missing) = true

ispositive(::Missing) = false
ispositive(x::AbstractFloat) = isnan(x) ? false : x > zero(x)
ispositive(x) = x > zero(x)

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

# STATISTICS

#CV2se
"""
    sdfromcv(cv::Real)::AbstractFloat

LnSD from CV.
"""
function sdfromcv(cv)
    return sqrt(varfromcv(cv))
end
"""
    varfromcv(cv::Real)::AbstractFloat

LnVariance from CV.
"""
function varfromcv(cv)
     return log(1 + cv^2)
end
"""
    cvfromvar(σ²::Real)::AbstractFloat

CV from variance.
"""
function cvfromvar(σ²)
    return sqrt(exp(σ²) - 1)
end
#CV2se
"""
    cvfromsd(σ::Real)::AbstractFloat

CV from variance.
"""
function cvfromsd(σ)
    return sqrt(exp(σ^2) - 1)
end


text:
using NonlinearEigenproblemsTest
using NonlinearEigenproblems
using NonlinearEigenproblems.RKHelper
using Random
using SparseArrays
using Test

@bench @testset "cached_lin_solver" begin
    dep = nep_gallery("dep0")
    Random.seed!(0)
    y = randn(dep.n)

    cache = nothing

    function solve_and_verify(v, add_to_cache = true)
        x = solve(cache, v, y, add_to_cache)
        @test compute_Mder(dep, v) * x ≈ y
    end

    test_size(cache_size) = @test length(cache.solver) == cache_size

    function run_tests(cache_creator)
        cache = cache_creator(); test_size(0)

        solve_and_verify(2.5, false); test_size(0)
        solve_and_verify(2.5); test_size(1)
        solve_and_verify(2.5); test_size(1)
        solve_and_verify(3.5); test_size(2)

        cache = cache_creator(); test_size(0)

        solve_and_verify(2.5); test_size(1)
    end

    run_tests(() -> LinSolverCache(Float64, dep, DefaultLinSolverCreator()))
    run_tests(() -> LinSolverCache(Float64, dep, BackslashLinSolverCreator()))
end


text:
# ---
# title: 1569. Number of Ways to Reorder Array to Get Same BST
# id: problem1569
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/number-of-ways-to-reorder-array-to-get-same-bst/description/>
# hidden: true
# ---
#
# Given an array `nums` that represents a permutation of integers from `1` to
# `n`. We are going to construct a binary search tree (BST) by inserting the
# elements of `nums` in order into an initially empty BST. Find the number of
# different ways to reorder `nums` so that the constructed BST is identical to
# that formed from the original array `nums`.
#
# For example, given `nums = [2,1,3]`, we will have 2 as the root, 1 as a left
# child, and 3 as a right child. The array `[2,3,1]` also yields the same BST
# but `[3,2,1]` yields a different BST.
#
# Return _the number of ways to reorder_  `nums`  _such that the BST formed is
# identical to the original BST formed from_  `nums`.
#
# Since the answer may be very large,  **return it modulo  **`10^9 + 7`.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/08/12/bb.png)
#
#
#
#     Input: nums = [2,1,3]
#     Output: 1
#     Explanation: We can reorder nums to be [2,3,1] which will yield the same BST. There are no other ways to reorder nums which will yield the same BST.
#
#
# **Example 2:**
#
# **![](https://assets.leetcode.com/uploads/2020/08/12/ex1.png)**
#
#
#
#     Input: nums = [3,4,5,1,2]
#     Output: 5
#     Explanation: The following 5 arrays will yield the same BST:
#     [3,1,2,4,5]
#     [3,1,4,2,5]
#     [3,1,4,5,2]
#     [3,4,1,2,5]
#     [3,4,1,5,2]
#
#
# **Example 3:**
#
# **![](https://assets.leetcode.com/uploads/2020/08/12/ex4.png)**
#
#
#
#     Input: nums = [1,2,3]
#     Output: 0
#     Explanation: There are no other orderings of nums that will yield the same BST.
#
#
# **Example 4:**
#
# **![](https://assets.leetcode.com/uploads/2020/08/12/abc.png)**
#
#
#
#     Input: nums = [3,1,2,5,4,6]
#     Output: 19
#
#
# **Example 5:**
#
#
#
#     Input: nums = [9,4,2,1,3,6,5,7,8,14,11,10,12,13,16,15,17,18]
#     Output: 216212978
#     Explanation: The number of ways to reorder nums to get the same BST is 3216212999. Taking this number modulo 10^9 + 7 gives 216212978.
#
#
#
#
# **Constraints:**
#
#   * `1 <= nums.length <= 1000`
#   * `1 <= nums[i] <= nums.length`
#   * All integers in `nums` are  **distinct**.
#
#
## @lc code=start
using LeetCode

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