text:
# This file contains common definitions for the three RLWE based SHE
# schemes we support (BGV, BFV, CKKS). See the cryptographic background section
# in the manual for an overview.

################################################################################
#                        SHEShemeParams
################################################################################

abstract type SHEShemeParams end

# These four functions (π, π⁻¹, 𝒩, 𝒢) determine the differences between the
# three schemes we support.
function π end
function π⁻¹ end
function 𝒩 end
function 𝒢 end

"""
    ℛ_cipher(params::SHEShemeParams)

Given the parameters of the SHE scheme, return the ciphertext ring.
"""
function ℛ_cipher end
const ciphertext_space = ℛ_cipher

ℛ_key(p::SHEShemeParams) = ℛ_cipher(p)

"""
    ℛ_plain(params::SHEShemeParams)

Given the parameters of the SHE scheme, return the plaintext ring.
"""
function ℛ_plain end
const plaintext_space = ℛ_plain

relin_window(params::SHEShemeParams) = params.relin_window

# These optional function change multiplication. By default they do nothing
mul_expand(params::SHEShemeParams, x) = x
mul_contract(params::SHEShemeParams, x) = x

# Some utilities
Base.show(io::IO, params::SHEShemeParams) = print(io, scheme_name(params), " parameters")
scheme_name(params::SHEShemeParams) = scheme_name(typeof(params))
Broadcast.broadcastable(params::SHEShemeParams) = Ref(params)

################################################################################
#                        PassthroughParams
################################################################################

abstract type PassthroughParams{P} <: SHEShemeParams; end
parent_params(p::PassthroughParams) = p.params
ℛ_plain(p::PassthroughParams) = ℛ_plain(parent_params(p))
ℛ_cipher(p::PassthroughParams) = ℛ_cipher(parent_params(p))
ℛ_key(p::PassthroughParams) = ℛ_key(parent_params(p))
relin_window(p::PassthroughParams) = relin_window(parent_params(p))
π⁻¹(p::PassthroughParams, plaintext) = π⁻¹(parent_params(p), plaintext)
π(p::PassthroughParams, b) = π(parent_params(p), b)
𝒩(p::PassthroughParams) = 𝒩(parent_params(p))
𝒢(p::PassthroughParams) = 𝒢(parent_params(p))

################################################################################
#                        The various FHE key types
################################################################################

struct PrivKey{P <: SHEShemeParams}
    params::P
    secret
end
Base.show(io::IO, kp::PrivKey{P}) where {P} = print(io, scheme_name(P), " private key")

struct KeyComponent
    mask
    masked
end

struct PubKey{P <: SHEShemeParams}
    params::P
    # One can think of a public key as a keyswitching key from 0 to the secret
    # secret key of the scheme. Of course, for 0, we have no noise growth
    # concerns, so we can get away with a single KeyComponent
    key::KeyComponent
end
Base.show(io::IO, kp::PubKey{P}) where {P} = print(io, scheme_name(P), " public key")

struct KeySwitchKey{P <: SHEShemeParams}
    params::P
    key::Vector{KeyComponent}
end
Base.show(io::IO, kp::KeySwitchKey{P}) where {P} = print(io, scheme_name(P), " key-switching key")

abstract type EvalKey end
struct EvalMultKey{P <: SHEShemeParams} <: EvalKey
    key::KeySwitchKey{P}
end
Base.show(io::IO, gk::EvalMultKey{P}) where {P} = print(io, scheme_name(P), " multiplication key")

struct GaloisKey{P <: SHEShemeParams} <: EvalKey
    galois_element::Int
    key::KeySwitchKey{P}
end
Base.show(io::IO, gk::GaloisKey{P}) where {P} = print(io, scheme_name(P), " galois key (element ", gk.galois_element, ")")

struct GaloisKeys{P <: SHEShemeParams}
    # Sorted collection of galois keys for various step sizes
    keys::Vector{GaloisKey{P}}
end
Base.show(io::IO, gk::GaloisKeys{P}) where {P} = print(io, scheme_name(P), " galois keys (elements ", join(map(k->k.galois_element, gk.keys), ", ", ")"))

struct KeyPair{P <: SHEShemeParams}
    priv::PrivKey{P}
    pub::PubKey{P}
end
Base.show(io::IO, kp::KeyPair{P}) where {P} = print(io, scheme_name(P), " key pair")

Base.broadcastable(k::Union{PrivKey, PubKey, EvalKey, KeyPair}) = Ref(k)

################################################################################
#                        CipherText
################################################################################

"""
    CipherText{P <: SHEShemeParams, Enc, T, N}

The CipherText for an RLWE-based FHE scheme (the scheme is determined by the
parameter `p`). Optionally, a corresponding plaintext encoding `Enc` may be
specified that will be applied on decryption (`Any` to obtain the raw
decryption result in the scheme's native plaintext space).
"""
struct CipherText{Plain, P <: SHEShemeParams, T, N}
    params::P
    cs::NTuple{N, T}
end
CipherText(params::P, cs::NTuple{N,T}) where {P <: SHEShemeParams,T,N} =
    CipherText{Any,P,T,N}(params, cs)
CipherText{Plain}(params::P, cs::NTuple{N,T}) where {Plain, P <: SHEShemeParams,T,N} =
    CipherText{Plain,P,T,N}(params, cs)
ring(c::CipherText) = ring(c.cs[1])

Base.length(c::CipherText) = length(c.cs)
Base.getindex(c::CipherText, i::Integer) = c.cs[i]
Base.lastindex(c::CipherText) = length(c)
Base.broadcastable(c::CipherText) = Ref(c)
function Base.show(io::IO, kp::CipherText{Enc, P, <:Any, N}) where {P, Enc, N}
    print(io, scheme_name(P), " ciphertext (length ", N)
    Enc != Any && print(io, ", encoding $Enc")
    print(io, ")")
end

################################################################################
#                        Key generation
################################################################################

function keygen(rng::AbstractRNG, params::SHEShemeParams)
    𝒰 = RingSampler(ℛ_key(params), DiscreteUniform(coefftype(ℛ_key(params))))

    mask = rand(rng, 𝒰)
    secret = rand(rng, 𝒢(params))
    error = rand(rng, 𝒩(params))

    masked = -(mask*secret + error)

    KeyPair(
        PrivKey(params, secret),
        PubKey(params, KeyComponent(mask, masked)))
end

# TODO: CSPRNG here
keygen(params::SHEShemeParams) = keygen(Random.GLOBAL_RNG, params)

################################################################################
#                        encryption/decryption
################################################################################

function encrypt(rng::AbstractRNG, pk::PubKey, ::Zero)
    @fields_as_locals pk::PubKey
    @fields_as_locals key::KeyComponent

    u = rand(rng, 𝒢(params))
    e₁, e₂ = rand(rng, 𝒩(params), 2)

    c₁ = masked*u + e₁
    c₂ = mask*u + e₂
    CipherText{Zero}(params, (c₁, c₂))
end

function encrypt(rng::AbstractRNG, key::PubKey, plaintext)
    c = encrypt(rng, key, Zero())
    c += oftype(c.cs[1], π⁻¹(key.params, plaintext))

    EncT = typeof(plaintext)
    EncT <: NTT.RingElement && (EncT = Any)
    return CipherText{EncT}(key.params, c.cs)
end
encrypt(rng::AbstractRNG, kp::KeyPair, plaintext) = encrypt(rng, kp.pub, plaintext)
encrypt(key::Union{KeyPair, PubKey}, plaintext) = encrypt(Random.GLOBAL_RNG, key, plaintext)

function decrypt(key::PrivKey, c::CipherText{T}) where T
    @fields_as_locals key::PrivKey

    while eltype(secret) != eltype(c[1])
        secret = modswitch_drop(secret)
    end

    b = c[1]
    spow = secret

    for i = 2:length(c)
        b += spow*c[i]
        spow *= secret
    end

    dec = π(params, b)
    T === Any ? dec : T(dec)
end
decrypt(key::KeyPair, plaintext) = decrypt(key.priv, plaintext)

################################################################################
#                            error handling
################################################################################

struct UsageError
    msg::AbstractString
end

################################################################################
#                        Homomorphic arithmetic
################################################################################

for f in (:+, :-)
    @eval function $f(c1::CipherText{P, Enc, T, N1}, c2::CipherText{P, Enc, T,N2}) where {P, Enc, T, N1, N2}
        if c1.params !== c2.params
            throw(UsageError("Attempting to add ciphertexts with differing parameters"))
        end
        CipherText{P, Enc, T, max(N1, N2)}(c1.params, tuple((
            i > length(c1) ? c2[i] :
            i > length(c2) ? c1[i] :
            $f(c1[i], c2[i]) for i in 1:max(N1, N2))...))
    end
end

function +(c::CipherText{Enc, P, T, N}, b::T) where {T, Enc, P, N}
    CipherText{Enc, P, T, N}(c.params, (c.cs[1] + b, c.cs[2:end]...))
end

function enc_mul(c1, c2)
    if c1.params !== c2.params
        throw(UsageError("Attempting to multiply ciphertexts with differing parameters"))
    end
    params = c1.params

    (c1, c2) = mul_expand.(params, (c1, c2))

    c = [zero(c1[1]) for i = 1:(length(c1) + length(c2) - 1)]
    for i = 1:length(c1), j = 1:length(c2)
        c[i+j-1] += c1[i] * c2[j]
    end

    c = mul_contract(params, c)
    (c...,)
end

function *(c1::CipherText{P, Enc, T}, c2::CipherText{P, Enc, T}) where {P, Enc, T}
    CipherText(c1.params, enc_mul(c1, c2))
end

################################################################################
#                        Key switching
################################################################################
default_relin_window(ℛ) = 1

function make_eval_key(rng::AbstractRNG, (old, new)::Pair{<:Any, <:PrivKey})
    @fields_as_locals new::PrivKey

    ℛ = ring(old)

    𝒰 = RingSampler(ℛ, DiscreteUniform(coefftype(ℛ)))
    𝒩gen = 𝒩(params)

    if relin_window(params) != 0
        nwindows = ndigits(modulus(coefftype(ℛ)), base=2^relin_window(params))
        evala = [old * coefftype(ℛ)(2)^(i*relin_window(params)) for i = 0:nwindows-1]
        evalb = eltype(evala)[]
    else
        # CRT basis decomposition
        evala = [typeof(old)(OffsetArray(map(x->eltype(ℛ)(CRTResidual(x)), arr), axes(old)), nothing) for arr in StructArrays.fieldarrays(parent(NTT.coeffs_primal(old)))]
        evalb = eltype(evala)[]
    end

    for i = 1:length(evala)
        mask = rand(rng, 𝒰)
        e = rand(rng, 𝒩gen)
        push!(evalb, mask)
        evala[i] -= mask*new.secret + e
    end
    KeySwitchKey(new.params, map(KeyComponent, evalb, evala))
end
keygen(rng::AbstractRNG, ::Type{EvalMultKey}, priv::PrivKey) = EvalMultKey(make_eval_key(rng, priv.secret^2=>priv))
function keygen(rng::AbstractRNG, ::Type{GaloisKey}, priv::PrivKey; galois_element = nothing, steps=nothing)
    @assert galois_element === nothing || steps === nothing && !(galois_element === nothing && steps === nothing)
    if galois_element === nothing
        @assert steps !== nothing
        galois_element = mod(steps > 0 ? 3^(2degree(priv.secret)-steps) : 3^-steps, 2degree(priv.secret))
    else
        @assert steps === nothing
    end
    GaloisKey(galois_element, make_eval_key(rng, NTT.apply_galois_element(priv.secret, galois_element)=>priv))
end
keygen(T::Type{<:EvalKey}, priv::PrivKey; kwargs...) = keygen(Random.GLOBAL_RNG, T, priv; kwargs...)

keyswitch_expand(ek, re) = re
keyswitch_contract(ek, c) = c
downswitch_keyelement(params, key::KeyComponent, elt) = key
function keyswitch(ek::KeySwitchKey, c::CipherText{Enc}) where {Enc}
    @fields_as_locals ek::KeySwitchKey

    @assert length(c.cs) in (2,3)
    ℛ = ring(c.cs[1])

    any_mask = ek.key[1].mask

    c1 = keyswitch_expand(ek, c[1])
    c2 = length(c) == 2 ? zero(c1) : keyswitch_expand(ek, c[2])

    if relin_window(params) == 0
        # CRT basis decomposition
        cendcoeffs = NTT.coeffs_primal(c[end])
        ps = [typeof(c1)(eltype(c1).(convert.(Integer, SignedMod.(OffsetArray(arr, axes(cendcoeffs))))), nothing) for arr in StructArrays.fieldarrays(parent(cendcoeffs))]
    else
        # digit-wise decomposition
        cendcoeffs = NTT.coeffs_primal(c[end])
        nwindows = ndigits(modulus(coefftype(ℛ)), base=2^relin_window(params))
        ds = Any[digits(convert(Integer, x), base=2^relin_window(params), pad=nwindows) for x in cendcoeffs]
        ps = map(1:nwindows) do i
            typeof(c1)([coefftype(ℛ)(ds[j][i]) for j in eachindex(cendcoeffs)], nothing)
        end
    end

    for i in eachindex(ps)
        key = downswitch_keyelement(ek.params, ek.key[i], typeof(ps[i]))
        c2 += key.mask * ps[i]
        c1 += key.masked * ps[i]
    end

    CipherText{Enc}(c.params, (keyswitch_contract(ek, c1), keyswitch_contract(ek, c2)))
end
keyswitch(k::GaloisKey, c::CipherText) = keyswitch(k.key, c)
keyswitch(k::EvalMultKey, c::CipherText) = keyswitch(k.key, c)

################################################################################
#     Rotations
################################################################################

function NTT.apply_galois_element(c::CipherText{Enc}, galois_element) where {Enc}
    CipherText{Enc}(c.params, map(re->NTT.apply_galois_element(re, galois_element), c.cs))
end

rotate(gk::GaloisKey, c::CipherText) = keyswitch(gk, NTT.apply_galois_element(c, gk.galois_element))
Base.circshift(c::CipherText, gk::GaloisKey) = rotate(gk, c)

################################################################################
#     Modulus switching
################################################################################

function modswitch(c::CipherText{Any}, new_modulus)
    error("Not implemented")
end

################################################################################
#     Plaintext spaces for a given ciphertext space an plain modulus
################################################################################

function plaintext_space(r::ResRing, p)
    ℤp = ResidueRing(Nemo.ZZ, p)
    ℤpx = PolynomialRing(ℤp, "x")[1]
    ResidueRing(ℤpx, Nemo.lift(Nemo.ZZ["x"][1], modulus(r)))
end

function plaintext_space(r::NegacyclicRing, p)
    coefft = Primes.isprime(p) ? GaloisField(p) :
        p == 256 ? UInt8 :
        Mod(p)
    if Primes.isprime(p) && p > 2degree(modulus(r))
        # TODO: Also needs to check here if the prime admits 2n-th roots of
        # unities.
        NegacyclicRing{coefft, degree(modulus(r))}(
            GaloisFields.minimal_primitive_root(coefft, 2degree(modulus(r))))
    else
        NegacyclicRing{coefft, degree(modulus(r))}()
    end
end


text:
using ITensors
using Symbolics

using ITensors: contraction_cost, optimal_contraction_sequence

function tensor_network(; m, k, d)
  l = Index(m, "l")
  r = Index(m, "r")
  h₁ = Index(k, "h₁")
  h₂ = Index(k, "h₂")
  h₃ = Index(k, "h₃")
  s₁ = Index(d, "s₁")
  s₂ = Index(d, "s₂")

  ψ = emptyITensor(l, s₁, s₂, r)
  L = emptyITensor(dag(l), l', h₁)
  H₁ = emptyITensor(dag(s₁), s₁', dag(h₁), h₂)
  H₂ = emptyITensor(dag(s₂), s₂', dag(h₂), h₃)
  R = emptyITensor(dag(r), r', h₃)
  return [ψ, L, H₁, H₂, R]
end

function main()
  mrange = 50:10:80
  krange = 50:10:80
  sequence_costs = Matrix{Any}(undef, length(mrange), length(krange))
  for iₘ in eachindex(mrange), iₖ in eachindex(krange)
    m_val = mrange[iₘ]
    k_val = krange[iₖ]
    d_val = 4

    TN = tensor_network(; m = m_val, k = k_val, d = d_val)
    sequence = optimal_contraction_sequence(TN)
    cost = contraction_cost(TN; sequence = sequence)

    @variables m, k, d
    TN_symbolic = tensor_network(; m = m, k = k, d = d)
    cost_symbolic = contraction_cost(TN_symbolic; sequence = sequence)
    sequence_cost = (dims = (m = m_val, k = k_val, d = d_val), sequence = sequence, cost = cost, symbolic_cost = cost_symbolic)
    sequence_costs[iₘ, iₖ] = sequence_cost
  end
  return sequence_costs
end

sequence_costs = main()

# Analyze the results.
println("Index dimensions")
display(getindex.(sequence_costs, :dims))

println("\nContraction sequences")
display(getindex.(sequence_costs, :sequence))

println("\nSymbolic contraction cost with d = 4")
# Fix d to a certain value (such as 4 for a Hubbard site)
@variables d
var_sub = Dict(d => 4)
display(substitute.(sum.(getindex.(sequence_costs, :symbolic_cost)), (var_sub,)))

ext: jl


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------
# TODO: review
"""
    CookieCutter(master, others)

A cookie-cutter simulation solver.

## Parameters

* `master` - Master simulation solver (a.k.a. facies solver)
* `others` - A list of pairs mapping categories to solvers

## Examples

Simulate lithology facies with image quilting and fill property
with direct Gaussian simulation:

```julia
julia> fsolver  = ImgQuilt(:facies => (TI=Strebelle, template=(30,30,1)))
julia> psolver₀ = DirectGaussSim(:property => (variogram=SphericalVariogram(range=10.),))
julia> psolver₁ = DirectGaussSim(:property => (variogram=SphericalVariogram(range=20.),))
julia> solver   = CookieCutter(fsolver, [0 => psolver₀, 1 => psolver₁])
```
"""
struct CookieCutter <: AbstractSimulationSolver
  master::AbstractSimulationSolver
  others::Dict{Number,AbstractSimulationSolver}
end

function solve(problem::SimulationProblem, solver::CookieCutter)
  # retrieve problem info
  pdata   = data(problem)
  pdomain = domain(problem)
  pvars   = variables(problem)
  preals  = nreals(problem)
  pmapper = mapper(problem)

  # master variable
  mvars = variables(solver.master)
  @assert length(mvars) == 1 "one single variable must be specified in master solver"
  mvar = mvars[1]
  @assert mvar ∈ keys(pvars) "invalid variable in master solver"
  mtype = pvars[mvar]

  # other variables
  ovars = Tuple(var => V for (var, V) in pvars if var ≠ mvar)
  @assert length(ovars) > 0 "cookie-cutter requires problem with more than one target variable"

  # define master and others problem
  if hasdata(problem)
    mproblem = SimulationProblem(pdata, pdomain, mvar, preals, mapper=pmapper)
    oproblem = SimulationProblem(pdata, pdomain, first.(ovars), preals, mapper=pmapper)
  else
    mproblem = SimulationProblem(pdomain, mvar => mtype, preals, mapper=pmapper)
    oproblem = SimulationProblem(pdomain, ovars, preals, mapper=pmapper)
  end

  # realizations of master variable
  msol = solve(mproblem, solver.master)
  mreals = msol[mvar]

  # pre-allocate memory for realizations
  reals = Dict(var => [Vector{V}(undef, nelms(pdomain)) for i in 1:preals] for (var,V) in ovars)

  # solve other problems
  for (mval, osolver) in solver.others
    osol = solve(oproblem, osolver)

    # cookie-cutter step
    for (var, V) in ovars
      vreals = osol[var]
      for i in 1:preals
        mask = mreals[i] .== mval
        reals[var][i][mask] .= vreals[i][mask]
      end
    end
  end

  realizations = merge(reals, Dict(mvar => mreals))

  SimulationSolution(pdomain, realizations)
end

# ------------
# IO methods
# ------------
function Base.show(io::IO, solver::CookieCutter)
  print(io, "CookieCutter")
end

function Base.show(io::IO, ::MIME"text/plain", solver::CookieCutter)
  mvar = variables(solver.master)[1]
  println(io, solver)
  println(io, "  └─", mvar, " ⇨ ", solver.master)
  for (val, osolver) in solver.others
    println(io, "    └─", val, " ⇨ ", osolver)
  end
end


text:
#####
##### Terms in the turbulent kinetic energy equation, all at cell centers
#####

@inline ϕ²(i, j, k, grid, ϕ) = ϕ(i, j, k, grid)^2

@inline function shear_production(i, j, k, grid, closure, velocities, diffusivities)
    ∂z_u² = ℑxzᶜᵃᶜ(i, j, k, grid, ϕ², ∂zᵃᵃᶠ, velocities.u)
    ∂z_v² = ℑyzᵃᶜᶜ(i, j, k, grid, ϕ², ∂zᵃᵃᶠ, velocities.v)
    @inbounds Ku = diffusivities.Kᵘ[i, j, k]
    return Ku * (∂z_u² + ∂z_v²)
end

@inline function buoyancy_flux(i, j, k, grid, closure, velocities, tracers, buoyancy, diffusivities)
    @inbounds Kc = diffusivities.Kᶜ[i, j, k]
    N² = ℑzᵃᵃᶜ(i, j, k, grid, ∂z_b, buoyancy, tracers)
    return - Kc * N²
end

@inline function dissipation(i, j, k, grid, closure, velocities, tracers, buoyancy, clock, tracer_bcs)
    e = tracers.e
    FT = eltype(grid)
    @inbounds eⁱʲᵏ = e[i, j, k]
    ẽ³² = sqrt(abs(eⁱʲᵏ^3))

    ℓ = TKE_mixing_lengthᶜᶜᶜ(i, j, k, grid, closure, velocities, tracers, buoyancy, clock, tracer_bcs)
    Cᴰ = closure.Cᴰ

    return Cᴰ * ẽ³² / ℓ
end


text:
punctuation = ['!','"','#','$','%','&','\'','(',')','*','+',',',
'-','.','/',':',';','<','=','>','?','@','[','\\',']','^','_','`',
'{','|','}','~']

function clean(string)
  string = Base.lowercase(string)
  sarray = convert(Array{typeof(string)}, Base.split(string))
  x = 1

  while x <= length(sarray)
    sarray[x] = Base.strip(sarray[x], punctuation)
    x += 1
  end

  return sarray
end

function parse_xml(doc)
  xdoc = parse_file(doc)
  xroot = root(xdoc)
  ces = get_elements_by_tagname(xroot, "EEBO")
  children = collect(child_elements(ces[1]))
  if(name(children[2]) == "GROUP")
    ces = get_elements_by_tagname(ces[1], "GROUP")
  else
    ces = get_elements_by_tagname(ces[1], "TEXT")
  end
  body = content(ces[1])
  text = string(body)

  free(xdoc)

  return text
end

function get_metadata(doc)
  metadata=[]
  xroot = root(doc)
  ces = get_elements_by_tagname(xroot, "HEADER")
  ces = get_elements_by_tagname(ces[1], "FILEDESC")
  ces = get_elements_by_tagname(ces[1], "TITLESTMT")
  ces1 = get_elements_by_tagname(ces[1], "TITLE")
  body = content(ces1[1])
  title = string(body)
  push!(metadata,title)

  ces2 = get_elements_by_tagname(ces[3], "AUTHOR")
  body2 = content(ces2[1])
  author = string(body2)
  push!(metadata,author)

  return metadata
end

function load_featurevector(path)
  if isfile(path)
    fv = FeatureVector()
    parsed_file = parse_xml(path)
    clean_text = clean(parsed_file)
    for word in clean_text
      if word in keys(fv)
        fv[word] += 1
      else
        setindex!(fv,1,word)
      end
    end
    return fv
  else
    Base.warn("$path is not to a valid file. Please check the path and try again.")
  end
end

function load_cluster(path)
  if isdir(path)
    curdir = pwd()
    cd(path)
    files = readdir(pwd())
    cl = Cluster()
    for file in files
      if endswith(file,"xml")
        fv = load_featurevector(file)
        cl[file] = fv
      else
        Base.warn("$file is not a valid file. Please enter a path to a directory of .xml files.")
      end
    end
    cd(curdir)
    return cl
  else
    Base.warn("$path is not to a valid directory. Please check the path and try again.")
  end
end

function load_dataset(path)
  if isdir(path)
    curdir = pwd()
    cd(path)
    dirs = readdir(pwd())
    ds = DataSet()
    for dir in dirs
      cl = Cluster()
      if isdir(dir)
        cl = load_cluster(dir)
        ds[dir] = cl
      else
        Base.warn("$dir is not a valid Directory. This will not be added.")
      end
    end
    cd(curdir)
    return ds
  else
    Base.warn("$path is not to a valid directory. Please check the path and try again.")
  end
end

max_stars_repo_name: drawkward/SLU-TextMining

max_issues_repo_name: drawkward/SLU-TextMining

max_forks_repo_name: drawkward/SLU-TextMining


text:
# Notations from Walker & Ni, "Anderson acceleration for fixed-point iterations", SINUM 2011
# Attempts to accelerate the iteration xₙ₊₁ = xₙ + beta*f(xₙ)

struct Anderson{Tm, Tb}
    m::Tm
    beta::Tb
end
struct AndersonCache{Txs, Tx, Tr, Ta, Tf} <: AbstractSolverCache
    xs::Txs
    gs::Txs
    old_x::Tx
    residuals::Tr
    alphas::Ta
    fx::Tf
end
function AndersonCache(df, method::Anderson)
    m = method.m
    N = length(df.x_f)
    T = eltype(df.x_f)

    xs = zeros(T, N, m+1) #ring buffer storing the iterates, from newest to oldest
    gs = zeros(T, N, m+1) #ring buffer storing the g of the iterates, from newest to oldest
    old_x = xs[:,1]
    if m > 0
        residuals = zeros(T, N, m) #matrix of residuals used for the least-squares problem
        alphas = zeros(T, m) #coefficients obtained by least-squares
    else
        residuals = nothing
        alphas = nothing
    end
    fx = similar(df.x_f, N) # temp variable to store f!

    AndersonCache(xs, gs, old_x, residuals, alphas, fx)
end

@views function anderson_(df::Union{NonDifferentiable, OnceDifferentiable},
                             x0::AbstractArray{T},
                             xtol::T,
                             ftol::T,
                             iterations::Integer,
                             store_trace::Bool,
                             show_trace::Bool,
                             extended_trace::Bool,
                             m::Integer,
                             beta::Real,
                             cache = AndersonCache(df, Anderson(m, beta))) where T
    picard_iteration = cache.alphas == nothing
    copyto!(cache.xs[:,1], x0)
    iters = 0
    tr = SolverTrace()
    tracing = store_trace || show_trace || extended_trace
    x_converged, f_converged, converged = false, false, false

    errs = zeros(iterations)

    for n = 1:iterations
        iters += 1
        # fixed-point iteration
        value!!(df, cache.fx, cache.xs[:,1])

        cache.gs[:,1] .= cache.xs[:,1] .+ beta.*cache.fx

        x_converged, f_converged, converged = assess_convergence(cache.gs[:,1], cache.old_x, cache.fx, xtol, ftol)

        if tracing
            dt = Dict()
            if extended_trace
                dt["x"] = copy(cache.xs[:,1])
                dt["f(x)"] = copy(cache.fx)
            end
            update!(tr,
                    n,
                    maximum(abs,cache.fx),
                    n > 1 ? sqeuclidean(cache.xs[:,1],cache.old_x) : convert(T,NaN),
                    dt,
                    store_trace,
                    show_trace)
        end

        if converged
            break
        end

        new_x = copy(cache.gs[:,1])

        if !picard_iteration
            #update of new_x
            m_eff = min(n-1,m)
            if m_eff > 0
                cache.residuals[:, 1:m_eff] .= (cache.gs[:,2:m_eff+1] .- cache.xs[:,2:m_eff+1]) .- (cache.gs[:,1] .- cache.xs[:,1])
                cache.alphas[1:m_eff] .= cache.residuals[:,1:m_eff] \ (cache.xs[:,1] .- cache.gs[:,1])
                for i = 1:m_eff
                    new_x .+= cache.alphas[i].*(cache.gs[:,i+1] .- cache.gs[:,1])
                end
            end

            cache.xs .= circshift(cache.xs,(0,1)) # no in-place circshift, unfortunately...
            cache.gs .= circshift(cache.gs,(0,1))

            if m > 1
                copyto!(cache.old_x, cache.xs[:,2])
            else
                copyto!(cache.old_x, cache.xs[:,1])
            end
            copyto!(cache.xs[:,1], new_x)
        else
            copyto!(cache.old_x, cache.xs[:,1])
            copyto!(cache.xs[:,1], cache.gs[:,1])
        end
    end

    # returning gs[:,1] rather than xs[:,1] would be better here if
    # xₙ₊₁ = xₙ + beta*f(xₙ) is convergent, but the convergence
    # criterion is not guaranteed

    x = similar(x0) # this is done to ensure that the final x is oftype(x0)
    copyto!(x, cache.xs[:,1])
    return SolverResults("Anderson m=$m beta=$beta",
                         x0, x, maximum(abs,cache.fx),
                         iters, x_converged, xtol, f_converged, ftol, tr,
                         first(df.f_calls), 0)
end

function anderson(df::Union{NonDifferentiable, OnceDifferentiable},
                     initial_x::AbstractArray{T},
                     xtol::Real,
                     ftol::Real,
                     iterations::Integer,
                     store_trace::Bool,
                     show_trace::Bool,
                     extended_trace::Bool,
                     m::Integer,
                     beta::Real,
                     cache = AndersonCache(df, Anderson(m, beta))) where T
    anderson_(df, initial_x, convert(T, xtol), convert(T, ftol), iterations, store_trace, show_trace, extended_trace, m, beta, cache)
end


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

module RouteTableClient

using Random
using Dates
using HTTP
using Swagger
import Swagger: field_name, property_type, hasproperty, validate_property, SwaggerApi, SwaggerModel
import Base: convert, propertynames

include("modelincludes.jl")

include("api_RouteTablesApi.jl")
include("api_RoutesApi.jl")

# export models
export convert, Resource
export convert, RouteListResult
export convert, RoutePropertiesFormat
export convert, RouteTableListResult
export convert, RouteTablePropertiesFormat
export convert, SubResource
export convert, Route
export convert, RouteTable
export convert, Subnet

# export operations
export convert, RouteTablesApi, RoutesApi

export check_required, field_name, property_type, hasproperty, propertynames, validate_property, convert

end


text:
@testset "WordPiece" begin
  wp1 = WordPiece(["[UNK]", "[CLS]", "[SEP]", "want", "##want", "##ed", "wa", "un", "runn", "##ing", ","])

  @test wp1(bert_uncased_tokenizer("UNwant\u00E9d,running")) == ["un", "##want", "##ed", ",", "runn", "##ing"]
  @test wp1(Int, bert_uncased_tokenizer("UNwant\u00E9d,running")) == [8, 5, 6, 11, 9, 10]

  wp2 = WordPiece(["[UNK]", "[CLS]", "[SEP]", "want", "##want", "##ed", "wa", "un", "runn", "##ing"])

  @test wp2(bert_uncased_tokenizer("")) == []
  @test wp2(bert_uncased_tokenizer("unwanted running")) == ["un", "##want", "##ed", "runn", "##ing"]
  @test wp2(bert_uncased_tokenizer("unwantedX running")) == ["[UNK]", "runn", "##ing"]

  vocab = Vocabulary(wp2)
  @test vocab(["un", "##want", "##ed", "runn", "##ing"]) == [8, 5, 6, 9, 10]
end


text:
@testset "Cerner" begin
    # TODO: figure out how to test this
end


text:
#The MySQL Driver must be previously installed

using DataFrames
using ODBC
listdsns()
listdrivers()
advancedconnect("Driver={MySQL ODBC 3.51 Driver};Server=ensembldb.ensembl.org;Port=5306;User=anonymous;Option=3;",ODBC.SQL_DRIVER_NOPROMPT)

query("use homo_sapiens_vega_69_37")
query("select count(*) from exon")
query("show columns from exon")
query("select * from exon where phase = 2")
#DBMS may not support batch statements
query("select count(*) from exon; select count(*) from exon;")
#Write results to a file (delimiter can be specified as Char arg after filename; default is ',')
query("select count(*) from exon","test.csv")
query("select * from exon where phase = 2","test.csv")

querymeta("select count(*) from exon")

disconnect(conn)


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

mutable struct IoK8sApimachineryPkgApisMetaV1APIVersions <: SwaggerModel
    apiVersion::Any # spec type: Union{ Nothing, String } # spec name: apiVersion
    kind::Any # spec type: Union{ Nothing, String } # spec name: kind
    serverAddressByClientCIDRs::Any # spec type: Union{ Nothing, Vector{IoK8sApimachineryPkgApisMetaV1ServerAddressByClientCIDR} } # spec name: serverAddressByClientCIDRs
    versions::Any # spec type: Union{ Nothing, Vector{String} } # spec name: versions

    function IoK8sApimachineryPkgApisMetaV1APIVersions(;apiVersion=nothing, kind=nothing, serverAddressByClientCIDRs=nothing, versions=nothing)
        o = new()
        validate_property(IoK8sApimachineryPkgApisMetaV1APIVersions, Symbol("apiVersion"), apiVersion)
        setfield!(o, Symbol("apiVersion"), apiVersion)
        validate_property(IoK8sApimachineryPkgApisMetaV1APIVersions, Symbol("kind"), kind)
        setfield!(o, Symbol("kind"), kind)
        validate_property(IoK8sApimachineryPkgApisMetaV1APIVersions, Symbol("serverAddressByClientCIDRs"), serverAddressByClientCIDRs)
        setfield!(o, Symbol("serverAddressByClientCIDRs"), serverAddressByClientCIDRs)
        validate_property(IoK8sApimachineryPkgApisMetaV1APIVersions, Symbol("versions"), versions)
        setfield!(o, Symbol("versions"), versions)
        o
    end
end # type IoK8sApimachineryPkgApisMetaV1APIVersions

const _property_map_IoK8sApimachineryPkgApisMetaV1APIVersions = Dict{Symbol,Symbol}(Symbol("apiVersion")=>Symbol("apiVersion"), Symbol("kind")=>Symbol("kind"), Symbol("serverAddressByClientCIDRs")=>Symbol("serverAddressByClientCIDRs"), Symbol("versions")=>Symbol("versions"))
const _property_types_IoK8sApimachineryPkgApisMetaV1APIVersions = Dict{Symbol,String}(Symbol("apiVersion")=>"String", Symbol("kind")=>"String", Symbol("serverAddressByClientCIDRs")=>"Vector{IoK8sApimachineryPkgApisMetaV1ServerAddressByClientCIDR}", Symbol("versions")=>"Vector{String}")
Base.propertynames(::Type{ IoK8sApimachineryPkgApisMetaV1APIVersions }) = collect(keys(_property_map_IoK8sApimachineryPkgApisMetaV1APIVersions))
Swagger.property_type(::Type{ IoK8sApimachineryPkgApisMetaV1APIVersions }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApimachineryPkgApisMetaV1APIVersions[name]))}
Swagger.field_name(::Type{ IoK8sApimachineryPkgApisMetaV1APIVersions }, property_name::Symbol) =  _property_map_IoK8sApimachineryPkgApisMetaV1APIVersions[property_name]

function check_required(o::IoK8sApimachineryPkgApisMetaV1APIVersions)
    (getproperty(o, Symbol("serverAddressByClientCIDRs")) === nothing) && (return false)
    (getproperty(o, Symbol("versions")) === nothing) && (return false)
    true
end

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


text:
using Test
# @time @testset "Debug" begin include("test_debug.jl") end
@time @testset "Miscellaneous" begin include("test_miscellaneous.jl") end
@time @testset "Acoustics" begin include("test_acoustics.jl") end
@time @testset "Heat diffusion" begin include("test_heat.jl") end
@time @testset "Linear deformation" begin include("test_linear_deformation.jl") end
@time @testset "Meshing" begin include("test_meshing.jl") end
@time @testset "Voxel box" begin include("test_voxel_box.jl") end
# @time @testset "Failing tests" begin include("test_failing.jl") end
true


text:
import Dates
import DateShifting
import Distributions
import Random
import TimeZones
import Test

Test.@testset "DateShifting.jl" begin
    Test.@testset "assert.jl" begin
        Test.@test nothing == Test.@test_nowarn DateShifting.always_assert(1 == 1)
        Test.@test_throws DateShifting.AlwaysAssertionError DateShifting.always_assert(1 == 2)
    end

    Test.@testset "compute_interval_nonrounded.jl" begin
        Test.@test DateShifting._compute_interval_nonrounded(Dates.DateTime("2000-01-01"), Dates.DateTime("2000-01-02")) == Dates.Day(1)
        Test.@test_throws ArgumentError DateShifting._compute_interval_nonrounded(Dates.DateTime("2000-01-02"), Dates.DateTime("2000-01-01"))
    end

    Test.@testset "sample.jl" begin
        Test.@test DateShifting._draw_sample(Random.GLOBAL_RNG, Dates.Day, nothing) == Dates.Day(0)
    end

    Test.@testset "intervals.jl" begin
        datetime_list = [
            Dates.DateTime("1900-01-01T00:00:00"),
            Dates.DateTime("1901-01-01T00:00:00"),
            Dates.DateTime("1900-02-01T00:00:00"),
            Dates.DateTime("1900-01-01T00:00:00"),
            Dates.DateTime("1900-01-03T00:00:00"),
            Dates.DateTime("1900-01-05T00:00:00"),
            Dates.DateTime("1900-01-02T00:00:00"),
            Dates.DateTime("1900-01-01T00:00:00.001"),
            Dates.DateTime("1900-01-01T00:20:00"),
            Dates.DateTime("1900-01-01T00:30:00"),
        ]

        zoneddatetime_list = [TimeZones.ZonedDateTime(x, TimeZones.tz"America/New_York") for x in datetime_list]

        expected_sequence = [1, 4, 8, 9, 10, 7, 5, 6, 3, 2]

        expected_intervals_round_to = Dict()
        expected_intervals_round_to[Dates.Day(1)] = [Dates.Day(0), Dates.Day(365), Dates.Day(31), Dates.Day(0), Dates.Day(2), Dates.Day(4), Dates.Day(1), Dates.Day(0), Dates.Day(0), Dates.Day(0)]
        expected_intervals_round_to[Dates.Day(2)] = [Dates.Day(0), Dates.Day(366), Dates.Day(32), Dates.Day(0), Dates.Day(2), Dates.Day(4), Dates.Day(2), Dates.Day(0), Dates.Day(0), Dates.Day(0)]
        expected_intervals_round_to[Dates.Hour(1)] = [Dates.Hour(0), Dates.Hour(8760), Dates.Hour(744), Dates.Hour(0), Dates.Hour(48), Dates.Hour(96), Dates.Hour(24), Dates.Hour(0), Dates.Hour(0), Dates.Hour(1)]
        expected_intervals_round_to[Dates.Hour(2)] = [Dates.Hour(0), Dates.Hour(8760), Dates.Hour(744), Dates.Hour(0), Dates.Hour(48), Dates.Hour(96), Dates.Hour(24), Dates.Hour(0), Dates.Hour(0), Dates.Hour(0)]
        expected_intervals_round_to[Dates.Minute(1)] = [Dates.Minute(0), Dates.Minute(525600), Dates.Minute(44640), Dates.Minute(0), Dates.Minute(2880), Dates.Minute(5760), Dates.Minute(1440), Dates.Minute(0), Dates.Minute(20), Dates.Minute(30)]
        expected_intervals_round_to[Dates.Minute(2)] = [Dates.Minute(0), Dates.Minute(525600), Dates.Minute(44640), Dates.Minute(0), Dates.Minute(2880), Dates.Minute(5760), Dates.Minute(1440), Dates.Minute(0), Dates.Minute(20), Dates.Minute(30)]
        expected_intervals_round_to[Dates.Second(1)] = [Dates.Second(0), Dates.Second(31536000), Dates.Second(2678400), Dates.Second(0), Dates.Second(172800), Dates.Second(345600), Dates.Second(86400), Dates.Second(0), Dates.Second(1200), Dates.Second(1800)]
        expected_intervals_round_to[Dates.Second(2)] = [Dates.Second(0), Dates.Second(31536000), Dates.Second(2678400), Dates.Second(0), Dates.Second(172800), Dates.Second(345600), Dates.Second(86400), Dates.Second(0), Dates.Second(1200), Dates.Second(1800)]

        for round_to_period in keys(expected_intervals_round_to)
            for input in Any[datetime_list, zoneddatetime_list]
                observed_sequence, observed_intervals = DateShifting.datetime_intervals(input;
                                                                                round_to = round_to_period)
                Test.@test observed_sequence == expected_sequence
                Test.@test observed_intervals == expected_intervals_round_to[round_to_period]
            end
        end
    end

    Test.@testset "public_sequence_and_random_date_shift.jl" begin
        datetime_list = [
            Dates.DateTime("1900-01-01T00:00:00"),
            Dates.DateTime("1901-01-01T00:00:00"),
            Dates.DateTime("1900-02-01T00:00:00"),
            Dates.DateTime("1900-01-01T00:00:00"),
            Dates.DateTime("1900-01-03T00:00:00"),
            Dates.DateTime("1900-01-05T00:00:00"),
            Dates.DateTime("1900-01-02T00:00:00"),
            Dates.DateTime("1900-01-01T00:20:00"),
            Dates.DateTime("1900-01-01T00:30:00"),
        ]

        zoneddatetime_list = [TimeZones.ZonedDateTime(x, TimeZones.tz"America/Los_Angeles") for x in datetime_list]

        expected_sequence = [1, 4, 8, 9, 7, 5, 6, 3, 2]

        round_to_periods = Dates.Period[
            Dates.Second(1),
            Dates.Second(2),
            Dates.Second(3),
            Dates.Second(5),
            Dates.Second(10),
            Dates.Millisecond(1),
            Dates.Millisecond(2),
            Dates.Millisecond(3),
            Dates.Millisecond(5),
            Dates.Millisecond(10),
            Dates.Millisecond(100),
            Dates.Millisecond(500),
            Dates.Millisecond(1000),
        ]

        for round_to_period in round_to_periods
            for input in Any[datetime_list, zoneddatetime_list]
                observed_sequence, observed_shifted_dates = DateShifting.sequence_and_random_date_shift(
                    Random.GLOBAL_RNG,
                    input;
                    round_to = round_to_period,
                    time_zone = TimeZones.TimeZone("America/Chicago"),
                    day = Distributions.DiscreteUniform(-31, 31),
                    hour = Distributions.DiscreteUniform(-24, 24),
                    minute = Distributions.DiscreteUniform(-60, 60),
                    second = Distributions.DiscreteUniform(-60, 60),
                    millisecond = Distributions.DiscreteUniform(-1000, 1000),
                    microsecond = Distributions.DiscreteUniform(-1000, 1000),
                    nanosecond = nothing,
                )
                Test.@test observed_sequence == expected_sequence
                Test.@test sortperm(input) == observed_sequence
                Test.@test sortperm(input) == expected_sequence
                Test.@test sortperm(observed_shifted_dates) == observed_sequence
                Test.@test sortperm(observed_shifted_dates) == expected_sequence
            end
        end
    end
end


text:
function advection(  m:: BurgerModel, x :: BurgerData , AD :: ResolveLL)
  xx= x.x
  n_off = AD.n_off
  result = zeros(xx)
  for ii = 1:n_off
    for kk = 1:n_off
      cross_iikk = -0.5*xx[ii,:] .*xx[kk,:]*kk
      if (ii + kk <= x.n )
        result[ii+kk,:] += cross_iikk
      end
      if (ii - kk > 0 )
        result[ii-kk,:] += cross_iikk
      elseif (ii - kk < 0 )
        result[kk-ii,:] -= cross_iikk
      end
    end
  end
  return result
end

function advection(  m:: BurgerModel, x :: BurgerData , AD :: ResolveSS)
  xx= x.x
  n_off = AD.n_off
  result = zeros(xx)
  for ii = (1+n_off):x.n
    for kk = (1+n_off):x.n
      cross_iikk = -0.5*xx[ii,:] .*xx[kk,:]*kk
      if (ii + kk <= x.n )
        result[ii+kk,:] += cross_iikk
      end
      if (ii - kk > 0 )
        result[ii-kk,:] += cross_iikk
      elseif (ii - kk < 0 )
        result[kk-ii,:] -= cross_iikk
      end
    end
  end
  return result
end

function advection(  m:: BurgerModel, x :: BurgerData , AD :: ResolveLS)
  xx= x.x
  n_off = AD.n_off
  result = zeros(xx)
  for ii = 1:n_off
    for kk = (1+n_off):x.n
      cross_iikk = -0.5*xx[ii,:] .*xx[kk,:]*kk
      if (ii + kk <= x.n )
        result[ii+kk,:] += cross_iikk
      end
      if (ii - kk > 0 )
        result[ii-kk,:] += cross_iikk
      elseif (ii - kk < 0 )
        result[kk-ii,:] -= cross_iikk
      end
    end
  end
  return result
end

function advection(  m:: BurgerModel, x :: BurgerData , AD :: ResolveSL)
  xx= x.x
  n_off = AD.n_off
  result = zeros(xx)
  for ii = 1:n_off
    for kk = (1+n_off):x.n
      cross_iikk = -0.5*xx[ii,:] .*xx[kk,:]*ii
      if (ii + kk <= x.n )
        result[ii+kk,:] += cross_iikk
      end
      if (ii - kk > 0 )
        result[ii-kk,:] -= cross_iikk
      elseif (ii - kk < 0 )
        result[kk-ii,:] += cross_iikk
      end
    end
  end
  return result
end


text:
#=

Implementations of cg and pcg.
Need to add conditions for termination when stagnate.
Look at two approaches: Matlab's, and
hypre: https://github.com/LLNL/hypre/blob/master/src/krylov/pcg.c
Started by Dan Spielman.
Contributors:

=#

"""
    x = cg(mat, b; tol, maxits, maxtime, verbose, pcgIts)

solves a symmetric linear system `mat x = b`.
# Arguments
* `tol` is set to 1e-6 by default,
* `maxits` defaults to Inf
* `maxtime` defaults to Inf.  It measures seconds.
* `verbose` defaults to false
* `pcgIts` is an array for returning the number of pcgIterations.  Default is length 0, in which case nothing is returned.
"""
function cg end

"""
    x = pcg(mat, b, pre; tol, maxits, maxtime, verbose, pcgIts, stag_test)`

solves a symmetric linear system using preconditioner `pre`.
# Arguments
* `pre` can be a function or a matrix.  If a matrix, a function to solve it is created with cholFact.
* `tol` is set to 1e-6 by default,
* `maxits` defaults to Inf
* `maxtime` defaults to Inf.  It measures seconds.
* `verbose` defaults to false
* `pcgIts` is an array for returning the number of pcgIterations.  Default is length 0, in which case nothing is returned.
* `stag_test=k` stops the code if rho[it] > (1-1/k) rho[it-k].  Set to 0 to deactivate.
"""
function pcg end

"""
    x = cgSolver(mat; tol, maxits, maxtime, verbose, pcgIts)

creates a solver for a PSD system `mat`.
The parameters are as described in cg.
"""
function cgSolver end

"""
    x = pcgSolver(mat, pre; tol, maxits, maxtime, verbose, pcgIts)

creates a solver for a PSD system using preconditioner `pre`.
The parameters are as described in pcg.
"""
function pcgSolver end

function cgSolver(mat; tol::Real=1e-6, maxits=Inf, maxtime=Inf, verbose=false, pcgIts=Int[])

    tol_=tol
    maxits_=maxits
    maxtime_=maxtime
    verbose_=verbose
    pcgIts_=pcgIts

    f(b; tol=tol_, maxits=maxits_, maxtime=maxtime_, verbose=verbose_, pcgIts=pcgIts_) =
       cg(mat, b, tol=tol, maxits=maxits, maxtime=maxtime, verbose=verbose, pcgIts=pcgIts)
end

"""
    x = cgLapSolver(A::AbstractMatrix; tol::Real=1e-6, maxits=Inf, maxtime=Inf, verbose=false, pcgIts=Int[])

Create a solver that uses cg to solve Laplacian systems in the laplacian of A.
This just exists to satisfy our interface.
It does nothing more than create the Laplacian and call cg on each connected component.
"""
function cgLapSolver(a::SparseMatrixCSC; tol::Real=1e-6, maxits=Inf, maxtime=Inf, verbose=false, pcgIts=Int[])

    return lapWrapComponents(cgLapSolver1, a, verbose=verbose, tol=tol, maxits=maxits, maxtime=maxtime, pcgIts=pcgIts)

end

function cgLapSolver1(A::AbstractMatrix; tol::Real=1e-6, maxits=Inf, maxtime=Inf, verbose=false, pcgIts=Int[])

    la = forceLap(A)

    tol_=tol
    maxits_=maxits
    maxtime_=maxtime
    verbose_=verbose
    pcgIts_=pcgIts

    f(b; tol=tol_, maxits=maxits_, maxtime=maxtime_, verbose=verbose_, pcgIts=pcgIts_) =
        cg(la, b .- mean(b), tol=tol, maxits=maxits, maxtime=maxtime, verbose=verbose, pcgIts=pcgIts)

end

function pcg(mat, b, pre::Union{AbstractArray,Matrix}; tol::Real=1e-6, maxits=Inf, maxtime=Inf, verbose=false, pcgIts=Int[])
    fact = cholesky(pre)
    F = x->(fact \ x)
    pcg(mat, b, F; tol=tol, maxits=maxits, maxtime=maxtime, verbose=verbose, pcgIts=pcgIts)
end

function pcgSolver(mat, pre; tol::Real=1e-6, maxits=Inf, maxtime=Inf, verbose=false, pcgIts=Int[])
    tol_=tol
    maxits_=maxits
    maxtime_=maxtime
    verbose_=verbose
    pcgIts_=pcgIts

    f(b; tol=tol_, maxits=maxits_, maxtime=maxtime_, verbose=verbose_, pcgIts=pcgIts_) =
        pcg(mat, b, pre, tol=tol, maxits=maxits, maxtime=maxtime, verbose=verbose, pcgIts=pcgIts)
end

"""
    x = pcgLapSolver(A, B; tol::Real=1e-6, maxits=Inf, maxtime=Inf, verbose=false, pcgIts=Int[])

Create a solver that uses pcg to solve Laplacian systems in `A`
Specialized for the case when the preconditioner the Laplacian matrix of `B`.
It solves the preconditioner by Cholesky Factorization.
"""
function pcgLapSolver(A::AbstractMatrix, B::AbstractMatrix; tol::Real=1e-6, maxits=Inf, maxtime=Inf, verbose=false, pcgIts=Int[])
    fact = chol_lap(B)

    tol_=tol
    maxits_=maxits
    maxtime_=maxtime
    verbose_=verbose
    pcgIts_=pcgIts

    la = forceLap(A)

    f(b; tol=tol_, maxits=maxits_, maxtime=maxtime_, verbose=verbose_, pcgIts=pcgIts_) =
        pcg(la, b, fact, tol=tol, maxits=maxits, maxtime=maxtime, verbose=verbose, pcgIts=pcgIts)

end

function cg(mat, b::Vector{Tval};
        tol::Real=1e-6, maxits=Inf, maxtime=Inf, verbose=false, pcgIts=Int[]) where Tval

    local al::Tval

    n = size(mat,2)

    nb = norm(b)

    # If input vector is zero, quit
    if nb == 0
      return zeros(size(b))
    end

    x = zeros(Tval,n)
    bestx = zeros(Tval,n)
    bestnr = one(Tval)

    r = copy(b)

    p = copy(r)

    rho = norm(r)^2

    t1 = time()

    itcnt = 0
    while itcnt < maxits
        itcnt = itcnt+1

        q = mat*p

        al = rho/dot(p, q)

        if (al < eps(Tval) || isinf(al))
          if verbose
            println("CG Stopped due to small or large alpha")
          end
          break
        end

        axpy2!(al,p,x)
        # x = x + al * p

        axpy2!(-al,q,r)
        #r .= r .- al.*q

        nr = norm(r)/nb

        if (nr < eps(Tval) || isinf(nr))
          if verbose
            println("CG Stopped due to small or large norm(r)")
          end
          break
        end

        if nr < bestnr
          bestnr = nr
          @inbounds @simd for i in 1:n
            bestx[i] = x[i]
          end
        end
        if nr < tol #Converged?
            break
        end

        # here is the top of the code in numerical templates

        oldrho = rho
        rho = norm(r)^2
        if (rho < eps(Tval) || isinf(rho))
          println("CG stopped for rho: ", rho)
          break
        end

        # the following would have higher accuracy
        #       rho = sum(r.^2)

        beta = rho/oldrho
        if (beta< eps(Tval) || isinf(beta))
          println("CG stopped for beta: ", beta)
          break
        end

        bzbeta!(beta,p,r)
        # p = z + beta*p

        if (time() - t1) > maxtime
            if verbose
                println("CG stopped at maxtime.")
            end
            break
        end

    end

    if verbose
        println("CG stopped after: ", round((time() - t1),digits=3), " seconds and ", itcnt, " iterations with relative error ", (norm(r)/norm(b)), ".")
    end

    if length(pcgIts) > 0
        pcgIts[1] = itcnt
    end

    return bestx
end

function pcg(mat, b::Vector{Tval}, pre::Function;
        tol::Real=1e-6, maxits=Inf, maxtime=Inf, verbose=false, pcgIts=Int[],
        stag_test::Integer=0) where Tval

    local al::Tval

    n = size(mat,2)

    nb = norm(b)

    # If input vector is zero, quit
    if nb == 0
      return zeros(size(b))
    end

    x = zeros(Tval,n)
    bestx = zeros(Tval,n)
    bestnr = 1.0

    r = copy(b)
    z = pre(r)
    p = copy(z)

    rho = dot(r, z)
    best_rho = rho
    stag_count = 0

    t1 = time()

    itcnt = 0
    while itcnt < maxits
        itcnt = itcnt+1

        q = mat*p

        pq = dot(p,q)

        if (pq < eps(Tval) || isinf(pq))
          if verbose
            println("PCG Stopped due to small or large pq")
          end
          break
        end

        al = rho/pq

        # the following line could cause slowdown
        if al*norm(p) < eps(Tval)*norm(x)
          if verbose
            println("PCG: Stopped due to stagnation.")
          end
          break
        end

        axpy2!(al,p,x)
        # x = x + al * p
        #=
        @inbounds @simd for i in 1:n
            x[i] += al*p[i]
        end
        =#
        #axpy

        axpy2!(-al,q,r)
        #r .= r .- al.*q
        #=
        @inbounds @simd for i in 1:n
            r[i] -= al*q[i]
        end
        =#

        nr = norm(r)/nb
        if nr < bestnr
          bestnr = nr
          @inbounds @simd for i in 1:n
            bestx[i] = x[i]
          end
        end
        if nr < tol #Converged?
            break
        end

        # here is the top of the code in numerical templates

        z = pre(r)

        oldrho = rho
        rho = dot(z, r) # this is gamma in hypre.

        if rho < best_rho*(1-1/stag_test)
          best_rho = rho
          stag_count = 0
        else
          if stag_test > 0
            if best_rho > (1-1/stag_test)*rho
              stag_count += 1
              if stag_count > stag_test
                println("PCG Stopped by stagnation test ", stag_test)
                break
              end
            end
          end
        end

        if (rho < eps(Tval) || isinf(rho))
          if verbose
            println("PCG Stopped due to small or large rho")
          end
          break
        end

        # the following would have higher accuracy
        #       rho = sum(r.^2)

        beta = rho/oldrho
        if (beta < eps(Tval) || isinf(beta))
          if verbose
            println("PCG Stopped due to small or large beta")
          end
          break
        end

        bzbeta!(beta,p,z)
        #=
        # p = z + beta*p
        @inbounds @simd for i in 1:n
            p[i] = z[i] + beta*p[i]
        end
        =#

        if (time() - t1) > maxtime
            if verbose
                println("PCG New stopped at maxtime.")
            end
            break
        end

    end

    if verbose
        println("PCG stopped after: ", round((time() - t1),digits=3), " seconds and ", itcnt, " iterations with relative error ", (norm(r)/norm(b)), ".")
    end

    if length(pcgIts) > 0
        pcgIts[1] = itcnt
    end

    return bestx
end

function axpy2!(al,p::Array,x::Array)
  n = length(x)
  @inbounds @simd for i in 1:n
      x[i] = x[i] + al*p[i]
  end
end

# p[i] = z[i] + beta*p[i]
function bzbeta!(beta,p::Array,z::Array)
  n = length(p)
  @inbounds @simd for i in 1:n
      p[i] = z[i] + beta*p[i]
  end
end

#=
function pcgDiagnose{Tval}(mat, b::Array{Tval,1}, pre::Function;
        tol::Real=1e-6, maxits=Inf, maxtime=Inf, verbose=false, pcgIts=Int[])

    local al::Tval

    n = size(mat,2)

    dic = Dict()
    dic["rho"] = []
    dic["alpha"] = []
    dic["beta"] = []
    dic["nr"] = []
    dic["pq"] = []

    nb = norm(b)

    # If input vector is zero, quit
    if nb == 0
      return bestx
    end

    x = zeros(Tval,n)
    bestx = zeros(Tval,n)
    bestnr = one(Tval)

    r = copy(b)
    z = pre(r)
    p = copy(z)

    rho = dot(r, z)
    push!(dic["rho"],rho)

    t1 = time()

    itcnt = 0
    while itcnt < maxits
        itcnt = itcnt+1

        q = mat*p

        pq = dot(p,q)
        push!(dic["pq"],pq)

        if (pq < eps(Tval) || isinf(pq))
          if verbose
            println("PCG Stopped due to small or large pq")
          end
          break
        end

        al = rho/pq
        push!(dic["alpha"],al)

        if al*norm(p) < eps(Tval)*norm(x)
          if verbose
            println("PCG: Stopped due to stagnation.")
          end
          break
        end

        axpy2!(al,p,x)
        # x = x + al * p

        #axpy

        axpy2!(-al,q,r)
        #r .= r .- al.*q

        nr = norm(r)/nb
        push!(dic["nr"],nr)

        if nr < bestnr
          bestnr = nr
          @inbounds @simd for i in 1:n
            bestx[i] = x[i]
          end
        end
        if nr < tol #Converged?
            break
        end

        # here is the top of the code in numerical templates

        z = pre(r)

        oldrho = rho
        rho = dot(z, r) # this is gamma in hypre.

        push!(dic["rho"],rho)
        if (rho < eps(Tval) || isinf(rho))
          if verbose
            println("PCG Stopped due to small or large rho")
          end
          break
        end

        # the following would have higher accuracy
        #       rho = sum(r.^2)

        beta = rho/oldrho

        push!(dic["beta"],beta)
        if (beta < eps(Tval) || isinf(beta))
          if verbose
            println("PCG Stopped due to small or large beta")
          end
          break
        end

        bzbeta!(beta,p,z)

        if (time() - t1) > maxtime
            if verbose
                println("PCG New stopped at maxtime.")
            end
            break
        end

    end

    if verbose
        println("PCG stopped after: ", round((time() - t1),3), " seconds and ", itcnt, " iterations with relative error ", (norm(r)/norm(b)), ".")
    end

    if length(pcgIts) > 0
        pcgIts[1] = itcnt
    end

    return bestx, dic
end
=#

#==========================================================
   Code we no longer use
===========================================================#

#==========================================================
   Code we no longer use
===========================================================#

# uses BLAS.  As fast as Matlab's pcg.
# set to use similar paramaters
function cgBLAS(mat, b::Vector{Tval};
        tol::Real=1e-6, maxits=Inf, maxtime=Inf, verbose=false, pcgIts=Int[]) where Tval

    n = size(mat,2)

    x = zeros(Tval,n)

    tol = tol * norm(b)
    r = copy(b)
    p = copy(r)

    # the following would have higher accuracy
    #       rho = sum(r.^2)
    rho = norm(r)^2

    t1 = time()

    itcnt = 0
    while itcnt < maxits
        itcnt = itcnt+1

        q = mat*p

        al = rho/dot(p, q)

        BLAS.axpy!(al,p,x)  # x = x + al * p

        BLAS.axpy!(-al,q,r)  # r -= al*q

        if norm(r) < tol #Converged?
            break
        end

        oldrho = rho
        rho = norm(r)^2

        # the following would have higher accuracy
        #       rho = sum(r.^2)

        beta = rho/oldrho

        # p = r + beta*p
        BLAS.scal!(n,beta,p,1) # p *= beta
        BLAS.axpy!(1.0,r,p) # p += r

        if (time() - t1) > maxtime
            if verbose
                println("CG BLAS stopped at maxtime.")
            end
            break
        end
      end

    if verbose
        println("CG BLAS stopped after: ", round((time() - t1),digits=3), " seconds and ", itcnt, " iterations with relative error ", (norm(r)/norm(b)), ".")

    end

    if length(pcgIts) > 0
        pcgIts[1] = itcnt
    end

    return x
end

# For when you can't use BLAS.
# set to use similar paramaters similar to MATLAB
function cgSlow(mat, b::Vector{Tval};
        tol::Real=1e-6, maxits=Inf, maxtime=Inf, verbose=false, pcgIts=Int[]) where Tval

    n = size(mat,2)

    x = zeros(Tval,n)

    tol = tol * norm(b)
    r = copy(b)
    p = copy(r)

    # the following would have higher accuracy
    #       rho = sum(r.^2)
    rho = norm(r)^2

    t1 = time()

    itcnt = 0
    while itcnt < maxits
        itcnt = itcnt+1

        q = mat*p

        al = rho/dot(p, q)

        # x = x + al * p
        for i in 1:n
            x[i] = x[i] + al*p[i]
        end

        # r = r - al*q
        for i in 1:n
            r[i] = r[i] - al*q[i]
        end

        if norm(r) < tol
            break
        end

        oldrho = rho
        rho = norm(r)^2

        # the following would have higher accuracy
        #       rho = sum(r.^2)

        beta = rho/oldrho

        # p = r + beta*p
        for i in 1:n
            p[i] = r[i] + beta*p[i]
        end

        if (time() - t1) > maxtime
            if verbose
                println("CG Slow stopped at maxtime.")
            end
            break
        end
      end

    if verbose
        println("CG Slow stopped after: ", round((time() - t1),digits=3), " seconds and ", itcnt, " iterations with relative error ", (norm(r)/norm(b)), ".")
    end

    if length(pcgIts) > 0
        pcgIts[1] = itcnt
    end

    return x
end

#====================================
   PCG Code
=====================================#

# uses BLAS.  As fast as Matlab's pcg.
# set to use similar parameters "
function pcgBLAS(mat, b::Vector{Tval}, pre;
        tol::Real=1e-6, maxits=Inf, maxtime=Inf, verbose=false, pcgIts=Int[]) where Tval

    n = size(mat,2)

    x = zeros(Tval,n)

    tol = tol * norm(b)
    r = copy(b)
    z = pre(r)
    p = copy(z)

    rho = dot(r, z) # BLAS.dot does not seem faster

    t1 = time()

    itcnt = 0
    while itcnt < maxits
        itcnt = itcnt+1

        q = mat*p

        if dot(p,q) == 0
            break
        end

        al = rho/dot(p, q)

        BLAS.axpy!(al,p,x)  # x = x + al * p

        BLAS.axpy!(-al,q,r)  # r -= al*q

        if norm(r) < tol #Converged?
            break
        end

        # here is the top of the code in numerical templates

        z = pre(r)

        oldrho = rho
        rho = dot(z, r)

        # the following would have higher accuracy
        #       rho = sum(r.^2)

        beta = rho/oldrho

        p = z + beta*p
        # BLAS.scal!(n,beta,p,1) # p *= beta
        # BLAS.axpy!(1.0,z,p) # p += z

        if (time() - t1) > maxtime
            if verbose
                println("PCG BLAS stopped at maxtime.")
            end
            break
        end

      end

    if verbose
        println("PCG BLAS stopped after: ", round((time() - t1),digits=3), " seconds and ", itcnt, " iterations with relative error ", (norm(r)/norm(b)), ".")
    end

    if length(pcgIts) > 0
        pcgIts[1] = itcnt
    end

    return x
end

function pcgSlow(mat, b::Vector{Tval}, pre;
        tol::Real=1e-6, maxits=Inf, maxtime=Inf, verbose=false, pcgIts=Int[]) where Tval

    n = size(mat,2)

    x = zeros(Tval,n)

    tol = tol * norm(b)
    r = copy(b)
    z = pre(r)
    p = copy(z)

    rho = dot(r, z)

    t1 = time()

    itcnt = 0
    while itcnt < maxits
        itcnt = itcnt+1

        q = mat*p

        al = rho/dot(p, q)

        # x = x + al * p
        for i in 1:n
            x[i] = x[i] + al*p[i]
        end

        # r = r - al*q
        for i in 1:n
            r[i] = r[i] - al*q[i]
        end

        if norm(r) < tol #Converged?
            break
        end

        # here is the top of the code in numerical templates

        z = pre(r)

        oldrho = rho
        rho = dot(z, r)

        # the following would have higher accuracy
        #       rho = sum(r.^2)

        beta = rho/oldrho

        # p = z + beta*p
        for i in 1:n
            p[i] = z[i] + beta*p[i]
        end

        if (time() - t1) > maxtime
            if verbose
                println("PCG Slow stopped at maxtime.")
            end
            break
        end

    end

    if verbose
        println("PCG Slow stopped after: ", round((time() - t1),digits=3), " seconds and ", itcnt, " iterations with relative error ", (norm(r)/norm(b)), ".")
    end

    if length(pcgIts) > 0
        pcgIts[1] = itcnt
    end

    return x
end


text:
using BrowseTables, Tables, Test

@test 1 == 1


text:
function binary_reduce(op::F, ::Type{T}, x, final) where {F,T}
	N = length(x)
	out = Vector{T}(undef, div(N+1,2))
	for i=1:div(N,2)
		out[i] = op(x[2i-1], x[2i])
	end
	isodd(N) && (out[end] = final)
	out
end
binary_reduce(op, x::AbstractVector{T}, final) where T = binary_reduce(op, T, x, final)

function pyramid(f, x, nbrLevels)
	v = [x]
	for i=2:nbrLevels
		push!(v, f(v[end]))
	end
	v
end

diagonalsum(x) = (N=length(x); sum(z->z[1]*z[2], zip(x,-N+1:2:N-1)))

# Function barrier to help compiler
diagonal_reduction(::Type{T},N,diagonalSums,intervalSums,intervalSize) where T =
	binary_reduce((i1,i2)->diagonalSums[i1]+diagonalSums[i2]+intervalSize*intervalSums[i2] - min(N-i1*intervalSize,intervalSize)*intervalSums[i1], T, 1:length(diagonalSums), diagonalSums[end])
function diagonal_reduction(t::Tuple{AbstractVector{T},Int}, allIntervalSums, N, leafSize) where T
	diagonalSums,level = t[1],t[2]
	diagonalSums = diagonal_reduction(T, N, diagonalSums, allIntervalSums[level], leafSize*(2^(level-1)))
	diagonalSums,level+1
end

nbrlevels(N,leafSize) = Int(log2(nextpow(2, div(N+leafSize-1,leafSize))))+1

abstract type Tree end
depth(t::Tree) = length(t.intervalSums)

struct SumTree{T} <: Tree
	leafSize::Int
	intervalSums::Vector{Vector{T}}
	diagonalSums::Vector{Vector{T}}
end
function SumTree(X::AbstractVector{T}, leafSize) where T
	N = length(X)
	nbrLevels = nbrlevels(N,leafSize)

	intervalSums = map(sum, Iterators.partition(X,leafSize))
	allIntervalSums = pyramid(z->binary_reduce(+,z,z[end]), intervalSums, nbrLevels)

	diagonalSums = map(diagonalsum, Iterators.partition(X,leafSize))
	allDiagonalSums = first.(pyramid(z->diagonal_reduction(z,allIntervalSums,N,leafSize), (diagonalSums,1), nbrLevels))

	SumTree(leafSize, reverse!(allIntervalSums), reverse!(allDiagonalSums))
end

struct MinMaxTree{T} <: Tree
	leafSize::Int
	intervalSums::Vector{Vector{T}}
	mins::Vector{Vector{T}}
	maxes::Vector{Vector{T}}
end
function MinMaxTree(X, Y, leafSize)
	@assert length(X)==length(Y)
	nbrLevels = nbrlevels(length(X),leafSize)
	intervalSums = map(sum, Iterators.partition(X,leafSize))
	allIntervalSums = pyramid(z->binary_reduce(+,z,z[end]), intervalSums, nbrLevels)
	mins = map(minimum, Iterators.partition(Y,leafSize))
	allMins = pyramid(z->binary_reduce(min,z,z[end]), mins, nbrLevels)
	maxes = map(maximum, Iterators.partition(Y,leafSize))
	allMaxes = pyramid(z->binary_reduce(max,z,z[end]), maxes, nbrLevels)
	MinMaxTree(leafSize,reverse!(allIntervalSums),reverse!(allMins),reverse!(allMaxes))
end

struct WeightTree{T} <: Tree
	leafSize::Int
	intervalSums::Vector{Vector{T}}
end
function WeightTree(X, leafSize)
	nbrLevels = nbrlevels(length(X),leafSize)
	intervalSums = map(sum, Iterators.partition(X,leafSize))
	allIntervalSums = pyramid(z->binary_reduce(+,z,z[end]), intervalSums, nbrLevels)
	WeightTree(leafSize,reverse!(allIntervalSums))
end
WeightTree(tree::MinMaxTree) = WeightTree(tree.leafSize, tree.intervalSums)


text:
struct TestSourceWithoutLength
end

function Base.eltype(iter::TestSourceWithoutLength)
    return NamedTuple{(:a, :b),Tuple{Int,Float64}}
end

Base.IteratorSize(::Type{T}) where {T <: TestSourceWithoutLength} = Base.SizeUnknown()

function Base.iterate(iter::TestSourceWithoutLength, state=1)
    if state == 1
        return (a = 1, b = 1.), 2
    elseif state == 2
        return (a = 2, b = 2.), 3
    else
        return nothing
    end
end


text:
@testset "Triangulate" begin
    @testset "Polygon - 2D" begin
        for T in Floats
            quad = setup_Quadrilateral2(T)
            triangles = triangulate(quad)
            @test length(triangles) == 2
            @test triangles[1].vertices[1] ≈ quad[1]
            @test triangles[1].vertices[2] ≈ quad[2]
            @test triangles[1].vertices[3] ≈ quad[3]
            @test triangles[2].vertices[1] ≈ quad[1]
            @test triangles[2].vertices[2] ≈ quad[3]
            @test triangles[2].vertices[3] ≈ quad[4]
        end
    end

    @testset "Quadrilateral - 3D" begin
        for T in Floats
            quad = setup_Quadrilateral3(T)
            triangles = triangulate(quad, Val(0))
            @test length(triangles) == 2
            @test triangles[1].vertices[1] ≈ quad[1]
            @test triangles[1].vertices[2] ≈ quad[2]
            @test triangles[1].vertices[3] ≈ quad[3]
            @test triangles[2].vertices[1] ≈ quad[3]
            @test triangles[2].vertices[2] ≈ quad[4]
            @test triangles[2].vertices[3] ≈ quad[1]
        end
    end

    @testset "QuadraticTriangle" begin
        for T in Floats
            tri6 = setup_QuadraticTriangle3(T)
            triangles = triangulate(tri6, Val(0))
            @test length(triangles) == 1
            @test triangles[1].vertices[1] ≈ tri6[1]
            @test triangles[1].vertices[2] ≈ tri6[2]
            @test triangles[1].vertices[3] ≈ tri6[3]

            triangles = triangulate(tri6, Val(1))
            @test length(triangles) == 4
            @test triangles[1].vertices[1] ≈ tri6[6]
            @test triangles[1].vertices[2] ≈ tri6[5]
            @test triangles[1].vertices[3] ≈ tri6[3]
            @test triangles[2].vertices[1] ≈ tri6[6]
            @test triangles[2].vertices[2] ≈ tri6[4]
            @test triangles[2].vertices[3] ≈ tri6[5]
            @test triangles[3].vertices[1] ≈ tri6[1]
            @test triangles[3].vertices[2] ≈ tri6[4]
            @test triangles[3].vertices[3] ≈ tri6[6]
            @test triangles[4].vertices[1] ≈ tri6[4]
            @test triangles[4].vertices[2] ≈ tri6[2]
            @test triangles[4].vertices[3] ≈ tri6[5]
        end
    end

    @testset "QuadraticQuadrilateral" begin
        for T in Floats
            quad8 = setup_QuadraticQuadrilateral3(T)
            triangles = triangulate(quad8, Val(0))
            @test length(triangles) == 2
            @test triangles[1].vertices[1] ≈ quad8[1]
            @test triangles[1].vertices[2] ≈ quad8[2]
            @test triangles[1].vertices[3] ≈ quad8[3]
            @test triangles[2].vertices[1] ≈ quad8[3]
            @test triangles[2].vertices[2] ≈ quad8[4]
            @test triangles[2].vertices[3] ≈ quad8[1]

            triangles = triangulate(quad8, Val(1))
            @test length(triangles) == 8
        end
    end
end

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
module HelpersTests

using Test
using InplaceArrays.Helpers

for D = 0:5
  @test tfill(2,Val(D)) == tuple(fill(2,D)...)
end

end # module


text:
N = 10
d = 2
g, weights, points = euclidean_graph(N, d)
@test nv(g) == N
@test ne(g) == N*(N-1) ÷ 2
@test (d,N) == size(points)
@test maximum(x->x[2], weights) <= sqrt(d)
@test minimum(x->x[2], weights) >= 0
@test maximum(points) <= 1
@test minimum(points) >= 0.

g, weights, points = euclidean_graph(N, d, bc=:periodic)
@test maximum(x->x[2], weights) <= sqrt(d/2)
@test minimum(x->x[2], weights) >= 0.
@test maximum(points) <= 1
@test minimum(points) >= 0.

g, weights = euclidean_graph(points, L=0.01,  bc=:periodic)
@test nv(g) == N
@test ne(g) == N*(N-1) ÷ 2
@test maximum(x->x[2], weights) <= sqrt(d/2)

@test_throws ErrorException euclidean_graph(points, bc=:ciao)


text:
using ArgParse
include("iPredict.jl")
using .iPredict

function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table s begin
        "--mode"
            help = "Set to 'global' to do global simulations, set to 'gene' to simulate a single gene"
            required = true
            arg_type = String
        "--gene"
            help = "The gene to simulate in gene-level simulations"
            required = false
            arg_type = String
        "--countsMatrix"
            help = "Path to isoform-level counts matrix"
            required = true
            arg_type = String
        "--GeneIsoformRelationships"
            help = "Path to file containing gene-isoform relationships"
            required = true
            arg_type = String
        "--NumSimulations"
            help = "Number of simulations"
            required = true
            arg_type = Int64
        "--output"
            help="Path to save output at"
            required = true
            arg_type = String
        "--NumIsoformsToSimulate"
            help = "Max number of isoforms to simulate"
            required = false
            arg_type = Int64
        "--filteringThreshold"
            help = "All genes must express at least this number of isoforms"
            required = false
            arg_type = Int64
        "--alpha"
            help = "Beta distribution parameter alpha"
            required = false
            arg_type = Float64
        "--beta"
            help = "Beta distribution parameter beta"
            required = false
            arg_type = Float64
        "--pFP"
            required = false
            default = 0.01
            arg_type = Float64
        "--pFN"
            required = false
            default = 0.04
            arg_type = Float64

    end

    return parse_args(s)
end

function main()

    #Get command line args. Very basic atm, should add some argument validity
    #checking.
    parsed_args = parse_commandline()
    mode = parsed_args["mode"]
    countsMatrix = parsed_args["countsMatrix"]
    GeneIsoformRelationships = parsed_args["GeneIsoformRelationships"]
    NumSimulations = parsed_args["NumSimulations"]
    output = parsed_args["output"]
    NumIsoformsToSimulate = parsed_args["NumIsoformsToSimulate"]
    filteringThreshold = parsed_args["filteringThreshold"]
    alpha::Union{Nothing,Float64} = parsed_args["alpha"]
    beta::Union{Nothing,Float64} = parsed_args["beta"]
    pFP::Float64 = parsed_args["pFP"]
    pFN::Float64 = parsed_args["pFN"]

    if mode == "global"
        globalPredict(countsMatrix, GeneIsoformRelationships,
        NumSimulations, output, NumIsoformsToSimulate, filteringThreshold,
        false, alpha, beta, pFN, pFP)

    elseif mode == "gene"
        gene = parsed_args["gene"]
        oneGenePrediction(countsMatrix, GeneIsoformRelationships, gene,
        NumSimulations, output, false)

    else
        throw(ArgumentError("Unsupported mode - valid modes are gene or global.
        See help manual for details."))
    end
end

main()

max_stars_repo_name: jenni-westoby/Obstacles

max_issues_repo_name: jenni-westoby/Obstacles

max_forks_repo_name: jenni-westoby/Obstacles


text:
function AbstractConvergenceChecker(arg0::jdouble, arg1::jdouble)
    return AbstractConvergenceChecker((jdouble, jdouble), arg0, arg1)
end

function converged(obj::AbstractConvergenceChecker, arg0::jint, arg1::Object, arg2::Object)
    return jcall(obj, "converged", jboolean, (jint, Object, Object), arg0, arg1, arg2)
end

function equals(obj::Object, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function get_absolute_threshold(obj::AbstractConvergenceChecker)
    return jcall(obj, "getAbsoluteThreshold", jdouble, ())
end

function get_class(obj::Object)
    return jcall(obj, "getClass", Class, ())
end

function get_relative_threshold(obj::AbstractConvergenceChecker)
    return jcall(obj, "getRelativeThreshold", jdouble, ())
end

function hash_code(obj::Object)
    return jcall(obj, "hashCode", jint, ())
end

function notify(obj::Object)
    return jcall(obj, "notify", void, ())
end

function notify_all(obj::Object)
    return jcall(obj, "notifyAll", void, ())
end

function to_string(obj::Object)
    return jcall(obj, "toString", JString, ())
end

function wait(obj::Object)
    return jcall(obj, "wait", void, ())
end

function wait(obj::Object, arg0::jlong)
    return jcall(obj, "wait", void, (jlong,), arg0)
end

function wait(obj::Object, arg0::jlong, arg1::jint)
    return jcall(obj, "wait", void, (jlong, jint), arg0, arg1)
end


text:
# ---
# title: 257. Binary Tree Paths
# id: problem257
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Tree, Depth-first Search
# link: <https://leetcode.com/problems/binary-tree-paths/description/>
# hidden: true
# ---
#
# Given a binary tree, return all root-to-leaf paths.
#
# **Note:**  A leaf is a node with no children.
#
# **Example:**
#
#
#
#     Input:
#
#        1
#      /   \
#     2     3
#      \
#       5
#
#     Output: ["1->2->5", "1->3"]
#
#     Explanation: All root-to-leaf paths are: 1->2->5, 1->3
#
#
#
## @lc code=start
using LeetCode

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


text:
function equals(obj::Object, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function get_class(obj::Object)
    return jcall(obj, "getClass", Class, ())
end

function get_footprint(obj::FieldOfView, arg0::Transform, arg1::OneAxisEllipsoid, arg2::jdouble)
    return jcall(obj, "getFootprint", List, (Transform, OneAxisEllipsoid, jdouble), arg0, arg1, arg2)
end

function get_margin(obj::AbstractFieldOfView)
    return jcall(obj, "getMargin", jdouble, ())
end

function hash_code(obj::Object)
    return jcall(obj, "hashCode", jint, ())
end

function notify(obj::Object)
    return jcall(obj, "notify", void, ())
end

function notify_all(obj::Object)
    return jcall(obj, "notifyAll", void, ())
end

function offset_from_boundary(obj::FieldOfView, arg0::Vector3D, arg1::jdouble, arg2::VisibilityTrigger)
    return jcall(obj, "offsetFromBoundary", jdouble, (Vector3D, jdouble, VisibilityTrigger), arg0, arg1, arg2)
end

function project_to_boundary(obj::FieldOfView, arg0::Vector3D)
    return jcall(obj, "projectToBoundary", Vector3D, (Vector3D,), arg0)
end

function to_string(obj::Object)
    return jcall(obj, "toString", JString, ())
end

function wait(obj::Object)
    return jcall(obj, "wait", void, ())
end

function wait(obj::Object, arg0::jlong)
    return jcall(obj, "wait", void, (jlong,), arg0)
end

function wait(obj::Object, arg0::jlong, arg1::jint)
    return jcall(obj, "wait", void, (jlong, jint), arg0, arg1)
end


text:
#= Structure of iteration:
for i in iter
    # body
end

next = iterate(iter)
while next !== nothing
    (i, state) = next
    # body
    next = iterate(iter, state)
end
=#

struct Fibonacci{T<:Real}
end

function Fibonacci(d::DataType)
    if d<:Real
        return Fibonacci{d}()
    else
        error("No Real type!")
    end
end

function Base.iterate(::Fibonacci{T}) where {T<:Real}
    return (zero(T), (one(T), one(T)))
end

function Base.iterate(::Fibonacci{T}, state::Tuple{T, T}) where {T<:Real}
    return (state[1], (state[2], state[1] + state[2]))
end

# Alternatively:
# Fibonacci(d::DataType) = d<:Real ? Fibonacci{d}() : error("No Real type!")
# Base.iterate(::Fibonacci{T}) where {T<:Real} = (zero(T), (one(T), one(T)))
# Base.iterate(::Fibonacci{T}, state::Tuple{T, T}) where {T<:Real} = (state[1], (state[2], state[1] + state[2]))

for e in Fibonacci(Int64)
    e > 100 && break
    print(e, " ")
end


text:
using JSServe, WGLMakie, AbstractPlotting, Colors
using JSServe.DOM

function styled_slider(slider, value)
    rows(slider, DOM.span(value, class="p-1"), class="w-64 p-2 items-center")
end
columns(args...; class="") = DOM.div(args..., class=class * " flex flex-col")
rows(args...; class="") = DOM.div(args..., class=class * " flex flex-row")
css_color(c) = "background-color: #$(hex(c))"
using MakieGallery, FileIO
# get a somewhat interesting point cloud:
cat = decompose(Point3f, MakieGallery.loadasset("cat.obj"))
# Create a little interactive app
function App(session, request)
    markersize = JSServe.Slider(range(1, stop=20, step=0.1))
    hue_slider = JSServe.Slider(1:360)
    color = map(hue_slider) do hue
        HSV(hue, 0.5, 0.5)
    end
    plot = scatter(cat, markersize=markersize, color=color, resolution=(500, 500), show_axis=false)
    m_slider = styled_slider(markersize, markersize.value)
    color_swatch = DOM.div(class="h-6 w-6 p-1 rounded dropshadow", style=map(css_color, color))
    h_slider = styled_slider(hue_slider, color_swatch)
    sliders = rows(m_slider, h_slider)
    return DOM.div(JSServe.Styling, JSServe.TailwindCSS, columns(sliders, plot))
end
AppInteraktive(s, r) = JSServe.record_state_map(s, App(s, r)).dom

s = JSServe.Session()
states = AppInteraktive(s, nothing)
sizeof(MsgPack.pack(x)) / 10^6
mkdir("simple")
JSServe.export_standalone(AppInteraktive, "simple")

server = JSServe.Server(app, "0.0.0.0", 8081)
close(server)


text:
include("../src/datastructs_module.jl")

module dstest
using ..datastructs

# import datastructs
# import .datastructs

# tree = datastructs.BinaryTree{Int}(100)
# n1 = datastructs.BinaryTreeNode{Int}(120)
# n2 = datastructs.BinaryTreeNode{Int}(110)
# n3 = datastructs.BinaryTreeNode{Int}(90)
# n4 = datastructs.BinaryTreeNode{Int}(80)

# insert!(tree, n1)
# datastructs.inorder_tree_walk(tree)

function test_linked_list()

    #=
    n1 = LinkedListNode();
    n2 = LinkedListNode(10);
    n3 = LinkedListNode(10.0);
    n4 = LinkedListNode([1,2,3]);
    =#

    ll = datastructs.LinkedList{Int}(0);

    n1 = datastructs.LinkedListNode{Int}(10);
    n2 = datastructs.LinkedListNode{Int}(20);
    n3 = datastructs.LinkedListNode{Int}(30);
    n4 = datastructs.LinkedListNode{Int}(40);
    n5 = datastructs.LinkedListNode{Int}(50);

    push!(ll, n1)
    push!(ll, n2)
    push!(ll, n3)
    push!(ll, n4)
    push!(ll, n5)
    #=
    =#
end

function test_basics_should_not_allocate()
    stack = datastructs.Stack{Int}(10)

    push!(stack, 10) # allocation
    push!(stack, 20) # allocation
    push!(stack, 30)
    push!(stack, 40) # allocation
    push!(stack, 50)
end

function test_performance_push()
    stack = datastructs.Stack{Float64}()

    N = 100
    rnb = rand(N) .+ 0.2

    for rr in rnb
        if rr >= 0.5
            push!(stack, rr)
        end
    end

end

function test_performance_push_pop()
    stack = datastructs.Stack{Float64}()

    N = 100
    rnb = rand(N) .+ 0.2

    for rr in rnb
        if rr >= 0.5
            push!(stack, rr)
        else
            bla = pop!(stack)
        end
    end

end

function test_stack_push(rnb::Vector{Float64})
    stack = datastructs.Stack{Float64}()
    #N = 1e3
    #rnb = rand(N) .+ 0.2
    for rr in rnb
        push!(stack, rr)
    end
end

function test_stack_allocation(rnb::Vector{Float64})
    # N = 1e3
    stack = datastructs.Stack{Float64}(length(rnb))

    #rnb = rand(N) .+ 0.2
    for rr in rnb
        push!(stack, rr)
    end

end

function test_array_push(rnb::Vector{Float64})
    stack = Array{Float64, 1}()

    for rr in rnb
        push!(stack, rr)
    end
end

using BenchmarkTools

function benchmark_it()

    global N = 2^20-1
    rnb = rand(Int(N)) .+ 0.2

    println("=== test_array_push ===")
    show(stdout, MIME("text/plain"), @benchmark dstest.test_array_push(rand(Int(N)) .+ 0.2))
    println("\n\n")

    println("=== test_stack_push ===")
    show(stdout, MIME("text/plain"), @benchmark dstest.test_stack_push(rand(Int(N)) .+ 0.2))
    println("\n\n")

    println("=== test_stack_allocation ===")
    show(stdout, MIME("text/plain"), @benchmark dstest.test_stack_allocation(rand(Int(N)) .+ 0.2))
    println("\n\n")

end

end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
# Let X be a Partition of N
#	X::Array{Int, 1}
#	X[i] > 0 for all i
#	sum(X) == N
#	X[i] >= X[j] when i < j

# Parameters:
#	N::Int
#	- the problem size
# Return Values:
#	P::Array{Array{Array{Int, 1}, 1}, 1} (Partitions)
#	- P[n][p] contains the pth Partition of n
#	WI::Array{Int, 2} (Width Information)
#	- WI[n, w] contains the number of Paritions of n whose first element is less than or equal to w
function partitions(N::Int)
	WI = zeros(Int, N, N)
	for n = 1:N
		WI[n, 1] = 1
	end
	for w = 2:N
		WI[w, w] = 1
		for n = (w + 1):N
			num = 0
			for i = 1:min(n - w, w)
				num += WI[n - w, i]
			end
			WI[n, w] = num
		end
	end
	for n = 2:N
		for w = 2:n
			WI[n, w] += WI[n, w - 1]
		end
	end
	P = Array(Array{Array{Int, 1},1}, N)
	Pn = Array(Array{Int, 1}, 1)
	Pn[1] = [1]
	P[1] = Pn
	for n = 2:N
		Pn = Array(Array{Int, 1}, WI[n, n])
		Pn[1] = ones(Int, n)
		i = 2
		for w = 2:(n - 1)
			for spi = 1:WI[n - w, min(n - w, w)]
				Pn[i] = [w; P[n - w][spi]]
				i += 1
			end
		end
		Pn[i] = [n]
		P[n] = Pn
	end
	return P, WI
end


text:
"""
    CalculusWithJulia

A package to accompany [notes](calculuswithjulia.github.io) on using Julia with calculus.

This package does two things: 1) it loads several other packages making it easier to use (and install) the functionality provided by them and 2) defines a handful of functions and plot recipes for convenience.

## Packages loaded by `CalculusWithJulia`

### Built in packages.

The `MathConstants` package is reexported, allowing `e` to be used instead of `\euler[tab]` for a value of Euler's constant.

The `SpecialFunctions` is loaded giving access to a few special functions used in these notes, e.g., `airyai`.

The `LinearAlgebra` package is loaded for access to several of its functions for working with vectors `norm`, `cdot` (`⋅`), `cross` (`×`), `det`.

### Plotting

The `Plots` package is loaded giving access to `plot`, `scatter`, `annotate`, etc. The backends  used in the notes are `plotly`, `gr`, and `pyplot`.

The `ImplicitEquations` package is loaded for plotting implicitly defined functions.

In addition, several plot recipes are provided to ease the creation of plots:
`plotif`, `trimplot`, and `signchart` are used for plotting univariate functions;
`plot_polar` and `plot_parametric_curve` are used to plot curves in 2 or 3 dimensions;
`plot_parametric_surface` makes the plotting os parameterically defined surfaces easier;
`vectorfieldplot` and `vectorfieldplot3d` can be used to plot vector fields; and
`arrow` is a simplified interface to `quiver` that also indicates 3D vectors.

### Symbolic math

The `Sympy` package is loaded for symbolic mathematics

### Zeros of functions

The `Roots` package is loaded for its `fzero` and `fzeros` functions.

### Derivatives

The `ForwardDiff` package is loaded giving access to its  `derivative`,  `gradient`, `jacobian`, and `hessian` functions for finding automatic derivatives of functions. In addition, this package defines `'` (for functions) to return a derivative (which commits [type piracy](https://docs.julialang.org/en/v1/manual/style-guide/index.html#Avoid-type-piracy-1)), `∇` to find the gradient (`∇(f)`), the divergence (`∇⋅F`). and the curl (`∇×F`), along with `divergence` and `curl`.

### Integration

This package reexports `QuadGK`, for one-dimensional integrals; `HCubature`, for multidimensional integrals; and provides `riemann`, for Riemann sums, and `fubini`, for iterated integrals over non-rectangular regions based on `QuadGK`.

"""
module CalculusWithJulia

using Reexport
@reexport using Plots
using RecipesBase
import ImplicitEquations # handle conflict with SymPy in plot-utils.jl
import ImplicitEquations: Pred
@reexport using LinearAlgebra
@reexport using Base.MathConstants
@reexport using SpecialFunctions
@reexport using Roots
@reexport using SymPy
import ForwardDiff
export ForwardDiff
@reexport using QuadGK
@reexport using HCubature

include("multidimensional.jl")
include("derivatives.jl")
include("integration.jl")
include("plot-utils.jl")

export unzip, divergence, gradient, curl, ∇
export tangent, secant, D
export riemann, fubini
export plotif, trimplot, signchart,
       plot_polar, plot_polar!,
       plot_parametric_curve,   plot_parametric_curve!,
       plot_parametric_surface, plot_parametric_surface!,
       vectorfieldplot,   vectorfieldplot!,
       vectorfieldplot3d, vectorfieldplot3d,
       arrow, arrow!

end # module


text:
module TestConstructors

using Test, DataFrames
using DataFrames: Index, _columns, index
const ≅ = isequal

#
# DataFrame
#
@testset "constructors" begin
    df = DataFrame()
    @inferred DataFrame()

    @test isempty(_columns(df))
    @test _columns(df) isa Vector{AbstractVector}
    @test index(df) == Index()
    @test size(DataFrame!()) == (0,0)

    df = DataFrame(Any[CategoricalVector{Union{Float64, Missing}}(zeros(3)),
                   CategoricalVector{Union{Float64, Missing}}(ones(3))],
                   Index([:x1, :x2]))
    @test size(df, 1) == 3
    @test size(df, 2) == 2

    df2 = DataFrame!(Any[CategoricalVector{Union{Float64, Missing}}(zeros(3)),
                     CategoricalVector{Union{Float64, Missing}}(ones(3))],
                     Index([:x1, :x2]))
    @test size(df2, 1) == 3
    @test size(df2, 2) == 2

    @test df == DataFrame([CategoricalVector{Union{Float64, Missing}}(zeros(3)),
                           CategoricalVector{Union{Float64, Missing}}(ones(3))])
    @test df == DataFrame!([CategoricalVector{Union{Float64, Missing}}(zeros(3)),
                            CategoricalVector{Union{Float64, Missing}}(ones(3))])
    @test df == DataFrame([CategoricalVector{Union{Float64, Missing}}(zeros(3)),
                           CategoricalVector{Union{Float64, Missing}}(ones(3))], [:x1, :x2])
    @test df == DataFrame(Any[CategoricalVector{Union{Float64, Missing}}(zeros(3)),
                              CategoricalVector{Union{Float64, Missing}}(ones(3))])
    @test df == DataFrame(Any[CategoricalVector{Union{Float64, Missing}}(zeros(3)),
                              CategoricalVector{Union{Float64, Missing}}(ones(3))], [:x1, :x2])
    @test df == DataFrame(AbstractVector[CategoricalVector{Union{Float64, Missing}}(zeros(3)),
                                         CategoricalVector{Union{Float64, Missing}}(ones(3))], [:x1, :x2])
    @test df == DataFrame((CategoricalVector{Union{Float64, Missing}}(zeros(3)),
                           CategoricalVector{Union{Float64, Missing}}(ones(3))))
    @test df == DataFrame((CategoricalVector{Union{Float64, Missing}}(zeros(3)),
                           CategoricalVector{Union{Float64, Missing}}(ones(3))), (:x1, :x2))
    @test df == DataFrame(x1 = Union{Int, Missing}[0.0, 0.0, 0.0],
                          x2 = Union{Int, Missing}[1.0, 1.0, 1.0])
    @test df == DataFrame!(x1 = Union{Int, Missing}[0.0, 0.0, 0.0],
                           x2 = Union{Int, Missing}[1.0, 1.0, 1.0])
    @test df == DataFrame([:x1=>Union{Int, Missing}[0.0, 0.0, 0.0],
                           :x2=>Union{Int, Missing}[1.0, 1.0, 1.0]])
    @test df == DataFrame!((:x1=>Union{Int, Missing}[0.0, 0.0, 0.0],
                            :x2=>Union{Int, Missing}[1.0, 1.0, 1.0]))

    @test DataFrame([1:3, 1:3]) == DataFrame(Any[1:3, 1:3]) ==
          DataFrame(UnitRange[1:3, 1:3]) == DataFrame(AbstractVector[1:3, 1:3]) ==
          DataFrame([[1,2,3], [1,2,3]]) == DataFrame(Any[[1,2,3], [1,2,3]]) ==
          DataFrame(([1,2,3], [1,2,3])) == DataFrame((1:3, 1:3)) ==
          DataFrame((1:3, [1,2,3])) == DataFrame([1:3, [1,2,3]])
          DataFrame((:x1=>1:3, :x2=>[1,2,3])) == DataFrame([:x1=>1:3, :x2=>[1,2,3]])

    @inferred DataFrame([1:3, 1:3])
    @inferred DataFrame((1:3, 1:3))
    @inferred DataFrame([1:3, 1:3], [:a, :b])
    @inferred DataFrame((1:3, 1:3), (:a, :b))

    @inferred DataFrame((:x1=>1:3, :x2=>[1,2,3]))
    @inferred DataFrame([:x1=>1:3, :x2=>[1,2,3]])

    @test df !== DataFrame(df)
    @test df == DataFrame(df)

    df2 = convert(DataFrame, Union{Float64, Missing}[0.0 1.0;
                                                     0.0 1.0;
                                                     0.0 1.0])
    names!(df2, [:x1, :x2])
    @test df[!, :x1] == df2[!, :x1]
    @test df[!, :x2] == df2[!, :x2]

    df2 = DataFrame([0.0 1.0;
                     0.0 1.0;
                     0.0 1.0])
    names!(df2, [:x1, :x2])
    @test df[!, :x1] == df2[!, :x1]
    @test df[!, :x2] == df2[!, :x2]

    @test_throws ArgumentError DataFrame!([0.0 1.0;
                                           0.0 1.0;
                                           0.0 1.0])

    df2 = DataFrame([0.0 1.0;
                     0.0 1.0;
                     0.0 1.0], [:a, :b])
    names!(df2, [:a, :b])
    @test df[!, :x1] == df2[!, :a]
    @test df[!, :x2] == df2[!, :b]

    @test_throws ArgumentError DataFrame!([0.0 1.0;
                                           0.0 1.0;
                                           0.0 1.0], [:a, :b])

    @test df == DataFrame(x1 = Union{Float64, Missing}[0.0, 0.0, 0.0],
                          x2 = Union{Float64, Missing}[1.0, 1.0, 1.0])
    @test df == DataFrame(x1 = Union{Float64, Missing}[0.0, 0.0, 0.0],
                          x2 = Union{Float64, Missing}[1.0, 1.0, 1.0],
                          x3 = Union{Float64, Missing}[2.0, 2.0, 2.0])[:, [:x1, :x2]]

    @test_throws BoundsError SubDataFrame(DataFrame(A=1), 0:0, :)
    @test_throws ArgumentError SubDataFrame(DataFrame(A=1), 0, :)
    @test_throws BoundsError DataFrame(A=1)[0, :]
    @test_throws BoundsError DataFrame(A=1)[0, 1:1]
    @test SubDataFrame(DataFrame(A=1), 1:1, :) == DataFrame(A=1)
    @test SubDataFrame(DataFrame(A=1:10), 1:4, :) == DataFrame(A=1:4)
    @test view(SubDataFrame(DataFrame(A=1:10), 1:4, :), 2:2, :) == DataFrame(A=2)
    @test view(SubDataFrame(DataFrame(A=1:10), 1:4, :), [true, true, false, false], :) == DataFrame(A=1:2)

    @test DataFrame(a=1, b=1:2) == DataFrame(a=[1,1], b=[1,2])
end

@testset "DataFrame keyword argument constructor" begin
    x = [1,2,3]
    y = [4,5,6]

    df = DataFrame(x=x, y=y)
    @test size(df) == (3, 2)
    @test names(df) == [:x, :y]
    @test df.x == x
    @test df.y == y
    @test df.x !== x
    @test df.y !== y
    df = DataFrame(x=x, y=y, copycols=true)
    @test size(df) == (3, 2)
    @test names(df) == [:x, :y]
    @test df.x == x
    @test df.y == y
    @test df.x !== x
    @test df.y !== y
    df = DataFrame(x=x, y=y, copycols=false)
    @test size(df) == (3, 2)
    @test names(df) == [:x, :y]
    @test df.x === x
    @test df.y === y
    @test_throws ArgumentError DataFrame(x=x, y=y, copycols=1)

    df = DataFrame!(x=x, y=y)
    @test size(df) == (3, 2)
    @test names(df) == [:x, :y]
    @test df.x === x
    @test df.y === y
    @test_throws ArgumentError DataFrame!(x=x, y=y, copycols=true)
end

@testset "DataFrame constructor" begin
    df1 = DataFrame(x=1:3, y=1:3)

    df2 = DataFrame(df1)
    df3 = copy(df1)
    @test df1 == df2 == df3
    @test df1.x !== df2.x
    @test df1.x !== df3.x
    @test df1.y !== df2.y
    @test df1.y !== df3.y

    df2 = DataFrame(df1, copycols=false)
    df3 = copy(df1, copycols=false)
    @test df1 == df2 == df3
    @test df1.x === df2.x
    @test df1.x === df3.x
    @test df1.y === df2.y
    @test df1.y === df3.y

    df1 = view(df1, :, :)
    df2 = DataFrame(df1)
    df3 = copy(df1)
    @test df1 == df2 == df3
    @test df1.x !== df2.x
    @test df1.x !== df3.x
    @test df1.y !== df2.y
    @test df1.y !== df3.y

    df2 = DataFrame(df1, copycols=false)
    @test df1 == df2
    @test df1.x === df2.x
    @test df1.y === df2.y

    df2 = DataFrame!(df1)
    @test df1 == df2
    @test df1.x === df2.x
    @test df1.y === df2.y
end

@testset "pair constructor" begin
    df = DataFrame(:x1 => zeros(3), :x2 => ones(3))
    @inferred DataFrame(:x1 => zeros(3), :x2 => ones(3))
    @test size(df, 1) == 3
    @test size(df, 2) == 2
    @test isequal(df, DataFrame(x1 = [0.0, 0.0, 0.0], x2 = [1.0, 1.0, 1.0]))

    df = DataFrame(:type => [], :begin => [])
    @test names(df) == [:type, :begin]

    a=[1,2,3]
    df = DataFrame(:a=>a, :b=>1, :c=>1:3)
    @test names(df) == [:a, :b, :c]
    @test df.a == a
    @test df.a !== a
    df = DataFrame(:a=>a, :b=>1, :c=>1:3, copycols=false)
    @test names(df) == [:a, :b, :c]
    @test df.a === a

    df = DataFrame!(:a=>a, :b=>1, :c=>1:3)
    @test names(df) == [:a, :b, :c]
    @test df.a === a
end

@testset "associative" begin
    df = DataFrame(Dict(:A => 1:3, :B => 4:6))
    @inferred DataFrame(Dict(:A => 1:3, :B => 4:6))
    @test df == DataFrame(A = 1:3, B = 4:6)
    @test eltypes(df) == [Int, Int]

    a=[1,2,3]
    df = DataFrame(Dict(:a=>a, :b=>1, :c=>1:3))
    @test names(df) == [:a, :b, :c]
    @test df.a == a
    @test df.a !== a
    df = DataFrame(Dict(:a=>a, :b=>1, :c=>1:3), copycols=false)
    @test names(df) == [:a, :b, :c]
    @test df.a === a

    df = DataFrame!(Dict(:a=>a, :b=>1, :c=>1:3))
    @test names(df) == [:a, :b, :c]
    @test df.a === a
end

@testset "vector constructors" begin
    x = [1,2,3]
    y = [1,2,3]

    df = DataFrame([x, y])
    @test names(df) == [:x1, :x2]
    @test df.x1 == x
    @test df.x2 == y
    @test df.x1 !== x
    @test df.x2 !== y
    df = DataFrame([x, y], copycols=true)
    @test names(df) == [:x1, :x2]
    @test df.x1 == x
    @test df.x2 == y
    @test df.x1 !== x
    @test df.x2 !== y
    df = DataFrame([x, y], copycols=false)
    @test names(df) == [:x1, :x2]
    @test df.x1 === x
    @test df.x2 === y

    df = DataFrame([x, y], [:x1, :x2])
    @test names(df) == [:x1, :x2]
    @test df.x1 == x
    @test df.x2 == y
    @test df.x1 !== x
    @test df.x2 !== y
    df = DataFrame([x, y], [:x1, :x2], copycols=true)
    @test names(df) == [:x1, :x2]
    @test df.x1 == x
    @test df.x2 == y
    @test df.x1 !== x
    @test df.x2 !== y
    df = DataFrame([x, y], [:x1, :x2], copycols=false)
    @test names(df) == [:x1, :x2]
    @test df.x1 === x
    @test df.x2 === y

    df = DataFrame((x, y))
    @test names(df) == [:x1, :x2]
    @test df.x1 == x
    @test df.x2 == y
    @test df.x1 !== x
    @test df.x2 !== y
    df = DataFrame((x, y), copycols=true)
    @test names(df) == [:x1, :x2]
    @test df.x1 == x
    @test df.x2 == y
    @test df.x1 !== x
    @test df.x2 !== y
    df = DataFrame((x, y), copycols=false)
    @test names(df) == [:x1, :x2]
    @test df.x1 === x
    @test df.x2 === y

    df = DataFrame!((x, y))
    @test names(df) == [:x1, :x2]
    @test df.x1 === x
    @test df.x2 === y

    df = DataFrame((x, y), (:x1, :x2))
    @test names(df) == [:x1, :x2]
    @test df.x1 == x
    @test df.x2 == y
    @test df.x1 !== x
    @test df.x2 !== y
    df = DataFrame((x, y), (:x1, :x2), copycols=true)
    @test names(df) == [:x1, :x2]
    @test df.x1 == x
    @test df.x2 == y
    @test df.x1 !== x
    @test df.x2 !== y
    df = DataFrame((x, y), (:x1, :x2), copycols=false)
    @test names(df) == [:x1, :x2]
    @test df.x1 === x
    @test df.x2 === y

    df = DataFrame!((x, y), (:x1, :x2))
    @test names(df) == [:x1, :x2]
    @test df.x1 === x
    @test df.x2 === y
end

@testset "recyclers" begin
    @test DataFrame(a = 1:5, b = 1) == DataFrame(a = collect(1:5), b = fill(1, 5))
    @test DataFrame(a = 1, b = 1:5) == DataFrame(a = fill(1, 5), b = collect(1:5))
end

@testset "constructor errors" begin
    @test_throws DimensionMismatch DataFrame(a=1, b=[])
    @test_throws DimensionMismatch DataFrame(Any[collect(1:10)], DataFrames.Index([:A, :B]))
    @test_throws DimensionMismatch DataFrame(A = rand(2,2))
    @test_throws DimensionMismatch DataFrame(A = rand(2,1))
    @test_throws ArgumentError DataFrame([1, 2, 3])

    @test_throws DimensionMismatch DataFrame!(a=1, b=[])
    @test_throws DimensionMismatch DataFrame!(Any[collect(1:10)], DataFrames.Index([:A, :B]))
    @test_throws DimensionMismatch DataFrame!(A = rand(2,2))
    @test_throws DimensionMismatch DataFrame!(A = rand(2,1))
    @test_throws ArgumentError DataFrame!([1, 2, 3])
    @test_throws MethodError DataFrame!([1 2; 3 4], copycols=false)
end

@testset "column types" begin
    df = DataFrame(A = 1:3, B = 2:4, C = 3:5)
    answer = [Array{Int,1}, Array{Int,1}, Array{Int,1}]
    @test map(typeof, eachcol(df)) == answer
    df[!, :D] = [4, 5, missing]
    push!(answer, Vector{Union{Int, Missing}})
    @test map(typeof, eachcol(df)) == answer
    df[!, :E] .= 'c'
    push!(answer, Vector{Char})
    @test map(typeof, eachcol(df)) == answer
end

@testset "Matrix constructor" begin
    df = DataFrame([1 2; 3 4])
    @test size(df) == (2, 2)
    @test df.x1 == [1, 3]
    @test df.x2 == [2, 4]
    @test_throws ArgumentError DataFrame!([1 2; 3 4])

end

@testset "constructor with types" begin
    df = DataFrame([Union{Int, Missing}, Union{Float64, Missing}, Union{String, Missing}],
                   [:A, :B, :C], 100)
    @test size(df, 1) == 100
    @test size(df, 2) == 3
    @test typeof(df[!, 1]) == Vector{Union{Int, Missing}}
    @test typeof(df[!, 2]) == Vector{Union{Float64, Missing}}
    @test typeof(df[!, 3]) == Vector{Union{String, Missing}}
    @test all(ismissing, df[!, 1])
    @test all(ismissing, df[!, 2])
    @test all(ismissing, df[!, 3])

    df = DataFrame([Union{Int, Missing}, Union{Float64, Missing}], [:x1, :x2], 2)
    @test size(df) == (2, 2)
    @test eltypes(df) == [Union{Int, Missing}, Union{Float64, Missing}]

    @test_throws ArgumentError DataFrame!([Union{Int, Missing}, Union{Float64, Missing}],
                                          [:x1, :x2], 2)
    @test size(df) == (2, 2)
    @test eltypes(df) == [Union{Int, Missing}, Union{Float64, Missing}]

    df = DataFrame([Union{Int, Missing}, Union{Float64, Missing}, Union{String, Missing}],
                   [:A, :B, :C])
    @test size(df, 1) == 0
    @test size(df, 2) == 3
    @test typeof(df[!, 1]) == Vector{Union{Int, Missing}}
    @test typeof(df[!, 2]) == Vector{Union{Float64, Missing}}
    @test typeof(df[!, 3]) == Vector{Union{String, Missing}}
    @test names(df) == [:A, :B, :C]
end

end # module