text:
module DeduplicateFiles

using CRC

export DeduplicationFile, same_file, identical_files,
       delete_duplicate_file, indexfiles, deduplicate_files

"A structure with data about files to keep/delete."
type DeduplicationFile
  start::String
  realstart::String
  relpath::String
  realpath::String
  dirname::String
  dirinode::UInt64
  basename::String
  stat::Base.Filesystem.StatStruct
end

"A hash that is unique for the file"
file_id(x::DeduplicationFile) = hash(x.basename, hash(x.dirinode, hash(x.stat.device)))

"Determine if two paths are hard-linked to the same inode"
function is_hardlink(a,b)
  sa = stat(a)
  sb = stat(b)
  return sa.device == sb.device && sa.inode == sb.inode
end
is_hardlink(a::DeduplicationFile,b::DeduplicationFile) =
  a.stat.device == b.stat.device && a.stat.inode == b.stat.inode

"""
Determine if two paths point to the same directory entry.

This function returns `true` iff `a` and `b` both refer to the same directory
entry. This can happen if directories on either path are symbolic or hard
links to directories on the other path.

If `a` and `b` point to the same directory entry, then deleting `a` will also
delete `b`, so it would be bad to mistake `b` for a backup copy of `a`.

If `a` and/or `b` are symbolic links, then this function checks whether the
*links* are stored in the same directory entry. It does *not* check if they
point to the same or different files.

NOTE: This has *never been tested* on a file-system with hard-linked
directories. (Why would you use hard-linked directories?)
"""
function same_file(a::AbstractString, b::AbstractString)
  basename(a) == basename(b) || return false
  dira = dirname(abspath(a))
  dirb = dirname(abspath(b))
  realpath(dira) != realpath(dirb) || return true
  sa = stat(dira)
  sb = stat(dirb)
  return sa.device == sb.device && sa.inode == sb.inode
end

"Determine if two files are byte-for-byte identical"
function identical_files(a::IO, b::IO; buflen::Int=1048576)
  stat(a).size == stat(b).size || return false
  # TODO: Option to pass buffers as arguments and avoid repeated allocation.
  buf_a = Array{UInt8}(buflen)
  buf_b = Array{UInt8}(buflen)
  while (n_a = readbytes!(a, buf_a)) > 0
    n_b = readbytes!(b, buf_b)
    n_a == n_b || return false
    view(buf_a, 1:n_a) == view(buf_b, 1:n_b) || return false
  end
  return true
end
function identical_files(a::AbstractString, b::AbstractString; kwargs...)
  !same_file(a,b) || throw(ArgumentError("Cannot compare a file to itself."))
  result = false
  open(a, "r") do af
    open(b, "r") do bf
      result = identical_files(af, bf; kwargs...)
    end
  end
  return result
end

"""
   `delete_duplicate_file(; delete=file_to_delete, keep=file_to_keep)`

Delete a file if, and only if, it is an identical copy of another file.

The files to delete and keep are compared byte-for-byte. If they are
identical, then `file_to_delete` is deleted, and the function returns `true`.

If the files differ, or if both file paths refer to the same file (e.g. via
symbolic links) then no file is removed, and the function returns `false`.

The `file_to_delete` may not be a symbolic link. (The `file_to_keep` may be a
symbolic link, but if it points to the `file_to_delete`, then, obviously, that
file will not be deleted.)

This function cannot be called as `delete_duplicate_file("file1", "file2")`
because there might be confusion as to which file should be kept.

`delete_duplicate(..., dry_run=true)` does not delete files, but only returns
`true` or `false` depending on if a file would have been deleted.

`delete_duplicate(..., replace_with=:symlink)` replaces the deleted file with
a symbolic link to the `file_to_keep`.

`delete_duplicate(..., replace_with=:hardlink)` replaces the deleted file with
a hard link. If the files are not on the same file system, then an error is
thrown and no file is deleted. (Requires `ln` to be available via system call.)

`delete_duplicate(..., delete_hardlinks=true)` deletes the `file-to-delete`
even if it is a hardlink to `file-to-keep` (and therefore didn't occupy any
extra disk space).
"""
function delete_duplicate_file(; delete::AbstractString="", keep::AbstractString="",
                               dry_run::Bool=false, replace_with::Symbol=:nothing,
                               delete_hardlinks::Bool=false, verbose::Bool=false)
  (!isempty(delete) && !isempty(keep)) || throw(ArgumentError("Missing file name."))
  !islink(delete) || error("The file to delete must not be as symbolic link.")
  replace_with in [:nothing, :symlink, :hardlink] ||
      throw(ArgumentError("Illegal value for replace_with: ", string(replace_with)))

  # Follow symlinks now, so that we're not caught by a symlink that changes
  # between comparing the files and deleting.
  delete = realpath(delete)
  keep = realpath(keep)

  replace_with != :hardlink || stat(delete).device == stat(keep).device ||
      error("Cannot create hard-links across file systems.")

  !same_file(delete, keep) || return false

  hl = is_hardlink(delete, keep)

  hl || identical_files(delete, keep) || return false

  verbose && println(dry_run?"Would delete ":"Deleting ", delete, " as duplicate of ", keep, " (", stat(delete).size, " byters)")

  if !dry_run && (!hl || ( hl && delete_hardlinks && replace_with != :hardlink ))
    rm(delete)
    if replace_with == :symlink
      symlink(keep, delete)
    elseif replace_with == :hardlink
      run(`ln $keep $delete`)
    end
  end

  # If deleting the `delete` file somehow removed the `keep` file,
  # then there's a serious problem. Stop immediately!
  # (This should never happen.)
  isfile(keep) || error("Deleting ", delete, " somehow removed ", keep)

  return true
end

"""
   `indexfiles(start; idx, follow_symlinks=false)`

Create an index of all files found when walking the directory tree.
Returns an iterator over `DeduplicationFile` objects.
"""
function indexfiles(start::AbstractString; idx=Dict{UInt64,DeduplicationFile}(),
      follow_symlinks=false)
  realstart = realpath(start)
  for (root, dirs, files) in walkdir(realstart, follow_symlinks=follow_symlinks)
    ds = stat(root)
    for file in files
      ff = joinpath(root, file)
      if !islink(ff)
        df = DeduplicationFile(start, realstart, relpath(ff, realstart),
            realpath(ff), dirname(realpath(ff)), ds.inode, file, stat(ff))
        idx[file_id(df)] = df
      end
    end
  end
  return values(idx)
end

function indexfiles(starts::AbstractArray; idx=Dict{UInt64,DeduplicationFile}(),
    follow_symlinks=false)
  for start in starts
    indexfiles(start, idx=idx, follow_symlinks=follow_symlinks)
  end
  return values(idx)
end

"""
   `processdups(by, proc, list)`

Apply the `by` function to each item in `list`, for every set of two or more
items that generated identical `by` values, call `proc` on that set.
The `proc` function should return a list, and processdups returns a
concatenation of all returned lists.
"""
function processdups(by::Function, proc::Function, list::AbstractArray; quickpair::Bool=false)
  if quickpair && length(list)==2
    return(proc(list))
  end

  ret = Array{Any,1}[]

  length(list) >= 2 || return ret

  bylist = map(i->(i, by(list[i])), 1:length(list));

  # First sort the list...
  sort!(bylist, by=i->i[2])

  # ...then go through the list. Duplicates will be next to each other.
  ci = 1;
  cmp = bylist[ci][2]
  for i=1:length(bylist)
    if i==length(bylist) || bylist[i+1][2] != cmp
      if i>ci
        sublist = map(j->list[bylist[j][1]], ci:i)
        ret = vcat(ret, proc(sublist))
      end
      ci = i+1;
      if i<length(bylist)
        cmp = bylist[ci][2];
      end
    end
  end
  return ret
end

"""
   `deldups(list, dfun; dry_run)`

Determine which files from a list of suspected duplicates that should be
deleted according to the decision function `dfun`, and delete them if they
actually are duplicates.
"""
function deldups(list, dfun; kwargs...)
  deleted = Array{Tuple{DeduplicationFile, DeduplicationFile}, 1}();
  sort!(list, lt=dfun)
  lmax = length(list)
  i = 1
  while i<lmax
    for j=length(list):-1:i+1
      if dfun(list[i], list[j]) && delete_duplicate_file(delete=list[i].realpath,
             keep=list[j].realpath; kwargs...)
        push!(deleted, (list[i], list[j]))
        lmax = j
        break
      end
    end
    i = i+1
  end
  return deleted
end

const crc32 = crc(CRC_32)

"Run crc32 on the first `buflen` bytes of a file"
function partialcrc32(a::IO; buflen::Int=1048576)
  buf_a = Array{UInt8}(buflen)
  n_a = readbytes!(a, buf_a)
  if n_a != length(buf_a)
    buf_a = buf_a[1:n_a]
  end
  return crc32(buf_a)
end

"Divide list by crc, and run deldups on each sub-list"
crcdeldups(list, dfun; kwargs...) =
  processdups(x->open(crc32, x.realpath), lst->deldups(lst, dfun; kwargs...), list, quickpair=true)

"Divide list by partial crc, and run crcdeldups on each sub-list"
pcrcdeldups(list, dfun; kwargs...) =
  processdups(x->open(partialcrc32, x.realpath), lst->crcdeldups(lst, dfun; kwargs...), list)

"Divide list by size, then run pcrcdeldups/crcdeldups on each sub-sub-list"
sizcrcdeldups(list, dfun; kwargs...) =
  processdups(x->x.stat.size, lst->lst[1].stat.size>1048576?pcrcdeldups(lst, dfun; kwargs...):crcdeldups(lst, dfun; kwargs...), list)

"""
   `list = deduplicate_files(startdirs, dfun)`

Scan `startdirs` recursively for files, call dfun on pairs of suspected
duplicates, and delete identical copies when dfun returns true.
"""
function deduplicate_files(startdirs, dfun; verbose::Bool=false, kwargs...)
  list = indexfiles(startdirs)
  if verbose
    println("Indexed ", length(list), " files.")
  end
  push!(kwargs, (:verbose, verbose))
  sizcrcdeldups(collect(list), dfun; kwargs...)
end

# Define a `show` method for tuples of deduplication files, to make it easier to
# read the returned array from `deduplicate_files`.
import Base.show
show(io::IO, t::Tuple{DeduplicationFile, DeduplicationFile}) = print(io,
  joinpath(t[1].start, t[1].relpath), " duplicate of ", joinpath(t[2].start, t[2].relpath))

"The total size (in bytes) of an array of `DeduplicationFile`s"
totalsize(x::DeduplicationFile) = x.stat.size
totalsize(x::Tuple{DeduplicationFile, DeduplicationFile}) = totalsize(x[1])
totalsize(x::AbstractArray) = sum(totalsize.(x))

# TODO: Check with dfun that at least one duplicate might be deleted
#       before doing CRC.

end # module


text:
using CSV, DataFrames, Feather, Gadfly, IJulia, JSON, RDatasets, StatsBase

# data = readdlm("iris_dataset.csv", ',')

iris = CSV.read("iris_dataset.csv", header=true, delim=',', copycols = true)

print(iris)

print(describe(iris))

# print(describe(iris, stats=[:q25, :q75, :first, :last]))

# print(describe(iris, stats=:all))

# Pkg.add("Statistics")

using Statistics

print(mean(iris[:sepallength]))

for x in names(iris)[1:end-1]
    for y in names(iris)[1:end-1]
        println("$x /t $y /t $(cor(iris[x], iris[y]))")
    end
    println("--------------------------------------------------")
end

# columns to array
sepallength = Array(iris[:sepallength])
# print(sepallength)

# DataFrame to matrix
irisarr = convert(Array, iris[:,:])

iris[1:3, [:petallength, :petalwidth]]

# CSV.write("test_iris.csv", iris)

ext: jl


text:
@testset "Kriging" begin
    d, n = 3, 10
    x = 2π * rand(d, n)
    y1 = [sum(sin, view(x, :, i)) / d for i in 1:n]
    y2 = [sum(cos, view(x, :, i)) / d for i in 1:n]
    y = vcat(transpose(y1), transpose(y2))
    mZero = MeanZero()
    kern = SE(0.0, 0.0)
    ntest = 5
    xtest = randn(d, ntest)
    yt1 = [sum(sin, view(xtest, :, i)) / d for i in 1:ntest]
    yt2 = [sum(cos, view(xtest, :, i)) / d for i in 1:ntest]
    ytest = vcat(transpose(yt1), transpose(yt2))
    krg = Kriging(x, y, mZero, kern)
    @testset "predict samples" begin
        y_pred, sig2 = predict_full(krg, x)
        @test maximum(abs, y - y_pred) ≈ 0.0 atol = 0.1
    end
    @testset "predict tests" begin
        y_pred, sig2 = predict_full(krg, xtest)
    end
    @testset "tune local" begin
        tune!(krg, ([-10.0, -10.0], [10.0, 10.0]), nothing; verbose = true)
    end
    @testset "tune" begin
        tune!(
            krg,
            ([-10.0, -10.0], [10.0, 10.0]),
            nothing;
            global_eval = 10000,
            verbose = true,
        )
    end
    @testset "append" begin
        append!(krg, xtest, ytest)
    end
    @testset "get" begin
        x1 = getx(krg, 1)
        y1 = gety(krg, 1)
        @test x1 ≈ x[:, 1] atol = 1e-6
        @test y1 ≈ y[:, 1] atol = 1e-6
    end
    @testset "get all" begin
        x1 = getx(krg)
        y1 = gety(krg)
        @test x1 ≈ hcat(x,xtest) atol = 1e-6
        @test y1 ≈ hcat(y,ytest) atol = 1e-6
    end
end


text:
module FastLevenshtein

    function fastlevenshtein(v1::Vector, v2::Vector, insert_cost::Real=1.0, delete_cost::Real=1.0, update_cost::Real=1.0)::Real
        """
        Efficiently Compute the Levenshtein (edit) distance between two vectors using dynamic programming.
        :param v1: The first vector is a Mx1 vector.
        :param v2: The second vector is a Nx1 vector.
        :return: The minimum number of edits required to make the vectors the same.
        """
        M::Int64, N::Int64 = length(v1), length(v2)

        # This table holds the subproblems we have to solve to get an answer
        table::Array{Real} = zeros(Real, (M+1, N+1))

        for i::Int64 ∈ 0:M
            for j::Int64 ∈ 0:N
                # if the length of the 1st vector is 0 then we need at least j edits
                if i == 0
                    table[i+1, j+1] = convert(AbstractFloat, j)
                # if the length of the 2nd vector is 0 then we need at least i edits
                elseif j == 0
                    table[i+1, j+1] = convert(AbstractFloat, i)
                # if the values are the same then the minumum number of edits hasn't changed
                elseif v1[i] == v2[j]
                    table[i+1, j+1] = table[i, j]
                # if the value are different then the minimum edits has increased by 1
                else
                    table[i+1, j+1] = min(
                        table[i+1, j] + insert_cost, # The minimum edits required if i > j (insert value @ current index)
                        table[i, j+1] + delete_cost, # The minimum edits required if i < j (remove value @ current index)
                        table[i, j]   + update_cost  # The minimum edits required if i = j (update value @ current index)
                    )
                end
            end
        end

        return table[M+1, N+1]
    end

    function fastlevenshtein(s1::AbstractString, s2::AbstractString, insert_cost::Real=1.0, delete_cost::Real=1.0, update_cost::Real=1.0)::Real
        """
        Efficiently Compute the Levenshtein (edit) distance between two strings using dynamic programming.
        :param s1: The first string contains M characters.
        :param s2: The second string contains N characters.
        :return: The minimum number of edits required to make the strings the same.
        """
        return fastlevenshtein(Vector{Char}(s1), Vector{Char}(s2), insert_cost, delete_cost, update_cost)
    end

    export fastlevenshtein
end


text:
export  fejer,
        fejer2,
        clenshaw_curtis,
        quadgp,
        golubwelsch,
        gauss,
        radau,
        lobatto
"""
    fejer(N::Int)
Fejer's first quadrature rule.
"""
function fejer(N::Int)
    @assert N >= 1 "N has to be positive"
    N == 1 && return zeros(1), [2.]
    θ = map(x->(2x-1)*pi/(2N),1:N)
    M = N ÷ 2
    return cos.(θ), map(x->2/N*(1-2*sum(cos(2*n*x)/(4*n^2-1) for n=1:M)),θ)
end

"""
    fejer2(n::Int)
Fejer's second quadrature rule according to [Waldvogel, J. Bit Numer Math (2006) 46: 195](https://doi.org/10.1007/s10543-006-0045-4).
"""
function fejer2(n::Int)
    @assert n >= 2
    N = 1:2:n-1
    m = n - length(N)
    v0 = push!(push!(map(x->2/x/(x-2.),N),1/last(N)), zeros(m)...)
    @inbounds v2 = -v0[1:n] - v0[n+1:-1:2];
    wf2 = real(ifft(v2))
    return map(x->cos(x*pi/n),0:n), push!(wf2, first(wf2))
end

"""
    clenshaw_curtis(n::Int)
Clenshaw-Curtis quadrature according to [Waldvogel, J. Bit Numer Math (2006) 46: 195](https://doi.org/10.1007/s10543-006-0045-4).
"""
function clenshaw_curtis(n::Int)
    @assert n >= 2
    N = 1:2:n-1
    l = length(N)
    m = n - l
    v0 = push!(push!(map(x->2/x/(x-2.),N),1/last(N)), zeros(m)...)
    @inbounds v2 = -v0[1:n] - v0[n+1:-1:2]
    g0= -ones(Float64,n)
    g0[1+l] += n
    g0[1+m] += n
    g = 1 / (n^2 - 1 + mod(n,2) ) * g0
    wcc = real(ifft(v2 + g));
    return map(i->cos(i/n*pi),0:n), push!(wcc,first(wcc))
end

"""
    quadgp(weight::Function,lb::Real,ub::Real,N::Int=10;quadrature::Function=clenshaw_curtis,bnd::Float64=Inf)
general purpose quadrature based on Gautschi, "Orthogonal Polynomials: Computation and Approximation", Section 2.2.2, pp. 93-95

Compute the `N`-point quadrature rule for `weight` with support (`lb`, `ub`).
The quadrature rule can be specified by the keyword `quadrature`.
The keyword `bnd` sets the numerical value for infinity.
"""
function quadgp(weight::Function,lb::Real,ub::Real,N::Int=10;quadrature::Function=clenshaw_curtis,bnd::Float64=Inf)
    @assert lb < ub "inconsistent interval bounds"
    t_fej,w_fej = quadrature(N)
    t, w, d = zeros(Float64,N),zeros(Float64,N),zeros(Float64,N)
    # transformation
    if (-bnd<lb<bnd) && (-bnd<ub<bnd)
        t = 0.5*((ub-lb)*t_fej .+ (ub+lb))
        d = 0.5*(ub-lb)*ones(Float64,size(t_fej))
    elseif -bnd>=lb && (-bnd<ub<bnd)
        t = ub .- (1 .- t_fej)./(1 .+ t_fej)
        d = 2 ./(1 .+ t_fej).^2
    elseif (-bnd<lb<bnd) && ub>=bnd
        t = lb .+ (1 .+ t_fej)./(1 .- t_fej)
        d = 2 ./(1 .- t_fej).^2
    elseif -bnd>=lb && ub>=bnd
        t = t_fej./(1 .- t_fej.^2)
        d = (1 .+ t_fej.^2)./( 1 .- t_fej.^2).^2
    end
    w = w_fej.*weight.(t).*d
    return t, w
end

####################################################
####################################################
####################################################
# Quadrature rules based on recurrence coefficients

function golubwelsch(α::AbstractVector{<:Real}, β::AbstractVector{<:Real}, maxiter::Int=30)
    N = length(α) - 1
    a, β0 = copy(α[1:N]), β[1]
    w = zero(a)
    special_eigenproblem!(a,copy(sqrt.(β)),w,maxiter)
    idx = sortperm(a)
    a[idx], (β0*w.^2)[idx]
end

golubwelsch(op::Union{OrthoPoly,AbstractCanonicalOrthoPoly}) = golubwelsch(op.α,op.β)

"""
    gauss(N::Int,α::AbstractVector{<:Real},β::AbstractVector{<:Real})
    gauss(α::AbstractVector{<:Real},β::AbstractVector{<:Real})
    gauss(N::Int,op::Union{OrthoPoly,AbstractCanonicalOrthoPoly})
    gauss(op::Union{OrthoPoly,AbstractCanonicalOrthoPoly})
Gauss quadrature rule, also known as Golub-Welsch algorithm

`gauss()` generates the `N` Gauss quadrature nodes and weights for a given weight function.
The weight function is represented by the `N` recurrence coefficients for the monic polynomials orthogonal
with respect to the weight function.

!!! note
    The function `gauss` accepts at most `N = length(α) - 1` quadrature points,
        hence providing at most an `(length(α) - 1)`-point quadrature rule.

!!! note
    If no `N` is provided, then `N = length(α) - 1`.
"""
function gauss(N::Int,α::AbstractVector{<:Real},β::AbstractVector{<:Real})
    N += 1
    @assert N > 0 "only positive N allowed"
    @assert length(α) == length(β) "inconsistent number of recurrence coefficients"
    N0 = length(α)
    @assert N0 >= N  "not enough recurrence coefficients"
    @inbounds golubwelsch(α[1:N],β[1:N])
end

gauss(α::AbstractVector{<:Real},β::AbstractVector{<:Real}) = gauss(length(α)-1,α,β)

gauss(N::Int,op::Union{OrthoPoly,AbstractCanonicalOrthoPoly}) = gauss(N::Int,op.α,op.β)

gauss(op::Union{OrthoPoly,AbstractCanonicalOrthoPoly}) = gauss(op.α,op.β)

"""
    radau(N::Int,α::AbstractVector{<:Real},β::AbstractVector{<:Real},end0::Real)
    radau(α::AbstractVector{<:Real},β::AbstractVector{<:Real},end0::Real)
    radau(N::Int,op::Union{OrthoPoly,AbstractCanonicalOrthoPoly},end0::Real)
    radau(op::Union{OrthoPoly,AbstractCanonicalOrthoPoly},end0::Real)
Gauss-Radau quadrature rule.
Given a weight function encoded by the recurrence coefficients `(α,β)`for the associated
orthogonal polynomials, the function generates the
nodes and weights `(N+1)`-point Gauss-Radau
quadrature rule for the weight function having a prescribed
node `end0` (typically at one of the end points of the support
interval of w, or outside thereof).

!!! note
    The function `radau` accepts at most `N = length(α) - 2` as an input,
        hence providing at most an `(length(α) - 1)`-point quadrature rule.

!!! note
    Reference: OPQ: A MATLAB SUITE OF PROGRAMS FOR GENERATING ORTHOGONAL POLYNOMIALS AND RELATED QUADRATURE RULES by Walter Gautschi
"""
function radau(N::Int,α::AbstractVector{<:Real},β::AbstractVector{<:Real},end0::Real)
    α_ = copy(α)
    @assert N > 0 "only positive N allowed"
    @assert length(α_) == length(β) > 0 "inconsistent number of recurrence coefficients"
    @assert length(α_) >= N + 1 "not enough recurrence coefficients"
    p0 = 0.
    p1 = 1.
    for n in Base.OneTo(N)
      pm1 = p0
      p0 = p1;
      @inbounds p1 = (end0 - α_[n])*p0 - β[n]*pm1;
    end
    @inbounds α_[N+1] = end0 - β[N+1]*p0/p1
    gauss(N+1 , α_, β)
end
radau(α::AbstractVector{<:Real},β::AbstractVector{<:Real},end0::Real) = radau(length(α)-2,α,β,end0)
radau(N::Int,op::Union{OrthoPoly,AbstractCanonicalOrthoPoly},end0::Real) = radau(N,op.α,op.β,end0)
radau(op::Union{OrthoPoly,AbstractCanonicalOrthoPoly},end0::Real) = radau(op.α,op.β,end0::Real)

"""
    lobatto(N::Int,α::AbstractVector{<:Real},β::AbstractVector{<:Real},endl::Real,endr::Real)
    lobatto(α::AbstractVector{<:Real},β::AbstractVector{<:Real},endl::Real,endr::Real)
    lobatto(N::Int,op::Union{OrthoPoly,AbstractCanonicalOrthoPoly},endl::Real,endr::Real)
    lobatto(op::Union{OrthoPoly,AbstractCanonicalOrthoPoly},endl::Real,endr::Real)
Gauss-Lobatto quadrature rule.
Given a weight function encoded by the recurrence coefficients for the associated
orthogonal polynomials, the function generates
the nodes and weights of the `(N+2)`-point Gauss-Lobatto
quadrature rule for the weight function, having two
prescribed nodes `endl`, `endr` (typically the left and right end
points of the support interval, or points to the left
resp. to the right thereof).

!!! note
    The function `radau` accepts at most `N = length(α) - 3` as an input,
    hence providing at most an `(length(α) - 1)`-point quadrature rule.

!!! note
    Reference: OPQ: A MATLAB SUITE OF PROGRAMS FOR GENERATING ORTHOGONAL POLYNOMIALS AND RELATED QUADRATURE RULES by Walter Gautschi
"""
function lobatto(N::Int,α_::AbstractVector{<:Real},β_::AbstractVector{<:Real},endl::Real,endr::Real)
    α, β = copy(α_), copy(β_)
    @assert N > 0 "only positive N allowed"
    @assert length(α) == length(β) > 0 "inconsistent number of recurrence coefficients"
    @assert endl < endr "inconsistent end points"
    @assert length(α) >= N+2 "not enough recurrence coefficients"
    p0l = 0.; p0r = 0.; p1l = 1.; p1r = 1.;
    for n=1:N+1
        pm1l=p0l
        p0l=p1l
        pm1r=p0r
        p0r=p1r
        @inbounds p1l=(endl-α[n])*p0l-β[n]*pm1l
        @inbounds p1r=(endr-α[n])*p0r-β[n]*pm1r
    end
    det = p1l * p0r - p1r * p0l
    @inbounds α[N+2] = (endl * p1l * p0r - endr * p1r * p0l) / det
    @inbounds β[N+2] = (endr - endl) * p1l * p1r / det;
    gauss(N+2,α,β)
end
lobatto(α::AbstractVector{<:Real},β::AbstractVector{<:Real},endl::Real,endr::Real) = lobatto(length(α)-3,α,β,endl,endr)
lobatto(N::Int,op::Union{OrthoPoly,AbstractCanonicalOrthoPoly},endl::Real,endr::Real) = lobatto(N,op.α,op.β,endl,endr)
lobatto(op::Union{OrthoPoly,AbstractCanonicalOrthoPoly},endl::Real,endr::Real) = lobatto(op.α,op.β,endl,endr)


text:
export rule

@rule BIFMHelper(:in, Marginalisation) (m_out::Any, ) = m_out


text:
using Pkg
using CUDA
using ScatterNNlib
using DataFrames
using CSV
using BenchmarkTools
using BenchmarkTools: Trial, TrialEstimate, median, mean

Pkg.status("ScatterNNlib")
Pkg.status("CUDA")

ENV["JULIA_NUM_THREADS"] = 1

d = 50
nbins = 20
getinfo(te::TrialEstimate) = te.time, te.gctime, te.memory
getstats(t::Trial) = [getinfo(minimum(t)), getinfo(mean(t)), getinfo(maximum(t))]

metadata = DataFrame(device=String[], dim=Int[], sample=Int[], bins=Int[])
mintime = DataFrame(min_time=Float64[], min_gc=Float64[], min_mem=Int[])
meantime = DataFrame(mean_time=Float64[], mean_gc=Float64[], mean_mem=Int[])
maxtime = DataFrame(max_time=Float64[], max_gc=Float64[], max_mem=Int[])

for l = [2^5, 2^10, 2^15, 2^20]
    hist = zeros(Float32, d, nbins)
    δ = rand(Float32, d, l)
    idx = rand(1:nbins, l)

    hist_gpu = CuArray(hist)
    δ_gpu = CuArray(δ)
    idx_gpu = CuArray(idx)

    scatter_add!(hist, δ, idx)
    scatter_add!(hist_gpu, δ_gpu, idx_gpu)

    b_cpu = @benchmark scatter_add!($hist, $δ, $idx)
    b_gpu = @benchmark scatter_add!($hist_gpu, $δ_gpu, $idx_gpu)
    s_cpu = getstats(b_cpu)
    s_gpu = getstats(b_gpu)

    push!(metadata, ("cpu", d, l, nbins))
    push!(mintime, s_cpu[1])
    push!(meantime, s_cpu[2])
    push!(maxtime, s_cpu[3])

    push!(metadata, ("gpu", d, l, nbins))
    push!(mintime, s_gpu[1])
    push!(meantime, s_gpu[2])
    push!(maxtime, s_gpu[3])
end

data = hcat(metadata, mintime, meantime, maxtime)
CSV.write("benchmark/scatter_julia.tsv", data; delim="\t")

## Benchmark
# @benchmark scatter_add!($hist, $δ, $idx)
# CuArrays.@time scatter_add!(hist_gpu, δ_gpu, idx_gpu)


text:
# ```@setup slurmgen
# # all examples are run on an example global.toml and downsized maps.
# ARGS = ["example.toml"]
# ```

configfile = ARGS[1]

# #
# # [SLURM Commands for Spectra (slurmgen.jl)](@id slurmgen)
# This command generates SLURM commands that executes [rawspectra.jl](@ref rawspectra)
# on all the pairs of maps in the config.
## this file just prints out the SLURM commands required to compute the spectra
using TOML
config = TOML.parsefile(configfile)

# Let's generate the commands we need for likelihood spectra.
mapids = [k for k in keys(config["map"])]
cmd = "sbatch scripts/8core2hr.cmd"
for i in 1:length(mapids)
    for j in i:length(mapids)
        println("$(cmd) \"julia src/rawspectra.jl $(configfile) $(mapids[i]) $(mapids[j])\"")
    end
end

# I just paste these in and wait a few hours. The resulting spectra are deposited in the
# `scratch` dir in the config.

# ## Fit Noise Model
# These SLURM commands fit the raw spectra with the camspec noise model.
## loop over freqs and noise channels
println()
for freq in ("100", "143", "217")
    for spec in ("TT", "EE")
        println("$(cmd) \"julia src/fitnoisemodel.jl $(configfile) $(freq) $(spec)\"")
    end
end

# ## Generating Signal-Only Simulations
# These SLURM commands correct the covariances for the insufficiently-apodized point source
# holes.
## loop over freqs and noise channels
println()
cmd = "sbatch scripts/8core2hr.cmd"
freqs = ("100", "143", "217")
for i in 1:3
    for j in i:3
        for spec in ("TT", "TE", "ET", "EE")
            freq1, freq2 = freqs[i], freqs[j]
            println("$(cmd) \"julia src/signalsim.jl $(configfile) $(freq1) $(freq2) $(spec) 5000\"")
        end
    end
end

# ## Fit Signal-Only Simulation Corrections
# These SLURM commands correct the covariances for the insufficiently-apodized point source
# holes.
## loop over freqs and noise channels
println()
freqs = ("100", "143", "217")
for i in 1:3
    for j in i:3
        for spec in ("TT", "TE", "ET", "EE")
            freq1, freq2 = freqs[i], freqs[j]
            println("$(cmd) \"julia src/corrections.jl $(configfile) $(freq1) $(freq2) $(spec)\"")
        end
    end
end

##
function unique_splits(X, Y, f1, f2)
    if (X == Y) && (f1 == f2)
        return [(1,2)]
    elseif (X==Y) && (f1 != f2)
        return [(1,2), (2,1)]
    else
        return [(1,1), (1,2), (2,1), (2,2)]
    end
end

using Test
using Base.Iterators

@test unique_splits("T", "T", "100", "100") == [(1,2)]
@test unique_splits("T", "E", "100", "100") == [(1,1), (1,2), (2,1), (2,2)]
@test unique_splits("T", "T", "100", "143") == [(1,2), (2,1)]
@test unique_splits("E", "E", "100", "143") == [(1,2), (2,1)]
@test unique_splits("T", "E", "100", "143") == [(1,1), (1,2), (2,1), (2,2)]

include("util.jl")
specs = (
    (:TT,"100","100"), (:TT,"143","143"), (:TT,"143","217"), (:TT,"217","217"),
    (:EE,"100","100"), (:EE,"100","143"), (:EE,"100","217"), (:EE,"143","143"),
    (:EE,"143","217"), (:EE,"217","217"),
    (:TE,"100","100"), (:TE,"100","143"), (:TE,"100","217"), (:TE,"143","143"),
    (:TE,"143","217"), (:TE,"217","217"),
    (:ET,"100","143"), (:ET,"100","217"), (:ET,"143","217")  # combined with TE in plic
)

constituents = []

for spec in specs
    AB, f1, f2 = spec
    A, B = string(AB)
    for (s1, s2) in unique_splits(A, B, f1, f2)
        push!(constituents, (A, B, f1, f2, s1, s2))
    end
end

run_name = config["general"]["name"]
cmd = "sbatch scripts/8core2hr.cmd"
scriptpath = joinpath(config["scratch"], "scripts")
mkpath(scriptpath)
open(joinpath(scriptpath, "$(run_name)_gen_covmats.sh"), "w") do f
    nspecs = length(constituents)

    ## generate the upper triangle
    for i1 in 1:nspecs
        for i2 in i1:nspecs
            A, B, f1, f2, s1, s2 = constituents[i1]
            C, D, f3, f4, s3, s4 = constituents[i2]
            write(f, "$(cmd) \"julia src/covmat.jl $(configfile) $f1 $f2 $f3 $f4 $A $B $C $D $s1 $s2 $s3 $s4\"\n")
        end
    end
end

max_forks_repo_name: sgiardie/PSpipe


text:
using PyPlot

vars=[]
for i in 1:length(emu_out.population)
    m, v = gp_regression(emu_out.population[i], emu_out.emulators[i])
    push!(vars, v)
end

ion()
for i in 1:length(emu_out.threshold_schedule)
    epsilon = emu_out.threshold_schedule[i]
    figure()
    scatter3D(emu_out.population[i][:, 3], emu_out.population[i][:, 2], vars[i])
    title("Variance #$(i); ϵ=$(epsilon)")
    zlabel("Variance")
    xlabel("Parameter 3")
    ylabel("Parameter 2")
    figure(figsize=(13,4.8))
    subplot(121)
    scatter(emu_out.population[i][:, 3], emu_out.population[i][:, 2], marker=".")
    xlabel("Parameter 3")
    ylabel("Parameter 2")
    title("Accepted population #$(i); ϵ=$(epsilon)")
    subplot(122)
    scatter(emu_out.emulators[i].gp_training_x[:, 3], emu_out.emulators[i].gp_training_x[:, 2], marker="x")
    xlabel("Parameter 3")
    ylabel("Parameter 2")
    title("Emulator #$(i) training data")
end


text:
#   Adapted from
#   https://github.com/dabeaz/bitey/blob/master/test/test_M.py
#
#   Copyright (C) 2012
#   David M. Beazley (Dabeaz LLC), http://www.dabeaz.com
#   BSD 3-clause "New" or "Revised" License

using Base.Test
using BitcodeRunner

M = bitcode_library("ctest.bc")

@test M.add_char(4,5) == Int8(9)
@test M.add_short(4,5) == Int16(9)
@test M.add_long(4,5) == Int32(9)
@test M.add_longlong(4,5) == Int64(9)

@test M.add_float(4,5) == 9f0
@test M.add_double(4,5) == 9.0

# Test passing Ref's and Arrays
a = Ref(Int16(2))
M.mutate_short(a)
@test a[] == Int16(4)
a = Ref(Int32(2))
M.mutate_long(a)
@test a[] == Int32(4)
a = Ref(Int64(2))
M.mutate_longlong(a)
@test a[] == Int64(4)
a = Ref(2f0)
M.mutate_float(a)
@test a[] == 4f0
a = Ref(2.0)
M.mutate_double(a)
@test a[] == 4.0
a = [2.0]
M.mutate_double(a)
@test a[1] == 4.0

# Array tests
@test M.arr_sum_int(collect(UInt32(1):UInt32(4))) == Int32(10)
@test M.arr_sum_double(collect(1.0:4.0)) == 10.0

# Struct test
p = Ref(M.Point(3, 4))
q = Ref(M.Point(6, 8))
@test M.distance(p, q) == 5.0


text:
for N in [Float64, Rational{Int}, Float32]
    B1 = BallInf(N[0, 0], N(1))
    B2 = BallInf(N[4, -4], N(1))
    B3 = BallInf(N[1, -1], N(1))
    C = Complement(B1)

    # double-complement is the identity
    @test Complement(C) == B1

    # dimension
    @test dim(C) == 2

    # membership
    @test N[1, 1] ∉ C && N[2, 2] ∈ C

    # emptiness
    @test !isempty(C) && isempty(Complement(Universe(2)))

    # inclusion
    subset, point = ⊆(B2, C, true)
    @test B2 ⊆ C && subset && point == N[]
    for X in [B1, B3]
        subset, point = ⊆(X, C, true)
        @test X ⊈ C && !subset && point ∈ X && point ∉ C
    end

    # isdisjoint
    for X in [B2, B3]
        res, w = isdisjoint(X, C, true)
        @test !isdisjoint(X, C) && !res && w ∈ X && w ∈ C
        res, w = isdisjoint(C, X, true)
        @test !isdisjoint(C, X) && !res && w ∈ X && w ∈ C
    end
    res, w = isdisjoint(B1, C, true)
    @test isdisjoint(B1, C) && res && w == N[]
    res, w = isdisjoint(C, B1, true)
    @test isdisjoint(C, B1) && res && w == N[]

    # test convexity from the type
    @test isconvextype(typeof(Complement(Universe{N}(2))))
    @test isconvextype(typeof(Complement(EmptySet{N}(2))))
    @test isconvextype(typeof(Complement(HalfSpace(N[1], N(0)))))

    # concrete complement
    H = HalfSpace(N[1, 0], N(1)) # x <= 1
    @test complement(H) == HalfSpace(N[-1, 0], N(-1))
    @test constraints_list(Complement(H)) == [HalfSpace(N[-1, 0], N(-1))]
    # test fallback for <: LazySet
    c = complement(N[0 1; 1 0] * H)
    @test c isa UnionSetArray && length(array(c)) == 1
    @test first(array(c)) == HalfSpace(N[0, -1], N(-1)) # complement of y <= 1 is y >= 1

    # boundedness
    @test isboundedtype(typeof(Complement(Universe{N}(2))))
    @test !isboundedtype(typeof(Complement(EmptySet{N}(2))))
end


text:
using GeometryTypes, Colors, GLVisualize, GLWindow, FileIO, GLAbstraction, ModernGL

function read_ifs(filename)
    open(filename) do io
        function str()
            n = read(io, UInt32)
            String(read(io, UInt8, n))
        end
        ifs = str()
        zero = read(io, Float32)
        name = str()
        vertices = str()
        nverts = read(io, UInt32)
        verts = Vector{Point3f0}(nverts)
        for i = 1:nverts
            verts[i] = Point3f0(
                read(io, Float32),
                read(io, Float32),
                read(io, Float32)
            )
        end
        tris = str()
        nfaces = read(io, UInt32)
        faces = Vector{GLTriangle}(nfaces)
        for i = 1:nfaces
            faces[i] = GLTriangle(
                read(io, UInt32),
                read(io, UInt32),
                read(io, UInt32)
            )
        end
        (verts, faces)
    end
end

window = glscreen()

function display_data(folder)
    meshpaths = filter(x-> endswith(x, ".ifs"), readdir(folder))[1:1024]
    v0 = (Point3f0(typemax(Float32)), Point3f0(typemin(Float32)))
    faces = GLTriangle[]
    verts = Point3f0[]
    norms = Normal{3, Float32}[]
    lastvidx = 0
    for (i, meshpath) in enumerate(meshpaths)
        vs, fs = read_ifs(joinpath(folder, meshpath))
        mini, maxi = extrema(vs)
        ns = normals(vs, fs)
        x, y = ind2sub((32, 32), i)
        s = maximum(maxi .- mini)
        append!(verts, vs .+ Point3f0(x, y, 0f0))
        append!(norms, ns)
        append!(faces, fs .+ lastvidx)
        lastvidx += length(fs) * 3
    end
    _view(visualize(
        GLNormalMesh(faces = faces, vertices = verts, normals = norms),
        color = RGBA{Float32}(rand(RGB{Float32}), 1f0)
    ))
end
display_data(homedir() * "/3dstuff/models")

function rendloop(window, frame_times)
    while isopen(window)
        tic()
        glFinish()
        render_frame(window)
        GLWindow.swapbuffers(window)
        GLWindow.poll_glfw()
        GLWindow.poll_reactive()
        push!(frame_times, toq())
    end
    GLFW.DestroyWindow(window)
    window.handle = C_NULL
    frame_times
end
timsies = Float64[]
rendloop(window, timsies)

open("/home/s/.julia/v0.6/Visualize/test2.csv", "w") do io
    for t in timsies
        println(io, t)
    end
end
println(pwd())


text:
@testset "memory" begin

let
    a,b = Mem.info()
    # NOTE: actually testing this is pretty fragile on CI
    #=@test a == =# CUDAdrv.available_memory()
    #=@test b == =# CUDAdrv.total_memory()
end

# dummy data
T = UInt32
N = 5
data = rand(T, N)
nb = sizeof(data)

# buffers are untyped, so we use a convenience function to get a typed pointer
# we prefer to return a device pointer (for managed buffers) to maximize CUDAdrv coverage
typed_pointer(buf::Union{Mem.Device, Mem.Unified}, T) = convert(CuPtr{T}, buf)
typed_pointer(buf::Mem.Host, T)                       = convert(Ptr{T},   buf)

# allocations and copies
for srcTy in [Mem.Device, Mem.Host, Mem.Unified],
    dstTy in [Mem.Device, Mem.Host, Mem.Unified]

    dummy = Mem.alloc(srcTy, 0)
    Mem.free(dummy)

    src = Mem.alloc(srcTy, nb)
    unsafe_copyto!(typed_pointer(src, T), pointer(data), N)

    dst = Mem.alloc(dstTy, nb)
    unsafe_copyto!(typed_pointer(dst, T), typed_pointer(src, T), N)

    ref = Array{T}(undef, N)
    unsafe_copyto!(pointer(ref), typed_pointer(dst, T), N)

    @test data == ref

    if isa(src, Mem.Device) || isa(src, Mem.Unified)
        Mem.set!(typed_pointer(src, T), zero(T), N)
    end

    Mem.free(src)
    Mem.free(dst)
end

# asynchronous operations
let
    src = Mem.alloc(Mem.Device, nb)

    @test_throws ArgumentError unsafe_copyto!(typed_pointer(src, T), pointer(data), N; async=true)
    unsafe_copyto!(typed_pointer(src, T), pointer(data), N; async=true, stream=CuDefaultStream())

    Mem.set!(typed_pointer(src, T), zero(T), N; async=true, stream=CuDefaultStream())

    Mem.free(src)
end

# pinned memory
let
    # can only get GPU pointer if the pinned buffer is mapped
    src = Mem.alloc(Mem.Host, nb)
    @test_throws ArgumentError convert(CuPtr{T}, src)
    Mem.free(src)

    # create a pinned and mapped buffer
    src = Mem.alloc(Mem.Host, nb, Mem.HOSTALLOC_DEVICEMAP)

    # get the CPU address and copy some data
    cpu_ptr = convert(Ptr{T}, src)
    unsafe_copyto!(cpu_ptr, pointer(data), N)

    # get the GPU address and construct a fake device buffer
    gpu_ptr = convert(CuPtr{Cvoid}, src)
    gpu_obj = Mem.alloc(Mem.Device, nb)
    dst = similar(gpu_obj, gpu_ptr)
    Mem.free(gpu_obj)

    # copy data back from the GPU and compare
    ref = Array{T}(undef, N)
    unsafe_copyto!(pointer(ref), typed_pointer(dst, T), N)
    @test ref == data

    Mem.free(src)
    # NOTE: don't free dst, it's just a mapped pointer
end

# pinned memory with existing memory
if attribute(dev, CUDAdrv.DEVICE_ATTRIBUTE_HOST_REGISTER_SUPPORTED) != 0
    # can only get GPU pointer if the pinned buffer is mapped
    src = Mem.register(Mem.Host, pointer(data), nb)
    @test_throws ArgumentError convert(CuPtr{T}, src)
    Mem.unregister(src)

    # register a pinned and mapped buffer
    src = Mem.register(Mem.Host, pointer(data), nb, Mem.HOSTREGISTER_DEVICEMAP)

    # get the GPU address and construct a fake device buffer
    gpu_ptr = convert(CuPtr{Cvoid}, src)
    gpu_obj = Mem.alloc(Mem.Device, nb)
    dst = similar(gpu_obj, gpu_ptr)
    Mem.free(gpu_obj)

    # copy data back from the GPU and compare
    ref = Array{T}(undef, N)
    unsafe_copyto!(pointer(ref), typed_pointer(dst, T), N)
    @test ref == data

    Mem.unregister(src)
    # NOTE: don't unregister dst, it's just a mapped pointer
end

# unified memory
let
    src = Mem.alloc(Mem.Unified, nb)

    #Mem.prefetch(src, nb; device=CUDAdrv.DEVICE_CPU)
    Mem.advise(src, Mem.ADVISE_SET_READ_MOSTLY)

    # get the CPU address and copy some data
    cpu_ptr = convert(Ptr{T}, src)
    unsafe_copyto!(cpu_ptr, pointer(data), N)

    # get the GPU address and construct a fake device buffer
    gpu_ptr = convert(CuPtr{Cvoid}, src)
    gpu_obj = Mem.alloc(Mem.Device, nb)
    dst = similar(gpu_obj, gpu_ptr)
    Mem.free(gpu_obj)

    # copy data back from the GPU and compare
    ref = Array{T}(undef, N)
    unsafe_copyto!(pointer(ref), typed_pointer(dst, T), N)
    @test ref == data

    Mem.free(src)
end

end


text:
export push, pop, shift, unshift, insert, deleteat

@generated function push(vec::StaticVector, x)
    newtype = similar_type(vec, (length(vec) + 1 ,))
    exprs = vcat([:(vec[$i]) for i = 1:length(vec)], :x)
    return quote
        $(Expr(:meta, :inline))
        @inbounds return $(Expr(:call, newtype, Expr(:tuple, exprs...)))
    end
end

@generated function unshift(vec::StaticVector, x)
    newtype = similar_type(vec, (length(vec) + 1 ,))
    exprs = vcat(:x, [:(vec[$i]) for i = 1:length(vec)])
    return quote
        $(Expr(:meta, :inline))
        @inbounds return $(Expr(:call, newtype, Expr(:tuple, exprs...)))
    end
end

# unfortunately index can't be used directly by a @generated function... this code is pretty slow
@generated function insert(vec::StaticVector, index, x)
    newtype = similar_type(vec, (length(vec) + 1 ,))
    expr = :((Tuple(vec[1:index-1])..., x, Tuple(vec[index:end])...))
    return quote
        $(Expr(:meta, :inline))
        @inbounds return $(Expr(:call, newtype, expr))
    end
end

@generated function pop(vec::StaticVector)
    newtype = similar_type(vec, (length(vec) - 1 ,))
    exprs = [:(vec[$i]) for i = 1:length(vec)-1]
    return quote
        $(Expr(:meta, :inline))
        @inbounds return $(Expr(:call, newtype, Expr(:tuple, exprs...)))
    end
end

@generated function shift(vec::StaticVector)
    newtype = similar_type(vec, (length(vec) - 1 ,))
    exprs = [:(vec[$i]) for i = 2:length(vec)]
    return quote
        $(Expr(:meta, :inline))
        @inbounds return $(Expr(:call, newtype, Expr(:tuple, exprs...)))
    end
end

# unfortunately index can't be used directly by a @generated function... this code is pretty slow
@generated function deleteat(vec::StaticVector, index)
    newtype = similar_type(vec, (length(vec) - 1 ,))
    expr = :((Tuple(vec[1:index-1])..., Tuple(vec[index+1:end])...))
    return quote
        $(Expr(:meta, :inline))
        @inbounds return $(Expr(:call, newtype, expr))
    end
end

# TODO consider prepend, append (can use vcat, but eltype might change), and
# maybe splice (a bit hard to get statically sized)


text:
###########################################################################
#
# This file contains following functions
#   glmxls - GLM regression models
#   bivariatexls - for two-way tables for discrete or continuous variables
#   univariatexls - for continuous variables
#   dfxls - for exporting a dataframe
#
###########################################################################

"""
"""
function hltest(glmout,q = 10)
    d = DataFrame(yhat = predict(glmout),y = glmout.model.rr.y)
    d[:group] = xtile(d[:yhat],nq = q)
    df3 = by(d,:group) do subdf
        n = size(subdf,1)
        o1 = sum(subdf[:y])
        e1 = sum(subdf[:yhat])
        DataFrame(
            observed1 = o1,
            observed0 = n - o1,
            expected1 = e1,
            expected0 = n - e1
        )
    end

    hlstat = sum((df3[:observed1] .- df3[:expected1]).^2 ./ df3[:expected1]
        .+ (df3[:observed0] .- df3[:expected0]).^2 ./ df3[:expected0])

    dof = q - 2

    pval = Distributions.ccdf(Distributions.Chisq(dof),hlstat)

    return (hlstat, dof, pval)
end

function nagelkerke(glmout)
    return StatsBase.r2(glmout,:Nagelkerke)
    # L = loglikelihood(glmout)
    # L0 = loglikelihood(glm(@eval(@formula($(glmout.mf.terms.eterms[1]) ~ 1)),glmout.mf.df,Bernoulli(),LogitLink()))
    # pow = 2/nobs(glmout)
    # return (1 - exp(-pow*(L - L0))) / (1 - exp(pow*L0))
end

function macfadden(glmout)
    return StatsBase.r2(glmout,:McFadden)
    # L = loglikelihood(glmout)
    # L0 = loglikelihood(glm(@eval(@formula($(glmout.mf.terms.eterms[1]) ~ 1)),glmout.mf.df,Bernoulli(),LogitLink()))
    # return 1 - L/L0
end

"""
    glmxls(glmout::DataFrames.DataFrameRegressionModel, workbook::PyObject, worksheet::AbstractString; labels::Union{Nothing,Label}=nothing,eform=false,ci=true, row = 0, col =0)

Outputs a GLM regression table to an excel spreadsheet.
To use this function, `PyCall` is required with a working version python and
a python package called `xlsxwriter` installed. If a label is found for a variable
or a value of a variable in a `Label`, the label will be output. Options are:

- `glmout`: returned value from a GLM regression model
- `workbook`: a returned value from xlsxwriter.Workbook() function (see an example below)
- `worksheet`: a string for the worksheet name
- `labels`: an option to specify a `label` dictionary (see an example below)
- `eform`: use `eform = true` to get exponentiated estimates, standard errors, or 95% confidence intervals
- `ci`: use `ci = true` (default) to get 95% confidence intervals. `ci = false` will produce standard errors and Z values instead.
- `row`: specify the row of the workbook to start the output table (default = 0 (for row 1))
- `col`: specify the column of the workbook to start the output table (default = 0 (for column A))

# Example 1
This example is useful when one wants to append a worksheet to an existing workbook.
It is responsibility of the user to open a workbook before the function call and close it
to actually create the physical file by close the workbook.

```
julia> using PyCall

julia> @pyimport xlsxwriter

julia> wb = xlsxwriter.Workbook("test_workbook.xlsx")
PyObject <xlsxwriter.workbook.Workbook object at 0x000000002A628E80>

julia> glmxls(ols1,wb,"OLS1",labels = label)

julia> bivairatexls(df,:incomecat,[:age,:race,:male,:bmicat],wb,"Bivariate",labels = label)

Julia> wb.close()
```

# Example 2
Alternatively, one can create a spreadsheet file directly. `PyCall` or `@pyimport`
does not need to be called before the function.

```
julia> glmxls(ols1,"test_workbook.xlsx","OLS1",labels = label)
```

"""
function glmxls(glmout,wbook::PyObject,wsheet::AbstractString;
    labels::Union{Nothing,Label} = nothing,
    eform::Bool = false, ci = true, row = 0, col = 0)

    if (typeof(glmout) <: StatsModels.RegressionModel) == false
        error("This is not a regression model output.")
    end

    if isa(glmout.model,GeneralizedLinearModel)
        distrib = glmout.model.rr.d
        linkfun = Link(glmout.model.rr)
    else
        distrib = nothing
        linkfun = nothing
    end

    # create a worksheet
    t = wbook.add_worksheet(wsheet)

    # attach formats to the workbook
    formats = ExcelTables.attach_formats(wbook)

    # starting location in the worksheet
    r = row
    c = col

    # set column widths
    t.set_column(c,c,40)
    t.set_column(c+1,c+4,9)

    # headings ---------------------------------------------------------
    # if ci == true, Estimate (95% CI) P-Value
    # if eform == true, Estimate is OR for logit, IRR for poisson
    otype = "Estimate"
    if eform == true
        otype = Stella.coeflab(distrib,linkfun)
    end

    t.write_string(r,c,"Variable",formats[:heading])
    if ci == true
        t.merge_range(r,c+1,r,c+3,string(otype," (95% CI)"),formats[:heading])
        t.write_string(r,c+4,"P-Value",formats[:heading])
    else
        t.write_string(r,c+1,otype,formats[:heading])
        t.write_string(r,c+2,"SE",formats[:heading])
        t.write_string(r,c+3,"Z Value",formats[:heading])
        t.write_string(r,c+4,"P-Value",formats[:heading])
    end

    #------------------------------------------------------------------
    r += 1
    c = col

    # go through each variable and construct variable name and value label arrays
    tdata = coeftable(glmout)
    nrows = length(tdata.rownms)
    varname = Vector{String}(undef,nrows)
    vals = Vector{String}(undef,nrows)
    nlev = zeros(Int,nrows)
    #nord = zeros(Int,nrow)
    for i = 1:nrows
    	# variable name
        # parse varname to separate variable name from value
        if occursin(":",tdata.rownms[i])
            (varname[i],vals[i]) = split(tdata.rownms[i],": ")
        elseif occursin(" - ",tdata.rownms[i])
            (varname[i],vals[i]) = split(tdata.rownms[i]," - ")
        else
            varname[i] = tdata.rownms[i]
            vals[i] = ""
        end

        # use labels if exist
        if labels != nothing
            if vals[i] != ""
                valn = vals[i] == "true" ? 1 : parse(Int,vals[i])
                vals[i] = vallab(labels,Symbol(varname[i]),valn)
            end
            varname[i] = varlab(labels,Symbol(varname[i]))
        end
    end
    for i = 1:nrows
        nlev[i] = countlev(varname[i],varname)
    end

    # write table
    tconfint = confint(glmout)
    lastvarname = ""

    for i = 1:nrows
        if otype == "OR" && varname[i] == "(Intercept)"
            continue
        end
        if varname[i] != lastvarname
            # output cell boundaries only and go to the next line
            if nlev[i] > 1
                t.write_string(r,c,varname[i],formats[:heading_left])

                if ci == true
                    t.write(r,c+1,"",formats[:empty_right])
                    t.write(r,c+2,"",formats[:empty_both])
                    t.write(r,c+3,"",formats[:empty_left])
                    t.write(r,c+4,"",formats[:p_fmt])
                else
                    t.write(r,c+1,"",formats[:empty_border])
                    t.write(r,c+2,"",formats[:empty_border])
                    t.write(r,c+3,"",formats[:empty_border])
                    t.write(r,c+4,"",formats[:p_fmt])
                end
                r += 1
                t.write_string(r,c,vals[i],formats[:varname_1indent])

            else
                if vals[i] != "" && vals[i] != "Yes"
                    t.write_string(r,c,string(varname[i]," - ",vals[i]),formats[:heading_left])
                else
                    t.write_string(r,c,varname[i],formats[:heading_left])
                end
            end
        else
            t.write_string(r,c,vals[i],formats[:varname_1indent])
        end

    	# estimates
        if eform == true
    	    t.write(r,c+1,exp(tdata.cols[1][i]),formats[:or_fmt])
        else
            t.write(r,c+1,tdata.cols[1][i],formats[:or_fmt])
        end

        if ci == true

            if eform == true
            	# 95% CI Lower
            	t.write(r,c+2,exp(tconfint[i,1]),formats[:cilb_fmt])

            	# 95% CI Upper
            	t.write(r,c+3,exp(tconfint[i,2]),formats[:ciub_fmt])
            else
                # 95% CI Lower
            	t.write(r,c+2,tconfint[i,1],formats[:cilb_fmt])

            	# 95% CI Upper
            	t.write(r,c+3,tconfint[i,2],formats[:ciub_fmt])
            end
        else
            # SE
            if eform == true
        	    t.write(r,c+2,exp(tdata.cols[1][i])*tdata.cols[2][i],formats[:or_fmt])
            else
                t.write(r,c+2,tdata.cols[1][i],formats[:or_fmt])
            end

            # Z value
            t.write(r,c+3,tdata.cols[3][i],formats[:or_fmt])

        end

        # P-Value
        t.write(r,c+4,tdata.cols[4][i] < 0.001 ? "< 0.001" : tdata.cols[4][i] ,formats[:p_fmt])

        lastvarname = varname[i]

        # update row
        r += 1
    end

    # Write model characteristics and goodness of fit statistics
    c = col

    # N
    t.write(r,c,"N",formats[:model_name])
    t.merge_range(r,c+1,r,c+4,nobs(glmout),formats[:n_fmt_center])

    # degress of freedom
    # r += 1
    # t.write(r,c,"DF",formats[:model_name])
    # t.merge_range(r,c+1,r,c+4,dof(glmout),formats[:n_fmt_center])

    # R² or pseudo R²
    r += 1
    if isa(linkfun,LogitLink)
        t.write(r,c,"Pseudo R² (MacFadden)",formats[:model_name])
        t.merge_range(r,c+1,r,c+4,macfadden(glmout),formats[:p_fmt_center])
        t.write(r+1,c,"Pseudo R² (Nagelkerke)",formats[:model_name])
        t.merge_range(r+1,c+1,r+1,c+4,nagelkerke(glmout),formats[:p_fmt_center])

        # -2 log-likelihood
        t.write(r+2,c,"-2 Log-Likelihood",formats[:model_name])
        t.merge_range(r+2,c+1,r+2,c+4,deviance(glmout),formats[:p_fmt_center])

        # Hosmer-Lemeshow GOF test
        t.write(r+3,c,"Hosmer-Lemeshow Chisq Test (df), p-value",formats[:model_name])
        hl = hltest(glmout)
        t.merge_range(r+3,c+1,r+3,c+4,string(round(hl[1],digits=4)," (",hl[2],"); p = ",round(hl[3],digits=4)),formats[:p_fmt_center])

        # ROC (c-statistic)
        t.write(r+4,c,"Area under the ROC Curve",formats[:model_name])
        roc = auc(roc(glmout.model.rr.y,predict(glmout)))
        n1 = sum(glmout.model.rr.y) # number of positive responses
        n2 = nobs(glmout) - n1
        q1 = roc / (2 - roc)
        q2 = (2*roc^2) / (1 + roc)
        rocse = sqrt((roc*(1-roc) + (n1-1)*(q1 - roc^2) + (n2 - 1)*(q2 - roc^2)) / (n1*n2))
        t.merge_range(r+4,c+1,r+4,c+4,string(round(roc,digits=4)," (95% CI, ", round(roc - 1.96*rocse,digits=4), " - ", round(roc + 1.96*rocse,digits=4),")"),formats[:p_fmt_center])

        r += 5
    else
        t.write(r,c,"R²",formats[:model_name])
        t.merge_range(r,c+1,r,c+4,r2(glmout),formats[:p_fmt_center])
        t.write(r+1,c,"Adjusted R²",formats[:model_name])
        t.merge_range(r+1,c+1,r+1,c+4,adjr2(glmout),formats[:p_fmt_center])

        r += 2
    end

    # AIC & BIC
    t.write(r,c,"AIC",formats[:model_name])
    t.merge_range(r,c+1,r,c+4,aic(glmout),formats[:p_fmt_center])

    r += 1
    t.write(r,c,"BIC",formats[:model_name])
    t.merge_range(r,c+1,r,c+4,bic(glmout),formats[:p_fmt_center])
end
function glmxls(glmout,
    wbook::AbstractString,
    wsheet::AbstractString;
    labels::Union{Nothing,Label} = nothing,
    eform::Bool = false,
    ci = true,
    row = 0,
    col = 0)

    xlsxwriter = pyimport("xlsxwriter")

    glmxls(glmout,xlsxwriter.Workbook(wbook),wsheet,labels=labels,eform=eform,ci=ci,row=row,col=col)
end

"""
    bivariatexls(df::DataFrame,colvar::Symbol,rowvars::Vector{Symbol},workbook::PyObject,worksheet::AbstractString; labels::Union{Nothing,Label}=nothing,row=0,col=0)

 Creates bivariate statistics and appends it in a nice tabular format to an existing workbook.
 To use this function, `PyCall` is required with a working version python and
 a python package called `xlsxwriter` installed.  If a label is found for a variable
 or a value of a variable in a `Label`, the label will be output. Options are:

- `df`: a DataFrame
- `colvar`: a categorical variable whose values will be displayed on the columns
- `rowvars`: a Vector of Symbols for variables to be displayed on the rows. Both continuous and categorical variables are allowed.
    For continuous variables, mean and standard deviations will be output and a p-value will be based on an ANOVA test. For categorical variables,
    a r x c table with cell counts and row percentages will be output with a p-value based on a chi-square test.
- `workbook`: a returned value from xlsxwriter.Workbook() function (see an example below)
- `worksheet`: a string for the worksheet name
- `labels`: an option to specify a `Label` object (see an example below)
- `row`: specify the row of the workbook to start the output table (default = 0 (for row 1))
- `col`: specify the column of the workbook to start the output table (default = 0 (for column A))
- `column_percent`: set this to `false` if you want row percentages in the output table (default = true)

### Example 1
This example is useful when one wants to append a worksheet to an existing workbook.
It is responsibility of the user to open a workbook before the function call and close it
to actually create the physical file by close the workbook.

```
julia> using PyCall

julia> xlsxwriter = pyimport("xlsxwriter")

julia> wb = xlsxwriter.Workbook("test_workbook.xlsx")
PyObject <xlsxwriter.workbook.Workbook object at 0x000000002A628E80>

julia> glmxls(ols1,wb,"OLS1",labels = label)

julia> bivairatexls(df,:incomecat,[:age,:race,:male,:bmicat],wb,"Bivariate",labels = label)

Julia> wb.close()
```

### Example 2
Alternatively, one can create a spreadsheet file directly. `PyCall` or `@pyimport`
does not need to be called before the function.

```
julia> bivariatexls(df,:incomecat,[:age,:race,:male,:bmicat],"test_workbook.xlsx","Bivariate",labels = label)
```
"""
function bivariatexls(df::DataFrame,
    colvar::Symbol,
    rowvars::Vector{Symbol},
    wbook::PyObject,
    wsheet::AbstractString;
    wt::Union{Nothing,Symbol} = nothing,
    labels::Union{Nothing,Label} = nothing,
    row::Int = 0,
    col::Int = 0,
    column_percent::Bool = true,
    verbose::Bool = false)

    # check data

    # colvar has to be a CategoricalArray and must have 2 or more categories
    if isa(df[!,colvar], CategoricalArray) == false || length(levels(df[!,colvar])) < 2
        error("`",colvar,"` is not a CategoricalArray or does not have two or more levels")
    end

    # create a worksheet
    t = wbook.add_worksheet(wsheet)

    # attach formats to the workbook
    formats = ExcelTables.attach_formats(wbook)

    # starting row and column
    r = row
    c = col

    # drop NAs in colvar
    df2 = df[completecases(df[!,[colvar]]),:]

    # number of columns
    # column values
    if wt == nothing
        collev = freqtable(df2,colvar,skipmissing=true)
    else
        collev = freqtable(df2,colvar,skipmissing=true,weights=df2[wt])
    end
    nlev = length(collev.array)
    colnms = names(collev,1)
    coltot = sum(collev.array,dims=1)

    # set column widths
    t.set_column(c,c,40)
    t.set_column(c+1,c+(nlev+1)*2+1,9)

    # create heading
    # column variable name
    # It uses three rows
    t.merge_range(r,c,r+2,c,"Variable",formats[:heading])

    # 1st row = variable name
    colvname = String(colvar)
    if labels != nothing
        colvname = varlab(labels,colvar)
    end
    t.merge_range(r,c+1,r,c+(nlev+1)*2+1,colvname,formats[:heading])

    # 2nd and 3rd rows
    r += 1

    t.merge_range(r,1,r,2,"All",formats[:heading])
    t.write_string(r+1,1,"N",formats[:n_fmt_right])
    t.write_string(r+1,2,"(%)",formats[:pct_fmt_parens])

    c = 3
    for i = 1:nlev

        # value label
        vals = string(colnms[i])
        if labels != nothing
            vals = vallab(labels,colvar,colnms[i])
        end

        t.merge_range(r,c+(i-1)*2,r,c+(i-1)*2+1,vals,formats[:heading])
        t.write_string(r+1,c+(i-1)*2,"N",formats[:n_fmt_right])
        t.write_string(r+1,c+(i-1)*2+1,"(%)",formats[:pct_fmt_parens])
    end

    # P-value
    t.merge_range(r,c+nlev*2,r+1,c+nlev*2,"P-Value",formats[:heading])

    # total
    c = 0
    r += 2
    t.write_string(r,c,"All, n (Row %)",formats[:model_name])
    if wt == nothing
        x = freqtable(df2,colvar,skipmissing=true)
    else
        x = freqtable(df2,colvar,skipmissing=true,weights=df2[wt])
    end
    tot = sum(x)
    t.write(r,c+1,tot,formats[:n_fmt_right])
    t.write(r,c+2,1.0,formats[:pct_fmt_parens])
    for i = 1:nlev
        t.write(r,c+i*2+1,x.array[i],formats[:n_fmt_right])
        t.write(r,c+i*2+2,x.array[i]/tot,formats[:pct_fmt_parens])
    end
    t.write(r,c+(nlev+1)*2+1,"",formats[:empty_border])

    # covariates
    c = 0
    r += 1
    for varname in rowvars

        if verbose == true
            println("Processing ",varname)
        end

        # sum of rowvars must be non-zero
        # if sum(length(skipmissing(df[varname]))) == 0
        #     warn("`",varname,"` in `rowvars` is empty.")
        #     continue
        # end

        # print the variable name
        vars = string(varname)
        if labels != nothing
            vars = varlab(labels,varname)
        end

        # determine if varname is categorical or continuous
        if isa(df2[!,varname], CategoricalArray) || eltype(df2[!,varname]) == String

            # categorial
            df3=df2[completecases(df2[:,[varname]]),[varname,colvar]]
            if wt == nothing
                x = freqtable(df3,varname,colvar,skipmissing=true)
            else
                x = freqtable(df3,varname,colvar,skipmissing=true,weights=df3[wt])
            end
            rowval = names(x,1)
            rowtot = sum(x.array,dims=2)
            coltot = sum(x.array,dims=1)

            # variable name
            # if there only two levels and one of the values is 1 or true
            # and the other values is 0 or false,
            # just output the frequency and percentage of the 1/true row

            # variable name
            t.write_string(r,c,vars,formats[:model_name])

            # two levels with [0,1] or [false,true]
            if length(rowval) <= 2 && rowval in ([1],[true],["Yes"],[0,1],[false,true],["No","Yes"])

                nrow = length(rowval)
                # row total
                t.write(r,c+1,rowtot[nrow],formats[:n_fmt_right])
                t.write(r,c+2,rowtot[nrow]/tot,formats[:pct_fmt_parens])

                for j = 1:nlev
                    t.write(r,c+j*2+1,x.array[nrow,j],formats[:n_fmt_right])
                    if column_percent
                        t.write(r,c+j*2+2, coltot[j] > 0 ? x.array[nrow,j]/coltot[j] : "",formats[:pct_fmt_parens])
                    elseif rowtot[2] > 0
                        t.write(r,c+j*2+2, rowtot[nrow] > 0 ? x.array[nrow,j]/rowtot[nrow] : "",formats[:pct_fmt_parens])
                    end
                end
                pval = pvalue(ChisqTest(x.array))
                if isnan(pval) || isinf(pval)
                    pval = ""
                elseif pval < 0.001
                    pval = "< 0.001"
                end
                t.write(r,c+(nlev+1)*2+1, pval,formats[:p_fmt])
                r += 1
            else
                for i = 1:nlev+1
                    t.write_string(r,c+(i-1)*2+1,"",formats[:empty_right])
                    t.write_string(r,c+(i-1)*2+2,"",formats[:empty_left])
                end
                t.write_string(r,c+(nlev+1)*2+1,"",formats[:empty_border])

                r += 1
                for i = 1:length(rowval)
                    # row value
                    vals = string(rowval[i])

                    if labels != nothing
                        vals = vallab(labels,varname,rowval[i])
                    end
                    t.write_string(r,c,vals,formats[:varname_1indent])

                    # row total
                    t.write(r,c+1,rowtot[i],formats[:n_fmt_right])
                    t.write(r,c+2,rowtot[i]/tot,formats[:pct_fmt_parens])

                    for j = 1:nlev
                        t.write(r,c+j*2+1,x.array[i,j],formats[:n_fmt_right])
                        if column_percent
                            t.write(r,c+j*2+2,coltot[j] > 0 ? x.array[i,j]/coltot[j] : "",formats[:pct_fmt_parens])
                        else
                            t.write(r,c+j*2+2,rowtot[i] > 0 ? x.array[i,j]/rowtot[i] : "",formats[:pct_fmt_parens])
                        end
                    end
                    # p-value - output only once
                    pval = pvalue(ChisqTest(x.array))

                    if isnan(pval) || isinf(pval)
                        pval = ""
                    elseif pval < 0.001
                        pval = "< 0.001"
                    end
                    if length(rowval) == 1
                        t.write(r,c+(nlev+1)*2,pval,formats[:p_fmt])
                    elseif i == 1
                        t.merge_range(r,c+(nlev+1)*2+1,r+length(rowval)-1,c+(nlev+1)*2+1,pval,formats[:p_fmt])
                    end
                    r += 1
                end
            end
        else
            # continuous variable
            df3=df2[completecases(df2[!,[varname]]),[varname,colvar]]
            y = tabstat(df3,varname,colvar) #,wt=df3[wt])

            # variable name
            t.write_string(r,c,string(vars,", mean (SD)"),formats[:model_name])

            # All
            t.write(r,c+1,mean(df3[!,varname]),formats[:f_fmt_right])
            t.write(r,c+2,std(df3[!,varname]),formats[:f_fmt_left_parens])

            # colvar levels
            for i = 1:nlev
                if i <= size(y,1) && y[i,:N] > 1
                    t.write(r,c+i*2+1,y[i,:mean],formats[:f_fmt_right])
                    t.write(r,c+i*2+2,y[i,:sd],formats[:f_fmt_left_parens])
                else
                    t.write(r,c+i*2+1,"",formats[:f_fmt_right])
                    t.write(r,c+i*2+2,"",formats[:f_fmt_left_parens])
                end
            end
            if size(y,1) > 1
                pval = anova(df3,varname,colvar).p[1]
                if isnan(pval) || isinf(pval)
                    pval = ""
                elseif pval < 0.001
                    pval = "< 0.001"
                end
                t.write(r,c+(nlev+1)*2+1,pval,formats[:p_fmt])
            else
                t.write(r,c+(nlev+1)*2+1,"",formats[:p_fmt])
            end

            r += 1
        end
    end
end
function bivariatexls(df::DataFrame,
    colvar::Symbol,
    rowvars::Vector{Symbol},
    wbook::AbstractString,
    wsheet::AbstractString;
    labels::Union{Nothing,Label} = nothing,
    row::Int = 0,
    col::Int = 0)

    xlsxwriter = pyimport("xlsxwriter")

    wb = xlsxwriter.Workbook(wbook)

    bivariatexls(df,colvar,rowvars,wb,wsheet,labels=labels,row=row,col=col)

    wb.close()
end

"""
    univariatexls(df::DataFrame,contvars::Vector{Symbol},workbook::PyObject,worksheet::AbstractString; labels::Union{Nothing,Label}=nothing,row=0,col=0)

Creates univariate statistics for a vector of continuous variable and
appends it to an existing workbook.
To use this function, `PyCall` is required with a working version python and
a python package called `xlsxwriter` installed.  If a label is found for a variable
in a `Label` object, the label will be output. Options are:

- `df`: a DataFrame
- `contvars`: a vector of continuous variables
- `workbook`: a returned value from xlsxwriter.Workbook() function (see an example below)
- `worksheet`: a string for the worksheet name
- `labels`: an option to specify a `label` dictionary (see an example below)
- `row`: specify the row of the workbook to start the output table (default = 0 (for row 1))
- `col`: specify the column of the workbook to start the output table (default = 0 (for column A))

# Example 1
This example is useful when one wants to append a worksheet to an existing workbook.
It is responsibility of the user to open a workbook before the function call and close it
to actually create the physical file by close the workbook.

```
julia> using PyCall

julia> @pyimport xlsxwriter

julia> wb = xlsxwriter.Workbook("test_workbook.xlsx")
PyObject <xlsxwriter.workbook.Workbook object at 0x000000002A628E80>

julia> glmxls(ols1,wb,"OLS1",labels = label)

julia> bivairatexls(df,:incomecat,[:age,:race,:male,:bmicat],wb,"Bivariate",labels = label)

julia> univariatexls(df,[:age,:income_amt,:bmi],wb,"Univariate",labels = label)

Julia> wb.close()
```

# Example 2
Alternatively, one can create a spreadsheet file directly. `PyCall` or `@pyimport`
does not need to be called before the function.

```jldoctest
julia> univariatexls(df,[:age,:income_amt,:bmi],"test_workbook.xlsx","Bivariate",labels = label)
```

"""
function univariatexls(df::DataFrame,
    contvars::Vector{Symbol},
    wbook::PyObject,
    wsheet::AbstractString;
    wt::Union{Nothing,Symbol} = nothing,
    labels::Union{Nothing,Label}=nothing,
    row = 0,
    col = 0)

    # create a worksheet
    t = wbook.add_worksheet(wsheet)

    # attach formats to the workbook
    formats = ExcelTables.attach_formats(wbook)

    # starting row and column
    r = row
    c = col

    # column width
    t.set_column(0,0,20)
    t.set_column(1,length(contvars),12)

    # output the row names
    t.write_string(r,c,"Statistic",formats[:heading])
    t.write_string(r+1,c,"N Total",formats[:heading_left])
    t.write_string(r+2,c,"N Miss",formats[:heading_left])
    t.write_string(r+3,c,"N Used",formats[:heading_left])
    t.write_string(r+4,c,"Sum",formats[:heading_left])
    t.write_string(r+5,c,"Mean",formats[:heading_left])
    t.write_string(r+6,c,"SD",formats[:heading_left])
    t.write_string(r+7,c,"Variance",formats[:heading_left])
    t.write_string(r+8,c,"Minimum",formats[:heading_left])
    t.write_string(r+9,c,"P25",formats[:heading_left])
    t.write_string(r+10,c,"Median",formats[:heading_left])
    t.write_string(r+11,c,"P75",formats[:heading_left])
    t.write_string(r+12,c,"Maximum",formats[:heading_left])
    t.write_string(r+13,c,"Skewness",formats[:heading_left])
    t.write_string(r+14,c,"Kurtosis",formats[:heading_left])
    t.write_string(r+15,c,"Smallest",formats[:heading_left])
    t.write_string(r+16,c,"",formats[:heading_left])
    t.write_string(r+17,c,"",formats[:heading_left])
    t.write_string(r+18,c,"",formats[:heading_left])
    t.write_string(r+19,c,"",formats[:heading_left])
    t.write_string(r+20,c,"Largest",formats[:heading_left])
    t.write_string(r+21,c,"",formats[:heading_left])
    t.write_string(r+22,c,"",formats[:heading_left])
    t.write_string(r+23,c,"",formats[:heading_left])
    t.write_string(r+24,c,"",formats[:heading_left])

    col = 1
    for vsym in contvars

        # if there is a label dictionary, pick up the variable label
        varstr = string(vsym)
        if labels != nothing
            varstr = varlab(labels,vsym)
        end

        t.write_string(0,col,varstr,formats[:heading])
        u = univariate(df,vsym) #,wt=df[wt])
        for j = 1:14
            if j<4
                fmttype = :n_fmt
            else
                fmttype = :p_fmt
            end
            if isnan(u[j,:Value]) || isinf(u[j,:Value])
                t.write(j,col,"",formats[fmttype])
            else
                t.write(j,col,u[j,:Value],formats[fmttype])
            end
        end
        smallest=Stella.smallest(df[vsym])
        if eltype(df) <: Integer
            fmttype = :n_fmt
        else
            fmttype = :p_fmt
        end
        for j = 1:5
            t.write(j+14,col,smallest[j],formats[fmttype])
        end
        largest=Stella.largest(df[vsym])
        for j = 1:5
            t.write(j+19,col,largest[j],formats[fmttype])
        end
        col += 1
    end
end
function univariatexls(df::DataFrame,contvars::Vector{Symbol},wbook::AbstractString,wsheet::AbstractString;
    wt::Union{Nothing,Symbol} = nothing,labels::Union{Nothing,Label}=nothing, row = 0, col = 0)

    xlsxwriter=pyimport("xlsxwriter")

    wb = wlsxwriter.Workbook(wbook)

    univariatexls(df,contvars,wb,wsheet,wt=wt,labels=labels,row=row,col=col)

    wb.close()
end

"""
    dfxls(df::DataFrame, workbook::PyObject, worksheet::AbstractString; nrows = 500, start = 1, row=0, col=0)

 To use this function, `PyCall` is required with a working version python and
 a python package called `xlsxwriter` installed. Options are:

- `df`: a DataFrame
- `workbook`: a returned value from xlsxwriter.Workbook() function (see an example below)
- `worksheet`: a string for the worksheet name (default: "Data1")
- `start`: specify the row number from which data will be output (default: 1)
- `nrows`: specify the number of rows to output (default: 500). If nows = 0, the entire dataframe will be output.
- `row`: specify the row of the workbook to start the output table (default = 0 (for row 1))
- `col`: specify the column of the workbook to start the output table (default = 0 (for column A))

# Example 1
This example is useful when one wants to append a worksheet to an existing workbook.
It is responsibility of the user to open a workbook before the function call and close it
to actually create the physical file by close the workbook.

```
julia> using PyCall

julia> xlsxwriter = pyimport("xlsxwriter")

julia> wb = xlsxwriter.Workbook("test_workbook.xlsx")
PyObject <xlsxwriter.workbook.Workbook object at 0x000000002A628E80>

julia> glmxls(ols1,wb,"OLS1",labels = label)

julia> bivairatexls(df,:incomecat,[:age,:race,:male,:bmicat],wb,"Bivariate",labels = label)

julia> univariatexls(df,[:age,:income_amt,:bmi],wb,"Univariate",labels = label)

julia> dfxls(df,wb,"dataframe",nrows = 0)

Julia> wb.close()
```

# Example 2
Alternatively, one can create a spreadsheet file directly. `PyCall` or `@pyimport`
does not need to be called before the function.

```
julia> dfxls(df,"test_workbook.xlsx","df",nrows = 0)
```

"""
function dfxls(df::DataFrame,
    wbook::PyObject,
    worksheet::AbstractString;
    nrows::Int64 = 0, start::Int64 = 1, col::Int64 = 0, row::Int64 = 0)

    # create a worksheet
    t = wbook.add_worksheet(worksheet)

    # attach formats to the workbook
    formats = ExcelTables.attach_formats(wbook)

    # starting row and column
    c = col

    # Eltype
    typ = Vector{DataType}(undef,size(df,2))
    for i=1:size(df,2)
        if typeof(df[i]) <: CategoricalArray
            typ[i] = eltype(df[i].pool.index)
        else
            typ[i] = nonmissingtype(eltype(df[i]))
        end
    end
    varnames = names(df)

    # if nrows = 0, output the full data set
    if nrows == 0
        start = 1
        nrows = size(df,1)
    else
        nrows = (start+nrows-1) < size(df,1) ? nrows : (size(df,1)-start+1)
    end

    for i = 1:size(df,2)

        r = row
        t.set_column(c,c,10)
        t.write_string(r,c,string(varnames[i]),formats[:heading])
        r += 1

        for j in start:(start+nrows-1)
            # println("i = ",i,"; j = ",j,". value = ",df[j,i])

            if ismissing(df[j,i])
                t.write_string(r,c," ",formats[:n_fmt])
            elseif typ[i] <: AbstractString
                if df[j,i] == ""
                    t.write_string(r,c," ",formats[:text])
                else
                    t.write(r,c,df[j,i],formats[:text])
                end
            elseif typ[i] <: Number
                if isnan(df[j,i]) || isinf(df[j,i])
                    t.write_string(r,c," ",formats[:text])
                elseif typ[i] <: Integer
                    t.write(r,c,df[j,i],formats[:n_fmt])
                elseif typ[i] <: AbstractFloat
                    t.write(r,c,df[j,i],formats[:f_fmt])
                end
            elseif typ[i] <: Date
                t.write(r,c,Dates.value(df[j,i] - Date(1899,12,30)),formats[:f_date])
            elseif typ[i] <: DateTime
                t.write(r,c,(Dates.value(df[j,i] - DateTime(1899,12,30,0,0,0)))/86400000,formats[:f_datetime])
            elseif typ[i] == Symbol || typ[i] == DataType
                t.write(r,c,string(df[j,i]),formats[:text])
            else
                # skip
            end

            r += 1
        end
        c += 1
    end

end
function dfxls(df::DataFrame,
    wbook::AbstractString,
    worksheet::AbstractString;
    nrows::Int64 = 0, start::Int64 = 1, col::Int64 = 0, row::Int64 = 0)

    # import xlsxwriter
    XlsxWriter = pyimport("xlsxwriter")

    # create a workbook
    wb = XlsxWriter.Workbook(wbook)

    dfxls(df,wb,worksheet, nrows = nrows, start = start, col = col, row = row)

    wb.close()
end

function newfilename(filen::AbstractString)
    while (isfile(filen))
        # separate the name into three parts, basename (number).extension
        (basename,ext) = splitext(filen)
        m = match(r"(.*)(\(.*\))$",basename)
        if m == nothing
            filen = string(filen," (1)")
        else
            m2 = parse(Int,replace(m[2],r"\((.*)\)",s"\1"))
            filen = string(m[1]," (",m2+1,")")
        end
    end
    return filen
end

function countlev(str::AbstractString,sarray::Vector{String})
    k = 0
    for i = 1:length(sarray)
        if str == sarray[i]
            k += 1
        end
    end
    return k
end


text:
using StaticArrays
using MeshCore: nshapes, VecAttrib
using MeshCore: ir_skeleton, ir_transpose, ir_bbyfacets
using MeshCore: attribute
using MeshCore: P1, L2, T3, ShapeColl, manifdim, nfacets, facetdesc, nshapes
using MeshCore: IncRel

"""
    orient_surface_mesh(t2v)

Orient surface mesh.

Return
- `orientedt2v`: oriented incidence relation
- `orientable`: Can the surface mesh be oriented? Bool flag.
    This would be false for a Mobius strip, for instance.

This is a purely topological operation. Creases in the surface are not
recognized as edges. The operation proceeds by flooding the surface across
manifold edges, stopping at sheet or non-manifold edges.

The incidence relation `orientedt2v` includes an attribute for each triangle
that provides the surface id to which the triangle belongs
(`orientedt2v.left.attributes["surfid"]`).
"""
function orient_surface_mesh(t2v)
    orientable = true
    all_oriented = false
    surfid = fill(0, nshapes(t2v.left))
    e2v = ir_skeleton(t2v)
    t2e = ir_bbyfacets(t2v, e2v)
    ot2ev = [SVector(t2e[i]) for i in 1:nshapes(t2e.left)]
    e2t = ir_transpose(t2e)
    while !all_oriented
        all_oriented = true
        for i in 1:nshapes(t2v.left)
            if surfid[i] == 0
                orientable = orientable && _orient_surface!(ot2ev, surfid, t2v, i, t2e, e2t);
                all_oriented = false;
                break;
            end
        end
    end
    elements = ShapeColl(T3, nshapes(t2v.left))
    vrts = ShapeColl(P1, nshapes(t2v.right))
    oc = [_vertfromedges(ot2ev[i], e2v) for i in 1:length(ot2ev)]
    orientedt2v = IncRel(elements, vrts, oc)
    orientedt2v.left.attributes["surfid"] = VecAttrib(surfid)
    orientedt2v.right.attributes["geom"] = attribute(t2v.right, "geom")
    return orientedt2v, orientable
end

function _vertfromedges(e, e2v)
    e1id = abs(e[1])
    e2id = abs(e[2])
    e3id = abs(e[3])
    _commonv(ea, eb) = (ea[1] == eb[1] ? ea[1] : (ea[2] == eb[1] ? ea[2] : (ea[2] == eb[2] ? ea[2] : ea[1])))
    if e[1] > 0
        return SVector(e2v[e1id][1], e2v[e1id][2], _commonv(e2v[e2id], e2v[e3id]))
    else
        return SVector(e2v[e1id][2], e2v[e1id][1], _commonv(e2v[e2id], e2v[e3id]))
    end
end

function _iscompatible(je, kes)
    for e in 1:3
        if je == -kes[e]
            return true
        end
    end
    return false
end

function _deduce_edge_order(je, kes)
    if _iscompatible(je, kes)
        return kes
    else
        return SVector(-kes[3], -kes[2], -kes[1])
    end
end

function _orient_surface!(ot2ev, surfid, t2v, i, t2e, e2t)
    orientable = true
    currsurfid = maximum(surfid) + 1
    ot2ev[i] = t2e[i] # First Triangle: Accept its default orientation
    surfid[i] = currsurfid
    stack = Int[]
    push!(stack, i)
    while (!isempty(stack))
        j = pop!(stack)
        for e in 1:3
            eid = abs(ot2ev[j][e])
            if length(e2t[eid]) == 2 # is this a manifold edge?
                k = e2t[eid][1] == j ? e2t[eid][2] : e2t[eid][1]
                if surfid[k] == 0
                    ot2ev[k] = _deduce_edge_order(ot2ev[j][e], ot2ev[k])
                    surfid[k] = currsurfid
                    push!(stack, k)
                else
                    orientable = orientable && _iscompatible(ot2ev[j][e], ot2ev[k])
                end
            end
        end
    end
    return orientable
end

function _normal(c, geom)
    n = StaticArrays.cross(geom[c[2]] - geom[c[1]], geom[c[3]] - geom[c[1]])
    n /= StaticArrays.norm(n)
    return n
end

"""
    _classify_mesh_edges(t2v, e2t, crease_ang)

Classify mesh edges as on/not on topological edge.

- `e2t`: the attribute `isonte`
"""
function _classify_mesh_edges(t2v, e2t, crease_ang)
    cada = cos(crease_ang);
    geom = attribute(t2v.right, "geom")
    surfid = attribute(t2v.left, "surfid")
    # e2v = ir_skeleton(t2v)
    # t2e = ir_bbyfacets(t2v, e2v)
    # ot2ev = [SVector(t2e[i]) for i in 1:nshapes(t2e.left)]
    # e2t = ir_transpose(t2e)
    isonte = fill(false, nshapes(e2t.left))
    for e in 1:nshapes(e2t.left)
        numt_at_e = length(e2t[e])
        @assert numt_at_e >= 1
        if numt_at_e == 2 # manifold edge
            # if the triangles belong to different surfaces, this is a TE
            if surfid[e2t[e][1]] != surfid[e2t[e][2]]
                println("different surfaces, $(surfid[e2t[e][1]]) vs $(surfid[e2t[e][2]])")
                isonte[e] = true
            else
                # check if it is a crease
                n1 = _normal(t2v[e2t[e][1]], geom)
                n2 = _normal(t2v[e2t[e][2]], geom)
                if StaticArrays.dot(n1, n2) < cada
                    isonte[e] = true
                end
            end
        elseif numt_at_e == 1 # sheet edge
            isonte[e] = true
        else # non manifold edge
            isonte[e] = true
        end
    end
    e2t.left.attributes["isonte"] = VecAttrib(isonte)
    return e2t
end

"""
    make_topo_faces(t2v, e2t)

Make topological faces.

Returns
- `t2v`: the attribute `tfid` lists the numbers of the topological faces.
"""
function make_topo_faces(t2v, crease_ang = 30/180*pi)
    e2v = ir_skeleton(t2v)
    e2v.right.attributes["geom"] = t2v.right.attributes["geom"]
    t2e = ir_bbyfacets(t2v, e2v)
    e2t = ir_transpose(t2e)
    e2t = _classify_mesh_edges(t2v, e2t, crease_ang)
    all_done = false
    tfid = fill(0, nshapes(t2v.left))
    while !all_done
        all_done = true
        for i in 1:nshapes(t2v.left)
            if tfid[i] == 0
                _collect_tf!(tfid, t2v, i, t2e, e2t)
                all_done = false;
                break;
            end
        end
    end
    t2v.left.attributes["tfid"] = VecAttrib(tfid)
    return t2v
end

function _collect_tf!(tfid, t2v, i, t2e, e2t)
    isonte = attribute(e2t.left, "isonte")
    currtfid = maximum(tfid) + 1
    tfid[i] = currtfid
    stack = Int[]
    push!(stack, i)
    while (!isempty(stack))
        j = pop!(stack)
        for e in 1:3
            eid = abs(t2e[j][e])
            if !isonte[eid]
                k = e2t[eid][1] == j ? e2t[eid][2] : e2t[eid][1]
                if tfid[k] == 0
                    tfid[k] = currtfid
                    push!(stack, k)
                end
            end
        end
    end
    return tfid
end


text:
module PDDL

using Base: @kwdef
using Julog
using AutoHashEquals

# Logical formulae and terms
export Term, Compound, Var, Const, Clause
# Abstract data types
export Domain, Problem, State, Action
# Generic concrete data types
export GenericDomain, GenericProblem, GenericState, GenericAction
# Interface methods
export satisfy, satisfiers, evaluate
export initstate, goalstate, transition, transition!
export available, execute, execute!, relevant, regress, regress!
# Parsing and writing
export parse_domain, parse_problem, parse_pddl, @pddl, @pddl_str
export write_domain, write_problem, write_pddl
export load_domain, load_problem
export save_domain, save_problem
# Simulators
export Simulator, EndStateSimulator, StateRecorder
# Abstract interpretation
export abstracted, abstractstate, lub, glb, widen, widen!
export AbstractedDomain
# Domain compilation
export compiled, compilestate
export CompiledDomain, CompiledAction, CompiledState
# Domain and action grounding
export ground, groundargs, groundactions
export GroundDomain, GroundAction, GroundActionGroup
# Extension interfaces
export attach!, register!
# Analysis tools
export infer_affected_fluents, infer_static_fluents
# Utilities
export find_matches, effect_diff, precond_diff

# PDDL requirement definitions and dependencies
include("requirements.jl")
# Abstract interface for PDDL-based planners and applications
include("interface/interface.jl")
# Generic concrete representations of PDDL types
include("generic/generic.jl")
# Parser for PDDL files
include("parser/parser.jl")
# Writer for PDDL files
include("writer/writer.jl")
# Simulators for executing action sequences and recording outputs
include("simulators/simulators.jl")
# Built-in functions and operators recognized by interpreters / compilers
include("builtins.jl")
# Methods for extending built-in functions and operators
include("extensions.jl")
# Theories for fluents with non-standard data types (e.g. sets, arrays)
include("theories/theories.jl")
# Abstractions for fluents of various types
include("abstractions/abstractions.jl")
# Interpreter-based interface implementation
include("interpreter/interpreter.jl")
# Compiler-based interface implementations
include("compiler/compiler.jl")
# Domain and action grounding
include("grounding/grounding.jl")
# Tools for analyzing domains
include("analysis/analysis.jl")

using .Parser, .Writer

end # module


text:
@symbol_func function v_te(cur_reactor::AbstractReactor, cur_rho::AbstractSymbol)
  cur_T = T_rho(cur_reactor, cur_rho)

  cur_v = 1.88e7

  cur_v *= sqrt(cur_T)

  cur_v
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 IoK8sApiAppsV1DeploymentStrategy <: SwaggerModel
    rollingUpdate::Any # spec type: Union{ Nothing, IoK8sApiAppsV1RollingUpdateDeployment } # spec name: rollingUpdate
    type::Any # spec type: Union{ Nothing, String } # spec name: type

    function IoK8sApiAppsV1DeploymentStrategy(;rollingUpdate=nothing, type=nothing)
        o = new()
        validate_property(IoK8sApiAppsV1DeploymentStrategy, Symbol("rollingUpdate"), rollingUpdate)
        setfield!(o, Symbol("rollingUpdate"), rollingUpdate)
        validate_property(IoK8sApiAppsV1DeploymentStrategy, Symbol("type"), type)
        setfield!(o, Symbol("type"), type)
        o
    end
end # type IoK8sApiAppsV1DeploymentStrategy

const _property_map_IoK8sApiAppsV1DeploymentStrategy = Dict{Symbol,Symbol}(Symbol("rollingUpdate")=>Symbol("rollingUpdate"), Symbol("type")=>Symbol("type"))
const _property_types_IoK8sApiAppsV1DeploymentStrategy = Dict{Symbol,String}(Symbol("rollingUpdate")=>"IoK8sApiAppsV1RollingUpdateDeployment", Symbol("type")=>"String")
Base.propertynames(::Type{ IoK8sApiAppsV1DeploymentStrategy }) = collect(keys(_property_map_IoK8sApiAppsV1DeploymentStrategy))
Swagger.property_type(::Type{ IoK8sApiAppsV1DeploymentStrategy }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiAppsV1DeploymentStrategy[name]))}
Swagger.field_name(::Type{ IoK8sApiAppsV1DeploymentStrategy }, property_name::Symbol) =  _property_map_IoK8sApiAppsV1DeploymentStrategy[property_name]

function check_required(o::IoK8sApiAppsV1DeploymentStrategy)
    true
end

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


text:
##############
# CONVERSION #
##############

function convert(::Type{RevisedAffineForm{N}}, x::T) where {N,T<:Real}
    TS = float(T)
    return RevisedAffineForm(x, SVector{N,TS}(zeros(TS, N)), zero(TS))
end

function convert(::Type{RevisedAffineForm{N}}, x::Interval{T}) where {N,T<:AbstractFloat}
    c, r = midpoint_radius(x)
    return RevisedAffineForm(c, SVector{N,T}(zeros(T, N)), r)
end

function convert(::Type{RevisedAffineForm{N,T}}, x::Real) where {N,T<:AbstractFloat}
    return RevisedAffineForm(convert(T, x), SVector{N,T}(zeros(T, N)), zero(T))
end

function convert(::Type{RevisedAffineForm{N,T}}, x::Interval) where {N,T<:AbstractFloat}
    c, r = midpoint_radius(x)
    return RevisedAffineForm(convert(T, c), SVector{N,T}(zeros(T, N)), convert(T, r))
end

#############
# PROMOTION #
#############

function promote_rule(
    ::Type{RevisedAffineForm{N,T}}, ::Type{Interval{S}}
) where {N,T<:AbstractFloat,S<:Real}
    return RevisedAffineForm{N,promote_type(T, S)}
end

function promote_rule(
    ::Type{RevisedAffineForm{N,T}}, ::Type{S}
) where {N,T<:AbstractFloat,S<:Real}
    return RevisedAffineForm{N,promote_type(T, S)}
end

function promote_rule(
    ::Type{RevisedAffineForm{N,T}}, ::Type{RevisedAffineForm{N,S}}
) where {N,T<:AbstractFloat,S<:AbstractFloat}
    return RevisedAffineForm{N,promote_type(T, S)}
end


text:
##
using DelimitedFiles
input_path = joinpath(@__DIR__(), "input.txt")
flip(x, y) = (y, x)
flip((x,y)) = (y, x)
flip(c::CartesianIndex) = CartesianIndex(flip(c.I))

##
input_data = readdlm(input_path,',',String)
encode_move(::Val{'U'}) = 0x00
encode_move(::Val{'D'}) = 0x11
encode_move(::Val{'R'}) = 0x01
encode_move(::Val{'L'}) = 0x10
data = map(input_data) do v
    k = encode_move(Val(v[1]))
    l = parse(Int, v[2:end])
    k, l
end |> permutedims
paths = [eachcol(data)...]

##
function decode_move(k, l)
    k == 0x00 && return CartesianIndex(0,l)
    k == 0x11 && return CartesianIndex(0,-l)
    k == 0x01 && return CartesianIndex(l,0)
    k == 0x10 && return CartesianIndex(-l,0)
    CartesianIndex(0,0)
end
decode_move(t::Tuple) = decode_move(t...)
function follow_path(paths::AbstractVector{<:AbstractVector{Tuple{UInt8,Int}}})
    map(paths) do mvs
        accumulate(mvs; init=CartesianIndex(0,0)) do cur, (k, l)
            cur + decode_move(k,l)
        end
    end
end
function intersections(paths::AbstractVector{<:AbstractVector{CartesianIndex{2}}})
    intersect(paths...)
end
dist(xs) = sum(abs,xs)
dist(c::CartesianIndex) = dist(c.I)

## TESTS
test_cases = (
    "R8,U5,L5,D3",
    "U7,R6,D4,L4",
    "R75,D30,R83,U83,L12,D49,R71,U7,L72,U62,R66,U55,R34,D71,R55,D58,R83",
    "R98,U47,R26,D63,R33,U87,L62,D20,R33,U53,R51,U98,R91,D20,R16,D67,R40,U7,R15,U6,R7",
)
res = map(test_cases) do x
    a = split(x,",")
    path = map(a) do v
        k = encode_move(Val(v[1]))
        l = parse(Int, v[2:end])
        k, l
    end
    follow_path([path])
end

##
test_dist = (
    "R8,U5,L5,D3",
    "U7,R6,D4,L4",
)
paths = map(test_dist) do x
    a = split(x,",")
    path = map(a) do v
        k = encode_move(Val(v[1]))
        l = parse(Int, v[2:end])
        k, l
    end
end
res = follow_path(paths|>collect)
display(res)

##
points = map(res) do ps
    prec = CartesianIndex(0,0)
    map(ps) do p
        x, y = p.I
        px, py = prec.I
        dx, dy = sign(x-px), sign(y-py)
        prec = p
        dx == 0 && return [CartesianIndex(px, j) for j in py+dy:dy:y]
        [CartesianIndex(j, py) for j in px+dx:dx:x]
    end
end
points = map(x -> reduce(vcat, x), points)
map(length, points)

##
res = follow_path(paths)
nothing

##
minimum(dist.(intersections(points)))

##
inters = intersections(points)

##
# steps2inters = map(points) do ps
#     findall(ps) do p
#         any(==(p), inters)
#     end
# end
steps2inters = map(inters) do i
    reduce(vcat, map(points) do ps
        findfirst(==(i), ps)
    end)
end

##
minimum(sum, steps2inters)


text:
module Models

using Parameters
using ACEtb.SlaterKoster: TwoCentreModel, SKH
using JuLIP.Potentials: fcut

import ACEtb.SlaterKoster: eval_bonds!, eval_diag!, getskh
import JuLIP: cutoff

export KwonTB

"""
`KwonTB <: SKModel`

Hamiltonian for an orthogonal sp TB model of Si developed by Kwon et al [1].

This implementation deviates  from [1] in how the cut-off is applied:
instead of "patching" a cubic spline between r1 and rcut, we simply multiply
with a quintic spline on the interval [0.5 (rcut + r0), rcut].

[1] I. Kwon, R. Biswas, C. Z. Wang, K. M. Ho and C. M. Soukoulis.
Transferable tight-binding models for silicon.
Phys Rev B 49 (11), 1994.
"""
@with_kw struct KwonTB <: TwoCentreModel
   skh::SKH = SKH("sp")
   # -------------------------------
   r0::Float64 = 2.360352   # Å
   Es::Float64 = -5.25      # eV
   Ep::Float64 = 1.2        # eV
   E0::Float64 = 8.7393204  # eV   8.7393204  is the original value; but it is just a constant shift
   # ------------------------------- Electronic Parameters
   # α   1     2    3     4
   #    ssσ   spσ  ppσ   ppπ
   hr0::NTuple{4, Float64} = (-2.038, 1.745, 2.75, -1.075)
   nc::NTuple{4, Float64} = (9.5, 8.5, 7.5, 7.5)
   rc::NTuple{4, Float64} = (3.4, 3.55, 3.7, 3.7)
   # ------------------------------- 2-body parameters
   m::Float64 = 6.8755
   mc::Float64 = 13.017
   dc::Float64 = 3.66995     # Å
   C::NTuple{4, Float64} = (2.1604385, -0.1384393, 5.8398423e-3, -8.0263577e-5)
   # --------------------------------
   r1::Float64 = 3.260176     # Å  (start of cut-off) >>> different r1 from Kwon paper
   rcut::Float64 = 4.16       # Å  (end of cut-off)
end

cutoff(model::KwonTB) = model.rcut

getskh(model::KwonTB) = model.skh

# isorthogonal(H::KwonTB) = true

kwon_bond(m::KwonTB, r::Real, α) = (
          m.hr0[α] * (m.r0 / r)^2 *
               exp( - 2 * (r/m.rc[α])^m.nc[α] + 2 * (m.r0/m.rc[α])^m.nc[α] ) )

function eval_bonds!(V, model::KwonTB, r)
   fc = fcut(r, model.r1, model.rcut)
   for α = 1:4
      V[α] = kwon_bond(model, r, α) * fc
   end
   return V
end

function eval_diag!(D, model::KwonTB, r, R)
   D[1] = model.Es
   D[2:4] .= model.Ep
   return D
end

end


text:
# # Fractional Zener
#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/notebooks/model_zener.ipynb)
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/notebooks/model_zener.ipynb)

using RHEOS
## include a helper function for plotting
include("assets/plothelper.jl");
#-

Fract_Zener

# #### Constitutive Equation
# ```math
# \sigma(t) + \frac{c_\alpha}{c_\beta} \frac{d^{\alpha-\beta} \sigma(t)}{dt^{\alpha-\beta}}= c_{\alpha} \frac{d^\alpha \epsilon(t)}{dt^\alpha} + c_\gamma \frac{d^\gamma \epsilon(t)}{dt^\gamma}+ \frac{c_\alpha c_\gamma}{c_\beta}\frac{d^{\alpha+\gamma-\beta} \epsilon(t)}{dt^{\alpha+\gamma-\beta}}
# ```

# ```math
# \text{for}\; \ 0 \leq \beta \leq \alpha \leq 1
# ```

# #### Relaxation Modulus
# ```math
# G(t) = c_\beta t^{-\beta} E_{\alpha-\beta,1-\beta}\left(-\frac{c_\beta}{c_\alpha} t^{\alpha-\beta}\right) + \frac{c_{\gamma} }{\Gamma(1-\gamma)} t^{-\gamma}
# ```

# #### Creep Modulus
# ```math
# \tilde{J}(s)= \frac{1}{s}\frac{c_\alpha s^{\alpha}+c_\beta s^{\beta}}{c_\alpha s^\alpha c_\beta s^{\beta} + c_\gamma s^\gamma (c_\alpha s^\alpha+c_\beta s^\beta)}
# ```

# #### Storage Modulus
# ```math
# G^{\prime}(\omega) = \frac{\left(c_\beta \omega^\beta\right)^2 \cdot c_\alpha \omega^\alpha \cos(\alpha \frac{\pi}{2}) + \left(c_\alpha \omega^\alpha\right)^2 \cdot c_\beta \omega^\beta \cos(\beta \frac{\pi}{2})}{\left(c_\alpha \omega^\alpha\right)^2+\left(c_\beta \omega^\beta\right)^2+2c_\alpha \omega^\alpha \cdot c_\beta \omega^\beta \cos((\alpha-\beta)\frac{\pi}{2})} + c_\gamma \omega^\gamma \cos\left(\gamma \frac{\pi}{2}\right)
# ```

# #### Loss Modulus
# ```math
# G^{\prime\prime}(\omega) = \frac{\left(c_\beta \omega^\beta\right)^2 \cdot c_\alpha \omega^\alpha \sin(\alpha \frac{\pi}{2}) + \left(c_\alpha \omega^\alpha\right)^2 \cdot c_\beta \omega^\beta \sin(\beta \frac{\pi}{2})}{\left(c_\alpha \omega^\alpha\right)^2+\left(c_\beta \omega^\beta\right)^2+2c_\alpha \omega^\alpha \cdot c_\beta \omega^\beta \cos((\alpha-\beta)\frac{\pi}{2})}+ c_\gamma \omega^\gamma \sin\left(\gamma \frac{\pi}{2}\right)
# ```

# # Fractional SLS (Zener)

FractSLS_Zener
#-

models = Vector{RheoModel}()

## plot moduli for varying α
for alpha in [0.1, 0.25, 0.5, 0.74, 0.9]

    push!(models, RheoModel(FractSLS_Zener, (cₐ = 1, a = alpha, kᵦ = 1, kᵧ = 1)))

end

plotmodel(models)

# # Fractional Jeffreys (Zener)

FractJeffreys_Zener
#-

models = Vector{RheoModel}()

## plot moduli for varying β
for beta in [0.1, 0.25, 0.5, 0.74, 0.9]

    push!(models, RheoModel(FractJeffreys_Zener, (ηₐ = 1, cᵦ = 1, β = beta, ηᵧ = 1)))

end

plotmodel(models, ymaxG = 2.0)

# # Standard Linear Solid (Zener)

SLS_Zener
#-

models = Vector{RheoModel}()

## plot moduli for varying kᵦ
for k in [1.0, 3.0, 5.0]

    push!(models, RheoModel(SLS_Zener, (η = 1, kᵦ = k, kᵧ = 1)))

end

plotmodel(models)

# # Jeffreys (Zener)

Jeffreys_Zener
#-

models = Vector{RheoModel}()

## plot moduli for varying ηₐ
for eta in [1.0, 5.0, 8.0]

    push!(models, RheoModel(Jeffreys_Zener, (ηₐ = eta, k = 3, ηᵧ = 1)))

end

plotmodel(models)

# # Fractional Solid

FractSolid
#-

models = Vector{RheoModel}()

## plot moduli for varying β
for beta in [0.1, 0.25, 0.5, 0.74, 0.9]

    push!(models, RheoModel(FractSolid, (η = 1, cᵦ = 1, β = beta, k = 0.5)))

end

plotmodel(models, ymaxG = 2)