text:
# ---
# title: 1577. Number of Ways Where Square of Number Is Equal to Product of Two Numbers
# id: problem1577
# author: Indigo
# date: 2021-01-27
# difficulty: Medium
# categories: Hash Table, Math
# link: <https://leetcode.com/problems/number-of-ways-where-square-of-number-is-equal-to-product-of-two-numbers/description/>
# hidden: true
# ---
#
# Given two arrays of integers `nums1` and `nums2`, return the number of
# triplets formed (type 1 and type 2) under the following rules:
#
#   * Type 1: Triplet (i, j, k) if `nums1[i]2 == nums2[j] * nums2[k]` where `0 <= i < nums1.length` and `0 <= j < k < nums2.length`.
#   * Type 2: Triplet (i, j, k) if `nums2[i]2 == nums1[j] * nums1[k]` where `0 <= i < nums2.length` and `0 <= j < k < nums1.length`.
#
#
#
# **Example 1:**
#
#
#
#     Input: nums1 = [7,4], nums2 = [5,2,8,9]
#     Output: 1
#     Explanation: Type 1: (1,1,2), nums1[1]^2 = nums2[1] * nums2[2]. (4^2 = 2 * 8).
#
#
# **Example 2:**
#
#
#
#     Input: nums1 = [1,1], nums2 = [1,1,1]
#     Output: 9
#     Explanation: All Triplets are valid, because 1^2 = 1 * 1.
#     Type 1: (0,0,1), (0,0,2), (0,1,2), (1,0,1), (1,0,2), (1,1,2).  nums1[i]^2 = nums2[j] * nums2[k].
#     Type 2: (0,0,1), (1,0,1), (2,0,1). nums2[i]^2 = nums1[j] * nums1[k].
#
#
# **Example 3:**
#
#
#
#     Input: nums1 = [7,7,8,3], nums2 = [1,2,9,7]
#     Output: 2
#     Explanation: There are 2 valid triplets.
#     Type 1: (3,0,2).  nums1[3]^2 = nums2[0] * nums2[2].
#     Type 2: (3,0,1).  nums2[3]^2 = nums1[0] * nums1[1].
#
#
# **Example 4:**
#
#
#
#     Input: nums1 = [4,7,9,11,23], nums2 = [3,5,1024,12,18]
#     Output: 0
#     Explanation: There are no valid triplets.
#
#
#
#
# **Constraints:**
#
#   * `1 <= nums1.length, nums2.length <= 1000`
#   * `1 <= nums1[i], nums2[i] <= 10^5`
#
#
## @lc code=start
using LeetCode

function num_triplets(nums1::Vector{Int}, nums2::Vector{Int})
    cnt1 = counter(nums1)
    cnt2 = counter(nums2)
    sqr_cnt1 = Dict(k ^ 2 => v for (k, v) in cnt1)
    sqr_cnt2 = Dict(k ^ 2 => v for (k, v) in cnt2)
    res = 0
    for (k1, v1) in cnt1
        for (k2, v2) in sqr_cnt2
            if k1 ^ 2 == k2
                res += (v1 - 1) * v1 * v2
            else
                res += cnt1[k2 ÷ k1] * v1 * v2
            end
        end
    end
    for (k1, v1) in cnt2
        for (k2, v2) in sqr_cnt1
            if k2 % k1 == 0
                if k1 ^ 2 == k2
                    res += (v1 - 1) * v1 * v2
                else
                    res += cnt2[k2 ÷ k1] * v1 * v2
                end
            end
        end
    end
    res ÷ 2
end
## @lc code=end


text:
"""
    Partition(limits::Matrix{Float64})

Construct an object capable of storing a hyperrectangular partition.
"""
mutable struct Partition
    limits::Matrix{Float64}
    regions::Vector{Matrix{Float64}}

    function Partition(limits::Matrix{Float64})
        check_limits(limits)
        return new(deepcopy(limits), deepcopy([limits]))
    end
end

"""
    split!(P::Partition, k::Int64, d::Int64)

Split the `k`'th subregion of `P` into equal halves in dimension `d`.
"""
function split!(P::Partition, k::Int64, d::Int64)
    if !(1 <= k <= length(P.regions))
        throw(ArgumentError("P does not contain a subregion at index $k"))
    elseif !(1 <= d <= size(P.limits, 1))
        throw(ArgumentError("P is not $d-dimensional"))
    end
    right = deepcopy(P.regions[k])
    loc = sum(right[d, :]) / 2
    P.regions[k][d, 2] = loc
    right[d, 1] = loc
    return push!(P.regions, right)
end

"""
    locate(P::Partition, X::AbstractMatrix)

Return a vector of integers giving the region index for each column of `X`.
"""
function locate(P::Partition, X::AbstractMatrix)
    d, n = size(X)
    if d != size(P.limits, 1)
        throw(DimensionMismatch("P does match the dimension of X"))
    end
    idx = Vector{Int64}(undef, n)
    for i in 1:n
        contained = true
        for k in 1:length(P.regions)
            for l in 1:d
                upper = X[l, i] == P.regions[k][l, 2] == P.limits[l, 2]
                below = X[l, i] >= P.regions[k][l, 1]
                above = X[l, i] < P.regions[k][l, 2]
                contained = below && (above || upper)
                if !contained
                    break
                end
            end
            if contained
                idx[i] = k
                break
            end
        end
        if !contained
            throw(ArgumentError("no region located for observation $i"))
        end
    end
    return idx
end

mutable struct PartitionedBayesPM <: AbstractBayesianLM
    P::Partition
    models::Vector{BayesPM}
    shape0::Float64
    scale0::Float64
    shape::Float64
    scale::Float64

    function PartitionedBayesPM(
        P::Partition, models::Vector{BayesPM}, shape0::Float64, scale0::Float64
    )
        shape = shape0
        scale = scale0
        for pm in models
            shape += pm.lm.shape - shape0
            scale += pm.lm.scale - scale0
        end
        return new(P, models, shape0, scale0, shape, scale)
    end
end

function shape_scale(ppm::PartitionedBayesPM)
    return ppm.shape, ppm.scale
end

"""
    PartitionedBayesPM(P::Partition, Js::Vector{Int64}; <keyword arguments>)

Contruct a partitioned polynomial model.

# Arguments

- `P::Partition`: A partition of the space.
- `λ::Float64=1.0`: Prior scaling.
- `shape0::Float64=1e-3`: Inverse-gamma prior shape hyperparameter.
- `scale0::Float64=1e-3`: Inverse-gamma prior scale hyperparameter.
"""
function PartitionedBayesPM(
    P::Partition,
    Js::Vector{Int64};
    λ::Float64=1.0,
    shape0::Float64=1e-3,
    scale0::Float64=1e-3,
)
    if length(P.regions) != length(Js)
        throw(ArgumentError("must supply a value of J for every region in P"))
    end
    d = size(P.limits, 1)
    basis = [tpbasis(d, J) for J in Js]
    models = [
        BayesPM(basis[i], P.regions[i]; λ=λ, shape0=shape0, scale0=scale0) for
        i in 1:length(Js)
    ]
    return PartitionedBayesPM(P, models, shape0, scale0)
end

function fit!(ppm::PartitionedBayesPM, X::AbstractMatrix, y::AbstractMatrix)
    check_regression_data(X, y)
    idx = locate(ppm.P, X)
    for k in unique(idx)
        region_mask = idx .== k
        ppm.scale -= ppm.models[k].lm.scale - ppm.scale0
        fit!(ppm.models[k], X[:, region_mask], y[:, region_mask])
        ppm.scale += ppm.models[k].lm.scale - ppm.scale0
    end
    return ppm.shape += size(X, 2) / 2
end

"""
    lasso_selection(X::AbstractMatrix, y::AbstractMatrix, Pmax::Int64, intercept::Bool)

Choose the first `Pmax` features introduced by a LASSO solution path.

# Arguments

- `X::AbstractMatrix`: A matrix with observations stored as columns.
- `y::AbstractMatrix`: A matrix with 1 row of response variables.
- `Pmax::Int64`: The maximum number of predictors.
- `intercept::Bool`: `true` if the first row of `X` are the intercept features
"""
function lasso_selection(X::AbstractMatrix, y::AbstractMatrix, Pmax::Int64, intercept::Bool)
    if size(X, 1) <= Pmax
        return [1:size(X, 1)...]
    end
    X, Pmax = intercept ? (X[2:end, :], Pmax - 1) : (X, Pmax)
    β = @suppress glmnet(X', y[1, :]; pmax=Pmax).betas[:, end]
    indices = intercept ? [1] : Int64[]
    @inbounds for j in 1:length(β)
        if β[j] != 0.0
            i = intercept ? j + 1 : j
            push!(indices, i)
        end
    end
    return indices
end

function evidence(models::Vector{BayesPM}, shape0::Float64, scale0::Float64)
    shape, scale, ev = shape0, scale0, 0.0
    for pm in models
        shape += pm.lm.shape - shape0
        scale += pm.lm.scale - scale0
        ev += 0.5 * logdet(pm.lm.Σ) - 0.5 * logdet(pm.lm.Σ0)
    end
    n = 2 * (shape - shape0)
    ev += -(n / 2) * log(2 * π)
    ev += shape0 * log(scale0) - shape * log(scale)
    ev += loggamma(shape) - loggamma(shape0)
    return ev
end

function _conditional_degree_selection!(
    X::AbstractMatrix, # Data to train new polynomial
    y::AbstractMatrix, # Data to train new polynomial
    k::Int64, # Subregion to be replaced
    d::Int64, # Dimension index of model cache
    lateral::Int64, # Lateral index of model cache
    sub_limits::Matrix{Float64}, # Limits of the new polynomial
    Jmax::Int64,
    Pmax::Int64,
    models::Vector{BayesPM},
    model_cache::Vector{Array{BayesPM,3}},
    basis_cache::Vector{Vector{Index}},
    λ::Float64,
    shape0::Float64,
    scale0::Float64,
    ratio::Float64,
)
    n = size(X, 2)
    best_ev = -Inf
    best_pm::Union{Nothing,BayesPM} = nothing
    models_cp = deepcopy(models)
    for J in 0:Jmax
        if J > 0 && n < ratio * length(basis_cache[J + 1])
            continue
        end
        if !isassigned(model_cache[k], d, lateral, J + 1)
            # Polynomial not found in cache, fit a new degree J polynmomial
            # and store it.
            basis = basis_cache[J + 1]
            Z = expand(X, basis, sub_limits; J=J)
            idx = lasso_selection(Z, y, Pmax, true)
            models_cp[k] = BayesPM(
                basis[idx], sub_limits; λ=λ, shape0=shape0, scale0=scale0
            )
            if n > 0
                fit!(models_cp[k].lm, Z[idx, :], y)
            end
            model_cache[k][d, lateral, J + 1] = models_cp[k]
        else
            # Polynomial found in cache.
            models_cp[k] = model_cache[k][d, lateral, J + 1]
        end
        ev = evidence(models_cp, shape0, scale0)
        if ev > best_ev
            best_ev = ev
            best_pm = models_cp[k]
        end
    end
    return best_pm
end

vol(limits::Matrix{Float64}) = prod(limits[:, 2] - limits[:, 1])

"""
    PartitionedBayesPM(X::AbstractMatrix, y::AbstractMatrix, limits::Matrix{Float64};
                       <keyword arguments>)

Perform a 1-step look ahead greedy search for a partitioned polynomial model.

# Keyword Arguments
- `Jmax::Int64=3`: The maximum degree of any polynomial model.
- `Pmax::Int64=500`: The maximum number of features in a particular regions.
- `Kmax::Int64=200`: The maximum number of regions
- `λ::Float64=1.0`: Prior scaling.
- `shape0::Float64=1e-3`: Inverse-gamma prior shape hyperparameter.
- `scale0::Float64=1e-3`: Inverse-gamma prior scale hyperparameter.
- `ratio::Float64=1.0`: Polynomial degrees are reduced until `size(X, 2) < ratio * length(tpbasis(d, J))`.
- `tol::Float64=1e-4`: The required increase in the model evidence to accept a split.
- `verbose::Bool=true`: Print details of the partition search.
"""
function PartitionedBayesPM(
    X::AbstractMatrix,
    y::AbstractMatrix,
    limits::Matrix{Float64};
    Jmax::Int64=3,
    Pmax::Int64=500,
    Kmax::Int64=200,
    λ::Float64=1.0,
    shape0::Float64=1e-3,
    scale0::Float64=1e-3,
    ratio::Float64=1.0,
    tol::Float64=1e-4,
    verbose::Bool=true,
)
    check_regression_data(X, y)
    check_limits(limits)
    if size(X, 1) != size(limits, 1)
        throw(ArgumentError("X and limits don't match in their first dimension"))
    elseif Jmax < 0
        throw(ArgumentError("Jmax must be non-negative"))
    elseif Kmax <= 0
        throw(ArgumentError("Kmax must be strictly positive"))
    elseif tol < 0.0
        throw(ArgumentError("tolerance must be non-negative"))
    end

    P = Partition(limits)
    idx = ones(Int64, size(X, 2))
    model_cache = [Array{BayesPM,3}(undef, size(X, 1), 2, Jmax + 1)]
    basis_cache = [tpbasis(size(X, 1), J) for J in 0:Jmax]
    models = Vector{BayesPM}(undef, 1)
    space_vol = vol(limits)

    # Set up the intial polynomial for the full space and clear the cache.
    # We don't need to cache the inital entry as we have already accepted it.
    # The cache indices supplied to _conditional_degree_selection! are just
    # placeholders
    models[1] = _conditional_degree_selection!(
        X,
        y,
        1,
        1,
        1,
        limits,
        Jmax,
        Pmax,
        models,
        model_cache,
        basis_cache,
        λ,
        shape0,
        scale0,
        ratio,
    )
    model_cache[1] = Array{BayesPM,3}(undef, size(X, 1), 2, Jmax + 1)
    ev = evidence(models, shape0, scale0)

    while length(models) < Kmax
        count = 0
        accepted = false
        K = length(models)
        for k in randperm(K)
            # Iterate over every subregion and test to see if any split
            # improves the model evidence beyond the tolerance.

            count += 1
            best = Dict{String,Any}("ev" => ev + tol)
            models_cp = deepcopy(models)
            push!(models_cp, models_cp[1])
            region_mask = idx .== k
            Xk, yk = X[:, region_mask], y[:, region_mask]
            for d in 1:size(X, 1)
                # Region k is split in every dimension and two polynomials are fitted
                # to the resulting subregions. The best details of the best dimensional
                # split are stored in the dictionary `best`

                if verbose
                    msg = "\rRegion = $count/$K; Dimension = $d/$(size(X, 1)); Evidence = $ev"
                    print(msg)
                end
                loc = sum(P.regions[k][d, :]) / 2
                left_lims, right_lims = deepcopy(P.regions[k]), deepcopy(P.regions[k])
                left_lims[d, 2], right_lims[d, 1] = loc, loc
                left_region_mask = Xk[d, :] .< loc
                left_pm = _conditional_degree_selection!(
                    Xk[:, left_region_mask],
                    yk[:, left_region_mask],
                    k,
                    d,
                    1,
                    left_lims,
                    Jmax,
                    Pmax,
                    models,
                    model_cache,
                    basis_cache,
                    λ * vol(left_lims) / space_vol,
                    shape0,
                    scale0,
                    ratio,
                )
                right_pm = _conditional_degree_selection!(
                    Xk[:, .!left_region_mask],
                    yk[:, .!left_region_mask],
                    k,
                    d,
                    2,
                    right_lims,
                    Jmax,
                    Pmax,
                    models,
                    model_cache,
                    basis_cache,
                    λ * vol(right_lims) / space_vol,
                    shape0,
                    scale0,
                    ratio,
                )
                models_cp[k] = left_pm
                models_cp[end] = right_pm
                tmp_ev = evidence(models_cp, shape0, scale0)
                if tmp_ev > best["ev"] >= ev + tol
                    accepted = true
                    best["ev"] = tmp_ev
                    best["d"] = d
                    best["left_region_mask"] = left_region_mask
                    best["left_pm"], best["right_pm"] = left_pm, right_pm
                end
            end
            if accepted
                # Update models and reset cache for k'th region. Add extra
                # cache entry for the additional region.
                models[k] = best["left_pm"]::BayesPM
                push!(models, best["right_pm"]::BayesPM)
                model_cache[k] = Array{BayesPM,3}(undef, size(X, 1), 2, Jmax + 1)
                push!(model_cache, Array{BayesPM,3}(undef, size(X, 1), 2, Jmax + 1))

                # Update the region indices
                region_idx = @view idx[region_mask]
                region_idx[.!best["left_region_mask"]::BitVector] .= length(models)

                split!(P, k, best["d"]::Int64)
                ev = best["ev"]::Float64
                break
            end
        end
        if !accepted
            break # Every split has been rejected. Terminate search.
        end
    end
    return PartitionedBayesPM(P, models, shape0, scale0)
end

function (ppm::PartitionedBayesPM)(X::AbstractMatrix)
    idx = locate(ppm.P, X)
    y = zeros(1, length(idx))
    for k in unique(idx)
        region_mask = idx .== k
        y[:, region_mask] = ppm.models[k](X[:, region_mask])
    end
    return y
end

ext: jl


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApimachineryPkgApisMetaV1GroupVersionForDiscovery <: SwaggerModel
    groupVersion::Any # spec type: Union{ Nothing, String } # spec name: groupVersion
    version::Any # spec type: Union{ Nothing, String } # spec name: version

    function IoK8sApimachineryPkgApisMetaV1GroupVersionForDiscovery(;groupVersion=nothing, version=nothing)
        o = new()
        validate_property(IoK8sApimachineryPkgApisMetaV1GroupVersionForDiscovery, Symbol("groupVersion"), groupVersion)
        setfield!(o, Symbol("groupVersion"), groupVersion)
        validate_property(IoK8sApimachineryPkgApisMetaV1GroupVersionForDiscovery, Symbol("version"), version)
        setfield!(o, Symbol("version"), version)
        o
    end
end # type IoK8sApimachineryPkgApisMetaV1GroupVersionForDiscovery

const _property_map_IoK8sApimachineryPkgApisMetaV1GroupVersionForDiscovery = Dict{Symbol,Symbol}(Symbol("groupVersion")=>Symbol("groupVersion"), Symbol("version")=>Symbol("version"))
const _property_types_IoK8sApimachineryPkgApisMetaV1GroupVersionForDiscovery = Dict{Symbol,String}(Symbol("groupVersion")=>"String", Symbol("version")=>"String")
Base.propertynames(::Type{ IoK8sApimachineryPkgApisMetaV1GroupVersionForDiscovery }) = collect(keys(_property_map_IoK8sApimachineryPkgApisMetaV1GroupVersionForDiscovery))
Swagger.property_type(::Type{ IoK8sApimachineryPkgApisMetaV1GroupVersionForDiscovery }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApimachineryPkgApisMetaV1GroupVersionForDiscovery[name]))}
Swagger.field_name(::Type{ IoK8sApimachineryPkgApisMetaV1GroupVersionForDiscovery }, property_name::Symbol) =  _property_map_IoK8sApimachineryPkgApisMetaV1GroupVersionForDiscovery[property_name]

function check_required(o::IoK8sApimachineryPkgApisMetaV1GroupVersionForDiscovery)
    (getproperty(o, Symbol("groupVersion")) === nothing) && (return false)
    (getproperty(o, Symbol("version")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApimachineryPkgApisMetaV1GroupVersionForDiscovery }, name::Symbol, val)
end


text:
# kernel execution

export AbstractGPUBackend, AbstractKernelContext, gpu_call

abstract type AbstractGPUBackend end

abstract type AbstractKernelContext end

"""
    backend(T::Type)
    backend(x)

Gets the GPUArrays back-end responsible for managing arrays of type `T`.
"""
backend(::Type) = error("This object is not a GPU array") # COV_EXCL_LINE
backend(x) = backend(typeof(x))

"""
    gpu_call(kernel::Function, arg0, args...; kwargs...)

Executes `kernel` on the device that backs `arg` (see [`backend`](@ref)), passing along any
arguments `args`. Additionally, the kernel will be passed the kernel execution context (see
[`AbstractKernelContext`]), so its signature should be `(ctx::AbstractKernelContext, arg0,
args...)`.

The keyword arguments `kwargs` are not passed to the function, but are interpreted on the
host to influence how the kernel is executed. The following keyword arguments are supported:

- `target::AbstractArray`: specify which array object to use for determining execution
  properties (defaults to the first argument `arg0`).
- `elements::Int`: how many elements will be processed by this kernel. In most
  circumstances, this will correspond to the total number of threads that needs to be
  launched, unless the kernel supports a variable number of elements to process per
  iteration. Defaults to the length of `arg0` if no other keyword arguments that influence
  the launch configuration are specified.
- `threads::Int` and `blocks::Int`: configure exactly how many threads and blocks are
  launched. This cannot be used in combination with the `elements` argument.
- `name::String`: inform the back end about the name of the kernel to be executed. This can
  be used to emit better diagnostics, and is useful with anonymous kernels.
"""
function gpu_call(kernel::F, args::Vararg{Any,N};
                  target::AbstractArray=first(args),
                  elements::Union{Int,Nothing}=nothing,
                  threads::Union{Int,Nothing}=nothing,
                  blocks::Union{Int,Nothing}=nothing,
                  name::Union{String,Nothing}=nothing) where {F,N}
    # non-trivial default values for launch configuration
    if elements===nothing && threads===nothing && blocks===nothing
        elements = length(target)
    elseif elements===nothing
        if threads === nothing
            threads = 1
        end
        if blocks === nothing
            blocks = 1
        end
    elseif threads!==nothing || blocks!==nothing
        error("Cannot specify both elements and threads/blocks configuration")
    end

    # the number of elements to process needs to be passed to the kernel somehow, so there's
    # no easy way to do this without passing additional arguments or changing the context.
    # both are expensive, so require manual use of `launch_heuristic` for those kernels.
    elements_per_thread = 1

    if elements !== nothing
        @assert elements > 0
        heuristic = launch_heuristic(backend(target), kernel, args...;
                                     elements, elements_per_thread)
        config = launch_configuration(backend(target), heuristic;
                                      elements, elements_per_thread)
        gpu_call(backend(target), kernel, args, config.threads, config.blocks; name=name)
    else
        @assert threads > 0
        @assert blocks > 0
        gpu_call(backend(target), kernel, args, threads, blocks; name=name)
    end
end

# how many threads and blocks `kernel` needs to be launched with, passing arguments `args`,
# to fully saturate the GPU. `elements` indicates the number of elements that needs to be
# processed, while `elements_per_threads` indicates the number of elements this kernel can
# process (i.e. if it's a grid-stride kernel, or 1 if otherwise).
#
# this heuristic should be specialized for the back-end, ideally using an API for maximizing
# the occupancy of the launch configuration (like CUDA's occupancy API).
function launch_heuristic(backend::AbstractGPUBackend, kernel, args...;
                          elements::Int, elements_per_thread::Int)
    return (threads=256, blocks=32)
end

# determine how many threads and blocks to actually launch given upper limits.
# returns a tuple of blocks, threads, and elements_per_thread (which is always 1
# unless specified that the kernel can handle a number of elements per thread)
function launch_configuration(backend::AbstractGPUBackend, heuristic;
                              elements::Int, elements_per_thread::Int)
    threads = clamp(elements, 1, heuristic.threads)
    blocks = max(cld(elements, threads), 1)

    if elements_per_thread > 1 && blocks > heuristic.blocks
        # we want to launch more blocks than required, so prefer a grid-stride loop instead
        ## try to stick to the number of blocks that the heuristic suggested
        blocks = heuristic.blocks
        nelem = cld(elements, blocks*threads)
        ## only bump the number of blocks if we really need to
        if nelem > elements_per_thread
            nelem = elements_per_thread
            blocks = cld(elements, nelem*threads)
        end
        (; threads, blocks, elements_per_thread=nelem)
    else
        (; threads, blocks, elements_per_thread=1)
    end
end

gpu_call(backend::AbstractGPUBackend, kernel, args, threads::Int, blocks::Int; kwargs...) =
    error("Not implemented") # COV_EXCL_LINE


text:
idxno = [rand(1:themis_nopkt.numParticles) for i in 1:10]
p1 = plot((5.2*themis_nopkt.allT[idxno]*Re)./(c), themis_nopkt.allPA[idxno], title = "Long packet (dPhi = 1)",
legend=false,xlabel="Time (ms)", ylabel="Energy (keV)",
xlim = [0,12], ylim = [0,90])
idx30 = [rand(1:themis_hilat.numParticles) for i in 1:10]
p2 = plot((5.2*themis_hilat.allT[idxno]*Re)./(c), themis_hilat.allPA[idxno], title = "Short Packet (dPhi = 30)",
legend=false,xlabel="Time (ms)", ylabel="Energy (keV)",
xlim = [0,12], ylim = [0,90])
idx100 = [rand(1:themis_100pkt.numParticles) for i in 1:10]
p3 = plot((5.2*themis_100pkt.allT[idxno]*Re)./(c), themis_100pkt.allPA[idxno], title = "Short Packet (dPhi = 100)",
legend=false,xlabel="Time (ms)", ylabel="Energy (keV)",
xlim = [0,12], ylim = [0,90])
bigplot = plot(p1,p2,p3, dpi = 300, layout = (3,1))

ylimmax = maximum([themis_nopkt.lostParticles[end,2], themis_hilat.lostParticles[end,2], themis_100pkt.lostParticles[end,2]])
p1 = plot(L*themis_nopkt.lostParticles[:,1]*Re/(c), themis_nopkt.lostParticles[:,2], legend=false,
    xlim = [0,L*maximum(maximum(themis_nopkt.allT))*Re/(c)], ylim = [0,ylimmax],
    title = "Number of Lost Particles (dphi = 1)", xlabel="Time (s)", ylabel="Number of Particles")
p2 = plot(L*themis_hilat.lostParticles[:,1]*Re/(c), themis_hilat.lostParticles[:,2], legend=false,
    xlim = [0,L*maximum(maximum(themis_hilat.allT))*Re/(c)], ylim = [0,ylimmax],
    title = "Number of Lost Particles (dphi = 30)", xlabel="Time (s)", ylabel="Number of Particles")
p3 = plot(L*themis_100pkt.lostParticles[:,1]*Re/(c), themis_100pkt.lostParticles[:,2], legend=false,
    xlim = [0,L*maximum(maximum(themis_100pkt.allT))*Re/(c)], ylim = [0,ylimmax],
    title = "Number of Lost Particles (dphi = 100)", xlabel="Time (s)", ylabel="Number of Particles")
bigplot = plot(p1,p2,p3, dpi = 300, layout = (3,1))

p1 = histogram(diff(one.lostParticles[:,2]),L*one.lostParticles[1:end-1,1]*Re/(c), bins= 100,ylim = (0,250))
p2 = histogram(diff(themis_hilat.lostParticles[:,2]),L*themis_hilat.lostParticles[1:end-1,1]*Re/(c), bins= 100,ylim = (0,250))
p3 = histogram(diff(themis_100pkt.lostParticles[:,2]),L*themis_100pkt.lostParticles[1:end-1,1]*Re/(c), bins= 100,ylim = (0,250))
bigplot = plot(p1,p2,p3, dpi = 300, layout = (3,1))
plot!(xlim  = [2.5,10], ylim = [0,75])

idxno = [rand(1:one.numParticles) for i in 1:10]
p1 = plot((5.2*one.allT*Re)./(c), one.allE, title = "Long packet (dPhi = 1)",
legend=false,xlabel="Time (ms)", ylabel="Energy (keV)",
xlim = [0,1], ylim = [0,1000])
idx30 = [rand(1:hund.numParticles) for i in 1:10]
p2 = plot((5.2*hund.allT*Re)./(c), hund.allE, title = "Short Packet (dPhi = 100)",
legend=false,xlabel="Time (ms)", ylabel="Energy (keV)",
xlim = [0,12], ylim = [0,1000])
Eplot = plot(p1,p2, dpi = 300, layout = (2,1))
p3 = plot((5.2*one.allT*Re)./(c), one.allPA, title = "Long packet (dPhi = 1)",
legend=false,xlabel="Time (ms)", ylabel="PA (deg)",
xlim = [0,12], ylim = [0,90])
idx30 = [rand(1:hund.numParticles) for i in 1:10]
p4 = plot((5.2*hund.allT*Re)./(c), hund.allPA, title = "Short Packet (dPhi = 100)",
legend=false,xlabel="Time (ms)", ylabel="PA (deg)",
xlim = [0,12], ylim = [0,90])
PAPlot = plot(p3,p4, dpi = 300, layout = (2,1))

p1 = plot((5.2*one.allT*Re)./(c), one.allPA, title = "Short packet (dPhi = 1)",
legend=false,xlabel="Time (ms)", ylabel="PA (deg)",
xlim = [0,0.5], ylim = [0,15])
p2 = plot((5.2*thi.allT*Re)./(c), thi.allPA, title = "Normal Packet (dPhi = 30)",
legend=false,xlabel="Time (ms)", ylabel="PA (deg)",
xlim = [0,0.5], ylim = [0,15])
p3 = plot((5.2*hun.allT*Re)./(c), hun.allPA, title = "Long Packet (dPhi = 100)",
legend=false,xlabel="Time (ms)", ylabel="PA (deg)",
xlim = [0,0.5], ylim = [0,15])
bigplot = plot(p1,p2,p3, dpi = 300, layout = (3,1))

one = load_resultant_matrix("test1", "results/jld2_211124_11", "PA_TEST_1_100", "setupasrun.conf", 1);
thi = load_resultant_matrix("test30", "results/jld2_211124_09", "PA_TEST_30_100", "setupasrun.conf", 1);
hun = load_resultant_matrix("test100", "results/jld2_211124_11", "PA_TEST_100_100", "setupasrun.conf", 1);

rm = hun
ipa = PAmatrix[1,:]
epa = Vector{Float64}()
for col in eachcol(PAmatrix)
    push!(epa,col[findall(!isnan, col)][end])
end
scatter!(ipa, epa)

directory="results/jld2_211124_22"
basename = "PA_TEST_100_n_100"
allZ = Vector{Vector{Float64}}();
allPZ = Vector{Vector{Float64}}();
allE = Vector{Vector{Float64}}();
allPA = Vector{Vector{Float64}}();
allT = Vector{Vector{Float64}}();
allG = Vector{Vector{Float64}}();
# allLambda = Vector{Vector{Float64}}();
# allPsinou = Vector{Vector{Float64}}();
for i in 1:1
    JLD2.@load directory*"/"*basename*"_$i.jld2" sol
    for traj in sol
        vars = Array(traj')
        timesteps = length(traj.t)
        b = zeros(timesteps)
        gamma = zeros(timesteps)
        Alpha = zeros(timesteps)

        @views calcb!(b,vars[:,5])
        @views calcGamma!(gamma,vars[:,2],vars[:,4],b)
        @views calcAlpha!(Alpha,vars[:,4],gamma)
        @views push!(allT, traj.t);
        @views push!(allZ, vars[:,1]);
        @views push!(allPZ, vars[:,2]);
        @views push!(allPA, Alpha);
        @views push!(allE, @. (511*(gamma - 1)));
        # @views push!(allG, vars[:,7])
        # @views push!(allLambda, vars[:,5])
        # @views push!(allPsinou, vars[:,8])
    end
end
@time tVec, Zmatrix, PZmatrix, Ematrix, PAmatrix = postProcessor(allT, allZ, allPZ, allE, allPA);
ipa = PAmatrix[1,:]
epa = Vector{Float64}()
for col in eachcol(PAmatrix)
    push!(epa,col[findall(!isnan, col)][end])
end
scatter!(ipa, epa, legend=false)

p1 = plot((5.2*one.allT*Re)./(c), one.allPA, title = "Short packet (dPhi = 1)",
legend=false,xlabel="Time (ms)", ylabel="PA (deg)",
xlim = [0,0.5], ylim = [0,15])
p2 = plot((5.2*thi.allT*Re)./(c), thi.allPA, title = "Normal Packet (dPhi = 30)",
legend=false,xlabel="Time (ms)", ylabel="PA (deg)",
xlim = [0,0.5], ylim = [0,15])
p3 = plot((allLambda, allPsi, title = "Long Packet (dPhi = 100)",
legend=false,xlabel="Time (ms)", ylabel="Psi (rad)",
xlim = [-pi/4, pi/4], ylim = [-5,5])
bigplot = plot(p1,p2,p3, dpi = 300, layout = (3,1))

max_stars_repo_name: ethantsai/nlwhistlers

max_issues_repo_name: ethantsai/nlwhistlers

max_forks_repo_name: ethantsai/nlwhistlers


text:
module QuantumOptics

using SparseArrays, LinearAlgebra

export bases, Basis, GenericBasis, CompositeBasis, basis,
        tensor, ⊗, permutesystems, @samebases,
        states, StateVector, Bra, Ket, basisstate, norm,
                dagger, normalize, normalize!,
        operators, AbstractOperator, DataOperator, expect, variance,
            identityoperator, ptrace, embed, dense, tr, sparse,
        operators_dense, DenseOperator, projector, dm,
        operators_sparse, SparseOperator, diagonaloperator,
        operators_lazysum, LazySum,
        operators_lazyproduct, LazyProduct,
        operators_lazytensor, LazyTensor,
        superoperators, SuperOperator, DenseSuperOperator, SparseSuperOperator,
                spre, spost, liouvillian,
        fock, FockBasis, number, destroy, create,
                fockstate, coherentstate, coherentstate!, displace,
        randstate, randoperator, thermalstate, coherentthermalstate, phase_average, passive_state,
        spin, SpinBasis, sigmax, sigmay, sigmaz, sigmap, sigmam, spinup, spindown,
        subspace, SubspaceBasis, projector,
        particle, PositionBasis, MomentumBasis, samplepoints, spacing, gaussianstate,
                position, momentum, potentialoperator, transform,
        nlevel, NLevelBasis, transition, nlevelstate,
        manybody, ManyBodyBasis, fermionstates, bosonstates,
                manybodyoperator, onebodyexpect, occupation,
        phasespace, qfunc, wigner, coherentspinstate, qfuncsu2, wignersu2, ylm,
        metrics, tracenorm, tracenorm_h, tracenorm_nh,
                tracedistance, tracedistance_h, tracedistance_nh,
                entropy_vn, fidelity, ptranspose, PPT,
                negativity, logarithmic_negativity,
        spectralanalysis, eigenstates, eigenenergies, simdiag,
        timeevolution, diagonaljumps, @skiptimechecks,
        steadystate,
        timecorrelations,
        semiclassical,
        stochastic

include("sortedindices.jl")
include("polynomials.jl")
include("bases.jl")
include("states.jl")
include("operators.jl")
include("operators_dense.jl")
include("sparsematrix.jl")
include("operators_sparse.jl")
include("operators_lazysum.jl")
include("operators_lazyproduct.jl")
include("operators_lazytensor.jl")
include("superoperators.jl")
include("spin.jl")
include("fock.jl")
include("state_definitions.jl")
include("subspace.jl")
include("particle.jl")
include("nlevel.jl")
include("manybody.jl")
include("transformations.jl")
include("phasespace.jl")
include("metrics.jl")
module timeevolution
    export diagonaljumps, @skiptimechecks
    include("timeevolution_base.jl")
    include("master.jl")
    include("schroedinger.jl")
    include("mcwf.jl")
    using .timeevolution_master
    using .timeevolution_schroedinger
    using .timeevolution_mcwf
end
include("steadystate.jl")
include("timecorrelations.jl")
include("spectralanalysis.jl")
include("semiclassical.jl")
module stochastic
    include("stochastic_definitions.jl")
    include("stochastic_schroedinger.jl")
    include("stochastic_master.jl")
    include("stochastic_semiclassical.jl")
    using .stochastic_schroedinger, .stochastic_master, .stochastic_semiclassical
    using .stochastic_definitions
end
include("printing.jl")

using .bases
using .states
using .operators
using .operators_dense
using .operators_sparse
using .operators_lazysum
using .operators_lazyproduct
using .operators_lazytensor
using .superoperators
using .spin
using .fock
using .state_definitions
using .subspace
using .particle
using .nlevel
using .manybody
using .phasespace
using .timeevolution
using .metrics
using .spectralanalysis
using .timecorrelations
using .printing

end # module


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

using Test
using REPL
using Random
import REPL.LineEdit
using Markdown

const BASE_TEST_PATH = joinpath(Sys.BINDIR, "..", "share", "julia", "test")
isdefined(Main, :FakePTYs) || @eval Main include(joinpath($(BASE_TEST_PATH), "testhelpers", "FakePTYs.jl"))
import .Main.FakePTYs: with_fake_pty

# For curmod_*
include(joinpath(BASE_TEST_PATH, "testenv.jl"))

include("FakeTerminals.jl")
import .FakeTerminals.FakeTerminal

function kill_timer(delay)
    # Give ourselves a generous timer here, just to prevent
    # this causing e.g. a CI hang when there's something unexpected in the output.
    # This is really messy and leaves the process in an undefined state.
    # the proper and correct way to do this in real code would be to destroy the
    # IO handles: `close(stdout_read); close(stdin_write)`
    test_task = current_task()
    function kill_test(t)
        # **DON'T COPY ME.**
        # The correct way to handle timeouts is to close the handle:
        # e.g. `close(stdout_read); close(stdin_write)`
        test_task.queue === nothing || Base.list_deletefirst!(test_task.queue, test_task)
        schedule(test_task, "hard kill repl test"; error=true)
        print(stderr, "WARNING: attempting hard kill of repl test after exceeding timeout\n")
    end
    return Timer(kill_test, delay)
end

# REPL tests
function fake_repl(@nospecialize(f); options::REPL.Options=REPL.Options(confirm_exit=false))
    # Use pipes so we can easily do blocking reads
    # In the future if we want we can add a test that the right object
    # gets displayed by intercepting the display
    input = Pipe()
    output = Pipe()
    err = Pipe()
    Base.link_pipe!(input, reader_supports_async=true, writer_supports_async=true)
    Base.link_pipe!(output, reader_supports_async=true, writer_supports_async=true)
    Base.link_pipe!(err, reader_supports_async=true, writer_supports_async=true)

    repl = REPL.LineEditREPL(FakeTerminal(input.out, output.in, err.in, options.hascolor), options.hascolor)
    repl.options = options

    hard_kill = kill_timer(900) # Your debugging session starts now. You have 15 minutes. Go.
    f(input.in, output.out, repl)
    t = @async begin
        close(input.in)
        close(output.in)
        close(err.in)
    end
    @test read(err.out, String) == ""
    #display(read(output.out, String))
    Base.wait(t)
    close(hard_kill)
    nothing
end

# Writing ^C to the repl will cause sigint, so let's not die on that
Base.exit_on_sigint(false)

# make sure `run_interface` can normally handle `eof`
# without any special handling by the user
fake_repl() do stdin_write, stdout_read, repl
    panel = LineEdit.Prompt("test";
        prompt_prefix = "",
        prompt_suffix = Base.text_colors[:white],
        on_enter = s -> true)
    panel.on_done = (s, buf, ok) -> begin
        @test !ok
        @test bytesavailable(buf) == position(buf) == 0
        nothing
    end
    repltask = @async REPL.run_interface(repl.t, LineEdit.ModalInterface(Any[panel]))
    close(stdin_write)
    Base.wait(repltask)
end

# These are integration tests. If you want to unit test test e.g. completion, or
# exact LineEdit behavior, put them in the appropriate test files.
# Furthermore since we are emulating an entire terminal, there may be control characters
# in the mix. If verification needs to be done, keep it to the bare minimum. Basically
# this should make sure nothing crashes without depending on how exactly the control
# characters are being used.
fake_repl(options = REPL.Options(confirm_exit=false,hascolor=true)) do stdin_write, stdout_read, repl
    repl.specialdisplay = REPL.REPLDisplay(repl)
    repl.history_file = false

    repltask = @async begin
        REPL.run_repl(repl)
    end

    global inc = false
    global b = Condition()
    global c = Condition()
    let cmd = "\"Hello REPL\""
        write(stdin_write, "$(curmod_prefix)inc || wait($(curmod_prefix)b); r = $cmd; notify($(curmod_prefix)c); r\r")
    end
    let t = @async begin
            inc = true
            notify(b)
            wait(c)
        end
        while (d = readline(stdout_read)) != ""
            # first line [optional]: until 80th char of input
            # second line: until end of input
            # third line: "Hello REPL"
            # last line: blank
            # last+1 line: next prompt
        end
        wait(t)
    end

    # Latex completions
    readuntil(stdout_read, "julia> ", keep=true)
    write(stdin_write, "\x32\\alpha\t")
    readuntil(stdout_read, "α")
    # Bracketed paste in search mode
    write(stdin_write, "\e[200~paste here ;)\e[201~")
    # Abort search (^C)
    write(stdin_write, '\x03')
    # Test basic completion in main mode
    write(stdin_write, "Base.REP\t")
    readuntil(stdout_read, "REPL")
    write(stdin_write, '\x03')
    write(stdin_write, "\\alpha\t")
    readuntil(stdout_read,"α")
    write(stdin_write, '\x03')
    # Test cd feature in shell mode.
    origpwd = pwd()
    mktempdir() do tmpdir
        try
            samefile = Base.Filesystem.samefile
            tmpdir_pwd = cd(pwd, tmpdir)
            homedir_pwd = cd(pwd, homedir())

            # Test `cd`'ing to an absolute path
            write(stdin_write, ";")
            readuntil(stdout_read, "shell> ")
            write(stdin_write, "cd $(escape_string(tmpdir))\n")
            readuntil(stdout_read, "cd $(escape_string(tmpdir))")
            readuntil(stdout_read, tmpdir_pwd)
            readuntil(stdout_read, "\n")
            readuntil(stdout_read, "\n")
            @test samefile(".", tmpdir)
            write(stdin_write, "\b")

            # Test using `cd` to move to the home directory
            write(stdin_write, ";")
            readuntil(stdout_read, "shell> ")
            write(stdin_write, "cd\n")
            readuntil(stdout_read, homedir_pwd)
            readuntil(stdout_read, "\n")
            readuntil(stdout_read, "\n")
            @test samefile(".", homedir_pwd)
            write(stdin_write, "\b")

            # Test using `-` to jump backward to tmpdir
            write(stdin_write, ";")
            readuntil(stdout_read, "shell> ")
            write(stdin_write, "cd -\n")
            readuntil(stdout_read, tmpdir_pwd)
            readuntil(stdout_read, "\n")
            readuntil(stdout_read, "\n")
            @test samefile(".", tmpdir)
            write(stdin_write, "\b")

            # Test using `~` (Base.expanduser) in `cd` commands
            if !Sys.iswindows()
                write(stdin_write, ";")
                readuntil(stdout_read, "shell> ")
                write(stdin_write, "cd ~\n")
                readuntil(stdout_read, homedir_pwd)
                readuntil(stdout_read, "\n")
                readuntil(stdout_read, "\n")
                @test samefile(".", homedir_pwd)
                write(stdin_write, "\b")
            end
        finally
            cd(origpwd)
        end
    end

    # issue #20482
    #if !Sys.iswindows()
    #    write(stdin_write, ";")
    #    readuntil(stdout_read, "shell> ")
    #    write(stdin_write, "echo hello >/dev/null\n")
    #    let s = readuntil(stdout_read, "\n", keep=true)
    #        @test occursin("shell> ", s) # make sure we echoed the prompt
    #        @test occursin("echo hello >/dev/null", s) # make sure we echoed the input
    #    end
    #    @test readuntil(stdout_read, "\n", keep=true) == "\e[0m\n"
    #end

    # issue #20771
    let s
        write(stdin_write, ";")
        readuntil(stdout_read, "shell> ")
        write(stdin_write, "'\n") # invalid input
        s = readuntil(stdout_read, "\n")
        @test occursin("shell> ", s) # check for the echo of the prompt
        @test occursin("'", s) # check for the echo of the input
        s = readuntil(stdout_read, "\n\n")
        @test startswith(s, "\e[0mERROR: unterminated single quote\nStacktrace:\n  [1] ") ||
              startswith(s, "\e[0m\e[1m\e[91mERROR: \e[39m\e[22m\e[91munterminated single quote\e[39m\nStacktrace:\n  [1] ")
        write(stdin_write, "\b")
    end

    # issue #27293
    if Sys.isunix()
        let s, old_stdout = stdout
            write(stdin_write, ";")
            readuntil(stdout_read, "shell> ")
            write(stdin_write, "echo ~")
            s = readuntil(stdout_read, "~")

            proc_stdout_read, proc_stdout = redirect_stdout()
            get_stdout = @async read(proc_stdout_read, String)
            try
                write(stdin_write, "\n")
                readuntil(stdout_read, "\n")
                s = readuntil(stdout_read, "\n")
            finally
                redirect_stdout(old_stdout)
            end
            @test s == "\e[0m" # the child has exited
            close(proc_stdout)
            # check for the correct, expanded response
            @test occursin(expanduser("~"), fetch(get_stdout))
            write(stdin_write, "\b")
        end
    end

    # issues #22176 & #20482
    # TODO: figure out how to test this on Windows
    #Sys.iswindows() || let tmp = tempname()
    #    try
    #        write(stdin_write, ";")
    #        readuntil(stdout_read, "shell> ")
    #        write(stdin_write, "echo \$123 >$tmp\n")
    #        let s = readuntil(stdout_read, "\n")
    #            @test occursin("shell> ", s) # make sure we echoed the prompt
    #            @test occursin("echo \$123 >$tmp", s) # make sure we echoed the input
    #        end
    #        @test readuntil(stdout_read, "\n", keep=true) == "\e[0m\n"
    #        @test read(tmp, String) == "123\n"
    #    finally
    #        rm(tmp, force=true)
    #    end
    #end

    # issue #10120
    # ensure that command quoting works correctly
    let s, old_stdout = stdout
        write(stdin_write, ";")
        readuntil(stdout_read, "shell> ")
        Base.print_shell_escaped(stdin_write, Base.julia_cmd().exec..., special=Base.shell_special)
        write(stdin_write, """ -e "println(\\"HI\\")\" """)
        readuntil(stdout_read, ")\"")
        proc_stdout_read, proc_stdout = redirect_stdout()
        get_stdout = @async read(proc_stdout_read, String)
        try
            write(stdin_write, '\n')
            s = readuntil(stdout_read, "\n", keep=true)
            if s == "\n"
                # if shell width is precisely the text width,
                # we may print some extra characters to fix the cursor state
                s = readuntil(stdout_read, "\n", keep=true)
                @test occursin("shell> ", s)
                s = readuntil(stdout_read, "\n", keep=true)
                @test s == "\r\r\n"
            else
                @test occursin("shell> ", s)
            end
            s = readuntil(stdout_read, "\n", keep=true)
            @test s == "\e[0m\n" # the child has exited
        finally
            redirect_stdout(old_stdout)
        end
        close(proc_stdout)
        @test fetch(get_stdout) == "HI\n"
        write(stdin_write, "\b")
    end

    # Issue #7001
    # Test ignoring '\0'
    let
        write(stdin_write, "\0\n")
        s = readuntil(stdout_read, "\n\n")
        @test !occursin("invalid character", s)
    end

    # Test that accepting a REPL result immediately shows up, not
    # just on the next keystroke
    write(stdin_write, "1+1\n") # populate history with a trivial input
    readline(stdout_read)
    write(stdin_write, "\e[A\n")
    let t = kill_timer(60)
        # yield make sure this got processed
        readuntil(stdout_read, "1+1")
        readuntil(stdout_read, "\n\n")
        close(t) # cancel timeout
    end

    # Issue #10222
    # Test ignoring insert key in standard and prefix search modes
    write(stdin_write, "\e[2h\e[2h\n") # insert (VT100-style)
    @test findfirst("[2h", readline(stdout_read)) === nothing
    readline(stdout_read)
    write(stdin_write, "\e[2~\e[2~\n") # insert (VT220-style)
    @test findfirst("[2~", readline(stdout_read)) === nothing
    readline(stdout_read)
    write(stdin_write, "1+1\n") # populate history with a trivial input
    readline(stdout_read)
    write(stdin_write, "\e[A\e[2h\n") # up arrow, insert (VT100-style)
    readline(stdout_read)
    readline(stdout_read)
    write(stdin_write, "\e[A\e[2~\n") # up arrow, insert (VT220-style)
    readline(stdout_read)
    readline(stdout_read)

    # Test down arrow to go back to history
    # populate history with a trivial input

    s1 = "12345678"; s2 = "23456789"
    write(stdin_write, s1, '\n')
    readuntil(stdout_read, s1)
    write(stdin_write, s2, '\n')
    readuntil(stdout_read, s2)
    # Two up arrow, enter, should get back to 1
    write(stdin_write, "\e[A\e[A\n")
    readuntil(stdout_read, s1)
    # Now, down arrow, enter, should get us back to 2
    write(stdin_write, "\e[B\n")
    readuntil(stdout_read, s2)

    # test that prefix history search "passes through" key bindings to parent mode
    write(stdin_write, "0x321\n")
    readuntil(stdout_read, "0x321")
    write(stdin_write, "\e[A\e[1;3C|||") # uparrow (go up history) and then Meta-rightarrow (indent right)
    s2 = readuntil(stdout_read, "|||", keep=true)
    @test endswith(s2, " 0x321\r\e[13C|||") # should have a space (from Meta-rightarrow) and not
                                            # have a spurious C before ||| (the one here is not spurious!)

    # "pass through" for ^x^x
    write(stdin_write, "\x030x4321\n") # \x03 == ^c
    readuntil(stdout_read, "0x4321")
    write(stdin_write, "\e[A\x18\x18||\x18\x18||||") # uparrow, ^x^x||^x^x||||
    s3 = readuntil(stdout_read, "||||", keep=true)
    @test endswith(s3, "||0x4321\r\e[15C||||")

    # Delete line (^U) and close REPL (^D)
    write(stdin_write, "\x15\x04")
    Base.wait(repltask)

    nothing
end

function buffercontents(buf::IOBuffer)
    p = position(buf)
    seek(buf,0)
    c = read(buf, String)
    seek(buf,p)
    c
end

function AddCustomMode(repl, prompt)
    # Custom REPL mode tests
    foobar_mode = LineEdit.Prompt(prompt;
        prompt_prefix="\e[38;5;166m",
        prompt_suffix=Base.text_colors[:white],
        on_enter = s->true,
        on_done = line->true)

    main_mode = repl.interface.modes[1]
    push!(repl.interface.modes,foobar_mode)

    hp = main_mode.hist
    hp.mode_mapping[:foobar] = foobar_mode
    foobar_mode.hist = hp

    foobar_keymap = Dict{Any,Any}(
        '<' => function (s,args...)
            if isempty(s)
                if !haskey(s.mode_state,foobar_mode)
                    s.mode_state[foobar_mode] = LineEdit.init_state(repl.t,foobar_mode)
                end
                LineEdit.transition(s,foobar_mode)
            else
                LineEdit.edit_insert(s,'<')
            end
        end
    )

    search_prompt, skeymap = LineEdit.setup_search_keymap(hp)
    mk = REPL.mode_keymap(main_mode)

    b = Dict{Any,Any}[skeymap, mk, LineEdit.history_keymap, LineEdit.default_keymap, LineEdit.escape_defaults]
    foobar_mode.keymap_dict = LineEdit.keymap(b)

    main_mode.keymap_dict = LineEdit.keymap_merge(main_mode.keymap_dict, foobar_keymap)
    foobar_mode, search_prompt
end

# Note: since the \t character matters for the REPL file history,
# it is important not to have the """ code reindent this line,
# possibly converting \t to spaces.
fakehistory = """
# time: 2014-06-29 20:44:29 EDT
# mode: julia
\té
# time: 2014-06-29 21:44:29 EDT
# mode: julia
\téé
# time: 2014-06-30 17:32:49 EDT
# mode: julia
\tshell
# time: 2014-06-30 17:32:59 EDT
# mode: shell
\tll
# time: 2014-06-30 99:99:99 EDT
# mode: julia
\tx ΔxΔ
# time: 2014-06-30 17:32:49 EDT
# mode: julia
\t1 + 1
# time: 2014-06-30 17:35:39 EDT
# mode: foobar
\tbarfoo
# time: 2014-06-30 18:44:29 EDT
# mode: shell
\tls
# time: 2014-06-30 19:44:29 EDT
# mode: foobar
\tls
# time: 2014-06-30 20:44:29 EDT
# mode: julia
\t2 + 2
"""

# Test various history related issues
for prompt = ["TestΠ", () -> randstring(rand(1:10))]
    fake_repl() do stdin_write, stdout_read, repl
        # In the future if we want we can add a test that the right object
        # gets displayed by intercepting the display
        repl.specialdisplay = REPL.REPLDisplay(repl)

        errormonitor(@async write(devnull, stdout_read)) # redirect stdout to devnull so we drain the output pipe

        repl.interface = REPL.setup_interface(repl)
        repl_mode = repl.interface.modes[1]
        shell_mode = repl.interface.modes[2]
        help_mode = repl.interface.modes[3]
        histp = repl.interface.modes[4]
        prefix_mode = repl.interface.modes[5]

        hp = REPL.REPLHistoryProvider(Dict{Symbol,Any}(:julia => repl_mode,
                                                       :shell => shell_mode,
                                                       :help  => help_mode))
        hist_path = tempname()
        write(hist_path, fakehistory)
        REPL.hist_from_file(hp, hist_path)
        f = open(hist_path, read=true, write=true, create=true)
        hp.history_file = f
        seekend(f)
        REPL.history_reset_state(hp)

        histp.hp = repl_mode.hist = shell_mode.hist = help_mode.hist = hp

        # Some manual setup
        s = LineEdit.init_state(repl.t, repl.interface)
        LineEdit.edit_insert(s, "wip")

        # LineEdit functions related to history
        LineEdit.edit_insert_last_word(s)
        @test buffercontents(LineEdit.buffer(s)) == "wip2"
        LineEdit.edit_backspace(s) # remove the "2"

        # Test that navigating history skips invalid modes
        # (in both directions)
        LineEdit.history_prev(s, hp)
        @test LineEdit.mode(s) == repl_mode
        @test buffercontents(LineEdit.buffer(s)) == "2 + 2"
        LineEdit.history_prev(s, hp)
        @test LineEdit.mode(s) == shell_mode
        @test buffercontents(LineEdit.buffer(s)) == "ls"
        LineEdit.history_prev(s, hp)
        @test LineEdit.mode(s) == repl_mode
        @test buffercontents(LineEdit.buffer(s)) == "1 + 1"
        LineEdit.history_next(s, hp)
        @test LineEdit.mode(s) == shell_mode
        @test buffercontents(LineEdit.buffer(s)) == "ls"
        LineEdit.history_next(s, hp)
        @test LineEdit.mode(s) == repl_mode
        @test buffercontents(LineEdit.buffer(s)) == "2 + 2"
        LineEdit.history_next(s, hp)
        @test LineEdit.mode(s) == repl_mode
        @test buffercontents(LineEdit.buffer(s)) == "wip"
        @test position(LineEdit.buffer(s)) == 3
        LineEdit.history_next(s, hp)
        @test buffercontents(LineEdit.buffer(s)) == "wip"
        LineEdit.history_prev(s, hp, 2)
        @test LineEdit.mode(s) == shell_mode
        @test buffercontents(LineEdit.buffer(s)) == "ls"
        LineEdit.history_prev(s, hp, -2) # equivalent to history_next(s, hp, 2)
        @test LineEdit.mode(s) == repl_mode
        @test buffercontents(LineEdit.buffer(s)) == "2 + 2"
        LineEdit.history_next(s, hp, -2) # equivalent to history_prev(s, hp, 2)
        @test LineEdit.mode(s) == shell_mode
        @test buffercontents(LineEdit.buffer(s)) == "ls"
        LineEdit.history_first(s, hp)
        @test LineEdit.mode(s) == repl_mode
        @test buffercontents(LineEdit.buffer(s)) == "é"
        LineEdit.history_next(s, hp, 6)
        @test LineEdit.mode(s) == shell_mode
        @test buffercontents(LineEdit.buffer(s)) == "ls"
        LineEdit.history_last(s, hp)
        @test buffercontents(LineEdit.buffer(s)) == "wip"
        @test position(LineEdit.buffer(s)) == 3
        # test that history_first jumps to beginning of current session's history
        hp.start_idx -= 5 # temporarily alter history
        LineEdit.history_first(s, hp)
        @test hp.cur_idx == 6
        # we are at the beginning of current session's history, so history_first
        # must now jump to the beginning of all history
        LineEdit.history_first(s, hp)
        @test hp.cur_idx == 1
        LineEdit.history_last(s, hp)
        @test hp.cur_idx-1 == length(hp.history)
        hp.start_idx += 5
        LineEdit.move_line_start(s)
        @test position(LineEdit.buffer(s)) == 0

        # Test that the same holds for prefix search
        ps = LineEdit.state(s, prefix_mode)::LineEdit.PrefixSearchState
        @test LineEdit.input_string(ps) == ""
        LineEdit.enter_prefix_search(s, prefix_mode, true)
        LineEdit.history_prev_prefix(ps, hp, "")
        @test ps.prefix == ""
        @test ps.parent == repl_mode
        @test LineEdit.input_string(ps) == "2 + 2"
        @test position(LineEdit.buffer(s)) == 5
        LineEdit.history_prev_prefix(ps, hp, "")
        @test ps.parent == shell_mode
        @test LineEdit.input_string(ps) == "ls"
        @test position(LineEdit.buffer(s)) == 2
        LineEdit.history_prev_prefix(ps, hp, "sh")
        @test ps.parent == repl_mode
        @test LineEdit.input_string(ps) == "shell"
        @test position(LineEdit.buffer(s)) == 2
        LineEdit.history_next_prefix(ps, hp, "sh")
        @test ps.parent == repl_mode
        @test LineEdit.input_string(ps) == "wip"
        @test position(LineEdit.buffer(s)) == 0
        LineEdit.move_input_end(s)
        LineEdit.history_prev_prefix(ps, hp, "é")
        @test ps.parent == repl_mode
        @test LineEdit.input_string(ps) == "éé"
        @test position(LineEdit.buffer(s)) == sizeof("é") > 1
        LineEdit.history_prev_prefix(ps, hp, "é")
        @test ps.parent == repl_mode
        @test LineEdit.input_string(ps) == "é"
        @test position(LineEdit.buffer(s)) == sizeof("é")
        LineEdit.history_next_prefix(ps, hp, "zzz")
        @test ps.parent == repl_mode
        @test LineEdit.input_string(ps) == "wip"
        @test position(LineEdit.buffer(s)) == 3
        LineEdit.accept_result(s, prefix_mode)

        # Test that searching backwards puts you into the correct mode and
        # skips invalid modes.
        LineEdit.enter_search(s, histp, true)
        ss = LineEdit.state(s, histp)
        write(ss.query_buffer, "l")
        LineEdit.update_display_buffer(ss, ss)
        LineEdit.accept_result(s, histp)
        @test LineEdit.mode(s) == shell_mode
        @test buffercontents(LineEdit.buffer(s)) == "ls"
        @test position(LineEdit.buffer(s)) == 0

        # Test that searching for `ll` actually matches `ll` after
        # both letters are types rather than jumping to `shell`
        LineEdit.history_prev(s, hp)
        LineEdit.enter_search(s, histp, true)
        write(ss.query_buffer, "l")
        LineEdit.update_display_buffer(ss, ss)
        @test buffercontents(ss.response_buffer) == "ll"
        @test position(ss.response_buffer) == 1
        write(ss.query_buffer, "l")
        LineEdit.update_display_buffer(ss, ss)
        LineEdit.accept_result(s, histp)
        @test LineEdit.mode(s) == shell_mode
        @test buffercontents(LineEdit.buffer(s)) == "ll"
        @test position(LineEdit.buffer(s)) == 0

        # Test that searching backwards with a one-letter query doesn't
        # return indefinitely the same match (#9352)
        LineEdit.enter_search(s, histp, true)
        write(ss.query_buffer, "l")
        LineEdit.update_display_buffer(ss, ss)
        LineEdit.history_next_result(s, ss)
        LineEdit.update_display_buffer(ss, ss)
        LineEdit.accept_result(s, histp)
        @test LineEdit.mode(s) == repl_mode
        @test buffercontents(LineEdit.buffer(s)) == "shell"
        @test position(LineEdit.buffer(s)) == 4

        # Test that searching backwards doesn't skip matches (#9352)
        # (for a search with multiple one-byte characters, or UTF-8 characters)
        LineEdit.enter_search(s, histp, true)
        write(ss.query_buffer, "é") # matches right-most "é" in "éé"
        LineEdit.update_display_buffer(ss, ss)
        @test position(ss.query_buffer) == sizeof("é")
        LineEdit.history_next_result(s, ss) # matches left-most "é" in "éé"
        LineEdit.update_display_buffer(ss, ss)
        LineEdit.accept_result(s, histp)
        @test buffercontents(LineEdit.buffer(s)) == "éé"
        @test position(LineEdit.buffer(s)) == 0

        # Issue #7551
        # Enter search mode and try accepting an empty result
        REPL.history_reset_state(hp)
        LineEdit.edit_clear(s)
        cur_mode = LineEdit.mode(s)
        LineEdit.enter_search(s, histp, true)
        LineEdit.accept_result(s, histp)
        @test LineEdit.mode(s) == cur_mode
        @test buffercontents(LineEdit.buffer(s)) == ""
        @test position(LineEdit.buffer(s)) == 0

        # Test that new modes can be dynamically added to the REPL and will
        # integrate nicely
        foobar_mode, custom_histp = AddCustomMode(repl, prompt)

        # ^R l, should now find `ls` in foobar mode
        LineEdit.enter_search(s, histp, true)
        ss = LineEdit.state(s, histp)
        write(ss.query_buffer, "l")
        LineEdit.update_display_buffer(ss, ss)
        LineEdit.accept_result(s, histp)
        @test LineEdit.mode(s) == foobar_mode
        @test buffercontents(LineEdit.buffer(s)) == "ls"
        @test position(LineEdit.buffer(s)) == 0

        # Try the same for prefix search
        LineEdit.history_next(s, hp)
        LineEdit.history_prev_prefix(ps, hp, "l")
        @test ps.parent == foobar_mode
        @test LineEdit.input_string(ps) == "ls"
        @test position(LineEdit.buffer(s)) == 1

        # Some Unicode handling testing
        LineEdit.history_prev(s, hp)
        LineEdit.enter_search(s, histp, true)
        write(ss.query_buffer, "x")
        LineEdit.update_display_buffer(ss, ss)
        @test buffercontents(ss.response_buffer) == "x ΔxΔ"
        @test position(ss.response_buffer) == 4
        write(ss.query_buffer, " ")
        LineEdit.update_display_buffer(ss, ss)
        LineEdit.accept_result(s, histp)
        @test LineEdit.mode(s) == repl_mode
        @test buffercontents(LineEdit.buffer(s)) == "x ΔxΔ"
        @test position(LineEdit.buffer(s)) == 0

        LineEdit.edit_clear(s)
        LineEdit.enter_search(s, histp, true)
        ss = LineEdit.state(s, histp)
        write(ss.query_buffer, "Å") # should not be in history
        LineEdit.update_display_buffer(ss, ss)
        @test buffercontents(ss.response_buffer) == ""
        @test position(ss.response_buffer) == 0
        LineEdit.history_next_result(s, ss) # should not throw BoundsError
        LineEdit.accept_result(s, histp)

        # Try entering search mode while in custom repl mode
        LineEdit.enter_search(s, custom_histp, true)
    end
end

# Test removal of prompt in bracket pasting
fake_repl() do stdin_write, stdout_read, repl
    repl.interface = REPL.setup_interface(repl)
    repl_mode = repl.interface.modes[1]
    shell_mode = repl.interface.modes[2]
    help_mode = repl.interface.modes[3]

    repltask = @async begin
        REPL.run_repl(repl)
    end

    global c = Condition()
    sendrepl2(cmd) = write(stdin_write, "$cmd\n notify($(curmod_prefix)c)\n")

    # Test removal of prefix in single statement paste
    sendrepl2("\e[200~julia> A = 2\e[201~\n")
    wait(c)
    @test Main.A == 2

    # Test removal of prefix in multiple statement paste
    sendrepl2("""\e[200~
            julia> mutable struct T17599; a::Int; end

            julia> function foo(julia)
            julia> 3
                end

                    julia> A = 3\e[201~
             """)
    wait(c)
    @test Main.A == 3
    @test Base.invokelatest(Main.foo, 4)
    @test Base.invokelatest(Main.T17599, 3).a == 3
    @test !Base.invokelatest(Main.foo, 2)

    sendrepl2("""\e[200~
            julia> goo(x) = x + 1
            error()

            julia> A = 4
            4\e[201~
             """)
    wait(c)
    @test Main.A == 4
    @test Base.invokelatest(Main.goo, 4) == 5

    # Test prefix removal only active in bracket paste mode
    sendrepl2("julia = 4\n julia> 3 && (A = 1)\n")
    wait(c)
    @test Main.A == 1

    # Test that indentation corresponding to the prompt is removed
    sendrepl2("""\e[200~julia> begin\n           α=1\n           β=2\n       end\n\e[201~""")
    wait(c)
    readuntil(stdout_read, "begin")
    @test readuntil(stdout_read, "end", keep=true) == "\n\r\e[7C    α=1\n\r\e[7C    β=2\n\r\e[7Cend"
    # for incomplete input (`end` below is added after the end of bracket paste)
    sendrepl2("""\e[200~julia> begin\n           α=1\n           β=2\n\e[201~end""")
    wait(c)
    readuntil(stdout_read, "begin")
    readuntil(stdout_read, "begin")
    @test readuntil(stdout_read, "end", keep=true) == "\n\r\e[7C    α=1\n\r\e[7C    β=2\n\r\e[7Cend"

    # Test switching repl modes
    sendrepl2("""\e[200~
            julia> A = 1
            1

            shell> echo foo
            foo

            shell> echo foo
                   foo
            foo foo

            help?> Int
            Dummy docstring

                Some text

                julia> error("If this error throws, the paste handler has failed to ignore this docstring example")

            julia> B = 2
            2\e[201~
             """)
    wait(c)
    @test Main.A == 1
    @test Main.B == 2

    # Close repl
    write(stdin_write, '\x04')
    Base.wait(repltask)
end

# Simple non-standard REPL tests
fake_repl() do stdin_write, stdout_read, repl
    panel = LineEdit.Prompt("testπ";
        prompt_prefix="\e[38;5;166m",
        prompt_suffix=Base.text_colors[:white],
        on_enter = s->true)

    hp = REPL.REPLHistoryProvider(Dict{Symbol,Any}(:parse => panel))
    search_prompt, skeymap = LineEdit.setup_prefix_keymap(hp, panel)
    REPL.history_reset_state(hp)

    panel.hist = hp
    panel.keymap_dict = LineEdit.keymap(Dict{Any,Any}[skeymap,
        LineEdit.default_keymap, LineEdit.escape_defaults])

    c = Condition()
    panel.on_done = (s, buf, ok) -> begin
        if !ok
            LineEdit.transition(s, :abort)
        end
        line = strip(String(take!(buf)))
        LineEdit.reset_state(s)
        notify(c, line)
        nothing
    end

    repltask = @async REPL.run_interface(repl.t, LineEdit.ModalInterface(Any[panel, search_prompt]))

    write(stdin_write,"a\n")
    @test wait(c) == "a"
    # Up arrow enter should recall history even at the start
    write(stdin_write,"\e[A\n")
    @test wait(c) == "a"
    # And again
    write(stdin_write,"\e[A\n")
    @test wait(c) == "a"
    # Close REPL ^D
    write(stdin_write, '\x04')
    Base.wait(repltask)
end

Base.exit_on_sigint(true)

let exename = `$(Base.julia_cmd()) --startup-file=no --color=no`
    # Test REPL in dumb mode
    with_fake_pty() do pts, ptm
        nENV = copy(ENV)
        nENV["TERM"] = "dumb"
        p = run(detach(setenv(`$exename -q`, nENV)), pts, pts, pts, wait=false)
        Base.close_stdio(pts)
        output = readuntil(ptm, "julia> ", keep=true)
        if ccall(:jl_running_on_valgrind, Cint,()) == 0
            # If --trace-children=yes is passed to valgrind, we will get a
            # valgrind banner here, not just the prompt.
            @test output == "julia> "
        end
        write(ptm, "1\nexit()\n")

        output = readuntil(ptm, ' ', keep=true)
        if Sys.iswindows()
        # Our fake pty is actually a pipe, and thus lacks the input echo feature of posix
            @test output == "1\n\njulia> "
        else
            @test output == "1\r\nexit()\r\n1\r\n\r\njulia> "
        end
        @test bytesavailable(ptm) == 0
        @test if Sys.iswindows() || Sys.isbsd()
                eof(ptm)
            else
                # Some platforms (such as linux) report EIO instead of EOF
                # possibly consume child-exited notification
                # for example, see discussion in https://bugs.python.org/issue5380
                try
                    eof(ptm) && !Sys.islinux()
                catch ex
                    (ex isa Base.IOError && ex.code == Base.UV_EIO) || rethrow()
                    @test_throws ex eof(ptm) # make sure the error is sticky
                    ptm.readerror = nothing
                    eof(ptm)
                end
            end
        @test read(ptm, String) == ""
        wait(p)
    end

    # Test stream mode
    p = open(`$exename -q`, "r+")
    write(p, "1\nexit()\n")
    @test read(p, String) == "1\n"
end # let exename

# issue #19864
mutable struct Error19864 <: Exception; end
function test19864()
    @eval Base.showerror(io::IO, e::Error19864) = print(io, "correct19864")
    buf = IOBuffer()
    fake_response = (Base.ExceptionStack([(exception=Error19864(),backtrace=Ptr{Cvoid}[])]),true)
    REPL.print_response(buf, fake_response, false, false, nothing)
    return String(take!(buf))
end
@test occursin("correct19864", test19864())

# Test containers in error messages are limited #18726
let io = IOBuffer()
    Base.display_error(io, Base.ExceptionStack(Any[(exception =
        (try
            [][trues(6000)]
            @assert false
        catch e
            e
        end), backtrace = [])]))
    @test length(String(take!(io))) < 1500
end

fake_repl() do stdin_write, stdout_read, repl
    # Relies on implementation detail to make sure we only have the single
    # replinit callback we want to test.
    saved_replinit = copy(Base.repl_hooks)
    slot = Ref(false)
    # Create a closure from a newer world to check if `_atreplinit`
    # can run it correctly
    atreplinit(@eval(repl::REPL.LineEditREPL -> ($slot[] = true)))
    Base._atreplinit(repl)
    @test slot[]
    @test_throws MethodError Base.repl_hooks[1](repl)
    copyto!(Base.repl_hooks, saved_replinit)
    nothing
end

let ends_with_semicolon = REPL.ends_with_semicolon
    @test !ends_with_semicolon("")
    @test ends_with_semicolon(";")
    @test !ends_with_semicolon("ä")
    @test !ends_with_semicolon("ä # äsdf ;")
    @test ends_with_semicolon("""a * "#ä" ;""")
    @test ends_with_semicolon("a; #=#=# =# =#\n")
    @test ends_with_semicolon("1;")
    @test ends_with_semicolon("1;\n")
    @test ends_with_semicolon("1;\r")
    @test ends_with_semicolon("1;\r\n   \t\f")
    @test ends_with_semicolon("1;#äsdf\n")
    @test ends_with_semicolon("""1;\n#äsdf\n""")
    @test !ends_with_semicolon("\"\\\";\"#\"")
    @test ends_with_semicolon("\"\\\\\";#\"")
    @test !ends_with_semicolon("begin\na;\nb;\nend")
    @test !ends_with_semicolon("begin\na; #=#=#\n=#b=#\nend")
    @test ends_with_semicolon("\na; #=#=#\n=#b=#\n# test\n#=\nfoobar\n=##bazbax\n")
    @test ends_with_semicolon("f()= 1; # é ; 2")
    @test ends_with_semicolon("f()= 1; # é")
    @test !ends_with_semicolon("f()= 1; \"é\"")
    @test !ends_with_semicolon("""("f()= 1; # é")""")
    @test !ends_with_semicolon(""" "f()= 1; # é" """)
    @test ends_with_semicolon("f()= 1;")
    # the next result does not matter because this is not legal syntax
    @test_nowarn ends_with_semicolon("1; #=# 2")
end

# PR #20794, TTYTerminal with other kinds of streams
let term = REPL.Terminals.TTYTerminal("dumb",IOBuffer("1+2\n"),IOContext(IOBuffer(),:foo=>true),IOBuffer())
    r = REPL.BasicREPL(term)
    REPL.run_repl(r)
    @test String(take!(term.out_stream.io)) == "julia> 3\n\njulia> \n"
    @test haskey(term, :foo) == true
    @test haskey(term, :bar) == false
    @test (:foo=>true) in term
    @test (:foo=>false) ∉ term
    @test term[:foo] == get(term, :foo, nothing) == true
    @test get(term, :bar, nothing) === nothing
    @test_throws KeyError term[:bar]
end

# Ensure even the dumb REPL elides content
let term = REPL.Terminals.TTYTerminal("dumb",IOBuffer("zeros(1000)\n"),IOBuffer(),IOBuffer())
    r = REPL.BasicREPL(term)
    REPL.run_repl(r)
    @test contains(String(take!(term.out_stream)), "⋮")
end

# a small module for alternative keymap tests
module AltLE
import REPL
import REPL.LineEdit

function history_move_prefix(s::LineEdit.MIState,
                             hist::REPL.REPLHistoryProvider,
                             backwards::Bool)
    buf = LineEdit.buffer(s)
    pos = position(buf)
    prefix = REPL.beforecursor(buf)
    allbuf = String(take!(copy(buf)))
    cur_idx = hist.cur_idx
    # when searching forward, start at last_idx
    if !backwards && hist.last_idx > 0
        cur_idx = hist.last_idx
    end
    hist.last_idx = -1
    idxs = backwards ? ((cur_idx-1):-1:1) : ((cur_idx+1):length(hist.history))
    for idx in idxs
        if startswith(hist.history[idx], prefix) && hist.history[idx] != allbuf
            REPL.history_move(s, hist, idx)
            seek(LineEdit.buffer(s), pos)
            LineEdit.refresh_line(s)
            return :ok
        end
    end
    REPL.Terminals.beep(LineEdit.terminal(s))
end
history_next_prefix(s::LineEdit.MIState, hist::REPL.REPLHistoryProvider) =
    history_move_prefix(s, hist, false)
history_prev_prefix(s::LineEdit.MIState, hist::REPL.REPLHistoryProvider) =
    history_move_prefix(s, hist, true)

end # module

# Test alternative keymaps and prompt
# (Alt. keymaps may be passed as a Vector{<:Dict} or as a Dict)

const altkeys = [Dict{Any,Any}("\e[A" => (s,o...)->(LineEdit.edit_move_up(s) || LineEdit.history_prev(s, LineEdit.mode(s).hist))), # Up Arrow
                 Dict{Any,Any}("\e[B" => (s,o...)->(LineEdit.edit_move_down(s) || LineEdit.history_next(s, LineEdit.mode(s).hist))), # Down Arrow
                 Dict{Any,Any}("\e[5~" => (s,o...)->(AltLE.history_prev_prefix(s, LineEdit.mode(s).hist))), # Page Up
                 Dict{Any,Any}("\e[6~" => (s,o...)->(AltLE.history_next_prefix(s, LineEdit.mode(s).hist))), # Page Down
                ]

for keys = [altkeys, merge(altkeys...)],
        altprompt = ["julia-$(VERSION.major).$(VERSION.minor)> ",
                     () -> "julia-$(Base.GIT_VERSION_INFO.commit_short)"]
    histfile = tempname()
    try
        fake_repl() do stdin_write, stdout_read, repl
            repl.specialdisplay = REPL.REPLDisplay(repl)
            repl.history_file = true
            withenv("JULIA_HISTORY" => histfile) do
                repl.interface = REPL.setup_interface(repl, extra_repl_keymap = altkeys)
            end
            repl.interface.modes[1].prompt = altprompt

            repltask = @async begin
                REPL.run_repl(repl)
            end

            sendrepl3(cmd) = write(stdin_write,"$cmd\n")

            sendrepl3("1 + 1;")                        # a simple line
            sendrepl3("multi=2;\e\nline=2;")           # a multiline input
            sendrepl3("ignoreme\e[A\b\b3;\e[B\b\b1;")  # edit the previous multiline input
            sendrepl3("1 +\e[5~\b*")                   # use prefix search to edit the 1st input

            # Close REPL ^D
            write(stdin_write, '\x04')
            Base.wait(repltask)

            # Close the history file
            # (otherwise trying to delete it fails on Windows)
            close(repl.interface.modes[1].hist.history_file)

            # Check that the correct prompt was displayed
            output = readuntil(stdout_read, "1 * 1;", keep=true)
            @test !occursin(output, LineEdit.prompt_string(altprompt))
            @test !occursin(output, "julia> ")

            # Check the history file
            history = read(histfile, String)
            @test occursin(r"""
                           ^\#\ time:\ .*\n
                            \#\ mode:\ julia\n
                            \t1\ \+\ 1;\n
                            \#\ time:\ .*\n
                            \#\ mode:\ julia\n
                            \tmulti=2;\n
                            \tline=2;\n
                            \#\ time:\ .*\n
                            \#\ mode:\ julia\n
                            \tmulti=3;\n
                            \tline=1;\n
                            \#\ time:\ .*\n
                            \#\ mode:\ julia\n
                            \t1\ \*\ 1;\n$
                           """xm, history)
        end
    finally
        rm(histfile, force=true)
    end
end

# Test that module prefix is omitted when type is reachable from Main (PR #23806)
fake_repl() do stdin_write, stdout_read, repl
    repl.specialdisplay = REPL.REPLDisplay(repl)
    repl.history_file = false

    repltask = @async begin
        REPL.run_repl(repl)
    end

    @eval Main module TestShowTypeREPL; export TypeA; struct TypeA end; end
    write(stdin_write, "TestShowTypeREPL.TypeA\n")
    @test endswith(readline(stdout_read), "\r\e[7CTestShowTypeREPL.TypeA\r\e[29C")
    readline(stdout_read)
    @test readline(stdout_read) == ""
    @eval Main using .TestShowTypeREPL
    readuntil(stdout_read, "julia> ", keep=true)
    write(stdin_write, "TypeA\n")
    @test endswith(readline(stdout_read), "\r\e[7CTypeA\r\e[12C")
    readline(stdout_read)
    @test readline(stdout_read) == ""

    # Close REPL ^D
    readuntil(stdout_read, "julia> ", keep=true)
    write(stdin_write, '\x04')
    Base.wait(repltask)
end

help_result(line) = Base.eval(REPL._helpmode(IOBuffer(), line))

# Docs.helpmode tests: we test whether the correct expressions are being generated here,
# rather than complete integration with Julia's REPL mode system.
for (line, expr) in Pair[
    "sin"          => :sin,
    "Base.sin"     => :(Base.sin),
    "@time(x)"     => Expr(:macrocall, Symbol("@time"), LineNumberNode(1, :none), :x),
    "@time"        => Expr(:macrocall, Symbol("@time"), LineNumberNode(1, :none)),
    ":@time"       => Expr(:quote, (Expr(:macrocall, Symbol("@time"), LineNumberNode(1, :none)))),
    "@time()"      => Expr(:macrocall, Symbol("@time"), LineNumberNode(1, :none)),
    "Base.@time()" => Expr(:macrocall, Expr(:., :Base, QuoteNode(Symbol("@time"))), LineNumberNode(1, :none)),
    "ccall"        => :ccall, # keyword
    "while       " => :while, # keyword, trailing spaces should be stripped.
    "0"            => 0,
    "\"...\""      => "...",
    "r\"...\""     => Expr(:macrocall, Symbol("@r_str"), LineNumberNode(1, :none), "..."),
    "using Foo"    => :using,
    "import Foo"   => :import,
    ]
    @test REPL._helpmode(line).args[4] == expr
    @test help_result(line) isa Union{Markdown.MD,Nothing}
end

# PR 30754, Issues #22013, #24871, #26933, #29282, #29361, #30348
for line in ["′", "type"]
    @test occursin("No documentation found.",
        sprint(show, help_result(line)::Union{Markdown.MD,Nothing}))
end

# PR 35154
@test occursin("|=", sprint(show, help_result("|=")))
@test occursin("broadcast", sprint(show, help_result(".=")))

# PR 35277
@test occursin("identical", sprint(show, help_result("===")))
@test occursin("broadcast", sprint(show, help_result(".<=")))

# Issue 39427
@test occursin("does not exist", sprint(show, help_result(":=")))

# Issue #40563
@test occursin("does not exist", sprint(show, help_result("..")))

# Issue #25930

# Brief and extended docs (issue #25930)
let text =
        """
            brief_extended()

        Short docs

        # Extended help

        Long docs
        """,
    md = Markdown.parse(text)
    @test md == REPL.trimdocs(md, false)
    @test !isa(md.content[end], REPL.Message)
    mdbrief = REPL.trimdocs(md, true)
    @test length(mdbrief.content) == 3
    @test isa(mdbrief.content[1], Markdown.Code)
    @test isa(mdbrief.content[2], Markdown.Paragraph)
    @test isa(mdbrief.content[3], REPL.Message)
    @test occursin("??", mdbrief.content[3].msg)
end

# issue #35216: empty and non-strings in H1 headers
let emptyH1 = Markdown.parse("# "),
    codeH1 = Markdown.parse("# `hello`")
    @test emptyH1 == REPL.trimdocs(emptyH1, false) == REPL.trimdocs(emptyH1, true)
    @test codeH1 == REPL.trimdocs(codeH1, false) == REPL.trimdocs(codeH1, true)
end

module BriefExtended
"""
    f()

Short docs

# Extended help

Long docs
"""
f() = nothing
@doc text"""
    f_plain()

Plain text docs
"""
f_plain() = nothing
@doc html"""
<h1><code>f_html()</code></h1>
<p>HTML docs.</p>
"""
f_html() = nothing
end # module BriefExtended

buf = IOBuffer()
md = Base.eval(REPL._helpmode(buf, "$(@__MODULE__).BriefExtended.f"))
@test length(md.content) == 2 && isa(md.content[2], REPL.Message)
buf = IOBuffer()
md = Base.eval(REPL._helpmode(buf, "?$(@__MODULE__).BriefExtended.f"))
@test length(md.content) == 1 && length(md.content[1].content[1].content) == 4
buf = IOBuffer()
txt = Base.eval(REPL._helpmode(buf, "$(@__MODULE__).BriefExtended.f_plain"))
@test !isempty(sprint(show, txt))
buf = IOBuffer()
html = Base.eval(REPL._helpmode(buf, "$(@__MODULE__).BriefExtended.f_html"))
@test !isempty(sprint(show, html))

# PR #27562
fake_repl() do stdin_write, stdout_read, repl
    repltask = @async begin
        REPL.run_repl(repl)
    end
    write(stdin_write, "Expr(:call, GlobalRef(Base.Math, :float), Core.SlotNumber(1))\n")
    readline(stdout_read)
    @test readline(stdout_read) == "\e[0m:(Base.Math.float(_1))"
    @test readline(stdout_read) == ""
    readuntil(stdout_read, "julia> ", keep=true)
    write(stdin_write, "ans\n")
    readline(stdout_read)
    @test readline(stdout_read) == "\e[0m:(Base.Math.float(_1))"
    @test readline(stdout_read) == ""
    readuntil(stdout_read, "julia> ", keep=true)
    write(stdin_write, '\x04')
    Base.wait(repltask)
end

# issue #31352
fake_repl() do stdin_write, stdout_read, repl
    repltask = @async begin
        REPL.run_repl(repl)
    end
    write(stdin_write, "struct Errs end\n")
    readline(stdout_read)
    readline(stdout_read)
    readuntil(stdout_read, "julia> ", keep=true)
    write(stdin_write, "Base.show(io::IO, ::Errs) = throw(Errs())\n")
    readline(stdout_read)
    readline(stdout_read)
    readuntil(stdout_read, "julia> ", keep=true)
    write(stdin_write, "Errs()\n")
    readline(stdout_read)
    readline(stdout_read)
    readuntil(stdout_read, "julia> ", keep=true)
    write(stdin_write, '\x04')
    wait(repltask)
    @test istaskdone(repltask)
end

# issue #34842
fake_repl() do stdin_write, stdout_read, repl
    repltask = @async begin
        REPL.run_repl(repl)
    end
    write(stdin_write, "?;\n")
    readline(stdout_read)
    @test endswith(readline(stdout_read), "search: ;")
    readuntil(stdout_read, "julia> ", keep=true)
    write(stdin_write, '\x04')
    Base.wait(repltask)
end

# issue #35771
fake_repl() do stdin_write, stdout_read, repl
    repltask = @async begin
        REPL.run_repl(repl)
    end
    write(stdin_write, "global x\n")
    readline(stdout_read)
    @test !occursin("ERROR", readline(stdout_read))
    readuntil(stdout_read, "julia> ", keep=true)
    write(stdin_write, '\x04')
    Base.wait(repltask)
end

fake_repl() do stdin_write, stdout_read, repl
    repltask = @async begin
        REPL.run_repl(repl)
    end
    write(stdin_write, "anything\x15\x19\x19") # ^u^y^y : kill line backwards + 2 yanks
    s1 = readuntil(stdout_read, "anything") # typed
    s2 = readuntil(stdout_read, "anything") # yanked (first ^y)
    s3 = readuntil(stdout_read, "anything") # previous yanked refreshed (from second ^y)
    s4 = readuntil(stdout_read, "anything", keep=true) # last yanked
    # necessary to read at least some part of the buffer,
    # for the "region_active" to have time to be updated

    @test LineEdit.state(repl.mistate).region_active == :off
    @test s4 == "anything" # no control characters between the last two occurrences of "anything"
    write(stdin_write, "\x15\x04")
    Base.wait(repltask)
end

# AST transformations (softscope, Revise, OhMyREPL, etc.)
@testset "AST Transformation" begin
    backend = REPL.REPLBackend()
    errormonitor(@async REPL.start_repl_backend(backend))
    put!(backend.repl_channel, (:(1+1), false))
    reply = take!(backend.response_channel)
    @test reply == Pair{Any, Bool}(2, false)
    twice(ex) = Expr(:tuple, ex, ex)
    push!(backend.ast_transforms, twice)
    put!(backend.repl_channel, (:(1+1), false))
    reply = take!(backend.response_channel)
    @test reply == Pair{Any, Bool}((2, 2), false)
    put!(backend.repl_channel, (nothing, -1))
    Base.wait(backend.backend_task)
end

backend = REPL.REPLBackend()
frontend_task = @async begin
    try
        @testset "AST Transformations Async" begin
            put!(backend.repl_channel, (:(1+1), false))
            reply = take!(backend.response_channel)
            @test reply == Pair{Any, Bool}(2, false)
            twice(ex) = Expr(:tuple, ex, ex)
            push!(backend.ast_transforms, twice)
            put!(backend.repl_channel, (:(1+1), false))
            reply = take!(backend.response_channel)
            @test reply == Pair{Any, Bool}((2, 2), false)
        end
    catch e
        Base.rethrow(e)
    finally
        put!(backend.repl_channel, (nothing, -1))
    end
end
REPL.start_repl_backend(backend)
Base.wait(frontend_task)

macro throw_with_linenumbernode(err)
    Expr(:block, LineNumberNode(42, Symbol("test.jl")), :(() -> throw($err)))
end

@testset "Install missing packages via hooks" begin
    @testset "Parse AST for packages" begin
        mods = REPL.modules_to_be_loaded(Base.parse_input_line("using Foo"))
        @test mods == [:Foo]
        mods = REPL.modules_to_be_loaded(Base.parse_input_line("import Foo"))
        @test mods == [:Foo]
        mods = REPL.modules_to_be_loaded(Base.parse_input_line("using Foo, Bar"))
        @test mods == [:Foo, :Bar]
        mods = REPL.modules_to_be_loaded(Base.parse_input_line("import Foo, Bar"))
        @test mods == [:Foo, :Bar]
        mods = REPL.modules_to_be_loaded(Base.parse_input_line("using Foo.bar, Foo.baz"))
        @test mods == [:Foo]

        mods = REPL.modules_to_be_loaded(Base.parse_input_line("if false using Foo end"))
        @test mods == [:Foo]
        mods = REPL.modules_to_be_loaded(Base.parse_input_line("if false if false using Foo end end"))
        @test mods == [:Foo]
        mods = REPL.modules_to_be_loaded(Base.parse_input_line("if false using Foo, Bar end"))
        @test mods == [:Foo, :Bar]
        mods = REPL.modules_to_be_loaded(Base.parse_input_line("if false using Foo: bar end"))
        @test mods == [:Foo]

        mods = REPL.modules_to_be_loaded(Base.parse_input_line("import Foo.bar as baz"))
        @test mods == [:Foo]
        mods = REPL.modules_to_be_loaded(Base.parse_input_line("using .Foo"))
        @test isempty(mods)
        mods = REPL.modules_to_be_loaded(Base.parse_input_line("using Base"))
        @test isempty(mods)
        mods = REPL.modules_to_be_loaded(Base.parse_input_line("using Base: nope"))
        @test isempty(mods)
        mods = REPL.modules_to_be_loaded(Base.parse_input_line("using Main"))
        @test isempty(mods)
        mods = REPL.modules_to_be_loaded(Base.parse_input_line("using Core"))
        @test isempty(mods)

        mods = REPL.modules_to_be_loaded(Base.parse_input_line(":(using Foo)"))
        @test isempty(mods)
        mods = REPL.modules_to_be_loaded(Base.parse_input_line("ex = :(using Foo)"))
        @test isempty(mods)

        mods = REPL.modules_to_be_loaded(Base.parse_input_line("Foo"))
        @test isempty(mods)

        mods = REPL.modules_to_be_loaded(Base.parse_input_line("@eval using Foo"))
        @test isempty(mods)
        mods = REPL.modules_to_be_loaded(Base.parse_input_line("begin using Foo; @eval using Bar end"))
        @test mods == [:Foo]
        mods = REPL.modules_to_be_loaded(Base.parse_input_line("Core.eval(Main,\"using Foo\")"))
        @test isempty(mods)
        mods = REPL.modules_to_be_loaded(Base.parse_input_line("begin using Foo; Core.eval(Main,\"using Foo\") end"))
        @test mods == [:Foo]
    end
end

# err should reprint error if deeper than top-level
fake_repl() do stdin_write, stdout_read, repl
    repltask = @async begin
        REPL.run_repl(repl)
    end
    # initialize `err` to `nothing`
    write(stdin_write, "global err = nothing\n")
    readline(stdout_read)
    readline(stdout_read) == "\e[0m"
    readuntil(stdout_read, "julia> ", keep=true)
    # generate top-level error
    write(stdin_write, "foobar\n")
    readline(stdout_read)
    @test readline(stdout_read) == "\e[0mERROR: UndefVarError: foobar not defined"
    @test readline(stdout_read) == ""
    readuntil(stdout_read, "julia> ", keep=true)
    # check that top-level error did not change `err`
    write(stdin_write, "err\n")
    readline(stdout_read)
    @test readline(stdout_read) == "\e[0m"
    readuntil(stdout_read, "julia> ", keep=true)
    # generate deeper error
    write(stdin_write, "foo() = foobar\n")
    readline(stdout_read)
    readuntil(stdout_read, "julia> ", keep=true)
    write(stdin_write, "foo()\n")
    readline(stdout_read)
    @test readline(stdout_read) == "\e[0mERROR: UndefVarError: foobar not defined"
    readuntil(stdout_read, "julia> ", keep=true)
    # check that deeper error did set `err`
    write(stdin_write, "err\n")
    readline(stdout_read)
    @test readline(stdout_read) == "\e[0m1-element ExceptionStack:"
    @test readline(stdout_read) == "UndefVarError: foobar not defined"
    @test readline(stdout_read) == "Stacktrace:"
    write(stdin_write, '\x04')
    Base.wait(repltask)
end

max_stars_repo_name: BSnelling/julia

max_issues_repo_name: BSnelling/julia

max_forks_repo_name: BSnelling/julia


text:
import ..Filter
import ..Wirreducible  #remove all polarization subdiagrams
import ..Girreducible  #remove all self-energy inseration
import ..NoHatree
import ..NoFock
import ..NoBubble  # true to remove all bubble subdiagram
import ..Proper  #ver4, ver3, and polarization diagrams may require to be irreducible along the transfer momentum/frequency

import ..DiagramType
import ..GreenDiag
import ..SigmaDiag
import ..PolarDiag
import ..Ver3Diag
import ..Ver4Diag

import ..Composite
import ..ChargeCharge
import ..SpinSpin
import ..UpUp
import ..UpDown
import ..Response

import ..Instant
import ..Dynamic
import ..D_Instant
import ..D_Dynamic
import ..AnalyticProperty

import ..symbol
import ..short

import ..Interaction
import ..GenericPara
import ..innerTauNum

import ..Diagram

import ..DiagramId
import ..Ver4Id
import ..Ver3Id
import ..GreenId
import ..SigmaId
import ..PolarId
import ..BareInteractionId
import ..BareGreenId

import ..TwoBodyChannel
import ..Alli
import ..PHr
import ..PHEr
import ..PPr
import ..AnyChan

import ..Permutation
import ..Di
import ..Ex
import ..DiEx

import ..uidreset
import ..toDataFrame
import ..mergeby

function build(para::GenericPara, extK = nothing, subdiagram = false)
    if para.diagType == Ver4Diag
        if isnothing(extK)
            extK = [DiagTree.getK(para.totalLoopNum, 1), DiagTree.getK(para.totalLoopNum, 2), DiagTree.getK(para.totalLoopNum, 3)]
        end
        return vertex4(para, extK, [PHr, PHEr, PPr, Alli], subdiagram)
    elseif para.diagType == SigmaDiag
        if isnothing(extK)
            extK = DiagTree.getK(para.totalLoopNum, 1)
        end
        return sigma(para, extK, subdiagram)
    elseif para.diagType == PolarDiag
        if isnothing(extK)
            extK = DiagTree.getK(para.totalLoopNum, 1)
        end
        return polarization(para, extK, subdiagram)
    elseif para.diagType == Ver3Diag
        if isnothing(extK)
            extK = [DiagTree.getK(para.totalLoopNum, 1), DiagTree.getK(para.totalLoopNum, 2)]
        end
        return vertex3(para, extK, subdiagram)
    else
        error("not implemented!")
    end
end

"""
    struct ParquetBlocks

    The channels of the left and right sub-vertex4 of a bubble diagram in the parquet equation

#Members
- `phi`   : channels of left sub-vertex for the particle-hole and particle-hole-exchange bubbles
- `ppi`   : channels of left sub-vertex for the particle-particle bubble
- `Γ4`   : channels of right sub-vertex of all channels
"""
struct ParquetBlocks
    phi::Vector{TwoBodyChannel}
    ppi::Vector{TwoBodyChannel}
    Γ4::Vector{TwoBodyChannel}
    function ParquetBlocks(; phi = [Alli, PHEr, PPr], ppi = [Alli, PHr, PHEr], Γ4 = union(phi, ppi))
        return new(phi, ppi, Γ4)
    end
end

function Base.isequal(a::ParquetBlocks, b::ParquetBlocks)
    if issetequal(a.phi, b.phi) && issetequal(a.ppi, b.ppi) && issetequal(a.Γ4, b.Γ4)
        return true
    else
        return false
    end
end
Base.:(==)(a::ParquetBlocks, b::ParquetBlocks) = Base.isequal(a, b)

function orderedPartition(_total, n, lowerbound = 1)
    @assert lowerbound >= 0
    total = _total - n * (lowerbound - 1)
    @assert total >= n
    unorderedPartition = collect(partitions(total, n))
    #e.g., loopNum =5, n =2 ==> unordered = [[4, 1], [3, 2]]
    orderedPartition = Vector{Vector{Int}}([])
    for p in unorderedPartition
        p = p .+ (lowerbound - 1)
        @assert sum(p) == _total
        for i in p
            @assert i >= lowerbound
        end
        append!(orderedPartition, Set(permutations(p)))
    end
    #e.g., loopNum =5, n =2 ==> ordered = [[4, 1], [1, 4], [3, 2], [2, 3]]
    return orderedPartition
end

function findFirstLoopIdx(partition, firstidx::Int)
    ## example: firstidx = 1
    # partition = [1, 1, 2, 1], then the loop partition = [1][2][34][5], thus firstTauIdx = [1, 2, 3, 5]
    # partition = [1, 0, 2, 0], then the loop partition = [1][][23][], thus firstTauIdx = [1, 2, 2, 4]
    # @assert length(partition) == length(isG)
    accumulated = accumulate(+, partition; init = firstidx) #  idx[i] = firstidx + p[1]+p[2]+...+p[i]
    firstLoopIdx = [firstidx,]
    append!(firstLoopIdx, accumulated[1:end-1])
    maxLoopIdx = accumulated[end] - 1
    return firstLoopIdx, maxLoopIdx
end

function findFirstTauIdx(partition::Vector{Int}, diagType::Vector{DiagramType}, firstidx::Int, _tauNum::Int)
    ## example: diagType =[Vertex4, GreenDiagram, Vertex4, GreenDiagram], firstidx = 1
    # n-loop G has n*_tauNum DOF, while n-loop ver4 has (n+1)*_tauNum DOF
    # partition = [1, 1, 2, 1], then the tau partition = [12][3][456][7], thus firstTauIdx = [1, 3, 4, 7]
    # partition = [1, 0, 2, 0], then the tau partition = [12][][345][], thus firstTauIdx = [1, 3, 3, 6]
    @assert length(partition) == length(diagType)
    @assert _tauNum >= 0
    taupartition = [innerTauNum(diagType[i], p, _tauNum) for (i, p) in enumerate(partition)]
    accumulated = accumulate(+, taupartition; init = firstidx) #  idx[i] = firstidx + p[1]+p[2]+...+p[i]
    firstTauidx = [firstidx,]
    append!(firstTauidx, accumulated[1:end-1])
    maxTauIdx = accumulated[end] - 1
    return firstTauidx, maxTauIdx
end

function allsame(df, name::Symbol)
    @assert all(x -> x == df[1, name], df[!, name]) "Not all rows of the $name field are the same.\n$df"
end
function allsame(df, names::Vector{Symbol})
    for name in names
        allsame(df, name)
    end
end
function allsametype(df, name::Symbol)
    @assert all(x -> typeof(x) == typeof(df[1, name]), df[!, name]) "Not all rows of the $name field are the same type.\n$df"
end
function allsametype(df, names::Vector{Symbol})
    for name in names
        allsametype(df, name)
    end
end


text:
using Pkg

Pkg.add("MultivariateStats")
Pkg.add("Arpack")
Pkg.add("StatsBase")
Pkg.add("CSV")
Pkg.add("DelimitedFiles")
Pkg.add("Statistics")
Pkg.add(Pkg.PackageSpec(;url="https://github.com/rikenbit/OnlinePCA.jl", version="0.2.5"))

max_stars_repo_name: rikenbit/onlinePCA-experiments

max_issues_repo_name: rikenbit/onlinePCA-experiments

max_forks_repo_name: rikenbit/onlinePCA-experiments


text:
# Geometry for vectors/arrows/segments

struct SegmentGeometry <: Gadfly.GeometryElement
    default_statistic::Gadfly.StatisticElement
     arrow::Bool
     filled::Bool
     tag::Symbol
end
SegmentGeometry(default_statistic=Gadfly.Stat.identity(); arrow=false, filled=false, tag=empty_tag) =
    SegmentGeometry(default_statistic, arrow, filled, tag)

const segment = SegmentGeometry

# Leave this as a function, pending extra arguments e.g. arrow attributes
vector(; filled::Bool=false) = SegmentGeometry(arrow=true, filled=filled)

hair(;intercept=0.0, orientation=:vertical) =
    SegmentGeometry(Gadfly.Stat.hair(intercept, orientation))

function vectorfield(;smoothness=1.0, scale=1.0, samples=20, filled::Bool=false)
    return SegmentGeometry(
        Gadfly.Stat.vectorfield(smoothness, scale, samples),
        arrow=true, filled=filled )
end

default_statistic(geom::SegmentGeometry) = geom.default_statistic
element_aesthetics(::SegmentGeometry) = [:x, :y, :xend, :yend, :color]

function render(geom::SegmentGeometry, theme::Gadfly.Theme, aes::Gadfly.Aesthetics)

    Gadfly.assert_aesthetics_defined("Geom.segment", aes, :x, :y, :xend, :yend)

    function arrow(x::T, y::T, xmax::T, ymax::T, xyrange::Vector{T}) where T<:Real
        dx = xmax-x
        dy = ymax-y
        vl = 0.225*hypot(dy/xyrange[2], dx/xyrange[1])
        θ =  atan2(dy/xyrange[2], dx/xyrange[1])
        ϕ = pi/15
        xr = -vl*xyrange[1]*[cos(θ+ϕ), cos(θ-ϕ)]
        yr = -vl*xyrange[2]*[sin(θ+ϕ), sin(θ-ϕ)]
        arr = [(xmax+xr[1],ymax+yr[1]), (xmax,ymax), (xmax+xr[2],ymax+yr[2]) ]
        return arr
    end

    n = length(aes.x)
    default_aes = Gadfly.Aesthetics()
    default_aes.color = fill(RGBA{Float32}(theme.default_color), n)

    aes = inherit(aes, default_aes)

    line_style = Gadfly.get_stroke_vector(theme.line_style)

    # Geom.vector requires information about scales

    if geom.arrow
        check = [aes.xviewmin, aes.xviewmax, aes.yviewmin, aes.yviewmax ]
        if any( map(x -> x === nothing, check) )
            error("For Geom.vector, Scale minvalue and maxvalue must be manually provided for both axes")
        end
         xyrange = [aes.xviewmax-aes.xviewmin, aes.yviewmax-aes.yviewmin]

         arrows = [ arrow(x, y, xend, yend, xyrange)
                for (x, y, xend, yend) in zip(aes.x, aes.y, aes.xend, aes.yend) ]
    end

    segments = [ [(x,y), (xend,yend)]
        for (x, y, xend, yend) in zip(aes.x, aes.y, aes.xend, aes.yend) ]

    classes = [svg_color_class_from_label( aes.color_label([c])[1] ) for c in aes.color ]

    ctx = context()

    compose!( ctx, (context(), Compose.line(segments, geom.tag),
                stroke(aes.color), linewidth(theme.line_width),
                strokedash(line_style), svgclass( classes )),
              svgclass("geometry")  )
    if geom.arrow
        if geom.filled
            compose!(ctx, (context(), Compose.polygon(arrows), fill(aes.color), strokedash([])) )
        else
            compose!(ctx, (context(), Compose.line(arrows), stroke(aes.color), linewidth(theme.line_width),
            strokedash([]))  )
        end
    end

    return ctx
end


text:
function barnes(x)

    a1 = 75.196
    a3 = 0.12694
    a5 = 1.0345e-5
    a7 = 0.030234
    a9 = 3.5256e-5
    a11 = 0.25645
    a13 = 1.3514e-5
    a15 = -5.2375e-6
    a17 = 7.0e-10
    a19 = -1.6638e-6
    a21 = 0.0005
    a2 = -3.8112
    a4 = -2.0567e-3
    a6 = -6.8306
    a8 = -1.28134e-3
    a10 = -2.266e-7
    a12 = -3.4604e-3
    a14 = -28.106
    a16 = -6.3e-8
    a18 = 3.4054e-4
    a20 = -2.8673

    x1 = x[1]
    x2 = x[2]
    y1 = x1*x2
    y2 = y1*x1
    y3 = x2^2
    y4 = x1^2

    # --- function value ---

    f = a1 + a2*x1 + a3*y4 + a4*y4*x1 + a5*y4^2 +
        a6*x2 + a7*y1 + a8*x1*y1 + a9*y1*y4 + a10*y2*y4 +
        a11*y3 + a12*x2*y3 + a13*y3^2 + a14/(x2+1) +
        a15*y3*y4 + a16*y1*y4*x2 + a17*y1*y3*y4 + a18*x1*y3 +
        a19*y1*y3 + a20*exp(a21*y1)

    # --- constraints ---

    c = zeros(3)
    c[1] = 1 - y1/700.0
    c[2] = y4/25.0^2 - x2/5.0
    c[3] = (x1/500.0- 0.11) - (x2/50.0-1)^2

    # --- derivatives of f ---

    dy1 = x2
    dy2 = y1 + x1*dy1
    dy4 = 2*x1
    dfdx1 = a2 + a3*dy4 + a4*y4 + a4*x1*dy4 + a5*2*y4*dy4 +
        a7*dy1 + a8*y1 + a8*x1*dy1 + a9*y1*dy4 + a9*y4*dy1 + a10*y2*dy4 + a10*y4*dy2 +
        a15*y3*dy4 + a16*x2*y1*dy4 + a16*x2*y4*dy1 + a17*y3*y1*dy4 + a17*y3*y4*dy1 + a18*y3 +
        a19*y3*dy1 + a20*exp(a21*y1)*a21*dy1

    dy1 = x1
    dy2 = x1*dy1
    dy3 = 2*x2
    dfdx2 = a6 + a7*dy1 + a8*x1*dy1 + a9*y4*dy1 + a10*y4*dy2 +
        a11*dy3 + a12*x2*dy3 + a12*y3 + a13*2*y3*dy3 + a14*-1/(x2+1)^2 +
        a15*y4*dy3 + a16*y4*y1 + a16*y4*x2*dy1 + a17*y4*y1*dy3 + a17*y4*y3*dy1 + a18*x1*dy3 +
        a19*y3*dy1 + a19*y1*dy3 + a20*exp(a21*y1)*a21*dy1

    dfdx = [dfdx1, dfdx2]

    # --- derivatives of c ---

    dcdx = zeros(3, 2)
    dcdx[1, 1] = -x2/700.0
    dcdx[1, 2] = -x1/700.0
    dcdx[2, 1] = 2*x1/25^2
    dcdx[2, 2] = -1.0/5
    dcdx[3, 1] = 1.0/500
    dcdx[3, 2] = -2*(x2/50.0-1)/50.0

    dcdx = dcdx'  # matlab format

    return f/30.0, c, Float64[], dfdx/30.0, dcdx, Float64[]
end