text:
#=
    Copyright (c) Yuki Ono.
    Licensed under the MIT License.
=#

module SendEML
    import Dates
    import TimeZones
    import Random
    import JSON
    import Sockets

    using Match

    const VERSION = 1.0

    const CR = UInt8('\r')
    const LF = UInt8('\n')
    const SPACE = UInt8(' ')
    const HTAB = UInt8('\t')
    const CRLF = "\r\n"

    const DATE_BYTES = Vector{UInt8}("Date:")
    const MESSAGE_ID_BYTES = Vector{UInt8}("Message-ID:")

    function find_cr(bytes::Vector{UInt8}, offset::Int)::Union{Int, Nothing}
        findnext(b -> b === CR, bytes, offset)
    end

    function find_lf(bytes::Vector{UInt8}, offset::Int)::Union{Int, Nothing}
        findnext(b -> b === LF, bytes, offset)
    end

    function find_all_lf(bytes::Vector{UInt8})::Vector{Int}
        indices = Int[]
        offset = 1
        while true
            @match find_lf(bytes, offset) begin
                _::Nothing => return indices
                idx::Int => begin
                    push!(indices, idx)
                    offset = idx + 1
                end
            end
        end
    end

    const Lines = Vector{Vector{UInt8}}

    function get_lines(bytes::Vector{UInt8})::Lines
        indices = find_all_lf(bytes)
        push!(indices, length(bytes))
        offset = 1
        map(i -> begin
            line = bytes[offset:i]
            offset = i + 1
            line
        end, indices)
    end

    function match_header(line::Vector{UInt8}, header::Vector{UInt8})::Bool
        if isempty(header)
            error("header is empty")
        end

        if length(line) < length(header)
            false
        else
            isequal(line[1:length(header)], header)
        end
    end

    function is_date_line(line::Vector{UInt8})::Bool
        match_header(line, DATE_BYTES)
    end

    function is_message_id_line(line::Vector{UInt8})::Bool
        match_header(line, MESSAGE_ID_BYTES)
    end

    function make_now_date_line()::String
        time = TimeZones.now(TimeZones.localzone())
        offset = replace(Dates.format(time, "zzzz"), ":" => "", count = 1)
        date_str = Dates.format(time, "eee, dd uuu yyyy HH:MM:SS ")
        "Date: " * date_str * offset * CRLF
    end

    function make_random_message_id_line()::String
        length = 62
        rand_str = Random.randstring(length)
        "Message-ID: <$rand_str>$CRLF"
    end

    function concat_bytes(lines::Lines)::Vector{UInt8}
        buf = Vector{UInt8}(undef, sum(l -> length(l), lines))
        offset = 1
        for l in lines
            copyto!(buf, offset, l)
            offset += length(l)
        end
        buf
    end

    function is_not_update(update_date::Bool, update_message_id::Bool)::Bool
        !update_date && !update_message_id
    end

    function is_wsp(b::UInt8)::Bool
        b === SPACE || b === HTAB
    end

    function first_byte(array::Vector{UInt8}, default::UInt8)::UInt8
        isempty(array) ? default : first(array)
    end

    function first_char(str::String, default::Char)::Char
        isempty(str) ? default : first(str)
    end

    function is_folded_line(bytes::Vector{UInt8})::Bool
        is_wsp(first_byte(bytes, UInt8(0)))
    end

    function replace_line(lines::Lines, match_line::Function, make_line::Function)::Lines
        idx = findnext(match_line, lines, 1)
        if isnothing(idx)
            return lines
        end

        p1 = Iterators.take(lines, idx - 1)
        p2 = [Vector{UInt8}(make_line())]
        p3 = Iterators.dropwhile(is_folded_line, Iterators.drop(lines, idx))

        collect(Iterators.flatten((p1, p2, p3)))
    end

    function replace_date_line(lines::Lines)::Lines
        replace_line(lines, is_date_line, make_now_date_line)
    end

    function replace_message_id_line(lines::Lines)::Lines
        replace_line(lines, is_message_id_line, make_random_message_id_line)
    end

    function replace_header(header::Vector{UInt8}, update_date::Bool, update_message_id::Bool)::Vector{UInt8}
        lines = get_lines(header)
        concat_bytes(@match (update_date, update_message_id) begin
            (true, true) => replace_message_id_line(replace_date_line(lines))
            (true, false) => replace_date_line(lines)
            (false, true) => replace_message_id_line(lines)
            (false, false) => lines
        end)
    end

    const EMPTY_LINE = [CR, LF, CR, LF]

    function combine_mail(header::Vector{UInt8}, body::Vector{UInt8})::Vector{UInt8}
        vcat(header, EMPTY_LINE, body)
    end

    function has_next_lf_cr_lf(bytes::Vector{UInt8}, idx::Int)::Bool
        if length(bytes) < (idx + 3)
            false
        else
            isequal(bytes[(idx + 1):(idx + 3)], [LF, CR, LF])
        end
    end

    function find_empty_line(bytes::Vector{UInt8})::Union{Int, Nothing}
        offset = 1
        while true
            @match find_cr(bytes, offset) begin
                _::Nothing => return nothing
                idx::Int, if has_next_lf_cr_lf(bytes, idx) end => return idx
                idx::Int => begin
                    offset = idx + 1
                end
            end
        end
    end

    function split_mail(bytes::Vector{UInt8})::Union{Tuple{Vector{UInt8}, Vector{UInt8}}, Nothing}
        idx = find_empty_line(bytes)
        if isnothing(idx)
            return nothing
        end

        header = bytes[1:(idx - 1)]
        body = bytes[(idx + length(EMPTY_LINE)):end]
        (header, body)
    end

    function replace_mail(bytes::Vector{UInt8}, update_date::Bool, update_message_id::Bool)::Union{Vector{UInt8}, Nothing}
        if is_not_update(update_date, update_message_id)
            return bytes
        end

        mail = split_mail(bytes)
        if isnothing(mail)
            return nothing
        end

        (header, body) = mail
        repl_header = replace_header(header, update_date, update_message_id)
        combine_mail(repl_header, body)
    end

    function make_json_sample()
        """{
            "smtpHost": "172.16.3.151",
            "smtpPort": 25,
            "fromAddress": "a001@ah62.example.jp",
            "toAddresses": [
                "a001@ah62.example.jp",
                "a002@ah62.example.jp",
                "a003@ah62.example.jp"
            ],
            "emlFiles": [
                "test1.eml",
                "test2.eml",
                "test3.eml"
            ],
            "updateDate": true,
            "updateMessageId": true,
            "useParallel": false
        }"""
    end

    function print_usage()
        println("Usage: {self} json_file ...")
        println("---")

        println("json_file sample:")
        println(make_json_sample());
    end

    function print_version()
        println("SendEML / Version: $VERSION")
    end

    const LAST_REPLY_REGEX = r"^\d{3} .+"

    function is_last_reply(line::String)::Bool
        !isnothing(match(LAST_REPLY_REGEX, line))
    end

    function is_positive_reply(line::String)::Bool
        return @match first_char(line, '0') begin
            '2' || '3' => true
            _ => false
        end
    end

    function make_id_prefix(use_parallel::Bool)::String
        use_parallel ? "id: $(Threads.threadid()), " : ""
    end

    function send_mail(sock::Sockets.TCPSocket, file::String, update_date::Bool, update_message_id::Bool, use_parallel::Bool = false)
        println(make_id_prefix(use_parallel) * "send: $file")

        mail = read(file)
        repl_mail = replace_mail(mail, update_date, update_message_id)
        if isnothing(repl_mail)
            println("error: Invalid mail: Disable updateDate, updateMessageId")
        end

        write(sock, something(repl_mail, mail))
        flush(sock)
    end

    function recv_line(sock::Sockets.TCPSocket, use_parallel::Bool = false)::String
        while true
            line = readline(sock)
            if isempty(line)
                error("Connection closed by foreign host")
            end

            println(make_id_prefix(use_parallel) * "recv: $line")

            if is_last_reply(line)
                if is_positive_reply(line)
                    return line
                end

                error(line)
            end
        end
    end

    function replace_crlf_dot(cmd::String)::String
        cmd === "$CRLF." ? "<CRLF>." : cmd
    end

    function send_line(sock::Sockets.TCPSocket, cmd::String, use_parallel::Bool = false)
        println(make_id_prefix(use_parallel) * "send: $(replace_crlf_dot(cmd))")

        write(sock, cmd * CRLF)
        flush(sock)
    end

    function make_send_cmd(sock::Sockets.TCPSocket, use_parallel::Bool)::Function
        cmd -> begin
            send_line(sock, cmd, use_parallel)
            recv_line(sock, use_parallel)
        end
    end

    function send_hello(send::Function)::String
        send("EHLO localhost")
    end

    function send_from(send::Function, from_addr::String)::String
        send("MAIL FROM: <$from_addr>")
    end

    function send_rcpt_to(send::Function, to_addrs::Vector{String})
        for addr in to_addrs
            send("RCPT TO: <$addr>")
        end
    end

    function send_data(send::Function)::String
        send("DATA")
    end

    function send_crlf_dot(send::Function)::String
        send("$CRLF.")
    end

    function send_rset(send::Function)::String
        send("RSET")
    end

    function send_quit(send::Function)::String
        send("QUIT")
    end

    struct Settings
        smtp_host::String
        smtp_port::Int
        from_address::String
        to_addresses::Vector{String}
        eml_files::Vector{String}
        update_date::Bool
        update_message_id::Bool
        use_parallel::Bool
    end

    function send_messages(settings::Settings, eml_files::Vector{String}, use_parallel::Bool)
        sock = Sockets.connect(settings.smtp_host, settings.smtp_port)
        send = make_send_cmd(sock, use_parallel)
        recv_line(sock, use_parallel)
        send_hello(send)

        try
            reset = false
            for file in eml_files
                if !isfile(file)
                    println("$file: EML file does not exist")
                    continue
                end

                if reset
                    println("---")
                    send_rset(send)
                end

                send_from(send, settings.from_address)
                send_rcpt_to(send, settings.to_addresses)
                send_data(send)
                send_mail(sock, file, settings.update_date, settings.update_message_id, use_parallel)

                send_crlf_dot(send)
                reset = true
            end
            send_quit(send)
        finally
            close(sock)
        end
    end

    function check_json_value(json::Dict{String, Any}, name::String, type::Type)
        if haskey(json, name)
            try
                convert(type, json[name])
            catch e
                error("$name: Invalid type: $(json[name])")
            end
        end
    end

    function check_settings(json::Dict{String, Any})
        names = ["smtpHost", "smtpPort", "fromAddress", "toAddresses", "emlFiles"];
        key_idx = findfirst(n -> !haskey(json, n), names);
        if !isnothing(key_idx)
            error("$(names[key_idx]) key does not exist")
        end

        check_json_value(json, "smtpHost", String)
        check_json_value(json, "smtpPort", Int)
        check_json_value(json, "fromAddress", String)
        check_json_value(json, "toAddresses", Vector{String})
        check_json_value(json, "emlFiles", Vector{String})
        check_json_value(json, "updateDate", Bool)
        check_json_value(json, "updateMessageId", Bool)
        check_json_value(json, "useParallel", Bool)
    end

    function get_settings_from_text(text::String)::Dict{String, Any}
        JSON.parse(text)
    end

    function get_settings(json_file::String)::Dict{String, Any}
        get_settings_from_text(read(json_file, String))
    end

    function map_settings(settings::Dict{String, Any})::Settings
        Settings(
            settings["smtpHost"],
            settings["smtpPort"],
            settings["fromAddress"],
            settings["toAddresses"],
            settings["emlFiles"],
            get(settings, "updateDate", true),
            get(settings, "updateMessageId", true),
            get(settings, "useParallel", false)
        )
    end

    function proc_json_file(json_file::String)
        if !isfile(json_file)
            error("JSON file does not exist")
        end

        # ! Avoid JSON.parsefile(): file not closed
        json = get_settings(json_file)
        check_settings(json)
        settings = map_settings(json)

        if settings.use_parallel && length(settings.eml_files) > 1
            if Threads.nthreads() === 1
                println("Threads.nthreads() === 1")
                println("    Windows: `set JULIA_NUM_THREADS=4` or `\$env:JULIA_NUM_THREADS=4`(PowerShell)")
                println("    Other: `export JULIA_NUM_THREADS=4`")
                println("---")
            end

            Threads.@threads for f in settings.eml_files
                send_messages(settings, Vector{String}([f]), true)
            end
        else
            send_messages(settings, settings.eml_files, false)
        end
    end

    function main()
        if isempty(ARGS)
            print_usage()
            exit(0)
        end

        if ARGS[1] === "--version"
            print_version()
            exit(0)
        end

        for json_file in ARGS
            try
                proc_json_file(json_file)
            catch e
                msg = isa(e, ErrorException) ? e.msg : e
                println("error: $json_file: $msg")
            end
        end
    end

    Base.@ccallable function julia_main()::Cint
        try
            main()
        catch
            Base.invokelatest(Base.display_error, Base.catch_stack())
            return 1
        end
        return 0
    end

    if abspath(PROGRAM_FILE) === @__FILE__
        main()
    end
end # module

max_stars_repo_name: ykon/sendeml-jl

max_issues_repo_name: ykon/sendeml-jl

max_forks_repo_name: ykon/sendeml-jl


text:
value(x) = x
cuify(x) = error("To use LinSolveGPUFactorize, you must do `using CuArrays`")
promote_u0(u0,p,t0) = u0
promote_tspan(u0,p,tspan,prob,kwargs) = tspan
get_tmp(x) = nothing
isdistribution(u0) = false

function SciMLBase.tmap(args...)
  error("Zygote must be added to differentiate Zygote? If you see this error, report it.")
end

if VERSION < v"1.4.0-DEV.635"
  # Piracy, should get upstreamed
  LinearAlgebra.ldiv!(Y::AbstractArray, A::AbstractArray, B::AbstractArray) = (copyto!(Y,B); ldiv!(A,Y))
  LinearAlgebra.ldiv!(Y::AbstractArray, A::AbstractMatrix, B::AbstractArray) = (copyto!(Y,B); ldiv!(A,Y))
  LinearAlgebra.ldiv!(x::AbstractArray,A::Diagonal,b::AbstractArray) = (x .= A.diag .\ b)
end

function __init__()
  @require ApproxFun="28f2ccd6-bb30-5033-b560-165f7b14dc2f" begin
    eval_u0(u0::ApproxFun.Fun) = false
  end

  @require Distributions="31c24e10-a181-5473-b8eb-7969acd0382f" begin
    handle_distribution_u0(_u0::Distributions.Sampleable) = rand(_u0)
    isdistribution(_u0::Distributions.Sampleable) = true
  end

  @require ForwardDiff="f6369f11-7733-5829-9624-2563aa707210" begin

    promote_u0(u0::AbstractArray{<:ForwardDiff.Dual},p::AbstractArray{<:ForwardDiff.Dual},t0) = u0
    promote_u0(u0,p::AbstractArray{<:ForwardDiff.Dual},t0) = eltype(p).(u0)
    promote_u0(u0,p::NTuple{N,<:ForwardDiff.Dual},t0) where N = eltype(p).(u0)
    promote_u0(u0,p::ForwardDiff.Dual,t0) where N = eltype(p).(u0)

    function promote_tspan(u0::AbstractArray{<:ForwardDiff.Dual},p,tspan::Tuple{<:ForwardDiff.Dual,<:ForwardDiff.Dual},prob,kwargs)
      return tspan
    end

    function promote_tspan(u0::AbstractArray{<:ForwardDiff.Dual},p,tspan,prob,kwargs)
      if (haskey(kwargs,:callback) && has_continuous_callback(kwargs[:callback])) ||
         (haskey(prob.kwargs,:callback) && has_continuous_callback(prob.kwargs[:callback]))

        return eltype(u0).(tspan)
      else
        return tspan
      end
    end

    value(x::Type{ForwardDiff.Dual{T,V,N}}) where {T,V,N} = V
    value(x::ForwardDiff.Dual) = value(ForwardDiff.value(x))

    @inline fastpow(x::ForwardDiff.Dual, y::ForwardDiff.Dual) = x^y

    sse(x::Number) = x^2
    sse(x::ForwardDiff.Dual) = sse(ForwardDiff.value(x)) + sum(sse, ForwardDiff.partials(x))
    totallength(x::Number) = 1
    totallength(x::ForwardDiff.Dual) = totallength(ForwardDiff.value(x)) + sum(totallength, ForwardDiff.partials(x))
    totallength(x::AbstractArray) = sum(totallength,x)

    @inline ODE_DEFAULT_NORM(u::ForwardDiff.Dual,::Any) = sqrt(sse(u))
    @inline ODE_DEFAULT_NORM(u::AbstractArray{<:ForwardDiff.Dual},t::Any) = sqrt(sum(sse,u) / totallength(u))
    @inline ODE_DEFAULT_NORM(u::ForwardDiff.Dual,::ForwardDiff.Dual) = sqrt(sse(u))
    @inline ODE_DEFAULT_NORM(u::AbstractArray{<:ForwardDiff.Dual},::ForwardDiff.Dual) = sqrt(sum(x->sse(x),u) / totallength(u))

    if !hasmethod(nextfloat, Tuple{ForwardDiff.Dual})
      # Type piracy. Should upstream
      Base.nextfloat(d::ForwardDiff.Dual{T,V,N}) where {T,V,N} = ForwardDiff.Dual{T}(nextfloat(d.value), d.partials)
      Base.prevfloat(d::ForwardDiff.Dual{T,V,N}) where {T,V,N} = ForwardDiff.Dual{T}(prevfloat(d.value), d.partials)
    end

    struct DiffCache{T<:AbstractArray, S<:AbstractArray}
        du::T
        dual_du::S
    end

    function DiffCache(u::AbstractArray{T}, siz, ::Type{Val{chunk_size}}) where {T, chunk_size}
        x = ArrayInterface.restructure(u,zeros(ForwardDiff.Dual{nothing,T,chunk_size}, siz...))
        DiffCache(u, x)
    end

    dualcache(u::AbstractArray, N=Val{ForwardDiff.pickchunksize(length(u))}) = DiffCache(u, size(u), N)

    function get_tmp(dc::DiffCache, u::T) where T<:ForwardDiff.Dual
      x = reinterpret(T, dc.dual_du)
    end

    function get_tmp(dc::DiffCache, u::AbstractArray{T}) where T<:ForwardDiff.Dual
      x = reinterpret(T, dc.dual_du)
    end

    function DiffEqBase.get_tmp(dc::DiffEqBase.DiffCache, u::LabelledArrays.LArray{T,N,D,Syms}) where {T,N,D,Syms}
      x = reinterpret(T, dc.dual_du.__x)
      LabelledArrays.LArray{T,N,D,Syms}(x)
    end

    get_tmp(dc::DiffCache, u::Number) = dc.du
    get_tmp(dc::DiffCache, u::AbstractArray) = dc.du

    # bisection(f, tup::Tuple{T,T}, t_forward::Bool) where {T<:ForwardDiff.Dual} = find_zero(f, tup, Roots.AlefeldPotraShi())
  end

  @require Measurements="eff96d63-e80a-5855-80a2-b1b0885c5ab7" begin

    promote_u0(u0::AbstractArray{<:Measurements.Measurement},p::AbstractArray{<:Measurements.Measurement},t0) = u0
    promote_u0(u0,p::AbstractArray{<:Measurements.Measurement},t0) = eltype(p).(u0)

    value(x::Type{Measurements.Measurement{T}}) where {T} = T
    value(x::Measurements.Measurement) = Measurements.value(x)

    @inline fastpow(x::Measurements.Measurement, y::Measurements.Measurement) = x^y

    # Support adaptive steps should be errorless
    @inline function ODE_DEFAULT_NORM(u::AbstractArray{<:Measurements.Measurement,N},t) where {N}
      sqrt(sum(x->ODE_DEFAULT_NORM(x[1],x[2]),zip((value(x) for x in u),Iterators.repeated(t))) / length(u))
    end
    @inline function ODE_DEFAULT_NORM(u::Array{<:Measurements.Measurement,N},t) where {N}
      sqrt(sum(x->ODE_DEFAULT_NORM(x[1],x[2]),zip((value(x) for x in u),Iterators.repeated(t))) / length(u))
    end
    @inline ODE_DEFAULT_NORM(u::Measurements.Measurement,t) = abs(Measurements.value(u))
  end

  @require MonteCarloMeasurements="0987c9cc-fe09-11e8-30f0-b96dd679fdca" begin

    promote_u0(u0::AbstractArray{<:MonteCarloMeasurements.AbstractParticles},p::AbstractArray{<:MonteCarloMeasurements.AbstractParticles},t0) = u0
    promote_u0(u0,p::AbstractArray{<:MonteCarloMeasurements.AbstractParticles},t0) = eltype(p).(u0)

    value(x::Type{MonteCarloMeasurements.AbstractParticles{T,N}}) where {T,N} = T
    value(x::MonteCarloMeasurements.AbstractParticles) = mean(x)

    @inline fastpow(x::MonteCarloMeasurements.AbstractParticles, y::MonteCarloMeasurements.AbstractParticles) = x^y

    # Support adaptive steps should be errorless
    @inline function ODE_DEFAULT_NORM(u::AbstractArray{<:MonteCarloMeasurements.AbstractParticles,N},t) where {N}
      sqrt(mean(x->ODE_DEFAULT_NORM(x[1],x[2]),zip((value(x) for x in u),Iterators.repeated(t))))
    end
    @inline function ODE_DEFAULT_NORM(u::AbstractArray{<:MonteCarloMeasurements.AbstractParticles,N},t::AbstractArray{<:MonteCarloMeasurements.AbstractParticles,N}) where {N}
      sqrt(mean(x->ODE_DEFAULT_NORM(x[1],x[2]),zip((value(x) for x in u),Iterators.repeated(value.(t)))))
    end
    @inline ODE_DEFAULT_NORM(u::MonteCarloMeasurements.AbstractParticles,t) = abs(value(u))
  end

  @require Unitful="1986cc42-f94f-5a68-af5c-568840ba703d" begin
    # Support adaptive errors should be errorless for exponentiation
    value(x::Type{Unitful.AbstractQuantity{T,D,U}}) where {T,D,U} = T
    value(x::Unitful.AbstractQuantity) = x.val
    @inline function ODE_DEFAULT_NORM(u::AbstractArray{<:Unitful.AbstractQuantity,N},t) where {N}
      sqrt(sum(x->ODE_DEFAULT_NORM(x[1],x[2]),zip((value(x) for x in u),Iterators.repeated(t))) / length(u))
    end
    @inline function ODE_DEFAULT_NORM(u::Array{<:Unitful.AbstractQuantity,N},t) where {N}
      sqrt(sum(x->ODE_DEFAULT_NORM(x[1],x[2]),zip((value(x) for x in u),Iterators.repeated(t))) / length(u))
    end
    @inline ODE_DEFAULT_NORM(u::Unitful.AbstractQuantity,t) = abs(value(u))
    @inline UNITLESS_ABS2(x::Unitful.AbstractQuantity) = real(abs2(x)/oneunit(x)*oneunit(x))
  end

  @require Tracker="9f7883ad-71c0-57eb-9f7f-b5c9e6d3789c" begin
    include("tracker.jl")
  end

  # Piracy, should get upstreamed
  @require CuArrays="3a865a2d-5b23-5a0f-bc46-62713ec82fae" begin
    cuify(x::AbstractArray) = CuArrays.CuArray(x)
    function LinearAlgebra.ldiv!(x::CuArrays.CuArray,_qr::CuArrays.CUSOLVER.CuQR,b::CuArrays.CuArray)
      _x = UpperTriangular(_qr.R) \ (_qr.Q' * reshape(b,length(b),1))
      x .= vec(_x)
      CuArrays.unsafe_free!(_x)
      return x
    end
    # make `\` work
    LinearAlgebra.ldiv!(F::CuArrays.CUSOLVER.CuQR, b::CuArrays.CuArray) = (x = similar(b); ldiv!(x, F, b); x)
    default_factorize(A::CuArrays.CuArray) = qr(A)
    function findall_events(affect!,affect_neg!,prev_sign::CuArrays.CuArray,next_sign::CuArrays.CuArray)
      hasaffect::Bool = affect! !== nothing
      hasaffectneg::Bool = affect_neg! !== nothing
      f = (p,n)-> ((p < 0 && hasaffect) || (p > 0 && hasaffectneg)) && p*n<=0
      A = map(f,prev_sign,next_sign)
      out = findall(A)
      CuArrays.unsafe_free!(A)
      out
    end

    ODE_DEFAULT_NORM(u::CuArrays.CuArray,t) = sqrt(real(sum(abs2,u))/length(u))

    @require ForwardDiff="f6369f11-7733-5829-9624-2563aa707210" begin
      @inline function ODE_DEFAULT_NORM(u::CuArrays.CuArray{<:ForwardDiff.Dual},t)
        sqrt(sum(abs2,value.(u)) / length(u))
      end

      @inline function ODE_DEFAULT_NORM(u::CuArrays.CuArray{<:ForwardDiff.Dual},t::ForwardDiff.Dual)
        sqrt(sum(abs2,u) / length(u))
      end
    end
  end

  @require CUDA="052768ef-5323-5732-b1bb-66c8b64840ba" begin
    cuify(x::AbstractArray) = CUDA.CuArray(x)
    function LinearAlgebra.ldiv!(x::CUDA.CuArray,_qr::CUDA.CUSOLVER.CuQR,b::CUDA.CuArray)
      _x = UpperTriangular(_qr.R) \ (_qr.Q' * reshape(b,length(b),1))
      x .= vec(_x)
      CUDA.unsafe_free!(_x)
      return x
    end
    # make `\` work
    LinearAlgebra.ldiv!(F::CUDA.CUSOLVER.CuQR, b::CUDA.CuArray) = (x = similar(b); ldiv!(x, F, b); x)
    default_factorize(A::CUDA.CuArray) = qr(A)
    function findall_events(affect!,affect_neg!,prev_sign::CUDA.CuArray,next_sign::CUDA.CuArray)
      hasaffect::Bool = affect! !== nothing
      hasaffectneg::Bool = affect_neg! !== nothing
      f = (p,n)-> ((p < 0 && hasaffect) || (p > 0 && hasaffectneg)) && p*n<=0
      A = map(f,prev_sign,next_sign)
      out = findall(A)
      CUDA.unsafe_free!(A)
      out
    end

    ODE_DEFAULT_NORM(u::CUDA.CuArray,t) = sqrt(real(sum(abs2,u))/length(u))

    @require ForwardDiff="f6369f11-7733-5829-9624-2563aa707210" begin
      @inline function ODE_DEFAULT_NORM(u::CUDA.CuArray{<:ForwardDiff.Dual},t)
        sqrt(sum(abs2,value.(u)) / length(u))
      end

      @inline function ODE_DEFAULT_NORM(u::CUDA.CuArray{<:ForwardDiff.Dual},t::ForwardDiff.Dual)
        sqrt(sum(abs2,u) / length(u))
      end
    end
  end

  @require ReverseDiff="37e2e3b7-166d-5795-8a7a-e32c996b4267" begin
    include("reversediff.jl")
  end

  @require Zygote="e88e6eb3-aa80-5325-afca-941959d7151f" begin
    function ∇tmap(cx, f, args...)
      ys_and_backs = SciMLBase.tmap((args...) -> Zygote._pullback(cx, f, args...), args...)
      if isempty(ys_and_backs)
        ys_and_backs, _ -> nothing
      else
        ys, backs = Zygote.unzip(ys_and_backs)
        function ∇tmap_internal(Δ)
          Δf_and_args_zipped = SciMLBase.tmap((f, δ) -> f(δ), backs, Δ)
          Δf_and_args = Zygote.unzip(Δf_and_args_zipped)
          Δf = reduce(Zygote.accum, Δf_and_args[1])
          (Δf, Δf_and_args[2:end]...)
        end
        ys,∇tmap_internal
      end
    end

    function ∇responsible_map(cx, f, args...)
      ys_and_backs = SciMLBase.responsible_map((args...) -> Zygote._pullback(cx, f, args...), args...)
      if isempty(ys_and_backs)
        ys_and_backs, _ -> nothing
      else
        ys, backs = Zygote.unzip(ys_and_backs)
        ys, function ∇responsible_map_internal(Δ)
          # Apply pullbacks in reverse order. Needed for correctness if `f` is stateful.
          Δf_and_args_zipped = SciMLBase.responsible_map((f, δ) -> f(δ), Zygote._tryreverse(SciMLBase.responsible_map, backs, Δ)...)
          Δf_and_args = Zygote.unzip(Zygote._tryreverse(SciMLBase.responsible_map, Δf_and_args_zipped))
          Δf = reduce(Zygote.accum, Δf_and_args[1])
          (Δf, Δf_and_args[2:end]...)
        end
      end
    end
  end

  @require GeneralizedGenerated="6b9d7cbe-bcb9-11e9-073f-15a7a543e2eb" begin
    SciMLBase.numargs(::GeneralizedGenerated.RuntimeFn{Args}) where Args = GeneralizedGenerated.from_type(Args) |> length
  end

  @require Pardiso="46dd5b70-b6fb-5a00-ae2d-e8fea33afaf2" begin
    mutable struct MKLPardisoFactorize
      A
      ps::Pardiso.MKLPardisoSolver
      verbose::Bool
      firsttime::Bool
    end
    MKLPardisoFactorize(;verbose=false) = MKLPardisoFactorize(nothing,Pardiso.MKLPardisoSolver(),verbose,true)
    function (p::MKLPardisoFactorize)(x,A,b,update_matrix=false;kwargs...)
      if p.firsttime
        Pardiso.set_phase!(p.ps, Pardiso.ANALYSIS)
        Pardiso.pardiso(p.ps, x, A, b)
        p.firsttime = false
      end

      if update_matrix
        Pardiso.set_phase!(p.ps, Pardiso.NUM_FACT)
        Pardiso.pardiso(p.ps, x, A, b)
        p.A = A
      end

      Pardiso.set_phase!(p.ps, Pardiso.SOLVE_ITERATIVE_REFINE)
      Pardiso.pardiso(p.ps, x, A, b)
    end
    function (p::MKLPardisoFactorize)(::Type{Val{:init}},f,u0_prototype)
      if eltype(u0_prototype) <: Complex
        mattype = Pardiso.COMPLEX_NONSYM
      else
        mattype = Pardiso.REAL_NONSYM
      end
      Pardiso.set_matrixtype!(p.ps, mattype)
      if p.verbose
          Pardiso.set_msglvl!(p.ps, Pardiso.MESSAGE_LEVEL_ON)
      end
      MKLPardisoFactorize(nothing,p.ps,p.verbose,true)
    end

    mutable struct PardisoFactorize
      A
      ps::Pardiso.PardisoSolver
      verbose::Bool
      firsttime::Bool
    end
    PardisoFactorize(;verbose=false) = PardisoFactorize(nothing,Pardiso.PardisoSolver(),verbose,true)
    function (p::PardisoFactorize)(x,A,b,update_matrix=false;kwargs...)
      if p.firsttime
        Pardiso.set_phase!(p.ps, Pardiso.ANALYSIS)
        Pardiso.pardiso(p.ps, x, A, b)
        p.firsttime = false
      end

      if update_matrix
        Pardiso.set_phase!(p.ps, Pardiso.NUM_FACT)
        Pardiso.pardiso(p.ps, x, A, b)
        p.A = A
      end

      Pardiso.set_phase!(p.ps, Pardiso.SOLVE_ITERATIVE_REFINE)
      Pardiso.pardiso(p.ps, x, A, b)
    end
    function (p::PardisoFactorize)(::Type{Val{:init}},f,u0_prototype)
      if eltype(u0_prototype) <: Complex
        mattype = Pardiso.COMPLEX_NONSYM
      else
        mattype = Pardiso.REAL_NONSYM
      end
      Pardiso.set_matrixtype!(p.ps, mattype)
      if p.verbose
          Pardiso.set_msglvl!(p.ps, Pardiso.MESSAGE_LEVEL_ON)
      end
      PardisoFactorize(nothing,p.ps,p.verbose,true)
    end

    mutable struct PardisoIterate
      A
      ps::Pardiso.PardisoSolver
      verbose::Bool
      firsttime::Bool
    end
    PardisoIterate(;verbose=false) = PardisoIterate(nothing,Pardiso.PardisoSolver(),verbose,true)
    function (p::PardisoIterate)(x,A,b,update_matrix=false;kwargs...)
      if p.firsttime
        Pardiso.set_phase!(p.ps, Pardiso.ANALYSIS)
        Pardiso.pardiso(p.ps, x, A, b)
        p.firsttime = false
      end

      #=
      # Pardiso iterative solver doesn't need a separate factorization phase
      # It will auto-update as it needs to.
      if update_matrix
        Pardiso.set_phase!(p.ps, Pardiso.NUM_FACT)
        Pardiso.pardiso(p.ps, x, A, b)
        p.A = A
      end
      =#

      Pardiso.set_phase!(ps, Pardiso.NUM_FACT_SOLVE_REFINE)
      Pardiso.pardiso(ps, X, A, B)
    end
    function (p::PardisoIterate)(::Type{Val{:init}},f,u0_prototype)
      if eltype(u0_prototype) <: Complex
        mattype = Pardiso.COMPLEX_NONSYM
      else
        mattype = Pardiso.REAL_NONSYM
      end
      Pardiso.set_matrixtype!(p.ps, mattype)
      Pardiso.set_solver!(ps, Pardiso.ITERATIVE_SOLVER)
      if p.verbose
          Pardiso.set_msglvl!(p.ps, Pardiso.MESSAGE_LEVEL_ON)
      end
      PardisoIterate(nothing,p.ps,p.verbose,true)
    end
  end
end

ext: jl


text:
module AndorInstaller

include("tools.jl")
target = "deps.jl"
rm(target; force=true)
AndorInstallTools.make_deps(target)

end # module


text:
using Pumas, LinearAlgebra, Plots

irm1 = @model begin
    @param begin
        θ ∈ VectorDomain(6)
        Ω ∈  PDiagDomain(5)
    end

    @random begin
        η ~ MvNormal(Ω)
    end

    @pre begin
        Ka      = θ[1]
        CL      = θ[2]*exp(η[1])
        Vc      = θ[3]*exp(η[2])
        Kin     = θ[4]*exp(η[3])
        Kout    = θ[5]*exp(η[4])
        IC50    = θ[6]*exp(η[5])
        IMAX    = 1
    end

    @init begin
        Resp = Kin/Kout
    end

    @vars begin
        cp = Cent/Vc
        inhibition = 1-(IMAX*cp/(IC50+cp))
    end

    @dynamics begin
        Gut'    = -Ka*Gut
        Cent'   =  Ka*Gut - (CL/Vc)*Cent
        Resp'   =  Kin*inhibition - Kout*Resp
    end

    @derived begin
        cp     = Cent / Vc
        resp   = Resp
    end
end

param = (θ = [
          0.5, # Ka  Absorption rate constant 1 (1/time)
          5, # CL   Clearance (volume/time)
          30, # Vc   Central volume (volume)
          10, # Kin  Response in rate constant (1/time)
          0.02, # Kout Response out rate constant (1/time)
          10, # IC50 Concentration for 50% of max inhibition (mass/volume)
          1, # IMAX Maximum inhibition
          ],
    Ω = Diagonal([0.04,0.04,0.04,0.04,0.04]));

regimen1 = DosageRegimen(1500, time=0,cmt=1)
subject1 = Subject(id=1,evs=regimen1)

sim = simobs(irm1,subject1,param,obstimes=0:0.1:100)
plot(sim, obsnames=[:cp,:resp])

irm2 = @model begin
    @param begin
        θ ∈ VectorDomain(6)
        Ω ∈  PDiagDomain(5)
    end

    @random begin
        η ~ MvNormal(Ω)
    end

    @pre begin
        Ka      = θ[1]
        CL      = θ[2]*exp(η[1])
        Vc      = θ[3]*exp(η[2])
        Kin     = θ[4]*exp(η[3])
        Kout    = θ[5]*exp(η[4])
        IC50    = θ[6]*exp(η[5])
        IMAX    = 1
    end

    @init begin
        Resp = Kin/Kout
    end

    @vars begin
        cp = Cent/Vc
        inhibition = 1-(IMAX*cp/(IC50+cp))
    end

    @dynamics begin
        Gut'    = -Ka*Gut
        Cent'   =  Ka*Gut - (CL/Vc)*Cent
        Resp'   =  Kin - Kout*inhibition*Resp
    end

    @derived begin
        cp     = Cent / Vc
        resp   = Resp
    end
end

sim = simobs(irm2,subject1,param,obstimes=0:0.1:100)
plot(sim, obsnames=[:cp,:resp])

irm3 = @model begin
    @param begin
        θ ∈ VectorDomain(6)
        Ω ∈  PDiagDomain(5)
    end

    @random begin
        η ~ MvNormal(Ω)
    end

    @pre begin
        Ka      = θ[1]
        CL      = θ[2]*exp(η[1])
        Vc      = θ[3]*exp(η[2])
        Kin     = θ[4]*exp(η[3])
        Kout    = θ[5]*exp(η[4])
        EC50    = θ[6]*exp(η[5])
        EMAX    = 1
    end

    @init begin
        Resp = Kin/Kout
    end

    @vars begin
        cp = Cent/Vc
        stimulation = 1 + (EMAX*cp/(EC50+cp))
    end

    @dynamics begin
        Gut'    = -Ka*Gut
        Cent'   =  Ka*Gut - (CL/Vc)*Cent
        Resp'   =  Kin*stimulation - Kout*Resp
    end

    @derived begin
        cp     = Cent / Vc
        resp   = Resp
    end
end

sim = simobs(irm3,subject1,param,obstimes=0:0.1:100)
plot(sim, obsnames=[:cp,:resp])

irm4 = @model begin
    @param begin
        θ ∈ VectorDomain(6)
        Ω ∈  PDiagDomain(5)
    end

    @random begin
        η ~ MvNormal(Ω)
    end

    @pre begin
        Ka      = θ[1]
        CL      = θ[2]*exp(η[1])
        Vc      = θ[3]*exp(η[2])
        Kin     = θ[4]*exp(η[3])
        Kout    = θ[5]*exp(η[4])
        EC50    = θ[6]*exp(η[5])
        EMAX    = 1
    end

    @init begin
        Resp = Kin/Kout
    end

    @vars begin
        cp = Cent/Vc
        stimulation = 1 + (EMAX*cp/(EC50+cp))
    end

    @dynamics begin
        Gut'    = -Ka*Gut
        Cent'   =  Ka*Gut - (CL/Vc)*Cent
        Resp'   =  Kin - Kout*stimulation*Resp
    end

    @derived begin
        cp     = Cent / Vc
        resp   = Resp
    end
end

sim = simobs(irm4,subject1,param,obstimes=0:0.1:100)
plot(sim, obsnames=[:cp,:resp])

using PumasTutorials
PumasTutorials.tutorial_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])


text:
struct NotExtensionOnOperation <:Rule
	enter::Function
	leave::Function
	function NotExtensionOnOperation()
		function enter(node::Node)
			if (node.kind== "TypeExtensionDefinition")
				throw(GraphQLError("{\"data\": null,\"errors\": [{\"message\": \"GraphQL cannot execute a request containing a TypeExtensionDefinition.\"}]}"))
			end
		end
		function leave(node::Node)

		end
		new(enter,leave)
	end
end

struct NotTypeOnOperation <:Rule
	enter::Function
	leave::Function
	function NotTypeOnOperation()
		function enter(node::Node)
			if (node.kind== "ObjectTypeDefinition")
				return throw(GraphQLError("{\"data\": null,\"errors\": [{\"message\": \"GraphQL cannot execute a request containing a ObjectTypeDefinition.\"}]}"))
			end
		end
		function leave(node::Node)

		end
		new(enter,leave)
	end
end

struct NotSchemaOnOperation <:Rule
	enter::Function
	leave::Function
	function NotSchemaOnOperation()
		function enter(node::Node)
			if (node.kind== "SchemaDefinition")
				return throw(GraphQLError("{\"data\": null,\"errors\": [{\"message\": \"GraphQL cannot execute a request containing a SchemaDefinition.\"}]}"))
			end
		end
		function leave(node::Node)

		end
		new(enter,leave)
	end
end

struct FragmentSubscription <:Rule
	enter::Function
	leave::Function
	function FragmentSubscription()
		function enter(node::Node)
			if (node.kind== "FragmentDefinition")
				if (node.typeCondition[2].name.value == "Subscription")
					if (length(node.selectionSet.selections)>1)
						return throw(GraphQLError("{\"data\": null,\"errors\": [{\"message\": \"subscription must select only one top level field\"}]}"))
					end
				end
			end
		end
		function leave(node::Node)

		end
		new(enter,leave)
	end
end

struct FragmentNames <:Rule
	enter::Function
	leave::Function
	function FragmentNames()
		nombres=[]
		function enter(node::Node)
			if (node.kind== "FragmentDefinition")
				valor= node.name.value
				if(valor in nombres)
					return throw(GraphQLError("{\"data\": null,\"errors\": [{\"message\": \"There can only be one fragment named \'$(valor)\'.\"}]}"))
				else
					push!(nombres,valor)
				end
			end
		end
		function leave(node::Node)

		end
		new(enter,leave)
	end
end

struct OperationNames <:Rule
	enter::Function
	leave::Function
	function OperationNames()
		nombres=[]
		function enter(node::Node)
			if (node.kind== "OperationDefinition")
				valor= node.name
				if (typeof(valor)<: Name)
					valor= valor.value
					if(valor in nombres)
						return throw(GraphQLError("{\"data\": null,\"errors\": [{\"message\": \"There can only be one operation named \'$(valor)\'.\"}]}"))
					else
						push!(nombres,valor)
					end
				end
			end
		end
		function leave(node::Node)

		end
		new(enter,leave)
	end
end

struct OperationAnonymous <:Rule
	enter::Function
	leave::Function
	function OperationAnonymous()
		n_operation = 0
		anonimo= false
		function enter(node::Node)
			if (node.kind== "OperationDefinition")
				valor = node.name
				n_operation=n_operation+1
				if !(typeof(valor)<: Name)
					anonimo=true
				end
				if ((anonimo== true) & (n_operation>1))
					return throw(GraphQLError("{\"data\": null,\"errors\": [{\"message\": \"This anonymous operation must be the only defined operation.\"}]}"))
				end
			end
		end
		function leave(node::Node)

		end
		new(enter,leave)
	end
end

struct SubscriptionFields <:Rule
	enter::Function
	leave::Function
	function SubscriptionFields()
		function enter(node::Node)
			if (node.kind== "OperationDefinition")
				if (node.operation == "subscription")
					if (length(node.selectionSet.selections)>1)
						return throw(GraphQLError("{\"data\": null,\"errors\": [{\"message\": \"subscription must select only one top level field.\"}]}"))
					end
				end
			end
		end
		function leave(node::Node)

		end
		new(enter,leave)
	end
end

struct FragmentUnknowNotUsed <:Rule
	enter::Function
	leave::Function
	function FragmentUnknowNotUsed()
		nombres=[]
		usados=[]
		yausados=[]
		texto = ""
		function enter(node::Node)
			if (node.kind== "FragmentSpread") #usados
				nombre = node.name.value
				if( nombre in nombres)
					deleteat!(nombres,findall(isequal(nombre), nombres))
				else
					if !(nombre in yausados)
					push!(usados,nombre)
					end
				end
			end
			if (node.kind== "FragmentDefinition")
				valor= node.name.value
				if(valor in usados)
					deleteat!(usados,findall(isequal(valor), usados))
					if(!(valor in yausados))
					    push!(yausados,valor)
					end
				else
					push!(nombres,valor)
				end
			end
		end
		function leave(node::Node)
			if (node.kind=="Document")
				if (length(nombres)>0)
				  for n in nombres
				      texto=texto*" "*n
				  end
				  return throw(GraphQLError("{\"data\": null,\"errors\": [{\"message\": \"Fragment $(texto) is not used.\"}]}"))
				end
				if (length(usados)>0)
				    for n in usados
				      texto=texto*" "*n
				    end
 				  return throw(GraphQLError("{\"data\": null,\"errors\": [{\"message\": \"Unknown fragment $(texto).\"}]}"))
				end
			end
		end
		new(enter,leave)
	end
end

struct FragmentCycles <:Rule
	enter::Function
	leave::Function
	function FragmentCycles()
		nombrescycles=[]
		usadoscycles=[]
		inicio = ""
		traza = ""
		leerspread =false
		function enter(node::Node)
			if (node.kind== "FragmentSpread") #usadoscycles
				if (leerspread == true)
				nombre = node.name.value
				push!(usadoscycles,nombre)
				end
			end

			if (node.kind== "FragmentDefinition")
				leerspread =true
				valor= node.name.value
				push!(nombrescycles,valor)
			end
		end

		function leave(node::Node)
			if (node.kind== "FragmentDefinition")
				leerspread =false
			end
			if (node.kind=="Document")
				if (length(nombrescycles)>0)
				if (length(nombrescycles)==length(usadoscycles))
					inicio = nombrescycles[1]
					for l in nombrescycles[2:end]
					    traza=traza*" "*l
					end
				    return throw(GraphQLError("{\"data\": null,\"errors\": [{\"message\": \"Cannot spread fragment $(inicio) within itself via $(traza).\"}]}"))
				end
				end
			end
		end
		new(enter,leave)
	end
end


text:
# Ref https://arxiv.org/pdf/1901.10739.pdf
# This term does not contain the kinetic energy, which must be added separately
# /!\ They have no 1/2 factor in front of the kinetic energy,
#     so for consistency the added kinetic energy must have a scaling of 2
# Energy = <u, ((-i hbar ∇ + βA)^2 + V) u>
# where ∇∧A = 2π ρ, ∇⋅A = 0 => A = x^⟂/|x|² ∗ ρ
# H = (-i hbar ∇ + βA)² + V - 2β x^⟂/|x|² ∗ (βAρ + J)
#   = - hbar^2 Δ + 2β hbar (-i∇)⋅A + β²|A|^2 - 2β x^⟂/|x|² ∗ (βAρ + J)
# where only the first three terms "count for the energy", and where
# J = 1/(2i) hbar (u* ∇u - u ∇u*)

# for numerical reasons, we solve the equation ∇∧A = 2π ρ in two parts: A = A_SR + A_LR, with
# ∇∧A_SR = 2π(ρ-ρref)
# ∇∧A_LR = 2πρref
# where ρref is a gaussian centered at the origin and with the same integral as ρ (=M)
# This way, ρ-ρref has zero mass, and A_SR is shorter range.
# The long-range part is computed explicitly

function ρref_real(x::T, y::T, M, σ) where {T <: Real}
    r = hypot(x, y)
    # gaussian normalized to have integral M
    M * exp(-T(1)/2 * (r/σ)^2) / (σ^2 * (2T(π))^(2/2))
end

function magnetic_field_produced_by_ρref(x::T, y::T, M, σ) where {T <: Real}
    # The solution of ∇∧A = 2π ρref is ϕ(r) [-y;x] where ϕ satisfies
    # ∇∧A = 2ϕ + r ϕ' => rϕ' + 2 ϕ = 2π ρref
    # Wolfram alpha (after a bit of coaxing) says the solution to
    # rϕ' + 2ϕ = C exp(-α x^2)
    # is
    # 1/x^2 (c1 - C exp(-α x^2)/2α), which coupled with smoothness at 0 gives
    # C/(2α)/x^2*(1 - exp(-α x^2))
    r = hypot(x, y)
    r == 0 && return magnetic_field_produced_by_ρref(1e-8, 0.0, M, σ) # hack

    α = 1 / (2*σ^2)
    C = 2T(π) * M / (σ^2 * (2T(π))^(2/2))
    r = hypot(x, y)
    ϕ = 1/2 * C / α / r^2 * (1 - exp(-α*r^2))
    ϕ * @SVector[-y, x]
end

function make_div_free(basis::PlaneWaveBasis{T}, A) where {T}
    out = [zeros(complex(T), basis.fft_size), zeros(complex(T), basis.fft_size)]
    A_fourier = [r_to_G(basis, A[α]) for α = 1:2]
    for (iG, G) in enumerate(G_vectors(basis))
        vec = [A_fourier[1][iG], A_fourier[2][iG]]
        G = [G[1], G[2]]
        # project on divergence free fields, ie in Fourier
        # project A(G) on the orthogonal of G
        if iG != 1
            out[1][iG], out[2][iG] = vec - (G'vec) * G / sum(abs2, G)
        else
            out[1][iG], out[2][iG] = vec
        end
    end
    [G_to_r(basis, out[α]) for α = 1:2]
end

struct Anyonic
    hbar
    β
end
function (A::Anyonic)(basis)
    @assert length(basis.kpoints) == 1
    @assert basis.kpoints[1].coordinate == [0, 0, 0]
    @assert basis.model.n_dim == 2
    # only square lattices allowed
    # (because I can't be bothered to think about the right formulas otherwise,
    # although they might already be correct)
    @assert basis.model.lattice[2, 1] == basis.model.lattice[1, 2] == 0
    @assert basis.model.lattice[1, 1] == basis.model.lattice[2, 2]
    @assert basis.model.n_spin_components == 1

    TermAnyonic(basis, eltype(basis)(A.hbar), eltype(basis)(A.β))
end

struct TermAnyonic{T <: Real, Tρ, TA} <: Term
    basis::PlaneWaveBasis{T}
    hbar::T
    β::T
    ρref::Tρ
    Aref::TA
end
function TermAnyonic(basis::PlaneWaveBasis{T}, hbar, β) where {T}
    # compute correction magnetic field
    ρref = zeros(T, basis.fft_size)
    Aref = [zeros(T, basis.fft_size), zeros(T, basis.fft_size)]
    M = basis.model.n_electrons  # put M to zero to disable
    σ = 2
    for (ir, r) in enumerate(r_vectors(basis))
        rcart = basis.model.lattice * (r - @SVector[.5, .5, .0])
        ρref[ir] = ρref_real(rcart[1], rcart[2], M, σ)
        Aref[1][ir], Aref[2][ir] = magnetic_field_produced_by_ρref(rcart[1], rcart[2], M, σ)
    end
    # Aref is not divergence-free in the finite basis, so we explicitly project it
    # This is because we assume elsewhere (eg to ensure self-adjointness of the Hamiltonian)
    Aref = make_div_free(basis, Aref)
    TermAnyonic(basis, hbar, β, ρref, Aref)
end

function ene_ops(term::TermAnyonic, ψ, occ; ρ, kwargs...)
    basis = term.basis
    T = eltype(basis)
    hbar = term.hbar
    β = term.β
    @assert ψ !== nothing # the hamiltonian depends explicitly on ψ

    # Compute A in Fourier domain
    # curl A = 2π ρ, ∇⋅A = 0
    # i [G1;G2] ∧ A(G1,G2) = 2π ρ(G1,G2), [G1;G2] ⋅ A(G1,G2) = 0
    # => A(G1, G2) = -2π i ρ(G1, G2) * [-G2;G1;0]/(G1^2 + G2^2)
    A1 = zeros(complex(T), basis.fft_size)
    A2 = zeros(complex(T), basis.fft_size)
    ρ_fourier = r_to_G(basis, ρ[:, :, :, 1])
    ρref_fourier = r_to_G(basis, term.ρref) # TODO optimize
    for (iG, G) in enumerate(G_vectors_cart(basis))
        G2 = sum(abs2, G)
        if G2 != 0
            A1[iG] = +2T(π) * G[2] / G2 * (ρ_fourier[iG] - ρref_fourier[iG]) * im
            A2[iG] = -2T(π) * G[1] / G2 * (ρ_fourier[iG] - ρref_fourier[iG]) * im
        end
    end
    Areal = [G_to_r(basis, A1) + term.Aref[1],
             G_to_r(basis, A2) + term.Aref[2],
             zeros(T, basis.fft_size)]

    # 2 hbar β (-i∇)⋅A + β^2 |A|^2
    ops_energy = [MagneticFieldOperator(basis, basis.kpoints[1],
                                        2 .* hbar .* β .* Areal),
                  RealSpaceMultiplication(basis, basis.kpoints[1],
                                          β^2 .* (abs2.(Areal[1]) .+ abs2.(Areal[2])))]

    # Now compute effective local potential - 2β x^⟂/|x|² ∗ (βAρ + J)
    J = compute_current(basis, ψ, occ)
    eff_current = [hbar .* J[α] .+
                   β .* ρ .* Areal[α] for α = 1:2]
    eff_current_fourier = [r_to_G(basis, eff_current[α]) for α = 1:2]
    # eff_pot = - 2β x^⟂/|x|² ∗ eff_current
    # => ∇∧eff_pot = -4πβ eff_current
    # => eff_pot(G1, G2) = 4πβ i eff_current(G1, G2) * [-G2;G1;0]/(G1^2 + G2^2)
    eff_pot_fourier = zeros(complex(T), basis.fft_size)
    for (iG, Gred) in enumerate(G_vectors(basis))
        G = basis.model.recip_lattice * Gred
        G2 = sum(abs2, G)
        if G2 != 0
            eff_pot_fourier[iG] += -4T(π)*β * im * G[2] / G2 * eff_current_fourier[1][iG]
            eff_pot_fourier[iG] += +4T(π)*β * im * G[1] / G2 * eff_current_fourier[2][iG]
        end
    end
    eff_pot_real = G_to_r(basis, eff_pot_fourier)
    ops_ham = [ops_energy..., RealSpaceMultiplication(basis, basis.kpoints[1], eff_pot_real)]

    E = zero(T)
    for iband = 1:size(ψ[1], 2)
        ψnk = @views ψ[1][:, iband]
        # TODO optimize this
        for op in ops_energy
            E += occ[1][iband] * real(dot(ψnk, op * ψnk))
        end
    end

    (E=E, ops=[ops_ham])
end


text:
using ProgressMeter

import DataFrames: DataFrame, colwise

export tapply_mean
function tapply_mean(subs, val, sz=(maximum(subs),))
    A = zeros(eltype(val), sz...)
    counter = zeros(Int64, sz...)
    for i = 1:length(val)
        @inbounds A[subs[i]] += val[i]
        @inbounds counter[subs[i]] += 1
    end
    A./counter
end

export update_tune_df
function update_tune_df(tune_df::Float64, acpt_rt::Float64, opt_rt::Float64, gamma1::Float64)
    max(exp(log(tune_df) - gamma1 * (acpt_rt - opt_rt)), 50)
end

export extend_mcmc
function extend_mcmc(dat::DataFrame,
                        chain::Array,
                        tune_df_chain::Array,
                        hyp::Tuple,
                        alpha::Array{Float64,1},
                        nstep::Int64,
                        labels::Array{String,1};
                        opt_rate::Float64 = 0.3)
    """
    Take old chain as input and add new iterations
    """
    nw = size(chain, 1)
    n, nd = size(dat)
    nm = size(chain[1,1][2], 1)
    labels = hcat([(parse.(Int64, split(x, "")) .-1) for x in labels])
    curLen = size(chain, 2)
    tune_df = tune_df_chain[:,curLen-1]

    # Initial estimates at NIW, prop, and z for each walker and each temperature
    param = Array{Tuple{Array{Dict{String,Array{Float64,N} where N},1},Array{Float64,1},Array{Int64,1}}}(undef, (nw, 1))
    for i in 1:nw
        dictionary = Dict{String,Array{Float64,N} where N}[]
        for m in 1:nm
                Sigma = chain[i,curLen][1][m]["Sigma"]
                mu = chain[i,curLen][1][m]["mu"]
                push!(dictionary, Dict("mu" => mu, "Sigma" => Sigma))
        end
        initprop = chain[i,curLen][2]
        z = chain[i,curLen][3]
        param[i,1] = (dictionary, initprop, z)
    end

    chain = Array{
                    Tuple{Array{Dict{String,Array{Float64,N} where N},1},
                    Array{Float64,1},
                    Array{Int64,1}}
                }(undef, nw, nstep + 1)

    for j in 1:nw
        chain[j,1] = deepcopy(param[j,1])
    end
    acpt_tracker = zeros(nm)

    # For tracking the acceptance rate, per cluster, for the adaptive tuning variance
    win_len = min(nstep, 50)
    acpt_win = zeros(nm, win_len)
    acpt_chain = zeros(nm, nstep)

    old_tune_df = deepcopy(tune_df)
    tune_df_chain = zeros(nm, nstep)

    @showprogress 1 "Running the MCMC..."  for i in 1:nstep
        # Proposes IW cluster at a time, accepts/rejects/returns new IW
        NIW, acpt = make_mcmc_move(dat, chain[:,i], hyp, alpha, labels, tune_df)
        acpt_tracker = acpt_tracker .+ acpt
        @inbounds acpt_win[:, (i-1) % win_len + 1] = acpt

        # Makes Gibbs updates of means, mixing weights, and class labels
        newz, newNIW, newprop = make_gibbs_update(dat, chain[:,i], hyp, alpha, labels)

        if i > 50
            # Update tuning parameter per cluster
            gamma1 = 10 / (i ^ 0.8)
            old_tune_df = copy(tune_df)
            @inbounds [tune_df[m] = update_tune_df(tune_df[m], mean(acpt_win[m,:]), opt_rate, gamma1) for m in 1:1:nm]
        end

        for j in 1:nw
            chain_link = deepcopy(param)

            @inbounds map!(x -> x, chain_link[j,1][1], newNIW[j,1])
            @inbounds map!(x -> x, chain_link[j,1][2], newprop[j,1])
            @inbounds map!(x -> x, chain_link[j,1][3], newz[j,1])

            @inbounds chain[j,i+1] = deepcopy(chain_link[j,1])
            @inbounds acpt_chain[:,i] = acpt_tracker ./ i
            @inbounds tune_df_chain[:,i] = tune_df
        end
    end
    (chain, acpt_chain, tune_df_chain)
end


text:
using ProbNumODE
using DifferentialEquations
using Plots
using LinearAlgebra
using LaTeXStrings

FILEDIR = dirname(@__FILE__)

function makeplot(sol, n=10)
    errs = sol.u - appxsol.(sol.t)
    p = plot(
        xlabel=L"$t$", ylabel=L"$(\hat{y}(t) - y^*(t))$",
    );
    stds = sqrt.(ProbNumODE.stack(diag.(sol.pu.Σ)))
    plot!(sol.t, [zero(sol.t) zero(sol.t)], ribbon=3stds,
          color=[1 2], linealpha=0, fillalpha=0.2, label="");

    if n > 0
        sp = ProbNumODE.sample(sol, n)
        sp_errs = ProbNumODE.stack(sol.u) .- sp
        for i in 1:size(sp)[3]
            plot!(p, sol.t, sp_errs[:, :, i], color=[1 2], label="", linewidth=0.5)
        end
        plot!(p, sol.t, ProbNumODE.stack(errs), color=:black, label="");
    else
        plot!(p, sol.t, ProbNumODE.stack(errs), color=[1 2], label="");
    end

    return p
end

prob = ProbNumODE.remake_prob_with_jac(fitzhugh_nagumo_iip())
appxsol = solve(remake(prob, u0=big.(prob.u0)), abstol=1e-30, reltol=1e-30);

sol = solve(prob, EKF1(), q=3, sigmarule=:dynamic, smooth=true, abstol=1e-10, reltol=1e-7)
p = makeplot(sol, 0);
savefig(joinpath(FILEDIR, "sigmas-tv-eks1_fitzhughnagumo.png"))

sol = solve(prob, EKF0(), q=3, sigmarule=:dynamicMV, smooth=true, abstol=1e-10, reltol=1e-7)
p = makeplot(sol, 0);
savefig(joinpath(FILEDIR, "sigmas-tvmv-eks0_fitzhughnagumo.png"))

sol = solve(prob, EKF0(), q=3, sigmarule=:dynamic, smooth=true, abstol=1e-10, reltol=1e-7)
p = makeplot(sol, 0);
savefig(joinpath(FILEDIR, "sigmas-tv-eks0_fitzhughnagumo_nosamples.png"))

println("Plots saved to $FILEDIR")

max_stars_repo_name: nathanaelbosch/capos

max_issues_repo_name: nathanaelbosch/capos

max_forks_repo_name: nathanaelbosch/capos


text:
accessType = Union{Colon,<:Integer,<:AbstractArray{<:Union{Integer,Bool}},Nothing}
"""
    struct pathLocs{m}
        indices
    end
    pathLocs(varargs...; m::Int=2, d::Int=1, exs=Colon())
some words
    pathLocs(s::scattered)
make a pathLocs that accesses every location in the output
    pathLocs(ii, s::scattered)
make a pathLocs that accesses every location in layer `ii` of the output.
"""
struct pathLocs{m}
    indices
end
# specifying on
function pathLocs(varargs...; m::Int=2, d::Int=1, exs=Colon())
    if length(varargs) == 0
        # default constructor that includes everything
        paired = [(i, :) for i = 0:m]
    elseif length(varargs) % 2 == 0 && (varargs[1] isa Integer)
        # specifically enumerated, so alternates between a number and a reference shape
        # pair each layer with the tuple
        paired = [(varargs[i], varargs[i+1]) for i = 1:2:length(varargs)]
    elseif length(varargs) == m + 1
        # every layer specified
        paired = [(i, varargs[i+1]) for i = 0:m]
    end
    present = [x[1] for x in paired]
    notPresent = [ii for ii = 0:m if !(ii in present)] # unspecified layers
    paired = [paired..., [(x, nothing) for x in notPresent]...]
    paired = sort(paired, by=(x) -> x[1])
    indices = map(x -> parseOne(x, d, exs), paired)
    fullList = map(ii -> paired[ii+1], 0:m)
    return pathLocs{m + 1}((indices...,))
end

function parseOne(x, d, exs)
    lay, select = x
    if select == nothing
        return nothing
    end
    if lay == 0
        if select == Colon()
            return (map(x -> Colon(), 1:d)..., Colon(), exs)
        elseif length(select) == d + 2
            return select
        elseif typeof(select) <: AbstractArray && ndims(select) == d + 2
            return select
        elseif typeof(select) <: AbstractArray
            return (select, 1, exs)
        else
            return (select..., 1, exs)
        end
    elseif lay == 1
        if select == Colon()
            return (map(x -> Colon(), 1:d)..., Colon(), exs)
        elseif length(select) == d + 2
            return select
        elseif typeof(select) <: AbstractArray && ndims(select) == d + 2
            return select
        elseif typeof(select) <: AbstractArray || typeof(select) <: Number
            return (map(x -> Colon(), 1:d)..., select, exs)
        elseif length(select) == 1
            return (map(x -> Colon(), 1:d)..., select[1], exs)
        else
            return (select..., exs)
        end
    elseif lay >= 2
        if select == Colon()
            return (map(x -> Colon(), 1:d)..., Colon(), Colon(), exs)
        elseif (select isa Tuple) && length(select) >= d + 3
            return select
        elseif typeof(select) <: AbstractArray && ndims(select) == d + 3
            return select
        elseif typeof(select) <: AbstractArray # only specify the first layer
            return (map(x -> Colon(), 1:d+1)..., select, exs)
        elseif length(select) == 2 # specify both
            return (map(x -> Colon(), 1:d)..., select..., exs)
        else
            return (select..., 1, exs)
        end
    end
end

import Base: size
function size(p::pathLocs)
    ł(x::Nothing) = Tuple{Vararg{<:Integer,0}}()
    ł(x::AbstractArray) = size(x)
    # all that's left is a tuple of things
    ł(x) = map(a -> length(a), x)
    map(ł, p.indices)
end


text:
@inline function get_interpolation(wpot::Interpolations.Extrapolation{T, 3}, pt::CartesianPoint{T}, ::Type{Cartesian})::T where {T <: SSDFloat}
    return wpot(pt.x, pt.y, pt.z)::T
end
@inline function get_interpolation(wpot::Interpolations.Extrapolation{T, 3}, pt::CartesianPoint{T}, ::Type{Cylindrical})::T where {T <: SSDFloat}
    pt_cyl::CylindricalPoint{T} = CylindricalPoint(pt)
    return wpot(pt_cyl.r, pt_cyl.φ, pt_cyl.z)::T
end
@inline function get_interpolation(wpot::Interpolations.Extrapolation{T, 3}, pt::CylindricalPoint{T}, ::Type{Cartesian})::T where {T <: SSDFloat}
    pt_car::CartesianPoint{T} = CartesianPoint(pt)
    return wpot(pt_car.x, pt_car.y, pt_car.z)::T
end
@inline function get_interpolation(wpot::Interpolations.Extrapolation{T, 3}, pt::CylindricalPoint{T}, ::Type{Cylindrical})::T where {T <: SSDFloat}
    return wpot(pt.r, pt.φ, pt.z)::T
end

function add_signal!(signal::AbstractVector{T}, timestamps::AbstractVector{T}, path::Vector{CartesianPoint{T}}, pathtimestamps::AbstractVector{T}, charge::T,
                        wpot::Interpolations.Extrapolation{T, 3}, S::CoordinateSystemType)::Nothing where {T <: SSDFloat}
    tmp_signal = Vector{T}(undef, length(pathtimestamps))
    @inbounds for i in eachindex(tmp_signal)
        tmp_signal[i] = get_interpolation(wpot, path[i], S)::T * charge
    end
    itp = interpolate( (pathtimestamps,), tmp_signal, Gridded(Linear()))
    t_max::T = last(pathtimestamps)
    i_max::Int = searchsortedlast(timestamps, t_max)
    signal[1:i_max] += itp(timestamps[1:i_max])
    if length(signal) > i_max
        signal[i_max+1:end] .+= tmp_signal[end]
    end
    nothing
end

function add_signal!(signal::AbstractVector{T}, timestamps::AbstractVector{TT}, path::EHDriftPath{T, TT}, charge::T, wpot::Interpolations.Extrapolation{T, 3}, S::CoordinateSystemType)::Nothing where {T <: SSDFloat, TT}
    add_signal!(signal, timestamps, path.e_path, path.timestamps_e, -charge, wpot, S) # electrons induce negative charge
    add_signal!(signal, timestamps, path.h_path, path.timestamps_h,  charge, wpot, S)
    nothing
end

function add_signal!(signal::AbstractVector{T}, timestamps::AbstractVector{<:RealQuantity}, paths::Vector{<:EHDriftPath{T}}, charges::Vector{T}, wpot::Interpolations.Extrapolation{T, 3}, S::CoordinateSystemType)::Nothing where {T <: SSDFloat}
    for ipath in eachindex(paths)
        add_signal!(signal, timestamps, paths[ipath], charges[ipath], wpot, S)
    end
    nothing
end


text:
#=

 Unit tests from libieeep1788 for interval numeric operations
 (Original author: Marco Nehmeier)
 converted into portable ITL format by Oliver Heimlich.

 Copyright 2013-2015 Marco Nehmeier (nehmeier@informatik.uni-wuerzburg.de)
 Copyright 2015-2017 Oliver Heimlich (oheim@posteo.de)

 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
 You may obtain a copy of the License at

     http://www.apache.org/licenses/LICENSE-2.0

 Unless required by applicable law or agreed to in writing, software
 distributed under the License is distributed on an "AS IS" BASIS,
 WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 See the License for the specific language governing permissions and
 limitations under the License.

=#
#Language imports

#Test library imports
using Test

#Arithmetic library imports
using IntervalArithmetic

#Preamble
setprecision(53)
setprecision(Interval, Float64)
setrounding(Interval, :tight)
# Set full format, and show decorations
@format full

@testset "minimal_inf_test" begin
    @test inf(∅) == Inf
    @test inf(Interval(-Inf, Inf)) == -Inf
    @test inf(Interval(1.0, 2.0)) == 1.0
    @test inf(Interval(-3.0, -2.0)) == -3.0
    @test inf(Interval(-Inf, 2.0)) == -Inf
    @test inf(Interval(-Inf, 0.0)) == -Inf
    @test inf(Interval(-Inf, -0.0)) == -Inf
    @test inf(Interval(-2.0, Inf)) == -2.0
    @test inf(Interval(0.0, Inf)) == -0.0
    @test inf(Interval(-0.0, Inf)) == -0.0
    @test inf(Interval(-0.0, 0.0)) == -0.0
    @test inf(Interval(0.0, -0.0)) == -0.0
    @test inf(Interval(0.0, 0.0)) == -0.0
    @test inf(Interval(-0.0, -0.0)) == -0.0
end

@testset "minimal_inf_dec_test" begin
    @test isnan(inf(nai()))
    @test inf(DecoratedInterval(∅, trv)) == Inf
    @test inf(DecoratedInterval(Interval(-Inf, Inf), def)) == -Inf
    @test inf(DecoratedInterval(Interval(1.0, 2.0), com)) == 1.0
    @test inf(DecoratedInterval(Interval(-3.0, -2.0), trv)) == -3.0
    @test inf(DecoratedInterval(Interval(-Inf, 2.0), dac)) == -Inf
    @test inf(DecoratedInterval(Interval(-Inf, 0.0), def)) == -Inf
    @test inf(DecoratedInterval(Interval(-Inf, -0.0), trv)) == -Inf
    @test inf(DecoratedInterval(Interval(-2.0, Inf), trv)) == -2.0
    @test inf(DecoratedInterval(Interval(0.0, Inf), def)) == -0.0
    @test inf(DecoratedInterval(Interval(-0.0, Inf), trv)) == -0.0
    @test inf(DecoratedInterval(Interval(-0.0, 0.0), dac)) == -0.0
    @test inf(DecoratedInterval(Interval(0.0, -0.0), trv)) == -0.0
    @test inf(DecoratedInterval(Interval(0.0, 0.0), trv)) == -0.0
    @test inf(DecoratedInterval(Interval(-0.0, -0.0), trv)) == -0.0
end

@testset "minimal_sup_test" begin
    @test sup(∅) == -Inf
    @test sup(Interval(-Inf, Inf)) == Inf
    @test sup(Interval(1.0, 2.0)) == 2.0
    @test sup(Interval(-3.0, -2.0)) == -2.0
    @test sup(Interval(-Inf, 2.0)) == 2.0
    @test sup(Interval(-Inf, 0.0)) == 0.0
    @test sup(Interval(-Inf, -0.0)) == 0.0
    @test sup(Interval(-2.0, Inf)) == Inf
    @test sup(Interval(0.0, Inf)) == Inf
    @test sup(Interval(-0.0, Inf)) == Inf
    @test sup(Interval(-0.0, 0.0)) == 0.0
    @test sup(Interval(0.0, -0.0)) == 0.0
    @test sup(Interval(0.0, 0.0)) == 0.0
    @test sup(Interval(-0.0, -0.0)) == 0.0
end

@testset "minimal_sup_dec_test" begin
    @test isnan(sup(nai()))
    @test sup(DecoratedInterval(∅, trv)) == -Inf
    @test sup(DecoratedInterval(Interval(-Inf, Inf), def)) == Inf
    @test sup(DecoratedInterval(Interval(1.0, 2.0), com)) == 2.0
    @test sup(DecoratedInterval(Interval(-3.0, -2.0), trv)) == -2.0
    @test sup(DecoratedInterval(Interval(-Inf, 2.0), dac)) == 2.0
    @test sup(DecoratedInterval(Interval(-Inf, 0.0), def)) == 0.0
    @test sup(DecoratedInterval(Interval(-Inf, -0.0), trv)) == 0.0
    @test sup(DecoratedInterval(Interval(-2.0, Inf), trv)) == Inf
    @test sup(DecoratedInterval(Interval(0.0, Inf), def)) == Inf
    @test sup(DecoratedInterval(Interval(-0.0, Inf), trv)) == Inf
    @test sup(DecoratedInterval(Interval(-0.0, 0.0), dac)) == +0.0
    @test sup(DecoratedInterval(Interval(0.0, -0.0), trv)) == +0.0
    @test sup(DecoratedInterval(Interval(0.0, 0.0), trv)) == +0.0
    @test sup(DecoratedInterval(Interval(-0.0, -0.0), trv)) == +0.0
end

@testset "minimal_mid_test" begin
    @test isnan(mid(∅))
    @test mid(Interval(-Inf, Inf)) == 0.0
    @test mid(Interval(-0x1.fffffffffffffP1023, +0x1.fffffffffffffP1023)) == 0.0
    @test mid(Interval(0.0, 2.0)) == 1.0
    @test mid(Interval(2.0, 2.0)) == 2.0
    @test mid(Interval(-2.0, 2.0)) == 0.0
    @test mid(Interval(0.0, Inf)) == 0x1.fffffffffffffP1023
    @test mid(Interval(-Inf, 1.2)) == -0x1.fffffffffffffP1023
    @test mid(Interval(-0x0.0000000000002p-1022, 0x0.0000000000001p-1022)) == 0.0
    @test mid(Interval(-0x0.0000000000001p-1022, 0x0.0000000000002p-1022)) == 0.0
    @test mid(Interval(0x1.fffffffffffffP+1022, 0x1.fffffffffffffP+1023)) == 0x1.7ffffffffffffp+1023
    @test mid(Interval(0x0.0000000000001p-1022, 0x0.0000000000003p-1022)) == 0x0.0000000000002p-1022
end

@testset "minimal_mid_dec_test" begin
    @test isnan(mid(DecoratedInterval(∅, trv)))
    @test isnan(mid(nai()))
    @test mid(DecoratedInterval(Interval(-Inf, Inf), def)) == 0.0
    @test mid(DecoratedInterval(Interval(-0x1.fffffffffffffP1023, +0x1.fffffffffffffP1023), trv)) == 0.0
    @test mid(DecoratedInterval(Interval(0.0, 2.0), com)) == 1.0
    @test mid(DecoratedInterval(Interval(2.0, 2.0), dac)) == 2.0
    @test mid(DecoratedInterval(Interval(-2.0, 2.0), trv)) == 0.0
    @test mid(DecoratedInterval(Interval(0.0, Inf), trv)) == 0x1.fffffffffffffP1023
    @test mid(DecoratedInterval(Interval(-Inf, 1.2), trv)) == -0x1.fffffffffffffP1023
    @test mid(DecoratedInterval(Interval(-0x0.0000000000002p-1022, 0x0.0000000000001p-1022), trv)) == 0.0
    @test mid(DecoratedInterval(Interval(-0x0.0000000000001p-1022, 0x0.0000000000002p-1022), trv)) == 0.0
    @test mid(DecoratedInterval(Interval(0x1.fffffffffffffP+1022, 0x1.fffffffffffffP+1023), trv)) == 0x1.7ffffffffffffp+1023
    @test mid(DecoratedInterval(Interval(0x0.0000000000001p-1022, 0x0.0000000000003p-1022), trv)) == 0x0.0000000000002p-1022
end

@testset "minimal_rad_test" begin
    @test radius(Interval(0.0, 2.0)) == 1.0
    @test radius(Interval(2.0, 2.0)) == 0.0
    @test isnan(radius(∅))
    @test radius(Interval(-Inf, Inf)) == Inf
    @test radius(Interval(0.0, Inf)) == Inf
    @test radius(Interval(-Inf, 1.2)) == Inf
    @test radius(Interval(-0x0.0000000000002p-1022, 0x0.0000000000001p-1022)) == 0x0.0000000000002p-1022
    @test radius(Interval(0x0.0000000000001p-1022, 0x0.0000000000002p-1022)) == 0x0.0000000000001p-1022
    @test radius(Interval(0x1p+0, 0x1.0000000000003p+0)) == 0x1p-51
end

@testset "minimal_rad_dec_test" begin
    @test radius(DecoratedInterval(Interval(0.0, 2.0), trv)) == 1.0
    @test radius(DecoratedInterval(Interval(2.0, 2.0), com)) == 0.0
    @test isnan(radius(DecoratedInterval(∅, trv)))
    @test isnan(radius(nai()))
    @test radius(DecoratedInterval(Interval(-Inf, Inf), trv)) == Inf
    @test radius(DecoratedInterval(Interval(0.0, Inf), def)) == Inf
    @test radius(DecoratedInterval(Interval(-Inf, 1.2), trv)) == Inf
    @test radius(DecoratedInterval(Interval(-0x0.0000000000002p-1022, 0x0.0000000000001p-1022), trv)) == 0x0.0000000000002p-1022
    @test radius(DecoratedInterval(Interval(0x0.0000000000001p-1022, 0x0.0000000000002p-1022), trv)) == 0x0.0000000000001p-1022
    @test radius(DecoratedInterval(Interval(0x1p+0, 0x1.0000000000003p+0), trv)) == 0x1p-51
end

@testset "minimal_mid_rad_test" begin
    @test midpoint_radius(∅) === (NaN, NaN)
    @test midpoint_radius(Interval(-Inf, Inf)) == (0.0, Inf)
    @test midpoint_radius(Interval(-0x1.fffffffffffffP1023, +0x1.fffffffffffffP1023)) == (0.0, 0x1.fffffffffffffP+1023)
    @test midpoint_radius(Interval(0.0, 2.0)) == (1.0, 1.0)
    @test midpoint_radius(Interval(2.0, 2.0)) == (2.0, 0.0)
    @test midpoint_radius(Interval(-2.0, 2.0)) == (0.0, 2.0)
    @test midpoint_radius(Interval(0, Inf)) == (0x1.fffffffffffffP+1023, Inf)
    @test midpoint_radius(Interval(-Inf, 1.2)) == (-0x1.fffffffffffffP+1023, Inf)
    @test midpoint_radius(Interval(-0x0.0000000000002p-1022, 0x0.0000000000001p-1022)) == (0.0, 0x0.0000000000002p-1022)
    @test midpoint_radius(Interval(0x1.fffffffffffffP+1022, 0x1.fffffffffffffP+1023)) == (0x1.7FFFFFFFFFFFFP+1023, 0x1.0p+1022)
    @test midpoint_radius(Interval(0x0.0000000000001p-1022, 0x0.0000000000003p-1022)) == (0x0.0000000000002P-1022, 0x0.0000000000001P-1022)
end

@testset "minimal_mid_rad_dec_test" begin
    @test midpoint_radius(DecoratedInterval(∅, trv)) === (NaN, NaN)
    @test midpoint_radius(DecoratedInterval(Interval(-Inf, Inf), dac)) == (0.0, Inf)
    @test midpoint_radius(DecoratedInterval(Interval(-0x1.fffffffffffffP+1023, 0x1.fffffffffffffP+1023), trv)) == (0.0, 0x1.fffffffffffffP+1023)
    @test midpoint_radius(DecoratedInterval(Interval(0.0, 2.0), com)) == (1.0, 1.0)
    @test midpoint_radius(DecoratedInterval(Interval(2.0, 2.0), dac)) == (2.0, 0.0)
    @test midpoint_radius(DecoratedInterval(Interval(-2.0, 2.0), trv)) == (0.0, 2.0)
    @test midpoint_radius(DecoratedInterval(Interval(0, Inf), trv)) == (0x1.fffffffffffffP+1023, Inf)
    @test midpoint_radius(DecoratedInterval(Interval(-Inf, 1.2), trv)) == (-0x1.fffffffffffffP+1023, Inf)
    @test midpoint_radius(DecoratedInterval(Interval(-0x0.0000000000002P-1022, 0x0.0000000000001P-1022), trv)) == (0.0, 0x0.0000000000002P-1022)
    @test midpoint_radius(DecoratedInterval(Interval(-0x0.0000000000001P-1022, 0x0.0000000000002P-1022), trv)) == (0.0, 0x0.0000000000002P-1022)
    @test midpoint_radius(DecoratedInterval(Interval(0x1.fffffffffffffP+1022, 0x1.fffffffffffffP+1023), trv)) == (0x1.7FFFFFFFFFFFFP+1023, 0x1.0p+1022)
    @test midpoint_radius(DecoratedInterval(Interval(0x0.0000000000001P-1022, 0x0.0000000000003P-1022), trv)) == (0x0.0000000000002P-1022, 0x0.0000000000001P-1022)
end

@testset "minimal_wid_test" begin
    @test diam(Interval(2.0, 2.0)) == 0.0
    @test diam(Interval(1.0, 2.0)) == 1.0
    @test diam(Interval(1.0, Inf)) == Inf
    @test diam(Interval(-Inf, 2.0)) == Inf
    @test diam(Interval(-Inf, Inf)) == Inf
    @test isnan(diam(∅))
    @test diam(Interval(0x1p+0, 0x1.0000000000001p+0)) == 0x1p-52
    @test diam(Interval(0x1p-1022, 0x1.0000000000001p-1022)) == 0x0.0000000000001p-1022
end

@testset "minimal_wid_dec_test" begin
    @test diam(DecoratedInterval(Interval(2.0, 2.0), com)) == 0.0
    @test diam(DecoratedInterval(Interval(1.0, 2.0), trv)) == 1.0
    @test diam(DecoratedInterval(Interval(1.0, Inf), trv)) == Inf
    @test diam(DecoratedInterval(Interval(-Inf, 2.0), def)) == Inf
    @test diam(DecoratedInterval(Interval(-Inf, Inf), trv)) == Inf
    @test isnan(diam(DecoratedInterval(∅, trv)))
    @test isnan(diam(nai()))
    @test diam(DecoratedInterval(Interval(0x1p+0, 0x1.0000000000001p+0), trv)) == 0x1p-52
    @test diam(DecoratedInterval(Interval(0x1p-1022, 0x1.0000000000001p-1022), trv)) == 0x0.0000000000001p-1022
end

@testset "minimal_mag_test" begin
    @test mag(Interval(1.0, 2.0)) == 2.0
    @test mag(Interval(-4.0, 2.0)) == 4.0
    @test mag(Interval(-Inf, 2.0)) == Inf
    @test mag(Interval(1.0, Inf)) == Inf
    @test mag(Interval(-Inf, Inf)) == Inf
    @test isnan(mag(∅))
    @test mag(Interval(-0.0, 0.0)) == 0.0
    @test mag(Interval(-0.0, -0.0)) == 0.0
end

@testset "minimal_mag_dec_test" begin
    @test mag(DecoratedInterval(Interval(1.0, 2.0), com)) == 2.0
    @test mag(DecoratedInterval(Interval(-4.0, 2.0), trv)) == 4.0
    @test mag(DecoratedInterval(Interval(-Inf, 2.0), trv)) == Inf
    @test mag(DecoratedInterval(Interval(1.0, Inf), def)) == Inf
    @test mag(DecoratedInterval(Interval(-Inf, Inf), trv)) == Inf
    @test isnan(mag(DecoratedInterval(∅, trv)))
    @test isnan(mag(nai()))
    @test mag(DecoratedInterval(Interval(-0.0, 0.0), trv)) == 0.0
    @test mag(DecoratedInterval(Interval(-0.0, -0.0), trv)) == 0.0
end

@testset "minimal_mig_test" begin
    @test mig(Interval(1.0, 2.0)) == 1.0
    @test mig(Interval(-4.0, 2.0)) == 0.0
    @test mig(Interval(-4.0, -2.0)) == 2.0
    @test mig(Interval(-Inf, 2.0)) == 0.0
    @test mig(Interval(-Inf, -2.0)) == 2.0
    @test mig(Interval(-1.0, Inf)) == 0.0
    @test mig(Interval(1.0, Inf)) == 1.0
    @test mig(Interval(-Inf, Inf)) == 0.0
    @test isnan(mig(∅))
    @test mig(Interval(-0.0, 0.0)) == 0.0
    @test mig(Interval(-0.0, -0.0)) == 0.0
end

@testset "minimal_mig_dec_test" begin
    @test mig(DecoratedInterval(Interval(1.0, 2.0), com)) == 1.0
    @test mig(DecoratedInterval(Interval(-4.0, 2.0), trv)) == 0.0
    @test mig(DecoratedInterval(Interval(-4.0, -2.0), trv)) == 2.0
    @test mig(DecoratedInterval(Interval(-Inf, 2.0), def)) == 0.0
    @test mig(DecoratedInterval(Interval(-Inf, -2.0), trv)) == 2.0
    @test mig(DecoratedInterval(Interval(-1.0, Inf), trv)) == 0.0
    @test mig(DecoratedInterval(Interval(1.0, Inf), trv)) == 1.0
    @test mig(DecoratedInterval(Interval(-Inf, Inf), trv)) == 0.0
    @test isnan(mig(DecoratedInterval(∅, trv)))
    @test isnan(mig(nai()))
    @test mig(DecoratedInterval(Interval(-0.0, 0.0), trv)) == 0.0
    @test mig(DecoratedInterval(Interval(-0.0, -0.0), trv)) == 0.0
end
# FactCheck.exitstatus()


text:
using Pkg
if ! ("PackageCompiler" ∈ keys(Pkg.project().dependencies))
    @info "Installing PackageCompiler ..."
    Pkg.add("PackageCompiler")
end
if ! ("Plots" ∈ keys(Pkg.project().dependencies))
    @info "Installing Plots ..."
    Pkg.add("Plots")
end
@info "Loading packages ..."
using KiteUtils, KitePodModels, KiteModels, Plots, AtmosphericModels, Reexport
using PackageCompiler

@info "Creating sysimage ..."
push!(LOAD_PATH,joinpath(pwd(),"src"))

PackageCompiler.create_sysimage(
    [:KiteUtils, :KitePodModels, :KiteModels, :Plots, :AtmosphericModels, :Reexport];
    sysimage_path="kps-image_tmp.so",
    precompile_execution_file=joinpath("test", "test_for_precompile.jl")
)


text:
import Pkg;
Pkg.add("ClimaComms");
import Pkg;
Pkg.add("ClimaCommsMPI");

using ClimaComms
using ClimaCommsMPI

is_distributed = get(ENV, "CLIMACORE_DISTRIBUTED", "") == "MPI"

using Logging
if is_distributed
    using ClimaComms
    if ENV["CLIMACORE_DISTRIBUTED"] == "MPI"
        using ClimaCommsMPI
        const comms_ctx = ClimaCommsMPI.MPICommsContext()
    else
        error("ENV[\"CLIMACORE_DISTRIBUTED\"] only supports the \"MPI\" option")
    end
    const pid, nprocs = ClimaComms.init(comms_ctx)
    logger_stream = ClimaComms.iamroot(comms_ctx) ? stderr : devnull
    prev_logger = global_logger(ConsoleLogger(logger_stream, Logging.Info))
    @info "Setting up distributed run on $nprocs \
        processor$(nprocs == 1 ? "" : "s")"
else
    const comms_ctx = nothing
    using TerminalLoggers: TerminalLogger
    prev_logger = global_logger(TerminalLogger())
end
atexit() do
    global_logger(prev_logger)
end

max_stars_repo_name: CliMA/ClimaCoupler


text:
x = [1.47, 1.50, 1.52, 1.55, 1.57, 1.60, 1.63, 1.65, 1.68, 1.70, 1.73, 1.75, 1.78, 1.80, 1.83]
y = [52.21, 53.12, 54.48, 55.84, 57.20, 58.57, 59.93, 61.29, 63.11, 64.47, 66.28, 68.10, 69.92, 72.19, 74.46]
X = [x.^0 x.^1 x.^2];
b = X \ y

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
#=
day12:
- Julia version: 1.7.0
- Author: Paul.Mealor
- Date: 2021-12-12
=#

using AoC, Test, Memoize

function parseLink!(links, link)
    (start,finish) = split(link, "-")
    if finish != "start"; links[start] = union(get(links,start,[]), [finish]); end
    if start != "start"; links[finish] = union(get(links,finish,[]), [start]); end
    return links
end
@test parseLink!(Dict(), "a-b") == Dict("a"=>["b"], "b"=>["a"])
@test parseLink!(parseLink!(Dict(), "a-b"), "a-c") == Dict("a"=>["b","c"], "b"=>["a"], "c"=>["a"])

function parseLinks(links)
    result = Dict{String, Vector{String}}()
    for link in links
        parseLink!(result, link)
    end
    return result
end

@memoize manyVisitsAllowed(cave) = occursin(r"^[A-Z]+$", cave)

visit(links, cave::T; allowTwoVisits=false) where T = visit!(links, cave, Set{Vector{T}}(), T[]; allowTwoVisits)
function visit!(links, cave, paths, visited; allowTwoVisits=false)
    push!(visited, cave)
    if cave == "end"
        if visited ∉ paths
            push!(paths, copy(visited))
        end
        pop!(visited)
        return paths
    end

    for target in links[cave]
        if manyVisitsAllowed(target)
            paths = visit!(links, target, paths, visited; allowTwoVisits=allowTwoVisits)
        elseif allowTwoVisits
            paths = visit!(links, target, paths, visited; allowTwoVisits=(target ∉ visited))
        elseif target ∉ visited
            paths = visit!(links, target, paths, visited; allowTwoVisits=allowTwoVisits)
        end
    end

    pop!(visited)
    return paths
end
@test visit(Dict("start"=>["end"]), "start") == Set([["start", "end"],])

function part1(lines)
    links = parseLinks(lines)
    return visit(links, "start") |> length
end
@test part1(exampleLines(12,1)) == 10
@test part1(exampleLines(12,2)) == 19
@test part1(exampleLines(12,3)) == 226

function part2(lines)
    links = parseLinks(lines)
    return visit(links, "start"; allowTwoVisits=true) |> length
end
@test part2(exampleLines(12,1)) == 36
@test part2(exampleLines(12,2)) == 103
@test part2(exampleLines(12,3)) == 3509

lines(12) |> ll -> @time part1(ll) |> show
lines(12) |> ll -> @time part2(ll) |> show


text:
using ArnoldiMethod
#computes normalized cut cost for partition `cut`
function _normalized_cut_cost(cut, W::AbstractMatrix, D)
    cut_cost = 0
    for j in axes(W, 2)
        for i in axes(W, 1)
            if cut[i] != cut[j]
                cut_cost += W[i, j]
            end
        end
    end
    cut_cost /= 2
    return cut_cost / sum(D * cut) + cut_cost / sum(D * (.~cut))
end

function _normalized_cut_cost(cut, W::SparseMatrixCSC, D)
    cut_cost = 0
    rows = rowvals(W)
    vals = nonzeros(W)
    n = size(W, 2)
    for i = 1:n
        for j in nzrange(W, i)
            row = rows[j]
            if cut[i] != cut[row]
                cut_cost += vals[j] / 2
            end
        end
    end
    return cut_cost / sum(D * cut) + cut_cost / sum(D * (.~cut))
end

function _partition_weightmx(cut, W::AbstractMatrix)
    nv = length(cut)
    nv2 = sum(cut)
    nv1 = nv - nv2
    newvid = Vector{Int}(undef, nv)
    vmap1 = Vector{Int}(undef, nv1)
    vmap2 = Vector{Int}(undef, nv2)
    j1 = 1
    j2 = 1
    for i in eachindex(cut)
        if cut[i] == false
            newvid[i] = j1
            vmap1[j1] = i
            j1 += 1
        else
            newvid[i] = j2
            vmap2[j2] = i
            j2 += 1
        end
    end

    W1 = similar(W, (nv1, nv1))
    W2 = similar(W, (nv2, nv2))

    for j in axes(W, 2)
        for i in axes(W, 1)
            if cut[i] == cut[j] == false
                W1[newvid[i], newvid[j]] = W[i, j]
            elseif cut[i] == cut[j] == true
                W2[newvid[i], newvid[j]] = W[i, j]
            end
        end
    end

    return (W1, W2, vmap1, vmap2)
end

function _partition_weightmx(cut, W::SparseMatrixCSC)
    nv = length(cut)
    nv2 = sum(cut)
    nv1 = nv - nv2
    newvid = Vector{Int}(undef, nv)
    vmap1 = Vector{Int}(undef, nv1)
    vmap2 = Vector{Int}(undef, nv2)
    j1 = 1
    j2 = 1
    for i in eachindex(cut)
        if cut[i] == false
            newvid[i] = j1
            vmap1[j1] = i
            j1 += 1
        else
            newvid[i] = j2
            vmap2[j2] = i
            j2 += 1
        end
    end

    rows = rowvals(W)
    vals = nonzeros(W)
    I1 = Vector{Int}(); I2 = Vector{Int}()
    J1 = Vector{Int}(); J2 = Vector{Int}()
    V1 = Vector{Float64}(); V2 = Vector{Float64}()
    for i = 1:nv
        for j in nzrange(W, i)
            row = rows[j]
            if cut[i] == cut[row] == false
                push!(I1, newvid[i])
                push!(J1, newvid[row])
                push!(V1, vals[j])
            elseif cut[i] == cut[row] == true
                push!(I2, newvid[i])
                push!(J2, newvid[row])
                push!(V2, vals[j])
            end
        end
    end
    W1 = sparse(I1, J1, V1)
    W2 = sparse(I2, J2, V2)
    return (W1, W2, vmap1, vmap2)
end

function _recursive_normalized_cut(W, thres=thres, num_cuts=num_cuts)
    m, n = size(W)
    D = Diagonal(vec(sum(W, dims=2)))

    m == 1 && return [1]

    #get eigenvector corresponding to second smallest eigenvalue
    # v = eigs(D-W, D, nev=2, which=SR())[2][:,2]
    # At least some versions of ARPACK have a bug, this is a workaround
    invDroot = sqrt.(inv(D)) # equal to Cholesky factorization for diagonal D
    if n > 12
        λ, Q = eigs(invDroot' * (D - W) * invDroot, nev=12, which=SR())
        ret = real(Q[:,2])
    else
        ret = eigen(Matrix(invDroot' * (D - W) * invDroot)).vectors[:,2]
    end
    v = invDroot * ret

    #perform n-cuts with different partitions of v and find best one
    min_cost = Inf
    best_thres = -1
    for t in range(minimum(v), stop=maximum(v), length=num_cuts)
        cut = v .> t
        cost = _normalized_cut_cost(cut, W, D)
        if cost < min_cost
            min_cost = cost
            best_thres = t
        end
    end

    if min_cost < thres
        #split graph, compute normalized_cut for each subgraph recursively and merge indices.
        cut = v .> best_thres
        W1, W2, vmap1, vmap2 = _partition_weightmx(cut, W)
        labels1 = _recursive_normalized_cut(W1, thres, num_cuts)
        labels2 = _recursive_normalized_cut(W2, thres, num_cuts)

        labels = Vector{Int}(undef, m)
        offset = maximum(labels1)

        for i in eachindex(labels1)
            labels[vmap1[i]] = labels1[i]
        end
        for i in eachindex(labels2)
            labels[vmap2[i]] = labels2[i] + offset
        end

        return labels
    else
        return ones(Int, m)
    end
end

"""
    normalized_cut(g, thres, distmx=weights(g), [num_cuts=10]);

Perform [recursive two-way normalized graph-cut](https://en.wikipedia.org/wiki/Segmentation-based_object_categorization#Normalized_cuts)
on a graph, partitioning the vertices into disjoint sets.
Return a vector that contains the set index for each vertex.

It is important to identify a good threshold for your application. A bisection search over the range (0,1) will help determine a good value of thres.

### Keyword Arguments
- `thres`: Subgraphs aren't split if best normalized cut is above this threshold
- `num_cuts`: Number of cuts performed to determine optimal cut

### References
"Normalized Cuts and Image Segmentation" - Jianbo Shi and Jitendra Malik
"""
function normalized_cut(g::AbstractGraph,
    thres::Real,
    W::AbstractMatrix{T}=adjacency_matrix(g),
    num_cuts::Int=10) where T <: Real

    return _recursive_normalized_cut(W, thres, num_cuts)
end


text:
print("base, ")

include("init-namedarrays.jl")

nn = @inferred similar(n)

@test copyto!(nn, 1, n, 1, length(n)) == n
@test copyto!(nn, 1, n.array, 1, length(n)) == n

@test_throws BoundsError copyto!(nn, 1, n.array, 1, length(n)+1)
@test_throws BoundsError copyto!(nn, 0, n.array, 1, length(n))
@test_throws BoundsError copyto!(nn, 2, n.array, 1, length(n))
@test_throws BoundsError copyto!(nn, length(n)+1, n.array, 1, 1)
@test_throws BoundsError copyto!(nn, 1, n.array, 0, length(n))
@test_throws BoundsError copyto!(nn, 1, n.array, length(n)+1, length(n))

nn = @inferred similar(n, Int)
@test namesanddim(nn) == namesanddim(n)
nn = @inferred similar(n, 3, 5)
@test hasdefaultnames(nn)
nn = @inferred similar(n, 2)
@test namesanddim(nn, 1) == namesanddim(n, 1)
nn = @inferred similar(n, 2, 10)
@test namesanddim(nn, 1) == namesanddim(n, 1)
@test hasdefaultnames(nn, 2)
@test hasdefaultdimnames(nn, 2)
nn = @inferred similar(n, 10, 4, 5)
@test namesanddim(nn, 2) == namesanddim(n, 2)
for i in [1,3]
	@test hasdefaultnames(nn, i)
	@test hasdefaultdimnames(nn, i)
end

# `ind2sub(dims, ind)` is deprecated, use `Tuple(CartesianIndices(dims)[ind])` for a
# direct replacement. In many cases, the conversion to `Tuple` is not necessary.
# @test ind2sub(n, 4) == ("two", "b")
# @test_throws BoundsError ind2sub(n, 0)
# @test_throws BoundsError ind2sub(n, 9)

DelimitedFiles.writedlm(stdout, n)
DelimitedFiles.writedlm(stdout, v)

## Issue #60
for func in [similar,
             a -> fill!(similar(a), zero(eltype(a))),
             a -> fill!(similar(a), one(eltype(a))),
             n -> hcat(n, n),
             n -> vcat(n, n)]
    fn = @inferred func(n)
    @test keytype.(fn.dicts) == (String, String)
end

## issue 61
d = Dict(v)
for key in keys(d)
	@test d[key] == v[key]
end


text:
using Cambrian
export MapElites
struct SparseArray{T,N} <: AbstractArray{T,N}
          data::Dict{NTuple{N,Int}, T}
          dims::NTuple{N,Int}
end

SparseArray(::Type{T}, dims::Int...) where {T} = SparseArray(T, dims)
SparseArray(::Type{T}, dims::NTuple{N,Int}) where {T,N} = SparseArray{T,N}(Dict{NTuple{N,Int}, T}(), dims)
Base.size(A::SparseArray) = A.dims
Base.similar(A::SparseArray, ::Type{T}, dims::Dims) where {T} = SparseArray(T, dims)
Base.getindex(A::SparseArray{T,N}, I::Vararg{Int,N}) where {T,N} = get(A.data, I, zero(T))
Base.setindex!(A::SparseArray{T,N}, v, I::Vararg{Int,N}) where {T,N} = (A.data[I] = v)
Base.zero(Ind::Type{Individual}) = missing
Base.zero(elt::Type{Union{Missing, Array{Float64,N} where N}}) = missing

struct MapElites
    feature_dimension::Int64
    grid_mesh::Int64
    solutions::SparseArray{Union{Missing,Individual}}
    performances::SparseArray{Union{Missing,Array{Float64}}}
end

function MapElites(f_dim::Int64,g_size::Int64)
    feature_dimension = f_dim
    grid_mesh = g_size
    size_map = Tuple([g_size for i in 1:f_dim])
    solutions = SparseArray(Union{Missing,Individual},size_map)
    performances = SparseArray(Union{Missing,Array{Float64}},size_map)
    MapElites(feature_dimension,grid_mesh,solutions,performances)
end