text:
"""
Fast L2-regularized Linear Regression

Kei Majima developed the original algorithm.
Soma Nonaka implemented the initial Julia version.
Shuntaro C. Aoki modified and refactored the original implementation.
"""

module FastL2LiR

export fit, predict, load_model

using Statistics
using LinearAlgebra
using Glob
using BrainDecoder

struct FastL2LiRModel{T<:Real}
    W::Array{T}
    b::Array{T}
end

function fit(X::Array{T, 2}, Y::Array{T}, alpha::T)::FastL2LiRModel where {T<:Real}
    return _fit(X, Y, alpha, false, 0)
end

function fit(X::Array{T, 2}, Y::Array{T}, alpha::T, num_features::Integer)::FastL2LiRModel where {T<:Real}
    return _fit(X, Y, alpha, true, num_features)
end

function _fit(X::Array{T, 2}, Y::Array{T}, alpha::T, feature_selection::Bool, num_features::Integer)::FastL2LiRModel where {T<:Real}
    """
    Fit L2-regularized linear regression model.
    """

    reshape_y::Bool = false
    if ndims(Y) > 2
        reshape_y = true
        y_shape = size(Y)[2:end]
        Y = reshape(Y, size(Y, 1), :)
    end

    n_sample = size(X, 1)
    n_feat = size(X, 2)
    n_trg = size(Y, 2)

    if n_feat < num_features
        feature_selection = false
    end

    # Append bias term in X
    Xb = [X ones(T, n_sample)]

    if feature_selection
        # With feature selection
        W = zeros(T, n_feat, n_trg)
        b = zeros(T, n_trg)

        C = cor(X, Y, dims=1)
        C = @. abs(C)
        C[C .== NaN] .= 0

        W0 = Xb' * Xb + alpha * I(n_feat + 1)
        W1 = Xb' * Y

        # @inbounds causes seg fault with MKL-linked Julia
        for i = 1:n_trg
            C0 = C[:, i]
            feature_index = sortperm(C0, rev=true)[1:num_features]
            push!(feature_index, size(Xb, 2))  # バイアス項のインデックスを追加する

            try
                Wb = W0[feature_index, feature_index] \ W1[feature_index, i]
            catch
                println("Singular matrix: use QR decomposition")
                Q, R = qr(W0)
                Wb = inv(R) * (Q' * W1[feature_index, i])
            end

            W[feature_index[1:end - 1], i] = Wb[1:end - 1, 1]
            b[i] = Wb[end, 1]
        end
    else
        # Without feature selection
        Wb = (Xb' * Xb + alpha * I(n_feat + 1)) \ (Xb' * Y)
        W = Wb[1:end - 1, :]
        b = Wb[end, :]
    end

    # Rerurns b as a row vector (compat. for PyFastL2LiR)
    b = collect(b')

    if reshape_y
        W = reshape(W, n_feat, y_shape...)
        b = reshape(b, 1, y_shape...)
    end

    return FastL2LiRModel(W, b)
end

function predict(model::FastL2LiRModel, X::Array{T, 2})::Array{T} where {T<:Real}
    """
    Predict y with given x and L2-regularized linear regression model.
    """

    reshape_y::Bool = false
    if ndims(model.W) > 2
        reshape_y = true
        w_shape = size(model.W)[2:end]
        W = reshape(model.W, size(model.W, 1), :)
        b = reshape(model.b, 1, :)
    else
        W = model.W
        b = model.b
    end

    Y_pred = X * W + repeat(b, size(X, 1), 1)

    if reshape_y
        Y_pred = reshape(Y_pred, size(Y_pred, 1), w_shape...)
    end

    return Y_pred
end

function load_model(model_dir::String; chunk_axis=0)::FastL2LiRModel
    """
    Load FastL2LiR model.
    """

    function _load_chunks(dir, name)
        files = sorted(glob(joinpath(dir, name, "*.mat")))
        return cat([load_array(f) for f in files]..., dims=chunk_axis)
    end

    # Load model
    if isfile(joinpath(model_dir, "W.mat"))
        W = load_array(joinpath(model_dir, "W.mat"), "W")
    else
        W = _load_chunks(model_dir, "W")
    end

    if isfile(joinpath(model_dir, "b.mat"))
        b = load_array(joinpath(model_dir, "b.mat"), "b")
    else
        b = _load_chunks(model_dir, "b")
    end

    return FastL2LiRModel(W, b)
end

end # module FastL2LiR


text:
@testset "Testing branch length distance LS score" begin

# on a tree:
net = readTopology("(A,(B,(C,D)));")
net2 = readTopology("(((A,B),C),D);")
divergenceTimes = Dict(1 => 0.0, 2 => 0.0, 3 => 0.0, 4 => 0.0, -4 => 10.0, -3 => 15.0, -2 => 25.0)
setBLfromDivergenceTimes!(net2, divergenceTimes)
dnaDistances = pairwiseTaxonDistanceMatrix(net2)
distanceMismatch = pairwiseDistanceLSscore(net, divergenceTimes, dnaDistances)
@test distanceMismatch==130.7669683062202

# on a network:
net = readTopology("(((A,(B)#H1:::0.9),(C,#H1:::0.1)),D);")
net2 = readTopology("(((D,(C)#H1:::0.9,(B,#H1:::0.1))),A);")
divergenceTimes = Dict(1 => 0.0, 2 => 0.0, 4 => 0.0, 5 => 0.0, -4 => 10.0, -3 => 15.0, -2 => 25.0, -6 => 7.0, 3 => 5.0)
setBLfromDivergenceTimes!(net2, divergenceTimes)
dnaDistances = pairwiseTaxonDistanceMatrix(net2)
distanceMismatch = pairwiseDistanceLSscore(net, divergenceTimes, dnaDistances)
@test distanceMismatch==63.75829357816911

end

ext: jl


text:
################################################################################
#  Copyright 2022, Tom Van Acker (BASF)                                        #
################################################################################
# TFT.jl                                                                       #
# A Julia package for Taylor-Fourier Transform.                                #
# See http://github.com/timmyfaraday/TFT.jl                                    #
################################################################################

"""
    TFT.sample_ospline(D::Int, K::Int, N::Int)

Function to obtain samples of the up-to-Dth-degree derivatives of Kth-degree
o-spline.

Input:
- `D::Int`  | maximum degree of the derivative [-]
- `K::Int`  | degree of the o-spline [-]
- `N::Int`  | number of samples of the fundamental cycle [-]

Output
- `Φ::Matrix{<:Real}`   | samples of the up-to-Dth-degree derivatives of the
                        | Kth-degree o-spline
"""
function sample_ospline(D::Int, K::Int, N::Int)
    # define the normalized time range `rU`
    ΔU  = 1 / N
    rU  = 0.0:ΔU:(1.0-ΔU)
    # define the knots range `rK`
    rK  = [-K:-1..., 1:K...]

    # initialize the Φ-matrix
    Φ = ones((K+1) * N, D+1)

    # fill the Φ-matrices
    for ctr in 0:K
        # one-based counter `cnt` wrt zero-based counter `ctr`
        ctn = ctr + 1
        # range of the first dimension of Φ wrt zero-based counter `ctr`
        rΦ  = (N * ctr + 1):(N * (ctr + 1))
        # base o-spline
        u   = -(K + 1) / 2 + ctr .+ rU
        ψ   = [1, -rK[ctn]] ./ -rK[ctn]
        for nk in 1:K
            Φ[rΦ, 1] .*= (u .- rK[ctr + nk]) ./ -rK[ctr + nk]
            if nk > 1
                ψ   = _DSP.conv(ψ, [1, -rK[ctr + nk]] ./ -rK[ctr + nk])
            end
        end
        # derivative by Horner scheme
        # NB: the Julia Polynomials pkg, takes its coefficients in reverse,
        # hence the need for the double reserve when determining ψ.
        for nd in 2:(D+1)
            ψ   = reverse(_POL.derivative(_POL.Polynomial(reverse(ψ))).coeffs)
            Φ[rΦ, nd] .*= ψ[1]
            for nk in 2:(K-nd+2)
                Φ[rΦ, nd] .*= u
                Φ[rΦ, nd] .+= ψ[nk]
            end
        end
    end

    # return the samples of the o-spline and its derivatives: Φ
    return Φ
end


text:
module M

export f

f(x) = x

g(x = 1, y = 2, z = 3; kwargs...) = x

h(x::Int, y::Int = 2, z::Int = 3; kwargs...) = x

const A{T} = Union{Array{T, 3}, Array{T, 4}}

h_1(x::A) = x
h_2(x::A{Int}) = x
h_3(x::A{T}) where {T} = x
h_4(x, ::Int, z) = x

i_1(x; y = x) = x * y
i_2(x::Int; y = x) = x * y
i_3(x::T; y = x) where {T} = x * y
i_4(x; y::T = zero(T), z::U = zero(U)) where {T, U} = x + y + z

j_1(x, y) = x * y # two arguments, no keyword arguments
j_1(x; y = x) = x * y # one argument, one keyword argument

k_1(x::String, y::T = 0, z::T = zero(T)) where T <: Number = x
k_2(x::String, y::U, z::T) where T <: Number where U <: Complex = x
k_3(x, y::T, z::U) where {T, U} = x + y + z
k_4(::String, ::Int = 0) = nothing
k_5(::Type{T}, x::String, func::Union{Nothing, Function} = nothing) where T <: Number = x
k_6(x::Vector{T}) where T <: Number = x
k_7(x::Union{T,Nothing}, y::T = zero(T)) where {T <: Integer} = x
k_8(x) = x
k_9(x::T where T<:Any) = x
k_10(x::T) where T = x
k_11(x::Int, xs...) = x
k_12(x::Int, xs::Real...) = x

mutable struct T
    a
    b
    c
end

struct K
    K(; a = 1) = new()
end

abstract type AbstractType1 <: Integer end
abstract type AbstractType2{S, T <: Integer} <: Integer end

struct CustomType{S, T <: Integer} <: Integer
end

primitive type BitType8 8 end

primitive type BitType32 <: Real 32 end

end


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

mutable struct IoK8sApiCoreV1NodeSpec <: SwaggerModel
    configSource::Any # spec type: Union{ Nothing, IoK8sApiCoreV1NodeConfigSource } # spec name: configSource
    externalID::Any # spec type: Union{ Nothing, String } # spec name: externalID
    podCIDR::Any # spec type: Union{ Nothing, String } # spec name: podCIDR
    providerID::Any # spec type: Union{ Nothing, String } # spec name: providerID
    taints::Any # spec type: Union{ Nothing, Vector{IoK8sApiCoreV1Taint} } # spec name: taints
    unschedulable::Any # spec type: Union{ Nothing, Bool } # spec name: unschedulable

    function IoK8sApiCoreV1NodeSpec(;configSource=nothing, externalID=nothing, podCIDR=nothing, providerID=nothing, taints=nothing, unschedulable=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1NodeSpec, Symbol("configSource"), configSource)
        setfield!(o, Symbol("configSource"), configSource)
        validate_property(IoK8sApiCoreV1NodeSpec, Symbol("externalID"), externalID)
        setfield!(o, Symbol("externalID"), externalID)
        validate_property(IoK8sApiCoreV1NodeSpec, Symbol("podCIDR"), podCIDR)
        setfield!(o, Symbol("podCIDR"), podCIDR)
        validate_property(IoK8sApiCoreV1NodeSpec, Symbol("providerID"), providerID)
        setfield!(o, Symbol("providerID"), providerID)
        validate_property(IoK8sApiCoreV1NodeSpec, Symbol("taints"), taints)
        setfield!(o, Symbol("taints"), taints)
        validate_property(IoK8sApiCoreV1NodeSpec, Symbol("unschedulable"), unschedulable)
        setfield!(o, Symbol("unschedulable"), unschedulable)
        o
    end
end # type IoK8sApiCoreV1NodeSpec

const _property_map_IoK8sApiCoreV1NodeSpec = Dict{Symbol,Symbol}(Symbol("configSource")=>Symbol("configSource"), Symbol("externalID")=>Symbol("externalID"), Symbol("podCIDR")=>Symbol("podCIDR"), Symbol("providerID")=>Symbol("providerID"), Symbol("taints")=>Symbol("taints"), Symbol("unschedulable")=>Symbol("unschedulable"))
const _property_types_IoK8sApiCoreV1NodeSpec = Dict{Symbol,String}(Symbol("configSource")=>"IoK8sApiCoreV1NodeConfigSource", Symbol("externalID")=>"String", Symbol("podCIDR")=>"String", Symbol("providerID")=>"String", Symbol("taints")=>"Vector{IoK8sApiCoreV1Taint}", Symbol("unschedulable")=>"Bool")
Base.propertynames(::Type{ IoK8sApiCoreV1NodeSpec }) = collect(keys(_property_map_IoK8sApiCoreV1NodeSpec))
Swagger.property_type(::Type{ IoK8sApiCoreV1NodeSpec }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApiCoreV1NodeSpec[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1NodeSpec }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1NodeSpec[property_name]

function check_required(o::IoK8sApiCoreV1NodeSpec)
    true
end

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


text:
const helix = let res_init=zeros(3), jac_init=zeros(3,3), x_init=[-1.0, 0.0, 0.0]

    function res(x, r)
        theta = atan(x[2], x[1]) / 2pi

        if x[1] <= 0.0 && x[2] <= 0.0
            theta = theta + 1.0
        end

        r[1] = 10.0 * (x[3] - 10.0 * theta )
        r[2] = 10.0 * (sqrt(x[1]^2 + x[2]^2) - 1.0)
        r[3] = x[3]
        return r
    end

    function jac(x, jac)
        t = x[1]^2 + x[2]^2
        ti = 100.0 / (2pi * t )
        jac[1, 1] = ti * x[2]
        t = 10.0 / sqrt(t)
        jac[2, 1] = x[1] * t
        jac[3, 1] = 0.0
        jac[1, 2] = -ti * x[1]
        jac[2, 2] = x[2] * t
        jac[3, 2] = 0.0
        jac[1, 3] = 10.0
        jac[2, 3] = 0.0
        jac[3, 3] = 1.0
        return jac
    end

    f(;scale=1, verbose=false, print_steps=false) = testone("helix", res, jac,
                                                            res_init, jac_init, x_init;
                                                            scale=scale, verbose=verbose,
                                                            print_steps=print_steps)
end


text:
export zoo_dc_file, zoo_discriminative, zoo_dc

struct  DiscriminativeFormat <: FileFormat end

const DiscriminativeVtreeFormat = Tuple{DiscriminativeFormat, VtreeFormat}
Tuple{DiscriminativeVtreeFormat,VtreeFormat}() = (DiscriminativeFormat(), VtreeFormat())

zoo_dc_file(name) =
    artifact"circuit_model_zoo" * zoo_version * "/lcs/$name"

"""
    zoo_dc(name)

Loads a discriminative circuit (logistic circuit or regression circuit) from the model zoo.
See https://github.com/UCLA-StarAI/Circuit-Model-Zoo.
"""
zoo_dc(name) =
    read(zoo_dc_file(name), DiscriminativeCircuit, DiscriminativeFormat())

const zoo_discriminative = zoo_dc

const discriminative_grammer = raw"""
    start: _header (_NL node)+ _NL?

    _header : "Logisitic" _WS "Circuit" | "Logistic" _WS "Circuit" | "Regression" _WS "Circuit"

    node: "T" _WS INT _WS INT _WS INT _WS params -> true_node
      | "F" _WS INT _WS INT _WS INT _WS params -> false_node
      | "D" _WS INT _WS INT _WS INT _WS elems -> or_node
      | "B" _WS params -> bias_node

    elem: "(" INT _WS INT _WS params ")"

    elems: elem (_WS elem)*

    params: LOGPROB (_WS LOGPROB)*

    %import common.INT
    %import common.SIGNED_INT
    %import common.SIGNED_NUMBER -> LOGPROB
    %import common.WS_INLINE -> _WS
    %import common.NEWLINE -> _NL
    """ * dimacs_comments

const discriminative_parser = Lark(discriminative_grammer)

abstract type DiscriminativeParse <: JuiceTransformer end

mutable struct DiscriminativeParsePlain <: DiscriminativeParse
    nodes::Dict{String, LogisticCircuit}
    root::LogisticCircuit
    DiscriminativeParsePlain() = new(Dict{String, LogisticCircuit}(), LogisticLiteralNode(0))
end

@rule start(t::DiscriminativeParsePlain, x) = begin
    x[end]
end

@rule elem(t::DiscriminativeParsePlain, x) = begin
    [t.nodes[x[1]], t.nodes[x[2]], x[3]]
end

@rule elems(t::DiscriminativeParsePlain, x) = begin
   Array(x)
end

@rule params(t::DiscriminativeParsePlain, x) = begin
    map(p -> Base.parse(Float32, p), Array(x))
end

@rule true_node(t::DiscriminativeParsePlain, x) = begin
    lit = Base.parse(Lit, x[3])
    # In Juice, LogisticCircuit leaves are an Or node with one child
    node = Logistic⋁Node([LogisticLiteralNode(lit)], length(x[4]))
    node.thetas[1,:] .= x[4]
    t.nodes[x[1]] = node
end

@rule false_node(t::DiscriminativeParsePlain, x) = begin
    lit = -Base.parse(Lit, x[3])
    # In Juice, LogisticCircuit leaves are an Or node with one child
    node = Logistic⋁Node([LogisticLiteralNode(lit)], length(x[4]))
    node.thetas[1,:] .= x[4]
    t.nodes[x[1]] = node
    t.root = node
end

@rule or_node(t::DiscriminativeParsePlain, x) = begin
    @assert length(x[4]) == Base.parse(Int, x[3])
    children = map(x[4]) do elem
        Logistic⋀Node(elem[1:2])
    end

    node = Logistic⋁Node(children, length(x[4][1][3]))
    for i = 1:length(children)
        node.thetas[i, :] .= x[4][i][3]
    end
    t.nodes[x[1]] = node
    t.root = node
end

@rule bias_node(t::DiscriminativeParsePlain, x) = begin
    node = Logistic⋁Node([t.root], length(x[1]))
    node.thetas[1, :] .= x[1]
    t.root = node
end

function Base.parse(::Type{DiscriminativeCircuit}, str::AbstractString, ::DiscriminativeFormat)
    ast = Lerche.parse(discriminative_parser, str)
    Lerche.transform(DiscriminativeParsePlain(), ast)
end

Base.read(io::IO, ::Type{LogisticCircuit}, ::DiscriminativeFormat) =
    parse(LogisticCircuit, read(io, String), DiscriminativeFormat())

# function Base

###########################################################
#  Write DiscriminativeCircuits
###########################################################

const DISCRIM_FORMAT = """c This file was saved by DiscriminativeCircuits.jl
c variables (from inputs) start from 1
c ids of logistic circuit nodes start from 0
c nodes appear bottom-up, children before parents
c the last line of the file records the bias parameter
c three types of nodes:
c	T (terminal nodes that correspond to true literals)
c	F (terminal nodes that correspond to false literals)
c	D (OR gates)
c
c file syntax:
c
c Logistic Circuit | Regression Circuit
c T id-of-true-literal-node id-of-vtree variable parameters
c F id-of-false-literal-node id-of-vtree variable parameters
c D id-of-or-gate id-of-vtree number-of-elements (id-of-prime id-of-sub parameters)s
c B bias-parameters
c
Logistic Circuit"""

function Base.write(io::IO, pc::DiscriminativeCircuit, ::DiscriminativeFormat, vtree2id::Function = (x -> 0))

    labeling = label_nodes(pc)
    map!(x -> x-1, values(labeling)) # vtree nodes are 0-based indexed

    println(io, DISCRIM_FORMAT)
    foreach(pc) do node
        if isliteralgate(node)
            # do nothing
        elseif isconstantgate(node)
            @assert false "Cannot have constant gates $(node)"
        elseif is⋁gate(node) && length(node.children) == 1 && isliteralgate(node.children[1])
            cc = node.children[1]
            t = ispositive(cc) ? "T" : "F"
            plit = ispositive(cc) ? literal(node.children[1]) : -literal(node.children[1])

            print(io, "$(t) $(labeling[node]) $(vtree2id(node)) $(plit) ")
            for p in node.thetas[1, :]
                print(io, " $p")
            end
            println(io)
        elseif is⋀gate(node)
            # do nothing
        else
            @assert is⋁gate(node)
            # Bias node
            if length(node.children) == 1 && is⋁gate(node.children[1])
                print(io, "B")
                for p in node.thetas[1,:]
                    print(io, " $p")
                end
                println(io)
            else
                print(io, "D $(labeling[node]) $(vtree2id(node)) $(length(node.children))")
                for (id, child) in enumerate(node.children)
                    print(io, " ($(labeling[child.children[1]]) $(labeling[child.children[2]])")
                    for p in node.thetas[id, :]
                        print(io, " $p")
                    end
                    print(io, ")")
                end
                println(io)
            end
        end
    end

end

# function Base.write(ios::Tuple{IO,IO}, pc::DiscriminativeCircuit, format::VtreeFormat)
#     vtree2id = write(ios[2], vtree(pc), format[2])
#     write(ios[1], pc, format[1], i -> vtree2id[vtree(i)])
# endl

max_stars_repo_name: Juice-jl/DiscriminativeCircuits

max_issues_repo_name: Juice-jl/DiscriminativeCircuits

max_forks_repo_name: Juice-jl/DiscriminativeCircuits


text:
# a Row "view" type for iterating `CSV.File`
struct Row <: Tables.AbstractRow
    names::Vector{Symbol}
    columns::Vector{AbstractVector}
    lookup::Dict{Symbol, AbstractVector}
    row::Int64
end

getnames(r::Row) = getfield(r, :names)
getcolumn(r::Row, col::Int) = getfield(r, :columns)[col]
getcolumn(r::Row, col::Symbol) = getfield(r, :lookup)[col]
getrow(r::Row) = getfield(r, :row)

Tables.columnnames(r::Row) = getnames(r)

@inline function Tables.getcolumn(row::Row, ::Type{T}, i::Int, nm::Symbol) where {T}
    column = getcolumn(row, i)
    @inbounds x = column[getrow(row)]
    return x
end

@inline function Tables.getcolumn(row::Row, col::Symbol)
    column = getcolumn(row, col)
    @inbounds x = column[getrow(row)]
    return x
end

@inline function Tables.getcolumn(row::Row, col::Int)
    column = getcolumn(row, col)
    @inbounds x = column[getrow(row)]
    return x
end

# main structure when parsing an entire file and inferring column types
struct File{threaded} <: AbstractVector{Row}
    name::String
    names::Vector{Symbol}
    types::Vector{Type}
    rows::Int64
    cols::Int64
    columns::Vector{AbstractVector}
    lookup::Dict{Symbol, AbstractVector}
end

getname(f::File) = getfield(f, :name)
getnames(f::File) = getfield(f, :names)
gettypes(f::File) = getfield(f, :types)
getrows(f::File) = getfield(f, :rows)
getcols(f::File) = getfield(f, :cols)
getcolumns(f::File) = getfield(f, :columns)
getlookup(f::File) = getfield(f, :lookup)
getcolumn(f::File, col::Int) = getfield(f, :columns)[col]
getcolumn(f::File, col::Symbol) = getfield(f, :lookup)[col]

function Base.show(io::IO, f::File)
    println(io, "CSV.File(\"$(getname(f))\"):")
    println(io, "Size: $(getrows(f)) x $(getcols(f))")
    show(io, Tables.schema(f))
end

Base.IndexStyle(::Type{<:File}) = Base.IndexLinear()
Base.eltype(f::File) = Row
Base.size(f::File) = (getrows(f),)

Tables.isrowtable(::Type{<:File}) = true
Tables.columnaccess(::Type{<:File}) = true
Tables.schema(f::File)  = Tables.Schema(getnames(f), gettypes(f))
Tables.columns(f::File) = f
Tables.columnnames(f::File) = getnames(f)
Base.propertynames(f::File) = getnames(f)

function Base.getproperty(f::File, col::Symbol)
    lookup = getfield(f, :lookup)
    return get(lookup, col) do
        getfield(f, col)
    end
end

Tables.getcolumn(f::File, nm::Symbol) = getcolumn(f, nm)
Tables.getcolumn(f::File, i::Int) = getcolumn(f, i)

Base.@propagate_inbounds function Base.getindex(f::File, row::Int)
    @boundscheck checkbounds(f, row)
    return Row(getnames(f), getcolumns(f), getlookup(f), row)
end

@inline function Base.iterate(f::File, st::Int=1)
    st > length(f) && return nothing
    return Row(getnames(f), getcolumns(f), getlookup(f), st), st + 1
end

"""
    CSV.File(source; kwargs...) => CSV.File

Read a UTF-8 CSV input and return a `CSV.File` object.

The `source` argument can be one of:
  * filename given as a string or FilePaths.jl type
  * an `AbstractVector{UInt8}` like a byte buffer or `codeunits(string)`
  * an `IOBuffer`

To read a csv file from a url, use the HTTP.jl package, where the `HTTP.Response` body can be passed like:
```julia
f = CSV.File(HTTP.get(url).body)
```

For other `IO` or `Cmd` inputs, you can pass them like: `f = CSV.File(read(obj))`.

Opens the file and uses passed arguments to detect the number of columns and column types, unless column types are provided
manually via the `types` keyword argument. Note that passing column types manually can slightly increase performance
for each column type provided (column types can be given as a `Vector` for all columns, or specified per column via
name or index in a `Dict`).

For text encodings other than UTF-8, load the [StringEncodings.jl](https://github.com/JuliaStrings/StringEncodings.jl)
package and call e.g. `CSV.File(open(read, source, enc"ISO-8859-1"))`.

The returned `CSV.File` object supports the [Tables.jl](https://github.com/JuliaData/Tables.jl) interface
and can iterate `CSV.Row`s. `CSV.Row` supports `propertynames` and `getproperty` to access individual row values. `CSV.File`
also supports entire column access like a `DataFrame` via direct property access on the file object, like `f = CSV.File(file); f.col1`.
Note that duplicate column names will be detected and adjusted to ensure uniqueness (duplicate column name `a` will become `a_1`).
For example, one could iterate over a csv file with column names `a`, `b`, and `c` by doing:

```julia
for row in CSV.File(file)
    println("a=\$(row.a), b=\$(row.b), c=\$(row.c)")
end
```

By supporting the Tables.jl interface, a `CSV.File` can also be a table input to any other table sink function. Like:

```julia
# materialize a csv file as a DataFrame, without copying columns from CSV.File
df = CSV.File(file) |> DataFrame

# load a csv file directly into an sqlite database table
db = SQLite.DB()
tbl = CSV.File(file) |> SQLite.load!(db, "sqlite_table")
```

Supported keyword arguments include:
* File layout options:
  * `header=1`: the `header` argument can be an `Int`, indicating the row to parse for column names; or a `Range`, indicating a span of rows to be concatenated together as column names; or an entire `Vector{Symbol}` or `Vector{String}` to use as column names; if a file doesn't have column names, either provide them as a `Vector`, or set `header=0` or `header=false` and column names will be auto-generated (`Column1`, `Column2`, etc.)
  * `normalizenames=false`: whether column names should be "normalized" into valid Julia identifier symbols; useful when iterating rows and accessing column values of a row via `getproperty` (e.g. `row.col1`)
  * `datarow`: an `Int` argument to specify the row where the data starts in the csv file; by default, the next row after the `header` row is used. If `header=0`, then the 1st row is assumed to be the start of data; providing a `datarow` or `skipto` argument does _not_ affect the `header` argument
  * `skipto::Int`: identical to `datarow`, specifies the number of rows to skip before starting to read data
  * `footerskip::Int`: number of rows at the end of a file to skip parsing
  * `limit`: an `Int` to indicate a limited number of rows to parse in a csv file; use in combination with `skipto` to read a specific, contiguous chunk within a file; note for large files when multiple threads are used for parsing, the `limit` argument may not result in exact an exact # of rows parsed; use `threaded=false` to ensure an exact limit if necessary
  * `transpose::Bool`: read a csv file "transposed", i.e. each column is parsed as a row
  * `comment`: rows that begin with this `String` will be skipped while parsing
  * `ignoreemptylines::Bool=true`: whether empty rows/lines in a file should be ignored (if `false`, each column will be assigned `missing` for that empty row)
  * `threaded::Bool`: whether parsing should utilize multiple threads; by default threads are used on large enough files, but isn't allowed when `transpose=true`; only available in Julia 1.3+
  * `tasks::Integer=Threads.nthreads()`: for multithreaded parsing, this controls the number of tasks spawned to read a file in chunks concurrently; defaults to the # of threads Julia was started with (i.e. `JULIA_NUM_THREADS` environment variable)
  * `select`: an `AbstractVector` of `Int`, `Symbol`, `String`, or `Bool`, or a "selector" function of the form `(i, name) -> keep::Bool`; only columns in the collection or for which the selector function returns `true` will be parsed and accessible in the resulting `CSV.File`. Invalid values in `select` are ignored.
  * `drop`: inverse of `select`; an `AbstractVector` of `Int`, `Symbol`, `String`, or `Bool`, or a "drop" function of the form `(i, name) -> drop::Bool`; columns in the collection or for which the drop function returns `true` will ignored in the resulting `CSV.File`. Invalid values in `drop` are ignored.
* Parsing options:
  * `missingstrings`, `missingstring`: either a `String`, or `Vector{String}` to use as sentinel values that will be parsed as `missing`; by default, only an empty field (two consecutive delimiters) is considered `missing`
  * `delim=','`: a `Char` or `String` that indicates how columns are delimited in a file; if no argument is provided, parsing will try to detect the most consistent delimiter on the first 10 rows of the file
  * `ignorerepeated::Bool=false`: whether repeated (consecutive) delimiters should be ignored while parsing; useful for fixed-width files with delimiter padding between cells
  * `quotechar='"'`, `openquotechar`, `closequotechar`: a `Char` (or different start and end characters) that indicate a quoted field which may contain textual delimiters or newline characters
  * `escapechar='"'`: the `Char` used to escape quote characters in a quoted field
  * `dateformat::Union{String, Dates.DateFormat, Nothing}`: a date format string to indicate how Date/DateTime columns are formatted for the entire file
  * `dateformats::Union{AbstractDict, Nothing}`: a Dict of date format strings to indicate how the Date/DateTime columns corresponding to the keys are formatted. The Dict can map column index `Int`, or name `Symbol` or `String` to the format string for that column.
  * `decimal='.'`: a `Char` indicating how decimals are separated in floats, i.e. `3.14` used '.', or `3,14` uses a comma ','
  * `truestrings`, `falsestrings`: `Vectors of Strings` that indicate how `true` or `false` values are represented; by default only `true` and `false` are treated as `Bool`
* Column Type Options:
  * `type`: a single type to use for parsing an entire file; i.e. all columns will be treated as the same type; useful for matrix-like data files
  * `types`: a Vector or Dict of types to be used for column types; a Dict can map column index `Int`, or name `Symbol` or `String` to type for a column, i.e. Dict(1=>Float64) will set the first column as a Float64, Dict(:column1=>Float64) will set the column named column1 to Float64 and, Dict("column1"=>Float64) will set the column1 to Float64; if a `Vector` if provided, it must match the # of columns provided or detected in `header`
  * `typemap::Dict{Type, Type}`: a mapping of a type that should be replaced in every instance with another type, i.e. `Dict(Float64=>String)` would change every detected `Float64` column to be parsed as `String`; only "standard" types are allowed to be mapped to another type, i.e. `Int64`, `Float64`, `Date`, `DateTime`, `Time`, and `Bool`. If a column of one of those types is "detected", it will be mapped to the specified type.
  * `pool::Union{Bool, Float64}=0.1`: if `true`, *all* columns detected as `String` will be internally pooled; alternatively, the proportion of unique values below which `String` columns should be pooled (by default 0.1, meaning that if the # of unique strings in a column is under 10%, it will be pooled)
  * `lazystrings::Bool=false`: avoid allocating full strings in string columns; returns a custom `LazyStringVector` array type that *does not* support mutable operations (e.g. `push!`, `append!`, or even `setindex!`). Calling `copy(x)` will materialize a full `Vector{String}`. Also note that each `LazyStringVector` holds a reference to the full input file buffer, so it won't be closed after parsing and trying to delete or modify the file may result in errors (particularly on windows) and generally has undefined behavior. Given these caveats, this setting can help avoid lots of string allocations in large files and lead to faster parsing times.
  * `strict::Bool=false`: whether invalid values should throw a parsing error or be replaced with `missing`
  * `silencewarnings::Bool=false`: if `strict=false`, whether invalid value warnings should be silenced
"""
function File(source;
    # file options
    # header can be a row number, range of rows, or actual string vector
    header::Union{Integer, Vector{Symbol}, Vector{String}, AbstractVector{<:Integer}}=1,
    normalizenames::Bool=false,
    # by default, data starts immediately after header or start of file
    datarow::Integer=-1,
    skipto::Union{Nothing, Integer}=nothing,
    footerskip::Integer=0,
    transpose::Bool=false,
    comment::Union{String, Nothing}=nothing,
    use_mmap=nothing,
    ignoreemptylines::Bool=true,
    select=nothing,
    drop=nothing,
    # parsing options
    missingstrings=String[],
    missingstring="",
    delim::Union{Nothing, Char, String}=nothing,
    ignorerepeated::Bool=false,
    quotechar::Union{UInt8, Char}='"',
    openquotechar::Union{UInt8, Char, Nothing}=nothing,
    closequotechar::Union{UInt8, Char, Nothing}=nothing,
    escapechar::Union{UInt8, Char}='"',
    dateformat::Union{String, Dates.DateFormat, Nothing}=nothing,
    dateformats::Union{AbstractDict, Nothing}=nothing,
    decimal::Union{UInt8, Char}=UInt8('.'),
    truestrings::Union{Vector{String}, Nothing}=["true", "True", "TRUE"],
    falsestrings::Union{Vector{String}, Nothing}=["false", "False", "FALSE"],
    # type options
    type=nothing,
    types=nothing,
    typemap::Dict=Dict{Type, Type}(),
    categorical=nothing,
    pool::Union{Bool, Real}=0.1,
    lazystrings::Bool=false,
    strict::Bool=false,
    silencewarnings::Bool=false,
    debug::Bool=false,
    parsingdebug::Bool=false,
    kw...)

    h = Header(source, header, normalizenames, datarow, skipto, footerskip, transpose, comment, use_mmap, ignoreemptylines, select, drop, missingstrings, missingstring, delim, ignorerepeated, quotechar, openquotechar, closequotechar, escapechar, dateformat, dateformats, decimal, truestrings, falsestrings, type, types, typemap, categorical, pool, lazystrings, strict, silencewarnings, debug, parsingdebug, false)
    return File(h; debug=debug, typemap=typemap, kw...)
end

function File(h::Header;
    startingbyteposition=nothing,
    endingbyteposition=nothing,
    limit::Union{Integer, Nothing}=nothing,
    threaded::Union{Bool, Nothing}=nothing,
    typemap::Dict=Dict{Type, Type}(),
    tasks::Integer=Threads.nthreads(),
    debug::Bool=false,
    )
    rowsguess, ncols, buf, len, datapos, datarow, options, coloptions, positions, types, flags, pool, categorical, customtypes = h.rowsguess, h.cols, h.buf, h.len, h.datapos, h.datarow, h.options, h.coloptions, h.positions, h.types, h.flags, h.pool, h.categorical, h.customtypes
    if startingbyteposition !== nothing
        datapos = startingbyteposition
    end
    if endingbyteposition !== nothing
        len = endingbyteposition
    end
    transpose = gettranspose(h)
    # determine if we can use threads while parsing
    minrows = min(something(limit, typemax(Int64)), rowsguess)
    if threaded === nothing && VERSION >= v"1.3-DEV" && tasks > 1 && !transpose && minrows > (tasks * 5) && (minrows * ncols) >= 5_000
        threaded = true
    elseif threaded === true
        if VERSION < v"1.3-DEV"
            @warn "incompatible julia version for `threaded=true`: $VERSION, requires >= v\"1.3\", setting `threaded=false`"
            threaded = false
        elseif transpose
            @warn "`threaded=true` not supported on transposed files"
            threaded = false
        elseif tasks == 1
            @warn "`threaded=true` but `tasks=1`; to support threaded parsing, pass `tasks=N` where `N > 1`; `tasks` defaults to `Threads.nthreads()`, so you may consider starting Julia with multiple threads"
            threaded = false
        end
    end

    # we now do our parsing pass over the file, starting at datapos
    refs = Vector{RefPool}(undef, ncols)
    if threaded === true
        # multithreaded
        finalrows, columns = multithreadparse(types, flags, buf, datapos, len, options, coloptions, rowsguess, datarow - 1, pool, refs, ncols, typemap, h.categorical, customtypes, limit, tasks, debug)
    else
        if limit !== nothing
            rowsguess = limit
        end
        columns = allocate(rowsguess, ncols, types, flags, refs)
        t = Base.time()
        finalrows, pos = parsefilechunk!(Val(transpose), ncols, typemap, columns, buf, datapos, len, something(limit, typemax(Int64)), positions, pool, refs, rowsguess, datarow - 1, types, flags, debug, options, coloptions, customtypes)
        debug && println("time for initial parsing: $(Base.time() - t)")
        # cleanup our columns if needed
        for i = 1:ncols
            column = columns[i]
            if column isa Vector{UInt32}
                # make a PooledArray given a columns' values (Vector{UInt32}) and refs (RefPool)
                makeandsetpooled!(columns, i, column, refs, flags, h.categorical)
            elseif column isa Vector{PosLen} || column isa ChainedVector{PosLen, Vector{PosLen}}
                # string column parsed lazily; return a LazyStringVector
                if anymissing(flags[i])
                    columns[i] = LazyStringVector{Union{String, Missing}}(buf, options.e, column)
                else
                    columns[i] = LazyStringVector{String}(buf, options.e, column)
                end
            else
                # if no missing values were parsed for a column, we want to "unwrap" it to a plain Vector{T}
                if !anymissing(flags[i])
                    if column isa Vector{Union{Missing, Bool}}
                        columns[i] = convert(Vector{Bool}, column)
                    elseif types[i] !== Union{} && types[i] <: SmallIntegers
                        columns[i] = convert(Vector{types[i]}, column)
                    else
                        columns[i] = parent(column)
                    end
                end
            end
            # if zero rows
            if types[i] === Union{}
                types[i] = Missing
                columns[i] = MissingVector(finalrows)
            end
        end
    end
    debug && println("types after parsing: $types, pool = $pool")
    columns = columns
    # delete any dropped columns from names, types, columns
    deleteat!(h.names, h.todrop)
    deleteat!(types, h.todrop)
    ncols -= length(h.todrop)
    deleteat!(columns, h.todrop)
    lookup = Dict(k => v for (k, v) in zip(h.names, columns))
    # for windows, it's particularly finicky about throwing errors when you try to modify an mmapped file
    # so we just want to make sure we finalize the input buffer so users don't run into surprises
    if Sys.iswindows() && ncols > 0 && !lazystrings(flags[1])
        finalize(buf)
    end
    return File{something(threaded, false)}(h.name, h.names, types, finalrows, ncols, columns, lookup)
end

const EMPTY_INT_ARRAY = Int64[]
const EMPTY_REFRECODE = UInt32[]

# after multithreaded parsing, we need to synchronize pooled refs from different chunks of the file
# we pick one chunk as "source of truth", then adjust other chunks as needed
function syncrefs!(refs, task_refs, col, task_rows, column)
    if !isassigned(refs, col)
        @inbounds refs[col] = task_refs[col]
    else
        @inbounds colrefs = refs[col]
        @inbounds tlcolrefs = task_refs[col]
        refrecodes = EMPTY_REFRECODE
        recode = false
        for (k, v) in tlcolrefs.refs
            refvalue = get(colrefs.refs, k, UInt32(0))
            if refvalue != v
                recode = true
                if isempty(refrecodes)
                    refrecodes = [UInt32(i) for i = 1:tlcolrefs.lastref]
                end
                if refvalue == 0
                    # "source of truth" chunk didn't know about this ref, so we add it to source of truth
                    refvalue = (colrefs.lastref += UInt32(1))
                end
                @inbounds colrefs.refs[k] = refvalue
                @inbounds refrecodes[v] = refvalue
            end
        end
        if recode
            for j = 1:task_rows
                @inbounds column[j] = refrecodes[column[j]]
            end
        end
    end
    return
end

vectype(::Type{A}) where {A <: SentinelVector{T}} where {T} = Vector{T}
vectype(T) = T

# while each chunk parses values into a SentinelVector, we invert after parsing is done
# so the final array type is SentinelVector wrapping a ChainedVector of the raw chunk vectors
function makechain(::Type{T}, column, N, col, pertaskcolumns, limit, sentref=nothing) where {T}
    if T === SVec{Int64}
        # synchronize int sentinels if necessary
        SentinelArrays.newsentinel!((pertaskcolumns[i][col]::SVec{Int64} for i = 1:N)...; force=false)
        sentref[] = column.sentinel
    end
    chain = Vector{vectype(T)}(undef, N)
    @inbounds chain[1] = parent(column)
    for i = 2:N
        @inbounds tp = pertaskcolumns[i][col]
        if tp isa T
            @inbounds chain[i] = parent(tp)
        else
            @inbounds chain[i] = convert(T, tp)
        end
    end
    x = ChainedVector(chain)
    if limit !== nothing && limit < length(x)
        resize!(x, limit)
    end
    return x
end

function makeandsetpooled!(columns, i, column, refs, flags, categorical)
    if categorical
        colrefs = isassigned(refs, i) ? refs[i].refs : Dict{String, UInt32}()
        if anymissing(flags[i]) && haskey(colrefs, missing)
            missingref = colrefs[missing]
            delete!(colrefs, missing)
            colrefs = convert(Dict{String, UInt32}, colrefs)
            for j = 1:length(column)
                @inbounds column[j] = ifelse(column[j] === missingref, UInt32(0), column[j])
            end
        else
            colrefs = convert(Dict{String, UInt32}, colrefs)
        end
        pool = CategoricalPool(colrefs)
        A = CategoricalArray{anymissing(flags[i]) ? Union{String, Missing} : String, 1}(column, pool)
        levels!(A, sort(levels(A)))
        columns[i] = A
    else
        r = isassigned(refs, i) ? (anymissing(flags[i]) ? refs[i].refs : convert(Dict{String, UInt32}, refs[i].refs)) : Dict{String, UInt32}()
        columns[i] = PooledArray(PooledArrays.RefArray(column), r)
    end
    return
end

function multithreadparse(types, flags, buf, datapos, len, options, coloptions, origrowsguess, datarow, pool, refs, ncols, typemap, categorical, customtypes, limit, N, debug)
    # when limiting w/ multithreaded parsing, we try to guess about where in the file the limit row # will be
    # then adjust our final file len to the end of that row
    # we add some cushion so we hopefully get the limit row correctly w/o shooting past too far and needing to resize! down
    # but we also don't guarantee limit will be exact w/ multithreaded parsing
    if limit !== nothing
        newlen = [0, ceil(Int64, (limit / (origrowsguess * 0.8)) * len), 0]
        findrowstarts!(buf, len, options, newlen, ncols, types, flags)
        len = newlen[2] - 1
        origrowsguess = limit
        debug && println("limiting, adjusting len to $len")
    end
    chunksize = div(len - datapos, N)
    ranges = [i == 0 ? datapos : (datapos + chunksize * i) for i = 0:N]
    ranges[end] = len
    debug && println("initial byte positions before adjusting for start of rows: $ranges")
    avgbytesperrow = findrowstarts!(buf, len, options, ranges, ncols, types, flags)
    origbytesperrow = ((len - datapos) / origrowsguess)
    weightedavgbytesperrow = ceil(Int64, avgbytesperrow * ((N - 2) / N) + origbytesperrow * (2 / N))
    rowsguess = ceil(Int64, (len - datapos) / weightedavgbytesperrow)
    debug && println("single-threaded estimated rows = $origrowsguess, multi-threaded estimated rows = $rowsguess")
    debug && println("multi-threaded column types sampled as: $types")
    rowchunkguess = cld(rowsguess, N)
    debug && println("parsing using $N tasks: $rowchunkguess rows chunked at positions: $ranges")
    pertaskcolumns = Vector{Vector{AbstractVector}}(undef, N)
    rows = zeros(Int64, N)
    locks = [ReentrantLock() for i = 1:ncols]
    @sync for i = 1:N
@static if VERSION >= v"1.3-DEV"
        Threads.@spawn begin
            tt = Base.time()
            task_flags = copy(flags)
            task_refs = Vector{RefPool}(undef, ncols)
            task_len = ranges[i + 1] - (i != N)
            task_pos = ranges[i]
            task_types = copy(types)
            task_columns = allocate(rowchunkguess, ncols, task_types, task_flags, task_refs)
            pertaskcolumns[i] = task_columns
            task_rows, task_pos = parsefilechunk!(Val(false), ncols, typemap, task_columns, buf, task_pos, task_len, typemax(Int64), EMPTY_INT_ARRAY, pool, task_refs, rowchunkguess, datarow + (rowchunkguess * (i - 1)), task_types, task_flags, debug, options, coloptions, customtypes)
            rows[i] = task_rows
            # promote column types/flags across task chunks
            for col = 1:ncols
                lock(locks[col])
                try
                    @inbounds T = types[col]
                    @inbounds types[col] = promote_types(T, task_types[col])
                    if T !== types[col]
                        debug && println("promoting col = $col from $T to $(types[col]), task chunk ($i) was type = $(task_types[col])")
                    end
                    @inbounds flags[col] |= task_flags[col]
                    # synchronize refs if needed
                    @inbounds column = task_columns[col]
                    if column isa Vector{UInt32}
                        syncrefs!(refs, task_refs, col, task_rows, column)
                    end
                finally
                    unlock(locks[col])
                end
            end
            debug && println("finished parsing $task_rows rows on task = $i: time for parsing: $(Base.time() - tt)")
        end
end # @static if VERSION >= v"1.3-DEV"
    end
    finalcolumns = Vector{AbstractVector}(undef, ncols)
    @sync for col = 1:ncols
@static if VERSION >= v"1.3-DEV"
        Threads.@spawn begin
            for i = 1:N
                task_columns = pertaskcolumns[i]
                task_rows = rows[i]
                # check if we need to promote a task-local column based on what other threads parsed
                @inbounds T = types[col]
                if (T === String || T === Union{String, Missing}) && !(task_columns[col] isa Vector{PosLen}) && !(task_columns[col] isa SVec2{String})
                    # promoting non-string to string column
                    debug && println("multithreaded promoting column $col to string from $(typeof(task_columns[col]))")
                    promotetostring!(col, Val(false), ncols, typemap, task_columns, buf, ranges[i], ranges[i + 1] - (i != N), task_rows, EMPTY_INT_ARRAY, pool, refs, task_rows, 0, types, flags, debug, options, coloptions, customtypes)
                elseif (T === Float64 || T === Union{Float64, Missing}) && task_columns[col] isa SVec{Int64}
                    # one chunk parsed as Int, another as Float64, promote to Float64
                    debug && println("multithreaded promoting column $col to float")
                    task_columns[col] = convert(SentinelVector{Float64}, task_columns[col])
                elseif T !== Union{} && T !== Missing && task_columns[col] isa MissingVector
                    # one chunk parsed all missing values, but another chunk had a typed value, promote to that
                    debug && println("multithreaded promoting column $col from missing on task $i")
                    task_columns[col] = allocate(T, task_rows)
                    if T == Union{PooledString, Missing}
                        colrefs = refs[col]
                        ref = getref!(colrefs, missing, Int16(0), options)
                        column = task_columns[col]
                        for j = 1:task_rows
                            @inbounds column[j] = ref
                        end
                    end
                end
            end
            @inbounds column = pertaskcolumns[1][col]
            if column isa SVec{Int64}
                sentref = Ref(column.sentinel)
                chain = makechain(SVec{Int64}, column, N, col, pertaskcolumns, limit, sentref)
                @inbounds finalcolumns[col] = anymissing(flags[col]) ? SentinelArray(chain, sentref[]) : chain
            elseif column isa SVec{Float64}
                chain = makechain(SVec{Float64}, column, N, col, pertaskcolumns, limit)
                @inbounds finalcolumns[col] = anymissing(flags[col]) ? SentinelArray(chain) : chain
            elseif column isa SVec2{String}
                chain = makechain(SVec2{String}, column, N, col, pertaskcolumns, limit)
                @inbounds finalcolumns[col] = anymissing(flags[col]) ? SentinelArray(chain) : chain
            elseif column isa SVec{Date}
                chain = makechain(SVec{Date}, column, N, col, pertaskcolumns, limit)
                @inbounds finalcolumns[col] = anymissing(flags[col]) ? SentinelArray(chain) : chain
            elseif column isa SVec{DateTime}
                chain = makechain(SVec{DateTime}, column, N, col, pertaskcolumns, limit)
                @inbounds finalcolumns[col] = anymissing(flags[col]) ? SentinelArray(chain) : chain
            elseif column isa SVec{Time}
                chain = makechain(SVec{Time}, column, N, col, pertaskcolumns, limit)
                @inbounds finalcolumns[col] = anymissing(flags[col]) ? SentinelArray(chain) : chain
            elseif column isa Vector{Union{Missing, Bool}}
                chain = makechain(Vector{anymissing(flags[col]) ? Union{Missing, Bool} : Bool}, column, N, col, pertaskcolumns, limit)
                @inbounds finalcolumns[col] = chain
            elseif column isa Vector{UInt32}
                chain = makechain(Vector{UInt32}, column, N, col, pertaskcolumns, limit)
                makeandsetpooled!(finalcolumns, col, chain, refs, flags, categorical)
            elseif column isa Vector{PosLen}
                chain = makechain(Vector{PosLen}, column, N, col, pertaskcolumns, limit)
                if anymissing(flags[col])
                    @inbounds finalcolumns[col] = LazyStringVector{Union{String, Missing}}(buf, options.e, chain)
                else
                    @inbounds finalcolumns[col] = LazyStringVector{String}(buf, options.e, chain)
                end
            elseif column isa MissingVector
                chain = makechain(MissingVector, column, N, col, pertaskcolumns, limit)
                @inbounds finalcolumns[col] = chain
                types[col] = Missing
            elseif nonmissingtype(types[col]) <: SmallIntegers
                T = types[col]
                chain = makechain(Vector{T}, column, N, col, pertaskcolumns, limit)
                @inbounds finalcolumns[col] = anymissing(flags[col]) ? chain : chain
            else
                chain = makechain(typeof(column), column, N, col, pertaskcolumns, limit)
                @inbounds finalcolumns[col] = anymissing(flags[col]) ? SentinelArray(chain) : chain
                # error("unhandled column type: $(typeof(column))")
            end
        end
end # @static if VERSION >= v"1.3-DEV"
    end
    finalrows = sum(rows)
    return limit !== nothing && limit < finalrows ? limit : finalrows, finalcolumns
end

function parsefilechunk!(TR::Val{transpose}, ncols, typemap, columns, buf, pos, len, limit, positions, pool, refs, rowsguess, rowoffset, types, flags, debug, options::Parsers.Options{ignorerepeated}, coloptions, ::Type{customtypes}) where {transpose, ignorerepeated, customtypes}
    row = 0
    startpos = pos
    if pos <= len && len > 0
        while row < limit
            row += 1
            pos = parserow(row, TR, ncols, typemap, columns, startpos, buf, pos, len, positions, pool, refs, rowsguess, rowoffset, types, flags, debug, options, coloptions, customtypes)
            (pos > len || row == limit) && break
            # if our initial row estimate was too few, we need to reallocate our columsn to read the rest of the file
            if row + 1 > rowsguess
                # (bytes left in file) / (avg bytes per row) == estimated rows left in file (+ 10 for kicks)
                estimated_rows_left = ceil(Int64, (len - pos) / ((pos - startpos) / row) + 10.0)
                newrowsguess = rowsguess + estimated_rows_left
                debug && reallocatetape(rowoffset + row, rowsguess, newrowsguess)
                for i = 1:ncols
                    reallocate!(columns[i], newrowsguess)
                end
                rowsguess = newrowsguess
            end
        end
    end
    # done parsing (at least this chunk), so resize columns to final row count
    for i = 1:ncols
        resize!(columns[i], row)
    end
    return row, pos
end

@noinline function promotetostring!(col, TR::Val{transpose}, ncols, typemap, columns, buf, pos, len, limit, positions, pool, refs, rowsguess, rowoffset, types, origflags, debug, options::Parsers.Options{ignorerepeated}, coloptions, ::Type{customtypes}) where {transpose, ignorerepeated, customtypes}
    if columns[col] isa Vector{UInt32} && !lazystrings(origflags[col])
        # optimize conversion from PooledString => String w/o reparsing
        types[col] = Union{String, anymissing(origflags[col]) ? Missing : Union{}}
        pooledvalues = columns[col]
        columns[col] = column = allocate(lazystrings(origflags[col]) ? PosLen : String, rowsguess)
        colrefs = refs[col].refs
        refvalues = anymissing(origflags[col]) ? Vector{Union{String, Missing}}(undef, length(colrefs)) : Vector{String}(undef, length(colrefs))
        for (k, v) in colrefs
            @inbounds refvalues[v] = k
        end
        for row = 1:limit
            @inbounds column[row] = refvalues[pooledvalues[row]]
        end
        return
    end
    flags = copy(origflags)
    for i = 1:ncols
        if i == col
            flags[i] |= TYPEDETECTED
            columns[i] = allocate(lazystrings(flags[i]) ? PosLen : String, rowsguess)
            types[i] = Union{String, anymissing(flags[i]) ? Missing : Union{}}
        else
            flags[i] |= WILLDROP
        end
    end
    row = 0
    startpos = pos
    if pos <= len && len > 0
        while row < limit
            row += 1
            pos = parserow(row, TR, ncols, typemap, columns, startpos, buf, pos, len, positions, pool, refs, rowsguess, rowoffset, types, flags, debug, options, coloptions, customtypes)
            pos > len && break
        end
    end
    return
end

@noinline reallocatetape(row, old, new) = println("thread = $(Threads.threadid()) warning: didn't pre-allocate enough column while parsing around row $row, re-allocating from $old to $new...")
@noinline notenoughcolumns(cols, ncols, row) = println("thread = $(Threads.threadid()) warning: only found $cols / $ncols columns around data row: $row. Filling remaining columns with `missing`")
@noinline toomanycolumns(cols, row) = println("thread = $(Threads.threadid()) warning: parsed expected $cols columns, but didn't reach end of line around data row: $row. Ignoring any extra columns on this row")
@noinline stricterror(T, buf, pos, len, code, row, col) = throw(Error("thread = $(Threads.threadid()) error parsing $T around row = $row, col = $col: \"$(String(buf[pos:pos+len-1]))\", error=$(Parsers.codes(code))"))
@noinline warning(T, buf, pos, len, code, row, col) = println("thread = $(Threads.threadid()) warning: error parsing $T around row = $row, col = $col: \"$(String(buf[pos:pos+len-1]))\", error=$(Parsers.codes(code))")
@noinline fatalerror(buf, pos, len, code, row, col) = throw(Error("thread = $(Threads.threadid()) fatal error, encountered an invalidly quoted field while parsing around row = $row, col = $col: \"$(String(buf[pos:pos+len-1]))\", error=$(Parsers.codes(code)), check your `quotechar` arguments or manually fix the field in the file itself"))

@inline function parsecustom!(::Type{T}, flag, columns, buf, pos, len, opts, row, rowoffset, col, types, flags) where {T}
    if @generated
        block = Expr(:block)
        push!(block.args, quote
            error("CSV.jl code-generation error, unexpected column type: $(typeof(column))")
        end)
        for i = 1:fieldcount(T)
            vec = fieldtype(T, i)
            pushfirst!(block.args, quote
                if column isa $(fieldtype(vec, 1))
                    return parsevalue!($(fieldtype(vec, 2)), flag, column, columns, buf, pos, len, opts, row, rowoffset, col, types, flags)
                end
            end)
        end
        pushfirst!(block.args, quote
            @inbounds column = columns[col]
        end)
        pushfirst!(block.args, Expr(:meta, :inline))
        # @show block
        return block
    else
        # println("generated function failed")
        @inbounds column = columns[col]
        return parsevalue!(eltype(parent(column)), flag, column, columns, buf, pos, len, opts, row, rowoffset, col, types, flags)
    end
end

@inline function parserow(row, TR::Val{transpose}, ncols, typemap, columns, startpos, buf, pos::A, len, positions, pool, refs, rowsguess, rowoffset, types, flags, debug, options::B, coloptions::C, ::Type{customtypes}) where {transpose, A, B, C, customtypes}
    for col = 1:ncols
        if transpose
            @inbounds pos = positions[col]
        end
        @inbounds flag = flags[col]
        @inbounds column = columns[col]
        @inbounds opts = coloptions === nothing ? options : coloptions[col]
        # @show typeof(column)
        if willdrop(flag) || (user(flag) && column isa MissingVector)
            pos, code = parsemissing!(buf, pos, len, opts, row, rowoffset, col)
        elseif !typedetected(flag)
            pos, code = detect(columns, buf, pos, len, opts, row, rowoffset, col, typemap, pool, refs, debug, types, flags, rowsguess)
        elseif column isa SVec{Int64}
            pos, code = parseint!(flag, column, columns, buf, pos, len, opts, row, rowoffset, col, types, flags)
        elseif column isa SVec{Float64}
            pos, code = parsevalue!(Float64, flag, column, columns, buf, pos, len, opts, row, rowoffset, col, types, flags)
        elseif column isa SVec2{String}
            pos, code = parsestring2!(flag, column, buf, pos, len, opts, row, rowoffset, col, types, flags)
        elseif column isa SVec{Date}
            pos, code = parsevalue!(Date, flag, column, columns, buf, pos, len, opts, row, rowoffset, col, types, flags)
        elseif column isa SVec{DateTime}
            pos, code = parsevalue!(DateTime, flag, column, columns, buf, pos, len, opts, row, rowoffset, col, types, flags)
        elseif column isa SVec{Time}
            pos, code = parsevalue!(Time, flag, column, columns, buf, pos, len, opts, row, rowoffset, col, types, flags)
        elseif column isa Vector{Union{Missing, Bool}}
            pos, code = parsevalue!(Bool, flag, column, columns, buf, pos, len, opts, row, rowoffset, col, types, flags)
        elseif column isa Vector{UInt32}
            pos, code = parsepooled!(flag, column, columns, buf, pos, len, opts, row, rowoffset, col, rowsguess, pool, refs, types, flags)
        elseif column isa Vector{PosLen}
            pos, code = parsestring!(flag, column, buf, pos, len, opts, row, rowoffset, col, types, flags)
        elseif customtypes !== Tuple{}
            pos, code = parsecustom!(customtypes, flag, columns, buf, pos, len, opts, row, rowoffset, col, types, flags)
        else
            error("bad array type: $(typeof(column))")
        end
        if promote_to_string(code)
            debug && println("promoting col = $col to string from $(typeof(column)) on chunk = $(Threads.threadid())")
            promotetostring!(col, TR, ncols, typemap, columns, buf, startpos, len, row, positions, pool, refs, rowsguess, rowoffset, types, flags, debug, options, coloptions, customtypes)
        end
        if transpose
            @inbounds positions[col] = pos
        else
            if col < ncols
                if Parsers.newline(code) || pos > len
                    options.silencewarnings || notenoughcolumns(col, ncols, rowoffset + row)
                    for j = (col + 1):ncols
                        @inbounds flags[j] |= ANYMISSING
                        @inbounds types[j] = Union{Missing, types[j]}
                        if columns[j] isa Vector{UInt32}
                            ref = getref!(refs[j], missing, code, options)
                            columns[j][row] = ref
                        end
                    end
                    break # from for col = 1:ncols
                end
            else
                if pos <= len && !Parsers.newline(code)
                    options.silencewarnings || toomanycolumns(ncols, rowoffset + row)
                    # ignore the rest of the line
                    pos = skiptorow(buf, pos, len, options.oq, options.e, options.cq, 1, 2)
                end
            end
        end
    end
    return pos
end

@inline function poslen(code, pos, len)
    pos = Core.bitcast(UInt64, pos) << 20
    pos |= ifelse(Parsers.sentinel(code), MISSING_BIT, UInt64(0))
    pos |= ifelse(Parsers.escapedstring(code), ESCAPE_BIT, UInt64(0))
    return pos | Core.bitcast(UInt64, len)
end

@inline function setposlen!(column, row, code, pos, len)
    @inbounds column[row] = poslen(code, pos, len)
    return
end

function detect(columns, buf, pos, len, options, row, rowoffset, col, typemap, pool, refs, debug, types, flags, rowsguess)
    # debug && println("detecting on task $(Threads.threadid())")
    int, code, vpos, vlen, tlen = Parsers.xparse(Int64, buf, pos, len, options)
    if Parsers.invalidquotedfield(code)
        fatalerror(buf, pos, tlen, code, rowoffset + row, col)
    end
    if Parsers.sentinel(code) && code > 0
        flags[col] |= ANYMISSING
        # return; parsing will continue to detect until a non-missing value is parsed
        @goto finaldone
    end
    if Parsers.ok(code)
        intT = get(typemap, Int64, Int64)
        if intT <: Integer
            # debug && println("detecting Int64 for column = $col on task = $(Threads.threadid())")
            column = allocate(intT, rowsguess)
            column[row] = int
            newT = intT
            @goto done
        elseif intT === String
            @goto stringdetect
        end
    end
    float, code, vpos, vlen, tlen = Parsers.xparse(Float64, buf, pos, len, options)
    if Parsers.ok(code)
        floatT = get(typemap, Float64, Float64)
        if floatT <: AbstractFloat
            # debug && println("detecting Float64 for column = $col on task = $(Threads.threadid())")
            column = allocate(floatT, rowsguess)
            column[row] = float
            newT = floatT
            @goto done
        elseif floatT === String
            @goto stringdetect
        end
    end
    if options.dateformat === nothing
        try
            date, code, vpos, vlen, tlen = Parsers.xparse(Date, buf, pos, len, options)
            if Parsers.ok(code)
                dtT = get(typemap, Date, Date)
                if dtT === Date
                    column = allocate(Date, rowsguess)
                    column[row] = date
                    newT = Date
                    @goto done
                elseif dtT === String
                    @goto stringdetect
                end
            end
        catch
        end
        try
            datetime, code, vpos, vlen, tlen = Parsers.xparse(DateTime, buf, pos, len, options)
            if Parsers.ok(code)
                dtmT = get(typemap, DateTime, DateTime)
                if dtmT === DateTime
                    column = allocate(DateTime, rowsguess)
                    column[row] = datetime
                    newT = DateTime
                    @goto done
                elseif dtmT === String
                    @goto stringdetect
                end
            end
        catch
        end
        try
            time, code, vpos, vlen, tlen = Parsers.xparse(Time, buf, pos, len, options)
            if Parsers.ok(code)
                tT = get(typemap, Time, Time)
                if tT === Time
                    column = allocate(Time, rowsguess)
                    column[row] = time
                    newT = Time
                    @goto done
                elseif tT === String
                    @goto stringdetect
                end
            end
        catch
        end
    else
        try
            # use user-provided dateformat
            DT = timetype(options.dateformat)
            dt, code, vpos, vlen, tlen = Parsers.xparse(DT, buf, pos, len, options)
            if Parsers.ok(code)
                column = allocate(DT, rowsguess)
                column[row] = dt
                newT = DT
                @goto done
            end
        catch
        end
    end
    bool, code, vpos, vlen, tlen = Parsers.xparse(Bool, buf, pos, len, options)
    if Parsers.ok(code) && get(typemap, Bool, Bool) === Bool
        column = allocate(Bool, rowsguess)
        column[row] = bool
        newT = Bool
        @goto done
    end
@label stringdetect
    _, code, vpos, vlen, tlen = Parsers.xparse(String, buf, pos, len, options)
    if pool > 0.0
        # debug && println("detecting PooledString for column = $col on task = $(Threads.threadid()); pool = $pool, row = $row, POOLSAMPLESIZE = $POOLSAMPLESIZE")
        r = RefPool()
        @inbounds refs[col] = r
        column = allocate(PooledString, pool < 1.0 && row < POOLSAMPLESIZE ? min(rowsguess, POOLSAMPLESIZE) : rowsguess)
        if pool < 1.0 && row < POOLSAMPLESIZE
            flags[col] |= MAYBEPOOLED
        end
        if anymissing(flags[col])
            ref = getref!(r, missing, code, options)
            for i = 1:(row - 1)
                @inbounds column[i] = ref
            end
        end
        ref = getref!(r, PointerString(pointer(buf, vpos), vlen), code, options)
        @inbounds column[row] = ref
        newT = PooledString
    else
        # debug && println("detecting String for column = $col on task = $(Threads.threadid()); pool = $pool, row = $row, POOLSAMPLESIZE = $POOLSAMPLESIZE")
        if lazystrings(flags[col])
            column = allocate(PosLen, rowsguess)
            @inbounds setposlen!(column, row, code, vpos, vlen)
        else
            column = allocate(String, rowsguess)
            @inbounds column[row] = str(buf, options.e, poslen(code, vpos, vlen))
        end
        newT = String
    end
@label done
    # if we're here, that means we found a non-missing value, so we need to update columns
    columns[col] = column
    flags[col] |= TYPEDETECTED
    types[col] = Union{newT, anymissing(flags[col]) ? Missing : Union{}}
@label finaldone
    return pos + tlen, code
end

function parseint!(flag, column, columns, buf, pos, len, options, row, rowoffset, col, types, flags)::Tuple{Int64, Int16}
    x, code, vpos, vlen, tlen = Parsers.xparse(Int64, buf, pos, len, options)
    if code > 0
        if !Parsers.sentinel(code)
            @inbounds column[row] = x
        else
            if !anymissing(flag)
                @inbounds flags[col] = flag | ANYMISSING
                @inbounds types[col] = Union{Int64, Missing}
            end
        end
    else
        if Parsers.invalidquotedfield(code)
            # this usually means parsing is borked because of an invalidly quoted field, hard error
            fatalerror(buf, pos, tlen, code, rowoffset + row, col)
        end
        if user(flag)
            if !options.strict
                options.silencewarnings || warning(Int64, buf, pos, tlen, code, rowoffset + row, col)
                flags[col] = flag | ANYMISSING
                types[col] = Union{Int64, Missing}
            else
                stricterror(Int64, buf, pos, tlen, code, rowoffset + row, col)
            end
        else
            # println("encountered an error parsing an Intger in auto-detect mode for column = $col, so let's try Float64")
            y, code, vpos, vlen, tlen = Parsers.xparse(Float64, buf, pos, len, options)
            if code > 0
                @inbounds columns[col] = convert(SentinelVector{Float64}, column)
                @inbounds columns[col][row] = y
                @inbounds types[col] = Union{Float64, anymissing(flag) ? Missing : Union{}}
            else
                # println("promoting to String instead")
                code |= PROMOTE_TO_STRING
            end
        end
    end
    return pos + tlen, code
end

function parsevalue!(::Type{type}, flag, column, columns, buf, pos, len, options, row, rowoffset, col, types, flags)::Tuple{Int64, Int16} where {type}
    x, code, vpos, vlen, tlen = Parsers.xparse(type, buf, pos, len, options)
    if code > 0
        if !Parsers.sentinel(code)
            @inbounds column[row] = x
        else
            if !anymissing(flag)
                @inbounds flags[col] = flag | ANYMISSING
                @inbounds types[col] = Union{type, Missing}
            end
        end
    else
        if Parsers.invalidquotedfield(code)
            # this usually means parsing is borked because of an invalidly quoted field, hard error
            fatalerror(buf, pos, tlen, code, rowoffset + row, col)
        end
        if user(flag)
            if !options.strict
                options.silencewarnings || warning(type, buf, pos, tlen, code, rowoffset + row, col)
                flags[col] = flag | ANYMISSING
                types[col] = Union{type, Missing}
            else
                stricterror(type, buf, pos, tlen, code, rowoffset + row, col)
            end
        else
            code |= PROMOTE_TO_STRING
        end
    end
    return pos + tlen, code
end

function parsestring!(flag, column, buf, pos, len, options, row, rowoffset, col, types, flags)::Tuple{Int64, Int16}
    x, code, vpos, vlen, tlen = Parsers.xparse(String, buf, pos, len, options)
    setposlen!(column, row, code, vpos, vlen)
    if Parsers.invalidquotedfield(code)
        # this usually means parsing is borked because of an invalidly quoted field, hard error
        fatalerror(buf, pos, tlen, code, rowoffset + row, col)
    end
    if Parsers.sentinel(code)
        if !anymissing(flag)
            @inbounds flags[col] |= ANYMISSING
            @inbounds types[col] = Union{String, Missing}
        end
    end
    return pos + tlen, code
end

function parsestring2!(flag, column, buf, pos, len, options, row, rowoffset, col, types, flags)::Tuple{Int64, Int16}
    x, code, vpos, vlen, tlen = Parsers.xparse(String, buf, pos, len, options)
    if Parsers.invalidquotedfield(code)
        # this usually means parsing is borked because of an invalidly quoted field, hard error
        fatalerror(buf, pos, tlen, code, rowoffset + row, col)
    end
    if Parsers.sentinel(code)
        if !anymissing(flag)
            @inbounds flags[col] |= ANYMISSING
            @inbounds types[col] = Union{String, Missing}
        end
    else
        @inbounds column[row] = str(buf, options.e, poslen(code, vpos, vlen))
    end
    return pos + tlen, code
end

function parsemissing!(buf, pos, len, options, row, rowoffset, col)::Tuple{Int64, Int16}
    x, code, vpos, vlen, tlen = Parsers.xparse(String, buf, pos, len, options)
    if Parsers.invalidquotedfield(code)
        # this usually means parsing is borked because of an invalidly quoted field, hard error
        fatalerror(buf, pos, tlen, code, rowoffset + row, col)
    end
    return pos + tlen, code
end

@inline function getref!(refpool, key::Missing, code, options)
    get!(refpool.refs, key) do
        refpool.lastref += UInt32(1)
    end
end

@inline function getref!(refpool, key::PointerString, code, options)
    x = refpool.refs
    if Parsers.escapedstring(code)
        key2 = unescape(key, options.e)
        index = Base.ht_keyindex2!(x, key2)
    else
        index = Base.ht_keyindex2!(x, key)
    end
    if index > 0
        @inbounds found_key = x.vals[index]
        ret = found_key::UInt32
    else
        @inbounds new = refpool.lastref += UInt32(1)
        @inbounds Base._setindex!(x, new, Parsers.escapedstring(code) ? key2 : String(key), -index)
        ret = new
    end
    return ret
end

function parsepooled!(flag, column, columns, buf, pos, len, options, row, rowoffset, col, rowsguess, pool, refs, types, flags)::Tuple{Int64, Int16}
    x, code, vpos, vlen, tlen = Parsers.xparse(String, buf, pos, len, options)
    if Parsers.invalidquotedfield(code)
        # this usually means parsing is borked because of an invalidly quoted field, hard error
        fatalerror(buf, pos, tlen, code, rowoffset + row, col)
    end
    @inbounds colrefs = refs[col]
    if Parsers.sentinel(code)
        @inbounds flags[col] = flag | ANYMISSING
        @inbounds types[col] = Union{PooledString, Missing}
        ref = getref!(colrefs, missing, code, options)
    else
        ref = getref!(colrefs, PointerString(pointer(buf, vpos), vlen), code, options)
    end
    @inbounds column[row] = ref
    if !user(flag) && maybepooled(flag)
        if rowsguess <= POOLSAMPLESIZE && ((length(colrefs.refs) - anymissing(flags[col])) / rowsguess) > pool
            # println("promoting col = $col to string from pooled on task = $(Threads.threadid())")
            code |= PROMOTE_TO_STRING
        elseif row == POOLSAMPLESIZE
            if ((length(colrefs.refs) - anymissing(flags[col])) / min(rowsguess, POOLSAMPLESIZE)) > pool
                numrefs = length(colrefs.refs) - anymissing(flags[col])
                # println("promoting col = $col to string from pooled on task = $(Threads.threadid()); numrefs = $numrefs, POOLSAMPLESIZE = $POOLSAMPLESIZE")
                code |= PROMOTE_TO_STRING
            else
                numrefs = length(colrefs.refs) - anymissing(flags[col])
                # println("confirming col = $col is pooled on task = $(Threads.threadid()); numrefs = $numrefs")
                resize!(column, rowsguess)
            end
            flags[col] &= ~MAYBEPOOLED
        end
    end
    return pos + tlen, code
end


text:
@testset ExtendedTestSet "find_circle" begin
    @testset ExtendedTestSet "find_circle" begin
        test = find_circle([309, 309], [312, 200], [155, 155])
        answer = [212, 251]
        @test test == answer
    end

    @testset ExtendedTestSet "find_circle" begin
        test = find_circle([20, 20], [30, 10], [10, 20])
        answer = [15, 5]
        @test test == answer
    end
end

@testset ExtendedTestSet "angle_calc" begin
    @testset ExtendedTestSet "angle_calc" begin
        test = angle_calc(4, 3)
        answer = 0.9272952180016122
        @test test ≈ answer
    end

    @testset ExtendedTestSet "angle_calc" begin
        test = angle_calc(1, 2)
        answer = 0.4636476090008061
        @test test ≈ answer
    end
end

@testset ExtendedTestSet "create_circular_mask" begin
    @testset ExtendedTestSet "create_circular_mask" begin
        test = create_circular_mask(4, 4, [2, 2], 1)
        answer = Bool.([
            0 1 0 0
            1 1 1 0
            0 1 0 0
            0 0 0 0
        ])
        @test test == answer
    end
end


text:
# using PyCall
# @pyimport ase.io as ase_io
# @pyimport ase.atoms as ase_atoms
# @pyimport ase.build as ase_build

module ASE

using Reexport

@reexport using JuLIP
@reexport using PyCall

# the functions to be implemented
import JuLIP:
      positions, set_positions!,
      cell, set_cell!,             # ✓
      pbc, set_pbc!,               # ✓
      calculator, set_calculator!, # ✓
      constraint, set_constraint!, # ✓
      neighbourlist,                # ✓
      energy, forces, virial, stress,
      momenta, set_momenta!,
      masses, set_masses!,
      set_transient!,
      atomic_numbers,
      Atoms, chemical_symbols,
      get_data, has_data, set_data!,
      bulk

import Base.length, Base.deleteat!, Base.deepcopy

# from arrayconversions:
using JuLIP: mat, vecs, JVecF, JVecs, JVecsF, JMatF,
      AbstractAtoms, AbstractConstraint, NullConstraint,
      AbstractCalculator, NullCalculator, maxdist, SVec,
      Dofs, set_dofs!

# extra ASE functionality:
import Base.repeat         # ✓

export ASEAtoms,      # ✓
      AbstractASECalculator, ASECalculator,
      extend!, get_info, set_info!, get_array, set_array!, has_array, has_info,
      get_transient, has_transient,
      velocities, set_velocities!,
      static_neighbourlist,
      read_xyz, write_xyz

@pyimport ase.build as ase_build
@pyimport ase.atoms as ase_atoms
@pyimport ase.io as ase_io

#################################################################
###  Wrapper for ASE Atoms object and its basic functionality
################################################################

"""
`TransientData`

some data which needs to be updated if the configuration (positions only!) has
changed too much.
"""
mutable struct TransientData
   max_change::Float64    # how much X may change before recomputing
   accum_change::Float64   # how much has it changed already
   data::Any
end

"""
`type ASEAtoms <: AbstractAtoms`

Julia wrapper for the ASE `Atoms` class.

## Constructors

* `ASEAtoms(s::AbstractString, X::JVecsF) -> at::ASEAtoms`

For internal usage there is also a constructor `ASEAtoms(po::PyObject)`
"""
type ASEAtoms <: AbstractAtoms
   po::PyObject       # ase.Atoms instance
   calc::AbstractCalculator
   cons::AbstractConstraint
   transient::Dict{Any, TransientData}
end

ASEAtoms(po::PyObject) = ASEAtoms(po, NullCalculator(), NullConstraint(),
                                 Dict{Any, TransientData}())

function set_calculator!(at::ASEAtoms, calc::AbstractCalculator)
   at.calc = calc
   return at
end
calculator(at::ASEAtoms) = at.calc
function set_constraint!(at::ASEAtoms, cons::AbstractConstraint)
   at.cons = cons
   return at
end
constraint(at::ASEAtoms) = at.cons

ASEAtoms(s::AbstractString, X::JVecsF) = ASEAtoms(s, mat(X))
ASEAtoms(s::AbstractString, X::Matrix) = ASEAtoms(ase_atoms.Atoms(s, X'))
ASEAtoms(s::AbstractString) = ASEAtoms(ase_atoms.Atoms(s))

"Return the PyObject associated with `a`"
pyobject(a::ASEAtoms) = a.po

length(at::ASEAtoms) = size(at.po[:get_positions](), 1)

"""
return a deep copy of this ASEAtoms object. Transient data is not copied.
TODO: deepcopy() should copy the constraints and calculator too.
"""
function deepcopy(at::ASEAtoms)
    new_at = ASEAtoms(at.po[:copy]())
    set_constraint!(new_at, constraint(at))
    set_calculator!(new_at, calculator(at))
    new_at.transient = deepcopy(at.transient)
    return new_at
end

JuLIP.Chemistry.rnn(s::AbstractString) = rnn(Symbol(s))

# ==========================================
#    some logic for storing permanent data

# temporarily reverted to the old `positions` implementation
# due to a bug in TightBinding.jl

positions(at::ASEAtoms) = at.po[:get_positions]()' |> vecs

function set_positions!(a::ASEAtoms, p::JVecsF)
   pold = positions(a)
   r = maxdist(pold, p)
   p_py = PyReverseDims(mat(p))
   a.po[:set_positions](p_py)
   update_transient_data!(a, r)
   return a
end

set_pbc!(at::ASEAtoms, val::Bool) = set_pbc!(at, (val,val,val))

function set_pbc!(at::ASEAtoms, val::AbstractVector{Bool})
   @assert length(val) == 3
   return set_pbc!(at, tuple(val...))
end

function set_pbc!(a::ASEAtoms, val::NTuple{3,Bool})
   a.po[:pbc] = val
   return a
end

pbc(a::ASEAtoms) = a.po[:pbc]

cell(at::ASEAtoms) = at.po[:get_cell]()

function set_cell!(a::ASEAtoms, p::Matrix)
   A = pinv(cell(a)) * p
   (a.po[:set_cell](p); a)
   update_transient_data!(a, vecnorm(A))
end

function deleteat!(at::ASEAtoms, n::Integer)
   at.po[:__delitem__](n-1) # delete in the actual array
   update_transient_data!(at, Inf)
   return at
end

function deleteat!(at::ASEAtoms, nn)
   deleteat!(at, collect(nn) - 1)
   return at
end

# ==========================================
#    some logic for storing transient data

function update_transient_data!(a::ASEAtoms, r::Real)
   for (key, t) in a.transient
      t.accum_change += r
      if t.accum_change + r >= t.max_change
         delete!(a.transient, key)
      end
   end
   return a
end

# Python arrays
has_array(a::ASEAtoms, name) = haskey(PyDict(a.po["arrays"]), name)
get_array(a::ASEAtoms, name) = a.po[:get_array](string(name))
function set_array!(a::ASEAtoms, name, value::Array)
   a.po[:set_array](string(name), value)
   return a
end

# Python info
has_info(a::ASEAtoms, name) = haskey(PyDict(a.po["info"]), name)
get_info(a::ASEAtoms, name) = PyDict(a.po["info"])[string(name)]
function set_info!(a::ASEAtoms, name, value::Any)
   PyDict(a.po["info"])[string(name)] = value
   return a
end

# Julia transient data
"""
`has_transient(a::ASEAtoms, name)`:
checks whether a transient datum with key `name` exists;
see also `set_transient!`.
"""
has_transient(a::ASEAtoms, name) = haskey(a.transient, name)

"""
`get_transient(a::ASEAtoms, name)`: retrieves a transient datum;
see also `set_transient!`.
"""
get_transient(a::ASEAtoms, name) = (a.transient[name]).data

"""
```
set_transient!(a::ASEAtoms, name, value, max_change=0.0)
```
Sets an arbitrary datum `value` under an arbitrary
key `name` in a special dictionary. The entry is paired with an
updatemeasure `chg`. Whenever atom positions change by a distance
`d`, the counter is incremented `chg += d`. When `chg > max_change`,
the entry is deleted from the dictionary.

The most common use of this functionality will be with `max_change = 0.0`;
in this case the entry will be deleted whenever atom positions
change.

Two examples where `max_change > 0`:
 * neighbourlist with buffer
 * Preconditioner is updated only when atom positions change significantly
"""
function set_transient!(a::ASEAtoms, name, value, max_change=0.0)
   a.transient[name] = TransientData(max_change, 0.0, value)
   return a
end

get_data(at::ASEAtoms, args...) = get_transient(at, args...)
has_data(at::ASEAtoms, args...) = has_transient(at::ASEAtoms, args...)
set_data!(at::ASEAtoms, args...) = set_transient!(at::ASEAtoms, args...)

# special arrays: momenta, velocities, masses, chemical_symbols
"get the momenta array"
momenta(at::ASEAtoms) = at.po[:get_momenta]()' |> vecs
"set the momenta array"
set_momenta!(at::ASEAtoms, p::JVecsF) = at.po[:set_momenta](p |> mat |> PyReverseDims)
"get the velocities array (convert from momenta)"
velocities(at::ASEAtoms) = at.po[:get_velocities]()' |> vecs
"convert to momenta, then set the momenta array"
set_velocities!(at::ASEAtoms, v::JVecsF) = at.po[:set_velocities](v |> mat |> PyReverseDims)
"get Vector of atom masses"
masses(at::ASEAtoms) = at.po[:get_masses]()
"set atom mass array as Vector{Float64}"
set_masses!(at::ASEAtoms, m::Vector{Float64}) = at.po[:set_masses](m)
"return vector of chemical symbols as strings"
chemical_symbols(at::ASEAtoms) = pyobject(at)[:get_chemical_symbols]()
"set the chemical symbols"
set_chemical_symbols!{T <: AbstractString}(at::ASEAtoms, s::Vector{T}) =
   pyobject(at)[:set_chemical_symbols](s)
"return vector of atomic numbers"
atomic_numbers(at::ASEAtoms) = pyobject(at)[:get_atomic_numbers]()

# ========================================================================
#     some nice Atoms generating and manipulating

"""
`repeat(a::ASEAtoms, n::(Int64, Int64, Int64)) -> ASEAtoms`

Takes an `ASEAtoms` configuration / cell and repeats is n_j times
into the j-th dimension. For example,
```
    atm = repeat( KeyError: key "z" not found"C"), (3,3,3) )
```
creates 3 x 3 x 3 unit cells of carbon.
"""
repeat(a::ASEAtoms, n::NTuple{3, Int64}) = ASEAtoms(a.po[:repeat](n))

repeat(a::ASEAtoms, n::AbstractArray) = repeat(a, tuple(n...))

import Base.*
*(at::ASEAtoms, n) = repeat(at, n)
*(n::NTuple{3, Int64}, at::ASEAtoms) = repeat(at, n)
*(at::ASEAtoms, n::Integer) = repeat(at, (n,n,n))
*(n::Integer, at::ASEAtoms) = repeat(at, (n,n,n))

export graphene_nanoribbon, nanotube, molecule

# @doc ase_build.bulk[:__doc__] ->
bulk(s::AbstractString, args...; pbc=true, kwargs...) =
   set_pbc!(ASEAtoms(ase_build.bulk(s, args...; kwargs...)), pbc)

# @doc ase_build.graphene_nanoribbon[:__doc__] ->
graphene_nanoribbon(args...; kwargs...) =
   ASEAtoms(ase_build.graphene_nanoribbon(args...; kwargs...))

"nanotube(n, m, length=1, bond=1.42, symbol=\"C\", verbose=False)"
function nanotube(args...; kwargs...)
   at = ASEAtoms(ase_build.nanotube(args...; kwargs...))
   C = Matrix(cell(at))
   if det(C) ≈ 0.0
      C[1,1] = 1.0
      C[2,2] = 1.0
      @assert !(det(C) ≈ 0.0)
   end
   set_cell!(at, C)
   return at
end

# @doc ase_build.molecule[:__doc__] ->
molecule(args...; kwargs...) =
   ASEAtoms(ase_build.molecule(args...; kwargs...))

# ############################################################
# # matscipy neighbourlist functionality
# ############################################################
#
# include("MatSciPy.jl")
#
# function neighbourlist(at::ASEAtoms, cutoff::Float64;
#                         recompute=false)::MatSciPy.NeighbourList
#    # TODO: also recompute if rcut is different !!!!!
#    # if no previous neighbourlist is available, compute a new one
#    if !has_transient(at, (:nlist, cutoff)) || recompute
#       # this nlist will be destroyed as soon as positions change
#       set_transient!(at, (:nlist, cutoff), MatSciPy.NeighbourList(at, cutoff))
#    end
#    return get_transient(at, (:nlist, cutoff))
# end
#
# """
# `static_neighbourlist(at::ASEAtoms, rcut::Float64)`
#
# This function first checks whether a static neighbourlist already exists
# with cutoff `rcut` and if it does then it returns the existing list.
# If it does not, then it computes a new neighbour list with the current
# configuration, stores it for later use and returns it.
# """
# function static_neighbourlist(at::ASEAtoms, rcut::Float64)
#    if !has_transient(at, (:snlist, rcut))
#       set_transient!( at, (:snlist, rcut),
#                           MatSciPy.NeighbourList(at, rcut),
#                           Inf )    # Inf means this is never deleted!
#    end
#    return get_transient(at, (:snlist, rcut))
# end
#

######################################################
#    Attaching an ASE-style calculator
#    ASE-style aliases
######################################################

"""
abstract type for all calculators that interface to ASE
"""
abstract type AbstractASECalculator <: AbstractCalculator end

"""
Concrete subtype of ASECalculator for classical potentials
"""
mutable struct ASECalculator <: AbstractASECalculator
   po::PyObject
end

function set_calculator!(at::ASEAtoms, calc::AbstractASECalculator)
   at.po[:set_calculator](calc.po)
   at.calc = calc
   return at
end

set_calculator!(at::ASEAtoms, po::PyObject) =
      set_calculator!(at, AbstractASECalculator(po))

forces(calc::AbstractASECalculator, at::ASEAtoms) = calc.po[:get_forces](at.po)' |> vecs
energy(calc::AbstractASECalculator, at::ASEAtoms) = calc.po[:get_potential_energy](at.po)

function virial(calc::AbstractASECalculator, at::ASEAtoms)
    s = calc.po[:get_stress](at.po)
    vol = det(cell(at))
    if size(s) == (6,)
      # unpack stress from compressed Voigt vector form
      s11, s22, s33, s23, s13, s12 = s
      return -JMatF([s11 s12 s13;
                     s12 s22 s23;
                     s13 s23 s33]) * vol
    elseif size(s) == (3,3)
      return -JMatF(s) * vol
    else
      error("got unxpected size(stress) $(size(stress)) from ASE")
    end
end

"""
Creates an `ASECalculator` that uses `ase.calculators.emt` to compute
energy and forces. This is very slow and is only included for
demonstration purposes.
"""
function EMTCalculator()
   @pyimport ase.calculators.emt as emt
   return ASECalculator(emt.EMT())
end

################### extra ASE functionality  ###################
# TODO: we probably want more of this
#       and a more structured way to translate
#

"""
* `extend!(at::ASEAtoms, atadd::ASEAtoms)::ASEAtoms`
* `extend!(at::ASEAtoms, S::AbstractString, x::JVecF)::ASEAtoms`

add `atadd` atoms to `at` and returns `at`; only `at` is modified.

A short variant is
```julia
#  extend!(at, (s, x))  <<<< deprecated
extend!(at, s, x)
```
where `s` is a string, `x::JVecF` a position
"""
function extend!(at::ASEAtoms, atadd::ASEAtoms)::ASEAtoms
   at.po[:extend](atadd.po)
   return at
end

function extend!{S <: AbstractString}(at::ASEAtoms, atnew::Tuple{S,JVecF})
   warn("`extend!(at, (s,x))` is deprecated; use `extend!(at, s, x)` instead")
   extend!(at, atnew[1], atnew[2])
end

extend!(at::ASEAtoms, S::AbstractString, x::JVecF) = extend!(at, ASEAtoms(S, [x]))

"""
* `write_xyz(filename, at, mode=:write)` : write atoms object to `filename`
* `write_xyz(filehandle, at)` : write atoms object as xyz file
* `write_xyz(filename, ats::Vector{ASEAtoms}, mode=:write)` : write a time series to a file
* `write_xyz(filename, at, x::Vector{Dofs}, mode=:write)` : write a time series to a file

to append to an existing file, use `:append` or `"a"` instead of `:write`.
"""
write_xyz(filename::AbstractString, at::ASEAtoms, mode=:write) =
   mode == :write ? ase_io.write(filename, at.po) : write_xyz(filename, [at], mode)

write_xyz(filehandle::PyObject, at::ASEAtoms) = ase_io.write(filehandle, at.po, format="xyz")

# open and close files from Python (to get a python filehandle)
pyopenf(filename::AbstractString, mode::AbstractString) = py"open($(filename), $(mode))"
pyclosef(filehandle) = filehandle[:close]()

function write_xyz(filename::AbstractString, at::ASEAtoms, xs::AbstractVector{Dofs}, mode=:write)
   x0 = dofs(at) # save the dofs
   filehandle = pyopenf(filename, string(mode)[1:1])
   for x in xs
     write_xyz(filehandle, set_dofs!(at, x))
   end
   pyclosef(filehandle)
   set_dofs!(at, x0)   # revert to original configuration
end

function write_xyz(filename::AbstractString, ats::AbstractVector{ASEAtoms}, mode=:write)
   filehandle = pyopenf(filename, string(mode)[1:1])
   for at in ats
      write_xyz(filehandle, at)
   end
   pyclosef(filehandle)
end

read_xyz(filename::AbstractString) = ASEAtoms(ase_io.read(filename))

# ------------------ Conversion to and from ASE Objects -----------------
#
# TODO: write tests for consistency of these conversions
#       for the new test_ase code  \\  X, P, M, Z, cell, pbc

Atoms(at_ase::ASE.ASEAtoms) =
   Atoms( positions(at_ase),
          momenta(at_ase),
          masses(at_ase),
          atomic_numbers(at_ase),
          JMat{Float64}(cell(at_ase)),
          pbc(at_ase);
          calc = calculator(at_ase),
          cons = constraint(at_ase) )

function ASEAtoms(at::Atoms)
   # simplify by assuming there is only one species
   @assert length(unique(chemical_symbols(at))) == 1
   sym = chemical_symbols(at)[1]
   at_ase = ASEAtoms("$sym$(length(at))")
   set_positions!(at_ase, positions(at))
   set_momenta!(at_ase, momenta(at))
   set_cell!(at_ase, Matrix(cell(at)))
   set_pbc!(at_ase, tuple(pbc(at)...))
   set_calculator!(at_ase, calculator(at))
   set_constraint!(at_ase, constraint(at))
end

# ----------- The Julia implementation of EMT (kind of for fun) ----------
# include("EMT.jl")

# # -------------- JuLIP NeighbourList Patch -------------
# using PyCall
# import NeighbourLists
# matscipy_neighbours = pyimport("matscipy.neighbours")
# function asenlist(at::Atoms, rcut)
#    pyat = ASEAtoms(at).po
#    return matscipy_neighbours[:neighbour_list]("ijdD", pyat, rcut)
# end
#
# function matscipy_nlist(at::Atoms{T}, rcut::T; recompute=false, kwargs...) where T <: AbstractFloat
#    i, j, r, R = asenlist(at, rcut)
#    i = copy(i)+1
#    j = copy(j)+1
#    r = copy(r)
#    R = vecs(copy(R'))
#    first = NeighbourLists.get_first(i, length(at))
#    NeighbourLists.sort_neigs!(j, r, R, first)
#    return NeighbourLists.PairList(positions(at), rcut, i, j, r, R, first)
# end

# -------------- Calling a JuLIP Calculator on ASEAtoms --------------

energy(V::AbstractCalculator, at::ASEAtoms) = energy(V, Atoms(at))
forces(V::AbstractCalculator, at::ASEAtoms) = forces(V, Atoms(at))
virial(V::AbstractCalculator, at::ASEAtoms) = virial(V, Atoms(at))
stress(V::AbstractCalculator, at::ASEAtoms) = stress(V, Atoms(at))

end # module


text:
using Test
using LogicCircuits

@testset "Circuit saver tests" begin

  circuit, vtree = load_struct_smooth_logical_circuit(zoo_lc_file("mnist-large.circuit"), zoo_vtree_file("balanced.vtree"))

  mktempdir() do tmp
    # the circuit above does not conform the the SDD requirements
    @test_throws Exception save_circuit("$tmp/temp.sdd", circuit, vtree)
  end

  # TODO add a test to load and save and load an .sdd file
  # currently we have no sdd+vtree in the model zoo to do this

end