text:
VERSION >= v"0.4.0-dev+6521" && __precompile__(true)

module IPNets
using Compat
import Base: IPAddr, IPv4, IPv6, parseipv4, parseipv6
import Base: length, size, endof, minimum, maximum, extrema, isless
import Base: in, contains, issubset, getindex
import Base: show, string, start, next

export
    # types
    IPNet, IPv4Net, IPv6Net, netmask

IPv4broadcast = typemax(UInt32)
IPv6broadcast = typemax(UInt128)

##################################################
# IPNet
##################################################

width(::Type{IPv4}) = @compat(UInt8(32))
width(::Type{IPv6}) = @compat(UInt8(128))

function _contiguousbitcount(n::Integer,t=UInt32)
    # takes an integer from 0 to 255 and a type, returns the number
    # of contiguous 1 bits in the number assuming it's of that type,
    # starting from the left.
    # cbc(240,UInt8) == 0x04 ("1111 0000")
    # cbc(252,UInt8) == 0x06 ("1111 1100")
    # cbc(127,UInt8) == error ("0111 1111")

    n = convert(t,n)
    invn = ~n
    bitct = log2(invn + 1)
    isinteger(bitct) || error("noncontiguous bits")

    bitct = @compat(floor(Int,bitct))
    return @compat(UInt8(sizeof(t)*8 - bitct))
end

function _mask2bits(t::Type, n::Unsigned)
    # takes a number of 1's bits in a
    # netmask and returns an integer representation
    maskbits = @compat(Int(width(t))) - @compat(Int(n))
    maskbits < 0 && throw(BoundsError())

    return (~(@compat(UInt128(2))^maskbits-1))
end

##################################################
# Network representations
##################################################
abstract IPNet

"""Returns the size of an IP network (# of hosts) as a tuple.
"""
function size(net::IPNet)
    numbits = width(typeof(net.netaddr)) - net.netmask
    return (big(2)^numbits, )
end

"""Returns the size of an IP network (# of hosts) as a tuple.
"""
length(net::IPNet) = size(net)[1]

"""String representation of an IP network"""
function string(net::IPNet)
    t = typeof(net)
    s = string("$t(\"")
    s = string(s, net.netaddr, "/", net.netmask, "\")")
    return s
end

function show(io::IO, net::IPNet)
    print(io, string(net))
end

# IP Networks are ordered first by starting network address
# and then by network mask. That is, smaller IP nets (with higher
# netmask values) are "less" than larger ones. This corresponds
# to secondary reordering by ending address.
function isless{T<:IPNet}(a::T, b::T)
    if a.netaddr == b.netaddr
        return isless(b.netmask, a.netmask)
    else
        return isless(a.netaddr, b.netaddr)
    end
end

function issubset{T<:IPNet}(a::T, b::T)
    astart, aend = extrema(a)
    bstart, bend = extrema(b)
    return (bstart <= astart <= aend <= bend)
end

"""Membership test for an IP address within an IP network"""
function in(ipaddr::IPAddr, net::IPNet)
    typeof(net.netaddr) == typeof(ipaddr) ||
        error("IPAddr is not the same type as IPNet")

    netstart = net.netaddr.host
    numbits = width(typeof(ipaddr)) - net.netmask
    netend = net.netaddr.host + big(2)^numbits - 1
    return netstart <= ipaddr.host <= netend
end

"""Membership test for an IP address within an IP network"""
function contains(net::IPNet, ipaddr::IPAddr)
    return in(ipaddr, net)
end

function getindex(net::IPNet, i::Integer)

    t = typeof(net.netaddr)
    ip = t(net.netaddr.host + i - 1)
    ip in net || throw(BoundsError())
    return ip
end

# Vector look-alikes
endof(net::IPNet) = @compat(UInt128(length(net)))
minimum(net::IPNet) = net[1]
maximum(net::IPNet) = net[end]
extrema(net::IPNet) = (minimum(net), maximum(net))
getindex(net::IPNet, r::Range) = [net[i] for i in r]
# getindex(net::IPNet, i::(Integer,)) = getindex(net,i[1])
start(net::IPNet) = net[1]
next{T<:IPAddr}(net::IPNet, s::T) = s, T(s.host + 1)
##################################################
# IPv4
##################################################

"""Type representing an IPv4 network"""
immutable IPv4Net <: IPNet
    netaddr::IPv4
    netmask::UInt8
    function IPv4Net(na::IPv4, nmi::Integer)
        (0 <= nmi <= width(IPv4)) || error("Invalid netmask")

        nm = @compat(UInt8(nmi))
        mask = _mask2bits(IPv4, nm)
        startip = @compat(UInt32(na.host & mask))
        new(IPv4(startip),nm)
    end
end

# "1.2.3.0/24"
function IPv4Net(ipmask::AbstractString)
    if search(ipmask,'/') > 0
        addrstr, netmaskstr = split(ipmask,"/")
        netmask = @compat(parse(UInt8,netmaskstr))
    else
        addrstr = ipmask
        netmask = width(IPv4)
    end
    netaddr = IPv4(addrstr)
    return IPv4Net(netaddr,netmask)
end

# "1.2.3.0", "255.255.255.0"
function IPv4Net(netaddr::AbstractString, netmask::AbstractString)
    netaddr = IPv4(netaddr).host
    netmask = _contiguousbitcount(IPv4(netmask).host)
    return IPv4Net(netaddr, netmask)
end

# 123872, 24
IPv4Net(ipaddr::Integer, netmask::Integer) = IPv4Net(IPv4(ipaddr), netmask)

# "(x,y)"
IPv4Net{A,M}(tuple::@compat(Tuple{A,M})) = IPv4Net(tuple[1],tuple[2])

# "1.2.3.0", 24
IPv4Net(netaddr::AbstractString, netmask::Integer) = IPv4Net(IPv4(netaddr), netmask)

"""Returns the netmask as an IPv4 address"""
netmask(n::IPv4Net) = IPv4(IPv4broadcast-2^(32-n.netmask)+1)

##################################################
# IPv6
##################################################

"""Type representing an IPv6 network"""
immutable IPv6Net <: IPNet
    # we treat the netmask as a potentially noncontiguous bitmask
    # for speed of calculation and consistency, but RFC2373, section
    # 2 provides for contiguous bitmasks only. We validate this
    # in the internal constructor. This wastes ~15 bytes per addr
    # for the benefit of rapid, consistent computation.
    netaddr::IPv6
    netmask::UInt32

    function IPv6Net(na::IPv6, nmi::Integer)
        (0 <= nmi <= width(IPv6)) || error("Invalid netmask")

        nm = @compat(UInt8(nmi))
        mask = _mask2bits(IPv6, nm)
        startip = @compat(UInt128(na.host & mask))
        return new(IPv6(startip), nm)
    end
end

# "2001::1/64"
function IPv6Net(ipmask::AbstractString)
    if search(ipmask,'/') > 0
        addrstr, netmaskbits = split(ipmask,"/")
        nmi = @compat(parse(Int,netmaskbits))
    else
        addrstr = ipmask
        nmi = width(IPv6)
    end
    netaddr = IPv6(addrstr)
    netmask = nmi
    return IPv6Net(netaddr,netmask)
end

# "2001::1", 64
function IPv6Net(netaddr::AbstractString, netmask::Integer)
    netaddr = IPv6(netaddr)
    return IPv6Net(netaddr, netmask)
end

# 123872, 128
IPv6Net(ipaddr::Integer, netmask::Integer) = IPv6Net(IPv6(ipaddr), netmask)

# (123872, 128)
IPv6Net{A,M}(tuple::@compat(Tuple{A,M})) = IPv6Net(tuple[1],tuple[2])

end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule libusb_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("libusb")
JLLWrappers.@generate_main_file("libusb", UUID("a877fdc9-fe69-5ed6-b93d-11ecd0dc2d49"))
end  # module libusb_jll

ext: jl


text:
# parameters
mass_bin_name = ARGS[1]
tslice = "t1"
path_wavelist = "src"
path_to_working_folder = "data"
Nbstrap_attempts = 100 # No. of bootstrap attempts

######################################################
using DelimitedFiles
push!(LOAD_PATH,"src")
using SDMHelper
using FittingPWALikelihood
using amplitudes_compass
using PWAHelper
using Random

###############################################################################
# Model setting

# read wavelist throw waves below threshol
wavelist = get_wavelist(joinpath(path_wavelist,"wavelist_formated.txt");
     path_to_thresholds=joinpath(path_wavelist,"thresholds_formated.txt"),
     M3pi=Meta.parse(mass_bin_name[1:4])/1000)

const Nwaves = size(wavelist,1)

# Blocks in wavelist
const noϵ =  [1] # flat wave
const posϵ = [i for (i,ϵ) in enumerate(wavelist[:,6]) if ϵ=="+"]
const negϵ = [i for (i,ϵ) in enumerate(wavelist[:,6]) if ϵ=="-"]

# Model description
const ModelBlocks = [noϵ, posϵ, negϵ, negϵ]
const Npar = get_npars(ModelBlocks)

# integrals from MC data
BmatMC = read_cmatrix(
    joinpath(path_to_working_folder,"integrmat_$(mass_bin_name)_$(tslice)_mc.txt"));

# get scale
const sqrtBscale = [sqrt(real(BmatMC[i,i])) for i in 1:size(BmatMC,1)];
const parscale = abs.(extnd(sqrtBscale, get_parameter_map(ModelBlocks, Nwaves)))
# normalize B
scaled_BmatMC = BmatMC ./ sqrtBscale ./ transpose(sqrtBscale)

# load precalculated data array
const PsiRD = read_precalc_basis(
    joinpath(path_to_working_folder,"functions_$(mass_bin_name)_$(tslice)_rd.bin"));

const Nd = size(PsiRD,1)
# normalize Psi
const scaled_PsiRD = PsiRD
scaled_PsiRD ./= transpose(sqrtBscale)
# for i in 1:size(PsiRD,2)
#     scaled_PsiRD[:,i] .*= 1.0/sqrtBscale[i]
# end

###############################################################################
# load the best minimum
path_and_filename = split(readline("data/llh_attmpts_$(mass_bin_name)_$(tslice).txt"),"\t")[1]
bestminpars = readdlm(path_and_filename)
best_minimum = vcat(bestminpars...)
best_minimum .*= parscale

###############################################################################
final_pars_all_att = Matrix{Float64}(undef, Npar+1, Nbstrap_attempts)
const scaled_pseudoPsiRD = copy(scaled_PsiRD)
#attempt at bootstrap
for b in 1:Nbstrap_attempts
    @show "Doing something"
    for e in 1:Nd
        scaled_pseudoPsiRD[e,:] .= scaled_PsiRD[rand(1:Nd),:]
    end
    # get functions to calculate LLH, derivative and hessian
    scaled_LLH, scaled_GRAD, scaled_LLH_and_GRAD!, scaled_HES =
        createLLHandGRAD(scaled_pseudoPsiRD, scaled_BmatMC, ModelBlocks);

    # fit
    @time minpars = minimize(scaled_LLH, scaled_LLH_and_GRAD!;
        algorithm = :LD_LBFGS, verbose=1, starting_pars=best_minimum,
        llhtolerance = 1e-4)
    # calculate LLH
    final_pars_all_att[1,b] = scaled_LLH(minpars)
    # scale parameters back
    final_pars = minpars ./ parscale
    final_pars_all_att[2:end,b] .= final_pars
end

# save results
output_name = joinpath(path_to_working_folder,"llhfit_bootstrap_$(mass_bin_name)_$(tslice).txt")
println("Done! Saving to $(output_name) ...")
writedlm(output_name, final_pars_all_att)


text:
using ADI
using HCIToolbox
using Test
using Statistics
using LinearAlgebra
using StableRNGs

rng = StableRNG(8799)

@testset "Median" begin include("median.jl") end
@testset "PCA" begin include("pca.jl") end
@testset "TPCA" begin include("tpca.jl") end
@testset "NMF" begin include("nmf.jl") end
@testset "GreeDS" begin include("greeds.jl") end

@testset "Metrics" begin include("metrics.jl") end

@testset "SDI" begin include("sdi.jl") end


text:
"""
Voxel corner and edge indexing conventions

        Z
        |

        5------5------6              Extra edges not drawn
       /|            /|              -----------
      8 |           6 |              - face diagonals
     /  9          /  10                - 13: 1 to 3
    8------7------7   |                 - 14: 1 to 8
    |   |         |   |                 - 15: 1 to 6
    |   1------1--|---2  -- Y           - 16: 5 to 7
    12 /          11 /                  - 17: 2 to 7
    | 4           | 2                   - 18: 4 to 7
    |/            |/                 - body diagonal
    4------3------3                     - 19: 1 to 7

  /
 X
"""
# this gets vectorized so we want to ensure it is the
# same type as out vertex
@inline function voxCrnrPos(::Type{PT}) where {PT}
    (PT(0, 0, 0),
    PT(0, 1, 0),
    PT(1, 1, 0),
    PT(1, 0, 0),
    PT(0, 0, 1),
    PT(0, 1, 1),
    PT(1, 1, 1),
    PT(1, 0, 1))
end

const voxCrnrPosInt = voxCrnrPos(SVector{3,UInt8})

# the voxel IDs at either end of the tetrahedra edges, by edge ID
const voxEdgeCrnrs = ((0x01, 0x02),
                      (0x02, 0x03),
                      (0x04, 0x03),
                      (0x01, 0x04),
                      (0x05, 0x06),
                      (0x06, 0x07),
                      (0x08, 0x07),
                      (0x05, 0x08),
                      (0x01, 0x05),
                      (0x02, 0x06),
                      (0x03, 0x07),
                      (0x04, 0x08),
                      (0x01, 0x03),
                      (0x01, 0x08),
                      (0x01, 0x06),
                      (0x05, 0x07),
                      (0x02, 0x07),
                      (0x04, 0x07),
                      (0x01, 0x07))

# direction codes:
# 0 => +x, 1 => +y, 2 => +z,
# 3 => +xy, 4 => +xz, 5 => +yz, 6 => +xyz
const voxEdgeDir = (0x01, 0x00, 0x01, 0x00, 0x01, 0x00, 0x01, 0x00, 0x02, 0x02, 0x02, 0x02, 0x03, 0x04, 0x05, 0x03, 0x04, 0x05, 0x06)

# For a pair of corner IDs, the edge ID joining them
# 0 denotes a pair with no edge
const voxEdgeIx = ((0x00, 0x01, 0x0d, 0x04, 0x09, 0x0f, 0x13, 0x0e),
                   (0x01, 0x00, 0x02, 0x00, 0x00, 0x0a, 0x11, 0x00),
                   (0x0d, 0x02, 0x00, 0x03, 0x00, 0x00, 0x0b, 0x00),
                   (0x04, 0x00, 0x03, 0x00, 0x00, 0x00, 0x12, 0x0c),
                   (0x09, 0x00, 0x00, 0x00, 0x00, 0x05, 0x10, 0x08),
                   (0x0f, 0x0a, 0x00, 0x00, 0x05, 0x00, 0x06, 0x00),
                   (0x13, 0x11, 0x0b, 0x12, 0x10, 0x06, 0x00, 0x07),
                   (0x0e, 0x00, 0x00, 0x0c, 0x08, 0x00, 0x07, 0x00))

# voxel corners that comprise each of the six tetrahedra
const subTets = ((0x01, 0x03, 0x02, 0x07),
                 (0x01, 0x08, 0x04, 0x07),
                 (0x01, 0x04, 0x03, 0x07),
                 (0x01, 0x02, 0x06, 0x07),
                 (0x01, 0x05, 0x08, 0x07),
                 (0x01, 0x06, 0x05, 0x07))
# tetrahedron corners for each edge (indices 1-4)
const tetEdgeCrnrs = ((0x01, 0x02),
                      (0x02, 0x03),
                      (0x01, 0x03),
                      (0x01, 0x04),
                      (0x02, 0x04),
                      (0x03, 0x04))

# triangle cases for a given tetrahedron edge code
const tetTri = ((0x01, 0x03, 0x04, 0x00, 0x00, 0x00),
                (0x01, 0x05, 0x02, 0x00, 0x00, 0x00),
                (0x03, 0x05, 0x02, 0x03, 0x04, 0x05),
                (0x02, 0x06, 0x03, 0x00, 0x00, 0x00),
                (0x01, 0x06, 0x04, 0x01, 0x02, 0x06),
                (0x01, 0x05, 0x06, 0x01, 0x06, 0x03),
                (0x04, 0x05, 0x06, 0x00, 0x00, 0x00),
                (0x04, 0x06, 0x05, 0x00, 0x00, 0x00),
                (0x01, 0x06, 0x05, 0x01, 0x03, 0x06),
                (0x01, 0x04, 0x06, 0x01, 0x06, 0x02),
                (0x02, 0x03, 0x06, 0x00, 0x00, 0x00),
                (0x03, 0x02, 0x05, 0x03, 0x05, 0x04),
                (0x01, 0x02, 0x05, 0x00, 0x00, 0x00),
                (0x01, 0x04, 0x03, 0x00, 0x00, 0x00))


text:
using DiskArrayTools
using Test
using DiskArrays: DiskArrays, ReshapedDiskArray, PermutedDiskArray, AbstractDiskArray

#Define a data structure that can be used for testing
struct _DiskArray{T,N,A<:AbstractArray{T,N}} <: AbstractDiskArray{T,N}
  getindex_list::Vector{NTuple{N,AbstractUnitRange}}
  setindex_list::Vector{NTuple{N,AbstractUnitRange}}
  parent::A
  chunksize::NTuple{N,Int}
end
_DiskArray(a;chunksize=size(a)) = _DiskArray(NTuple{ndims(a),AbstractUnitRange}[],NTuple{ndims(a),AbstractUnitRange}[],a,chunksize)
Base.size(a::_DiskArray) = size(a.parent)
DiskArrays.haschunks(::_DiskArray) = DiskArrays.Chunked()
DiskArrays.eachchunk(a::_DiskArray) = DiskArrays.GridChunks(a,a.chunksize)
getindex_count(a::_DiskArray) = length(a.getindex_list)
setindex_count(a::_DiskArray) = length(a.setindex_list)
trueparent(a::_DiskArray) = a.parent
getindex_count(a::ReshapedDiskArray) = getindex_count(a.parent)
setindex_count(a::ReshapedDiskArray) = setindex_count(a.parent)
trueparent(a::ReshapedDiskArray) = trueparent(a.parent)
getindex_count(a::PermutedDiskArray) = getindex_count(a.a.parent)
setindex_count(a::PermutedDiskArray) = setindex_count(a.a.parent)
trueparent(a::PermutedDiskArray{T,N,<:PermutedDimsArray{T,N,perm,iperm}}) where {T,N,perm,iperm} = permutedims(trueparent(a.a.parent),perm)
function DiskArrays.readblock!(a::_DiskArray,aout,i::AbstractUnitRange...)
  ndims(a) == length(i) || error("Number of indices is not correct")
  all(r->isa(r,AbstractUnitRange),i) || error("Not all indices are unit ranges")
  #println("reading from indices ", join(string.(i)," "))
  push!(a.getindex_list,i)
  aout .= a.parent[i...]
end
function DiskArrays.writeblock!(a::_DiskArray,v,i::AbstractUnitRange...)
  ndims(a) == length(i) || error("Number of indices is not correct")
  all(r->isa(r,AbstractUnitRange),i) || error("Not all indices are unit ranges")
  #println("Writing to indices ", join(string.(i)," "))
  push!(a.setindex_list,i)
  view(a.parent,i...) .= v
end

include("diskstack.jl")
include("interpolate.jl")
include("cfdiskarray.jl")
include("concatdiskarray.jl")


text:
#
# With minor changes from https://github.com/JuliaGaussianProcesses/AbstractGPs.jl/docs
#
### Process examples
# Always rerun examples
const EXAMPLES_OUT = joinpath(@__DIR__, "src", "examples")
ispath(EXAMPLES_OUT) && rm(EXAMPLES_OUT; recursive=true)
mkpath(EXAMPLES_OUT)

# Install and precompile all packages
# Workaround for https://github.com/JuliaLang/Pkg.jl/issues/2219
examples = filter!(isdir, readdir(joinpath(@__DIR__, "..", "examples"); join=true))
let script = "using Pkg; Pkg.activate(ARGS[1]); Pkg.instantiate()"
    for example in examples
        if !success(`$(Base.julia_cmd()) -e $script $example`)
            error(
                "project environment of example ",
                basename(example),
                " could not be instantiated",
            )
        end
    end
end
# Run examples asynchronously
processes = let literatejl = joinpath(@__DIR__, "literate.jl")
    map(examples) do example
        return run(
            pipeline(
                `$(Base.julia_cmd()) $literatejl $(basename(example)) $EXAMPLES_OUT`;
                stdin=devnull,
                stdout=devnull,
                stderr=stderr,
            );
            wait=false,
        )::Base.Process
    end
end

# Check that all examples were run successfully
isempty(processes) || success(processes) || error("some examples were not run successfully")

# Building Documenter
using Documenter
using AdvancedPS

DocMeta.setdocmeta!(AdvancedPS, :DocTestSetup, :(using AdvancedPS); recursive=true)

makedocs(;
    sitename="AdvancedPS",
    format=Documenter.HTML(),
    modules=[AdvancedPS],
    pages=[
        "Home" => "index.md",
        "api.md",
        "Examples" => [
            "example.md",
            map(
                (x) -> joinpath("examples", x),
                filter!(filename -> endswith(filename, ".md"), readdir(EXAMPLES_OUT)),
            )...,
        ],
    ],
    strict=true,
    checkdocs=:exports,
    doctestfilters=[
        # Older versions will show "0 element Array" instead of "Type[]".
        r"(Any\[\]|0-element Array{.+,[0-9]+})",
        # Older versions will show "Array{...,1}" instead of "Vector{...}".
        r"(Array{.+,\s?1}|Vector{.+})",
        # Older versions will show "Array{...,2}" instead of "Matrix{...}".
        r"(Array{.+,\s?2}|Matrix{.+})",
    ],
)

deploydocs(; repo="github.com/TuringLang/AdvancedPS.jl.git", push_preview=true)


text:
function io_next_token(io::IO, delim::Char, na::AbstractString,
                       line::Int, hasheader::Bool)
    isheader = line == 1 && hasheader
    eof(io) && return (val="", newline=false)
    buf = IOBuffer()
    c = read(io, Char)
    is_quoted = (c == QUOTE_CHAR)

    newline = false
    if is_quoted
        while true
            if eof(io)
                throw(ArgumentError("unexpected end of file in line $line"))
            end
            c = read(io, Char)
            if c == QUOTE_CHAR
                eof(io) && break
                c = read(io, Char)
                if c == QUOTE_CHAR
                    write(buf, QUOTE_CHAR)
                elseif c == delim
                    break
                elseif c == '\n'
                    newline = true
                    break
                elseif c == '\r'
                    newline = true
                    !eof(io) && Base.peek(io) == Int('\n') && read(io, Char)
                    break
                else
                    throw(ArgumentError("separator or newline expected "*
                                        "after closing quote in line $line"))
                end
            else
                write(buf, c)
            end
        end
    else
        if c == delim
            return (val = !isheader && na == "" ? missing : "", newline=false)
        end
        if c == '\n'
            return (val = !isheader && na == "" ? missing : "", newline=true)
        end
        if c == '\r'
            !eof(io) && Base.peek(io) == Int('\n') && read(io, Char)
            return (val = !isheader && na == "" ? missing : "", newline=true)
        end
        write(buf, c)
        while !eof(io)
            c = read(io, Char)
            c == delim && break
            if c == '\n'
                newline = true
                break
            end
            if c == '\r'
                newline = true
                eof(io) && break
                Base.peek(io) == Int('\n') && read(io, Char)
                break
            end
            write(buf, c)
        end
    end
    val = String(take!(buf))
    (val= !is_quoted && !isheader && val == na ? missing : val, newline=newline)
end

function ingest_csv(io::IO, delim::Char, na::AbstractString, hasheader::Bool)
    local line
    line_idx = 0
    data = Vector{Union{Missing,String}}[]
    newline = true
    while !eof(io)
        if newline
            line_idx += 1
            line = Union{Missing,String}[]
            push!(data, line)
        end
        val, newline = io_next_token(io, delim, na, line_idx, hasheader)
        push!(line, val)
    end
    data
end

lift_parse(val::Union{Missing, String}, parser::DataType) =
    ismissing(val) ? missing : tryparse(parser, val)

function try_parser(col::Vector{<:Union{Missing, String}}, parser::DataType)
    new_col = lift_parse.(col, parser)
    pos = findfirst(x -> x === nothing, new_col)
    pos === nothing && return new_col
    nothing
end

function rep_try_parser(col::Vector{<:Union{Missing, String}},
                        parsers::Vector)
    for parser in parsers
        new_col = try_parser(col, parser)
        new_col === nothing || return new_col
    end
    any(ismissing.(col)) ? col : Vector{String}(col)
end

"""
    read_csv(io::IO;
         delim::Char=',', header::Bool=true, na::AbstractString="",
         parsers::Vector = [Int, Float64])
    read_csv(filename::AbstractString;
             delim::Char=',', header::Bool=true, na::AbstractString="",
             parsers::Vector = [Int, Float64])

    Reads CSV from `io` stream or `filename` file to a `DataFrame` using
    `delim` separator and `na` as string for representing missing value.

    If `header` is `true` then first line of the file is assumed to contain
    column names. In such a case it is assumed that verbatim value `na` is allowed
    and is parsed as column name.

    `"` character is used for quoting fields. In quoted field use `""` to
    represent `"`.

    If `na` is read from CSV verbatim then `DataFrame` will contain `missing`
    in that entry. If a string equal to `na` should be read as `na` then then
    it should be quoted in source file as `"na"`.

    `parsers` controls what conversions `read_csv` tries to perform on
    data that is read in. Conversion is performed using `parse` function and
    `parsers` gives information in which sequence different types should be tried.

    Additionally you can control how much data is read in using the parameters:
    * `skiphead`: number of lines to skip on top of the file
    * `nrows`: if `nothing`: all lines are read in,
               if negative: number of parsed data lines to remove from the tail,
               if positive: number of parsed data lines to keep at head
"""
function read_csv(io::IO;
                  delim::Char=',', header::Bool=true, na::AbstractString="",
                  parsers::Vector = [Int, Float64],
                  skiphead::Int=0, nrows::Union{Int, Nothing}=nothing)
    if delim in [QUOTE_CHAR, '\n', '\r']
        throw(ArgumentError("delim is a quote char or a newline"))
    end
    if any(occursin.([QUOTE_CHAR, delim, '\n', '\r'], na))
        throw(ArgumentError("na contains quote, separator or a newline"))
    end
    for i in 1:skiphead
        readline(io)
    end

    data = ingest_csv(io, delim, na, header)
    if header
        if isempty(data)
            throw(ArgumentError("Source has zero rows and header" *
                                " was requested"))
        end
        if any(ismissing.(data[1]))
            @error "Unexpected error when parsing header of source"
        end
        names = Symbol.(popfirst!(data))
    end
    if !(nrows === nothing)
        if nrows < 0
            data = data[1:(end-nrows)]
        else
            data = data[1:nrows]
        end
    end
    if isempty(data)
        if header
            return DataFrame([[] for i in 1:length(names)], names)
        else
            return DataFrame()
        end
    end

    ref_length = length(data[1])
    if header && length(names) != ref_length
        throw(ArgumentError("data in line 2 has different number of " *
                            "columns than header"))
    end
    for (i, data_line) in enumerate(data)
        if length(data_line) != ref_length
            throw(ArgumentError("data in line $(i+header) has different " *
                                "number of columns than data in line $(i+header)"))
        end
    end
    cols = [rep_try_parser(getindex.(data, i), parsers) for i in 1:ref_length]
    header ? DataFrame(cols, names) : DataFrame(cols)
end

function read_csv(filename::AbstractString;
                  delim::Char=',', header::Bool=true, na::AbstractString="",
                  parsers::Vector = [Int, Float64],
                  skiphead::Int=0, nrows::Union{Int, Nothing}=nothing)
    open(filename) do io
        read_csv(io, delim=delim, header=header, na=na, parsers=parsers,
                 skiphead=skiphead, nrows=nrows)
    end
end


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

mutable struct EncryptionScopeKeyVaultProperties <: SwaggerModel
    keyUri::Any # spec type: Union{ Nothing, String } # spec name: keyUri

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

const _property_map_EncryptionScopeKeyVaultProperties = Dict{Symbol,Symbol}(Symbol("keyUri")=>Symbol("keyUri"))
const _property_types_EncryptionScopeKeyVaultProperties = Dict{Symbol,String}(Symbol("keyUri")=>"String")
Base.propertynames(::Type{ EncryptionScopeKeyVaultProperties }) = collect(keys(_property_map_EncryptionScopeKeyVaultProperties))
Swagger.property_type(::Type{ EncryptionScopeKeyVaultProperties }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_EncryptionScopeKeyVaultProperties[name]))}
Swagger.field_name(::Type{ EncryptionScopeKeyVaultProperties }, property_name::Symbol) =  _property_map_EncryptionScopeKeyVaultProperties[property_name]

function check_required(o::EncryptionScopeKeyVaultProperties)
    true
end

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


text:
module QuantumStateBase
    using DataDeps
    using LinearAlgebra
    using ClassicalOrthogonalPolynomials

    const DIM = 100

    function __init__()
        register(DataDep(
            "QuantumStateBase",
            """Data for QuantumStateBase.""",
            ""
        ))

        mkpath(joinpath(DataDeps.standard_loadpath[1], "QuantumStateBase"))
    end

    datadep_root() = datadep"QuantumStateBase"

    # polynomial
    include("polynomial.jl")

    # state
    include("representation.jl")
    include("basis.jl")
    include("operator.jl")
    include("state.jl")

    # wigner
    include("wigner_util.jl")
    include("wigner.jl")

    # pdf of quadrature
    # include("quadrature_pdf.jl")
    # include("sampler.jl")
end


text:
@testset "4.7 Miscellaneous" begin
    include("4.7.1 (c trig)^m (d trig)^n.jl")
    include("4.7.2 trig^m (a trig+b trig)^n.jl")
    include("4.7.3 (c+d x)^m trig^n trig^p.jl")
    include("4.7.4 x^m (a+b trig^n)^p.jl")
    include("4.7.5 x^m trig(a+b log(c x^n))^p.jl")
    include("4.7.6 f^(a+b x+c x^2) trig(d+e x+f x^2)^n.jl")
    include("4.7.7 Trig functions.jl")
end

max_stars_repo_path: test/4 Trig functions/4.7 Miscellaneous/.jl


text:
module ArrayIterators

export EachRow, EachCol

@static if VERSION < v"1.1"
    # added in https://github.com/JuliaLang/julia/pull/29749
    # once available in Compat, we should re-export here
    # https://github.com/JuliaLang/Compat.jl/issues/639
    error("Julia version $VERSION is not supported")
end

const EachRow_eachrow{A,I} = Base.Generator{I,typeof(eachrow(zeros(0,0)).f).name.wrapper{A}}
const EachCol_eachcol{A,I} = Base.Generator{I,typeof(eachcol(zeros(0,0)).f).name.wrapper{A}}

const EachRow_eachslice{A,I} = Base.Generator{I,typeof(eachslice(zeros(0,0),dims=1).f).name.wrapper{A,Tuple{},Tuple{Colon}}}
const EachCol_eachslice{A,I} = Base.Generator{I,typeof(eachslice(zeros(0,0),dims=2).f).name.wrapper{A,Tuple{Colon},Tuple{}}}

const EachRow{A,I} = Union{EachRow_eachrow{A,I}, EachRow_eachslice{A,I}}
const EachCol{A,I} = Union{EachCol_eachcol{A,I}, EachCol_eachslice{A,I}}

Base.parent(x::EachRow) = x.f.A
Base.parent(x::EachCol) = x.f.A

end # module


text:
using MCMCSamples
using Base.Test

@testset "MCMCDraws constructor and methods" begin
    va = 1:10
    vb = randn(length(va))
    draws = MCMCDraws(:b => vb, :a => va)
    @test keys(draws) == [:a, :b] # order important
    @test draws[:a] == va
    @test draws[:b] == vb
    @test length(draws) == length(va)
    @test MCMCDraws(:a => va, :b => vb).dict == draws.dict
    @test_throws ArgumentError MCMCDraws()
    @test_throws ArgumentError MCMCDraws(:a => va, :b => vb, :c => 1:8)
end

@testset "MCMCChains pooling" begin
    va1 = randn(10)
    vb1 = randn(length(va1))
    draws1 = MCMCDraws(:b => vb1, :a => va1)
    drawsc = MCMCDraws(:c => randn(length(va1)))
    @test MCMCSamples.all_same_keys([draws1, draws1])
    @test !MCMCSamples.all_same_keys([draws1, drawsc])
    va2 = randn(20)
    vb2 = randn(length(va2))
    draws2 = MCMCDraws(:a => va2, :b => vb2)
    chain1 = MCMCChain(draws1; discard = 5)
    chain2 = MCMCChain(draws2; discard = 2)
    drawsp = pool(chain1, chain2)
    @test length(drawsp) == length(draws1) - 5 + length(draws2) - 2
    @test keys(drawsp) == [:a, :b]
    @test drawsp[:a] == vcat(va1[6:end], va2[3:end])
    @test drawsp[:b] == vcat(vb1[6:end], vb2[3:end])
end


text:
#!/usr/bin/env julia

using Pkg
Pkg.add("HTTP")
Pkg.add("CodecZlib")
Pkg.add("EzXML")
Pkg.add("JuliaDB")
Pkg.add("Gumbo")
Pkg.add("Cascadia")


text:
import MsgPack
import UUIDs: UUID
import HTTP
import Sockets

import Base: endswith
function endswith(vec::Vector{T}, suffix::Vector{T}) where T
    local liv = lastindex(vec)
    local lis = lastindex(suffix)
    liv >= lis && (view(vec, (liv - lis + 1):liv) == suffix)
end

# to fix lots of false error messages from HTTP
# https://github.com/JuliaWeb/HTTP.jl/pull/546
# we do HTTP.Stream{HTTP.Messages.Request,S} instead of just HTTP.Stream to prevent the Julia warning about incremental compilation
function HTTP.closebody(http::HTTP.Stream{HTTP.Messages.Request,S}) where S <: IO
    if http.writechunked
        http.writechunked = false
        try
            write(http.stream, "0\r\n\r\n")
        catch end
    end
end

"""
    run([host,] port=1234)

Start Pluto! Are you excited? I am!

## Technobabble

The default `host` is `"127.0.0.1"`. For wild setups like Docker and heroku, you might need to change this to `"0.0.0.0"`.

This will start the static HTTP server and a WebSocket server. The server runs _synchronously_ (i.e. blocking call) on `http://[host]:[port]/`. Pluto notebooks can be started from the main menu in the web browser.
"""
function run(host, port::Union{Missing,Integer}; launchbrowser::Bool=false, session=ServerSession())
    pluto_router = http_router_for(session)

    hostIP = parse(Sockets.IPAddr, host)
    if ismissing(port)
        port, serversocket = Sockets.listenany(hostIP, UInt16(1234))
    else
        try
            serversocket = Sockets.listen(hostIP, UInt16(port))
        catch e
            @error "Port with number $port is already in use. Use Pluto.run() to automatically select an available port."
            return
        end
    end

    servertask = @async HTTP.serve(hostIP, UInt16(port), stream=true, server=serversocket) do http::HTTP.Stream
        # messy messy code so that we can use the websocket on the same port as the HTTP server

        if HTTP.WebSockets.is_upgrade(http.message)
            try
                requestURI = http.message.target |> HTTP.URIs.unescapeuri |> HTTP.URI
                @assert endswith(requestURI.path, string(session.secret))

                HTTP.WebSockets.upgrade(http) do clientstream
                    if !isopen(clientstream)
                        return
                    end
                    while !eof(clientstream)
                        # This stream contains data received over the WebSocket.
                        # It is formatted and MsgPack-encoded by send(...) in PlutoConnection.js
                        try
                            parentbody = let
                                # For some reason, long (>256*512 bytes) WS messages get split up - `readavailable` only gives the first 256*512
                                data = UInt8[]
                                while !endswith(data, MSG_DELIM)
                                    if eof(clientstream)
                                        if isempty(data)
                                            return
                                        end
                                        @warn "Unexpected eof after" data
                                        append!(data, MSG_DELIM)
                                        break
                                    end
                                    append!(data, readavailable(clientstream))
                                end
                                # TODO: view to avoid memory allocation
                                unpack(data[1:end - length(MSG_DELIM)])
                            end
                            process_ws_message(session, parentbody, clientstream)
                        catch ex
                            if ex isa InterruptException
                                rethrow(ex)
                                # TODO: this won't work, `upgrade` wraps our function in a try without catch
                            elseif ex isa HTTP.WebSockets.WebSocketError
                                # that's fine!
                            elseif ex isa InexactError
                                # that's fine! this is a (fixed) HTTP.jl bug: https://github.com/JuliaWeb/HTTP.jl/issues/471
                                # TODO: remove this switch
                            else
                                bt = stacktrace(catch_backtrace())
                                @warn "Reading WebSocket client stream failed for unknown reason:" exception = (ex, bt)
                            end
                        end
                    end
                end
            catch ex
                if ex isa InterruptException
                    rethrow(ex)
                elseif ex isa Base.IOError
                    # that's fine!
                elseif ex isa ArgumentError && occursin("stream is closed", ex.msg)
                    # that's fine!
                else
                    bt = stacktrace(catch_backtrace())
                    @warn "HTTP upgrade failed for unknown reason" exception = (ex, bt)
                end
            end
        else
            request::HTTP.Request = http.message
            request.body = read(http)
            HTTP.closeread(http)

            request_body = IOBuffer(HTTP.payload(request))
            if eof(request_body)
                # no request body
                response_body = HTTP.handle(pluto_router, request)
            else
                @warn "HTTP request contains a body, huh?" request_body
            end

            request.response::HTTP.Response = response_body
            request.response.request = request
            try
                HTTP.startwrite(http)
                write(http, request.response.body)
                HTTP.closewrite(http)
            catch e
                if isa(e, Base.IOError) || isa(e, ArgumentError)
                    # @warn "Attempted to write to a closed stream at $(request.target)"
                else
                    rethrow(e)
                end
            end
        end
    end

    address = let
        hostPretty = (hostStr = string(hostIP)) == "127.0.0.1" ? "localhost" : hostStr
        portPretty = Int(port)
        "http://$(hostPretty):$(portPretty)/"
    end
    Sys.set_process_title("Pluto server - $address")
    println("Go to $address to start writing ~ have fun!")
    println()
    println("Press Ctrl+C in this terminal to stop Pluto")
    println()

    launchbrowser && @warn "Not implemented yet"

    # create blocking call:
    try
        wait(servertask)
    catch e
        if isa(e, InterruptException)
            @sync begin
                println("\n\nClosing Pluto... Restart Julia for a fresh session. \n\nHave a nice day! 🎈")
                @async close(serversocket)
                # TODO: HTTP has a kill signal?
                # TODO: put do_work tokens back
                empty!(session.notebooks)
                for client in values(session.connected_clients)
                    @async close(client.stream)
                end
                empty!(session.connected_clients)
                for (notebook_id, ws) in WorkspaceManager.workspaces
                    @async WorkspaceManager.unmake_workspace(ws)
                end
            end
        else
            rethrow(e)
        end
    end
end

run(port::Union{Missing,Integer}=missing; kwargs...) = run("127.0.0.1", port; kwargs...)

"All messages sent over the WebSocket get decoded+deserialized and end up here."
function process_ws_message(session::ServerSession, parentbody::Dict, clientstream::IO)
    client_id = Symbol(parentbody["client_id"])
    client = get!(session.connected_clients, client_id, ClientSession(client_id, clientstream))

    messagetype = Symbol(parentbody["type"])
    request_id = Symbol(parentbody["request_id"])

    args = []
    if haskey(parentbody, "notebook_id")
        notebook = let
            notebook_id = UUID(parentbody["notebook_id"])
            get(session.notebooks, notebook_id, nothing)
        end

        if messagetype === :connect
            if notebook === nothing
                messagetype === :connect || @warn "Remote notebook not found locally!"
            else
                client.connected_notebook = notebook
            end
        end

        push!(args, notebook)

        if haskey(parentbody, "cell_id")
            cell_id = UUID(parentbody["cell_id"])
            index = cell_index_from_id(notebook, cell_id)
            if index === nothing
                @warn "Remote cell not found locally!"
            else
                push!(args, notebook.cells[index])
            end
        end
    end

    body = parentbody["body"]

    if haskey(responses, messagetype)
        responsefunc = responses[messagetype]
        try
            responsefunc(session, body, args..., initiator=Initiator(client.id, request_id))
        catch ex
            @warn "Response function to message of type $(messagetype) failed"
            rethrow(ex)
        end
    else
        @warn "Message of type $(messagetype) not recognised"
    end
end


text:
module MyInterpolations
export my_lin_interp
function my_lin_interp(grid,vals)
    function func(x::Real)
        if x <= grid[1]
            return (vals[2] - vals[1])/(grid[2] - grid[1])*(x - grid[1]) + vals[1]
        elseif x >= grid[end]
            return (vals[end] - vals[end-1])/(grid[end] - grid[end-1])*(x - grid[end]) + vals[end]
        else
            index = searchsortedlast(grid,x)
            return (vals[index+1] - vals[index])/(grid[index+1] - grid[index])*(x - grid[index]) + vals[index]
        end
    end
end
end


text:
"""Symbolic one-hot Boolean encoding of unknown integer value in the range `1:k`.
The value is `i` if `vars[i]` is true and `vars[j]` is false for all j ≠ i.
"""
struct BoundedInteger{k,N,V}
    name::N
    vars::V
end

BoundedInteger{k}(name::N, vars::V) where {k,N,V} = BoundedInteger{k,N,V}(name, vars)

Base.max(x::BoundedInteger{k}) where {k} = k

function BoundedInteger{k}(name) where {k}
    vars = [Variable(Symbol(name, subscript(i))) for i in 1:k]

    return BoundedInteger{k}(name, vars)
end

function decode(var_dict, x::BoundedInteger)
    values = [var_dict[v] for v in x.vars]

    num_true = count(values)

    if num_true ≠ 1
        error("Variable $x has not been successfully solved: values $(x.vars .=> values)")
    end

    return findfirst(values)
end

Base.show(io::IO, i::BoundedInteger) = print(io, i.name)

Base.getindex(i::BoundedInteger, n) = i.vars[n]

function clauses(i::BoundedInteger{k}) where {k}
    vars = i.vars

    # takes one of the possible values:
    has_value = [∨(vars...)]

    # but not two simultaneously:
    unique_value =
        [¬vars[i] ∨ ¬vars[j] for i in 1:k for j in i+1:k]

    return has_value ∪ unique_value
end

subscript_digit(i) = '₀' + i
subscript(i::Integer) = join(subscript_digit.(reverse(digits(i))))


text:
function vars_max(vars::Vector{CS.Variable})
    return MinMaxObjective(vars_max, [var.idx for var in vars])
end

function vars_max(com::CS.CoM, set_var_idx, var_val)
    objective = com.objective
    # return the minimal possible maximum value
    sense = com.sense
    search_space = com.search_space
    c_best_bound = minimum(values(search_space[1]))
    if sense == :Min
        for var_idx in objective.indices
            if var_idx == set_var_idx
                min_val = var_val
            else
                min_val = minimum(values(search_space[var_idx]))
            end
            if min_val > c_best_bound
                c_best_bound = min_val
            end
        end
    end
    @assert c_best_bound > 0
    return c_best_bound
end


text:
include("../utils.jl")
include("group_utils.jl")
using NLsolve

@testset "General Linear group" begin
    @testset "basic properties" begin
        G = GeneralLinear(3)
        @test G === GeneralLinear(3, ℝ)
        @test repr(G) == "GeneralLinear(3, ℝ)"
        @test base_manifold(G) === GeneralLinear(3)
        @test decorated_manifold(G) === Euclidean(3, 3)
        @test number_system(G) === ℝ
        @test manifold_dimension(G) == 9
        @test representation_size(G) == (3, 3)
        Gc = GeneralLinear(2, ℂ)
        @test decorated_manifold(Gc) === Euclidean(2, 2; field=ℂ)
        @test repr(Gc) == "GeneralLinear(2, ℂ)"
        @test number_system(Gc) == ℂ
        @test manifold_dimension(Gc) == 8
        @test representation_size(Gc) == (2, 2)
        Gh = GeneralLinear(4, ℍ)
        @test repr(Gh) == "GeneralLinear(4, ℍ)"
        @test number_system(Gh) == ℍ
        @test manifold_dimension(Gh) == 4 * 16
        @test representation_size(Gh) == (4, 4)

        @test (@inferred invariant_metric_dispatch(G, LeftAction())) === Val(true)
        @test (@inferred invariant_metric_dispatch(G, RightAction())) === Val(false)
        @test is_default_metric(
            MetricManifold(G, InvariantMetric(EuclideanMetric(), LeftAction())),
        ) === true
        @test @inferred(Manifolds.default_metric_dispatch(G, EuclideanMetric())) ===
              Val(true)
        @test @inferred(
            Manifolds.default_metric_dispatch(
                G,
                InvariantMetric(EuclideanMetric(), LeftAction()),
            )
        ) === Val(true)
        @test @inferred(
            Manifolds.default_metric_dispatch(
                MetricManifold(G, InvariantMetric(EuclideanMetric(), LeftAction())),
            )
        ) === Val(true)
        @test Manifolds.allocation_promotion_function(Gc, exp!, (1,)) === complex

        q = identity_element(G)
        @test is_identity(G, q)
        @test isapprox(G, q, Identity(G))
        @test isapprox(G, Identity(G), q)
    end

    @testset "GL(1,𝔽) special cases" begin
        @testset "real" begin
            G = GeneralLinear(1)
            p = 3.0 * ones(1, 1)
            X = 1.0 * ones(1, 1)
            @test exp(G, p, X) ≈ p * exp(X)' * exp(X - X')
            q = exp(G, p, X)
            Y = log(G, p, q)
            @test Y ≈ X
            @test exp_lie(G, X) ≈ exp(X)
            @test log_lie(G, exp(X)) ≈ X
            @test log_lie(G, [1.0]) == zeros(1) # vector to vector
            log_lie(G, Identity(G)) == zeros(1, 1) # Matrix to matrix
        end
        @testset "complex" begin
            G = GeneralLinear(1, ℂ)
            p = (1 + im) * ones(1, 1)
            X = (1 - im) * ones(1, 1)
            @test exp(G, p, X) ≈ p * exp(X)' * exp(X - X')
            q = exp(G, p, X)
            Y = log(G, p, q)
            @test Y ≈ X
            @test exp_lie(G, X) ≈ exp(X)
            @test log_lie(G, exp(X)) ≈ X
        end
    end

    @testset "Real" begin
        G = GeneralLinear(3)

        @test_throws DomainError is_point(G, randn(2, 3), true)
        @test_throws DomainError is_point(G, randn(2, 2), true)
        @test_throws DomainError is_point(G, randn(ComplexF64, 3, 3), true)
        @test_throws DomainError is_point(G, zeros(3, 3), true)
        @test_throws DomainError is_point(G, Float64[0 0 0; 0 1 1; 1 1 1], true)
        @test is_point(G, Float64[0 0 1; 0 1 1; 1 1 1], true)
        @test is_point(G, Identity(G), true)
        @test_throws DomainError is_vector(
            G,
            Float64[0 1 1; 0 1 1; 1 0 0],
            randn(3, 3),
            true,
        )
        @test is_vector(G, Float64[0 0 1; 0 1 1; 1 1 1], randn(3, 3), true)

        types = [Matrix{Float64}]
        pts = [
            Matrix(Diagonal([1, 2, 3])),
            [-2 5 -5; 0 2 -1; -3 -5 -2],
            [-5 1 0; 1 0 1; 0 1 3],
        ]
        vpts = [[-1 -2 0; -2 1 -2; 2 0 2], [1 1 1; 0 0 -2; 2 0 0]]

        retraction_methods = [
            Manifolds.GroupExponentialRetraction(LeftAction()),
            Manifolds.GroupExponentialRetraction(RightAction()),
        ]

        inverse_retraction_methods = [
            Manifolds.GroupLogarithmicInverseRetraction(LeftAction()),
            Manifolds.GroupLogarithmicInverseRetraction(RightAction()),
        ]

        basis_types = [DefaultOrthonormalBasis(), ProjectedOrthonormalBasis(:svd)]

        for T in types
            gpts = convert.(T, pts)
            vgpts = convert.(T, vpts)
            test_group(
                G,
                gpts,
                vgpts,
                vgpts;
                test_diff=true,
                test_invariance=true,
                test_lie_bracket=true,
                test_adjoint_action=true,
            )
            test_manifold(
                G,
                gpts;
                test_reverse_diff=false,
                test_forward_diff=false,
                test_project_point=true,
                test_injectivity_radius=false,
                test_project_tangent=true,
                test_musical_isomorphisms=true,
                test_default_vector_transport=true,
                vector_transport_methods=[
                    ParallelTransport(),
                    SchildsLadderTransport(),
                    PoleLadderTransport(),
                ],
                retraction_methods=retraction_methods,
                inverse_retraction_methods=inverse_retraction_methods,
                basis_types_vecs=basis_types,
                basis_types_to_from=basis_types,
                exp_log_atol_multiplier=1e7,
                retraction_atol_multiplier=1e7,
            )
        end
    end

    @testset "Complex" begin
        G = GeneralLinear(2, ℂ)

        @test_throws DomainError is_point(G, randn(ComplexF64, 2, 3), true)
        @test_throws DomainError is_point(G, randn(ComplexF64, 3, 3), true)
        @test_throws DomainError is_point(G, zeros(2, 2), true)
        @test_throws DomainError is_point(G, ComplexF64[1 im; 1 im], true)
        @test is_point(G, ComplexF64[1 1; im 1], true)
        @test is_point(G, Identity(G), true)
        @test_throws DomainError is_point(G, Float64[0 0 0; 0 1 1; 1 1 1], true)
        @test_throws DomainError is_vector(
            G,
            ComplexF64[im im; im im],
            randn(ComplexF64, 2, 2),
            true,
        )
        @test is_vector(G, ComplexF64[1 im; im im], randn(ComplexF64, 2, 2), true)

        types = [Matrix{ComplexF64}]
        pts = [
            [-1-5im -1+3im; -6-4im 4+6im],
            [1+3im -1-4im; -2-2im -3-1im],
            [-6+0im 1+1im; 1-1im -4+0im],
        ]
        vpts = [[1+0im -2-1im; -1-2im -4+1im], [-2+2im -1-1im; -1-1im -3+0im]]

        retraction_methods = [
            Manifolds.GroupExponentialRetraction(LeftAction()),
            Manifolds.GroupExponentialRetraction(RightAction()),
        ]

        inverse_retraction_methods = [
            Manifolds.GroupLogarithmicInverseRetraction(LeftAction()),
            Manifolds.GroupLogarithmicInverseRetraction(RightAction()),
        ]

        for T in types
            gpts = convert.(T, pts)
            vgpts = convert.(T, vpts)
            test_group(G, gpts, vgpts, vgpts; test_diff=true, test_invariance=true)
            test_manifold(
                G,
                gpts;
                test_reverse_diff=false,
                test_forward_diff=false,
                test_project_point=true,
                test_injectivity_radius=false,
                test_project_tangent=true,
                test_musical_isomorphisms=true,
                test_default_vector_transport=true,
                vector_transport_methods=[
                    ParallelTransport(),
                    SchildsLadderTransport(),
                    PoleLadderTransport(),
                ],
                retraction_methods=retraction_methods,
                inverse_retraction_methods=inverse_retraction_methods,
                exp_log_atol_multiplier=1e8,
                retraction_atol_multiplier=1e8,
            )
        end
    end
end


text:
# Three letters (for PDB)
# =======================

# download("ftp://ftp.wwpdb.org/pub/pdb/data/monomers/components.cif","components.cif")
#
# open("components.cif","r") do fh
#     m3 = nothing
#     m1 = nothing
#     peptide = false
#     for line in eachline(fh)
#         if startswith(line, "_chem_comp")
#             if m3 == nothing
#                 m3 = match(r"^_chem_comp.three_letter_code\s+(\w{3})\s*$",line) # It avoids: DA, DT, ...
#             end
#             if m1 == nothing
#                 m1 = match(r"^_chem_comp.one_letter_code\s+(\w{1})\s*$",line) # It avoids: ?
#             end
#             if startswith(line,"_chem_comp.type")
#                 peptide = occursin(r"peptide"i,line) # It avoids 0DA, 2DA, etc...
#             end
#             if m1 != nothing && m3 != nothing && parent != nothing && peptide
#                 res = collect(m1[1])[1]
#                 if (Residue(res) != XAA) # It avoids non standar residues, i.e. PYL
#                     println("\"$(m3[1])\" => '$(Char(Residue(res)))',")
#                 end
#                 m3 = nothing
#                 m1 = nothing
#                 peptide = false
#             end
#         else
#             m3 = nothing
#             m1 = nothing
#             peptide = false
#         end
#     end
# end

const _res2three = [ "ALA", "ARG", "ASN", "ASP", "CYS", "GLN", "GLU", "GLY", "HIS", "ILE",
                     "LEU", "LYS", "MET", "PHE", "PRO", "SER", "THR", "TRP", "TYR", "VAL",
                     "   ", "XAA" ]

"""
This function returns the three letter name of the `Residue`.

```jldoctest
julia> using MIToS.MSA

julia> residue2three(Residue('G'))
"GLY"

```
"""
function residue2three(res::Residue)
    int_res = Int(res)
    if int_res == 21  || !isvalid(res)
        #         GAP
        throw(ErrorException("Residue($(int_res)) has not three letter name."))
    end
    _res2three[int_res]
end

"""
It takes a three letter residue name and returns the corresponding `Residue`.
If the name isn't in the MIToS dictionary, a `XAA` is returned.

```jldoctest
julia> using MIToS.MSA

julia> three2residue("ALA")
A

```
"""
function three2residue(res::String)
    if length(res) == 3
        get(_three2res, uppercase(res), XAA)
    else
        throw(ErrorException("The residue name should have 3 letters."))
    end
end

const _three2res = Dict{String, Residue}(
"00C" => 'C',
"02K" => 'A',
"02L" => 'N',
"03Y" => 'C',
"07O" => 'C',
"08P" => 'C',
"0A0" => 'D',
"0A1" => 'Y',
"0A2" => 'K',
"0A8" => 'C',
"0AA" => 'V',
"0AB" => 'V',
"0AC" => 'G',
"0AF" => 'W',
"0AG" => 'L',
"0AH" => 'S',
"0AK" => 'D',
"0BN" => 'F',
"0CS" => 'A',
"0E5" => 'T',
"0EA" => 'Y',
"0FL" => 'A',
"0NC" => 'A',
"0WZ" => 'Y',
"0Y8" => 'P',
"143" => 'C',
"1OP" => 'Y',
"1PA" => 'F',
"1PI" => 'A',
"1TQ" => 'W',
"1TY" => 'Y',
"1X6" => 'S',
"200" => 'F',
"23F" => 'F',
"26B" => 'T',
"2AG" => 'A',
"2CO" => 'C',
"2FM" => 'M',
"2HF" => 'H',
"2KK" => 'K',
"2KP" => 'K',
"2LU" => 'L',
"2ML" => 'L',
"2MR" => 'R',
"2MT" => 'P',
"2OR" => 'R',
"2QZ" => 'T',
"2R3" => 'Y',
"2TL" => 'T',
"2TY" => 'Y',
"2VA" => 'V',
"2XA" => 'C',
"33X" => 'A',
"3AH" => 'H',
"3CF" => 'F',
"3GA" => 'A',
"3MD" => 'D',
"3NF" => 'Y',
"3QN" => 'K',
"3XH" => 'G',
"4BF" => 'Y',
"4CF" => 'F',
"4CY" => 'M',
"4DP" => 'W',
"4FB" => 'P',
"4FW" => 'W',
"4GJ" => 'C',
"4HT" => 'W',
"4IN" => 'W',
"4PH" => 'F',
"4U7" => 'A',
"56A" => 'H',
"5AB" => 'A',
"5CR" => 'F',
"5CS" => 'C',
"5CT" => 'K',
"5CW" => 'W',
"5HP" => 'E',
"5OW" => 'K',
"5VV" => 'N',
"6CL" => 'K',
"6CW" => 'W',
"6GL" => 'A',
"6HN" => 'K',
"6V1" => 'C',
"6WK" => 'C',
"7JA" => 'I',
"9NE" => 'E',
"9NF" => 'F',
"9NR" => 'R',
"9NV" => 'V',
"A5N" => 'N',
"AA3" => 'A',
"AA4" => 'A',
"AAR" => 'R',
"ABA" => 'A',
"ACB" => 'D',
"ACL" => 'R',
"AEI" => 'D',
"AFA" => 'N',
"AGM" => 'R',
"AGT" => 'C',
"AHB" => 'N',
"AHO" => 'A',
"AHP" => 'A',
"AIB" => 'A',
"AKL" => 'D',
"AKZ" => 'D',
"ALA" => 'A',
"ALC" => 'A',
"ALM" => 'A',
"ALN" => 'A',
"ALO" => 'T',
"ALS" => 'A',
"ALT" => 'A',
"ALV" => 'A',
"ALY" => 'K',
"AN8" => 'A',
"APH" => 'A',
"API" => 'K',
"APK" => 'K',
"AR2" => 'R',
"AR4" => 'E',
"AR7" => 'R',
"ARG" => 'R',
"ARM" => 'R',
"ARO" => 'R',
"AS2" => 'D',
"ASA" => 'D',
"ASB" => 'D',
"ASI" => 'D',
"ASK" => 'D',
"ASL" => 'D',
"ASN" => 'N',
"ASP" => 'D',
"ASQ" => 'D',
"AYA" => 'A',
"AZK" => 'K',
"AZS" => 'S',
"AZY" => 'Y',
"B1F" => 'F',
"B2A" => 'A',
"B2F" => 'F',
"B2I" => 'I',
"B2V" => 'V',
"B3A" => 'A',
"B3D" => 'D',
"B3E" => 'E',
"B3K" => 'K',
"B3S" => 'S',
"B3U" => 'H',
"B3X" => 'N',
"B3Y" => 'Y',
"BB6" => 'C',
"BB7" => 'C',
"BB8" => 'F',
"BB9" => 'C',
"BBC" => 'C',
"BCS" => 'C',
"BFD" => 'D',
"BG1" => 'S',
"BH2" => 'D',
"BHD" => 'D',
"BIF" => 'F',
"BIU" => 'I',
"BL2" => 'L',
"BLE" => 'L',
"BLY" => 'K',
"BMT" => 'T',
"BNN" => 'F',
"BOR" => 'R',
"BPE" => 'C',
"BSE" => 'S',
"BTA" => 'L',
"BTC" => 'C',
"BTR" => 'W',
"BUC" => 'C',
"BUG" => 'V',
"C1X" => 'K',
"C22" => 'A',
"C3Y" => 'C',
"C4R" => 'C',
"C5C" => 'C',
"C6C" => 'C',
"CAF" => 'C',
"CAS" => 'C',
"CAY" => 'C',
"CCL" => 'K',
"CCS" => 'C',
"CEA" => 'C',
"CG6" => 'C',
"CGA" => 'E',
"CGU" => 'E',
"CHP" => 'G',
"CIR" => 'R',
"CLE" => 'L',
"CLG" => 'K',
"CLH" => 'K',
"CME" => 'C',
"CMH" => 'C',
"CML" => 'C',
"CMT" => 'C',
"CR5" => 'G',
"CS0" => 'C',
"CS1" => 'C',
"CS3" => 'C',
"CS4" => 'C',
"CSA" => 'C',
"CSB" => 'C',
"CSD" => 'C',
"CSE" => 'C',
"CSJ" => 'C',
"CSO" => 'C',
"CSP" => 'C',
"CSR" => 'C',
"CSS" => 'C',
"CSU" => 'C',
"CSW" => 'C',
"CSX" => 'C',
"CSZ" => 'C',
"CTE" => 'W',
"CTH" => 'T',
"CWR" => 'S',
"CXM" => 'M',
"CY0" => 'C',
"CY1" => 'C',
"CY3" => 'C',
"CY4" => 'C',
"CYA" => 'C',
"CYD" => 'C',
"CYF" => 'C',
"CYG" => 'C',
"CYJ" => 'K',
"CYM" => 'C',
"CYQ" => 'C',
"CYR" => 'C',
"CYS" => 'C',
"CZ2" => 'C',
"CZZ" => 'C',
"D11" => 'T',
"D3P" => 'G',
"DAB" => 'A',
"DAH" => 'F',
"DAL" => 'A',
"DAR" => 'R',
"DAS" => 'D',
"DBB" => 'T',
"DBS" => 'S',
"DBU" => 'T',
"DBY" => 'Y',
"DBZ" => 'A',
"DC2" => 'C',
"DCY" => 'C',
"DDE" => 'H',
"DGH" => 'G',
"DGL" => 'E',
"DGN" => 'Q',
"DHA" => 'S',
"DHI" => 'H',
"DHN" => 'V',
"DHV" => 'V',
"DI7" => 'Y',
"DIL" => 'I',
"DIR" => 'R',
"DIV" => 'V',
"DLE" => 'L',
"DLS" => 'K',
"DLY" => 'K',
"DM0" => 'K',
"DMH" => 'N',
"DMK" => 'D',
"DNE" => 'L',
"DNL" => 'K',
"DNP" => 'A',
"DNS" => 'K',
"DOH" => 'D',
"DON" => 'L',
"DPL" => 'P',
"DPN" => 'F',
"DPP" => 'A',
"DPQ" => 'Y',
"DPR" => 'P',
"DSE" => 'S',
"DSG" => 'N',
"DSN" => 'S',
"DSP" => 'D',
"DTH" => 'T',
"DTR" => 'W',
"DTY" => 'Y',
"DVA" => 'V',
"DYS" => 'C',
"ECC" => 'Q',
"EFC" => 'C',
"EHP" => 'F',
"ESB" => 'Y',
"ESC" => 'M',
"EXY" => 'L',
"F2F" => 'F',
"FAK" => 'K',
"FB5" => 'A',
"FB6" => 'A',
"FCL" => 'F',
"FGA" => 'E',
"FGL" => 'G',
"FGP" => 'S',
"FH7" => 'K',
"FHL" => 'K',
"FHO" => 'K',
"FLA" => 'A',
"FLE" => 'L',
"FLT" => 'Y',
"FME" => 'M',
"FOE" => 'C',
"FP9" => 'P',
"FT6" => 'W',
"FTR" => 'W',
"FTY" => 'Y',
"FVA" => 'V',
"FZN" => 'K',
"GAU" => 'E',
"GFT" => 'S',
"GGL" => 'E',
"GHG" => 'Q',
"GHP" => 'G',
"GL3" => 'G',
"GLH" => 'Q',
"GLJ" => 'E',
"GLK" => 'E',
"GLN" => 'Q',
"GLQ" => 'E',
"GLU" => 'E',
"GLY" => 'G',
"GLZ" => 'G',
"GMA" => 'E',
"GPL" => 'K',
"GSC" => 'G',
"GSU" => 'E',
"GT9" => 'C',
"GVL" => 'S',
"H14" => 'F',
"H5M" => 'P',
"HAC" => 'A',
"HAR" => 'R',
"HBN" => 'H',
"HHI" => 'H',
"HIA" => 'H',
"HIC" => 'H',
"HIP" => 'H',
"HIQ" => 'H',
"HIS" => 'H',
"HL2" => 'L',
"HLU" => 'L',
"HMR" => 'R',
"HPC" => 'F',
"HPE" => 'F',
"HPH" => 'F',
"HPQ" => 'F',
"HQA" => 'A',
"HRG" => 'R',
"HRP" => 'W',
"HS8" => 'H',
"HS9" => 'H',
"HSE" => 'S',
"HSL" => 'S',
"HSO" => 'H',
"HSV" => 'H',
"HTI" => 'C',
"HTN" => 'N',
"HTR" => 'W',
"HV5" => 'A',
"HVA" => 'V',
"HY3" => 'P',
"HYP" => 'P',
"HZP" => 'P',
"I2M" => 'I',
"I58" => 'K',
"IAM" => 'A',
"IAR" => 'R',
"IAS" => 'D',
"IEL" => 'K',
"IGL" => 'G',
"IIL" => 'I',
"ILE" => 'I',
"ILG" => 'E',
"ILX" => 'I',
"IML" => 'I',
"IOY" => 'F',
"IPG" => 'G',
"IT1" => 'K',
"IYR" => 'Y',
"IYT" => 'T',
"IZO" => 'M',
"JJJ" => 'C',
"JJK" => 'C',
"JJL" => 'C',
"K1R" => 'C',
"KCX" => 'K',
"KGC" => 'K',
"KNB" => 'A',
"KOR" => 'M',
"KPF" => 'K',
"KPI" => 'K',
"KST" => 'K',
"KYN" => 'W',
"KYQ" => 'K',
"LA2" => 'K',
"LAA" => 'D',
"LAL" => 'A',
"LBY" => 'K',
"LCK" => 'K',
"LCX" => 'K',
"LDH" => 'K',
"LED" => 'L',
"LEF" => 'L',
"LEH" => 'L',
"LEI" => 'V',
"LEM" => 'L',
"LEN" => 'L',
"LET" => 'K',
"LEU" => 'L',
"LEX" => 'L',
"LLP" => 'K',
"LLY" => 'K',
"LME" => 'E',
"LMF" => 'K',
"LMQ" => 'Q',
"LP6" => 'K',
"LPD" => 'P',
"LPG" => 'G',
"LPS" => 'S',
"LSO" => 'K',
"LTR" => 'W',
"LVG" => 'G',
"LVN" => 'V',
"LYF" => 'K',
"LYK" => 'K',
"LYM" => 'K',
"LYN" => 'K',
"LYR" => 'K',
"LYS" => 'K',
"LYX" => 'K',
"LYZ" => 'K',
"M0H" => 'C',
"M2L" => 'K',
"M2S" => 'M',
"M30" => 'G',
"M3L" => 'K',
"MAA" => 'A',
"MAI" => 'R',
"MBQ" => 'Y',
"MC1" => 'S',
"MCL" => 'K',
"MCS" => 'C',
"MD3" => 'C',
"MD6" => 'G',
"MDF" => 'Y',
"MEA" => 'F',
"MED" => 'M',
"MEG" => 'E',
"MEN" => 'N',
"MEQ" => 'Q',
"MET" => 'M',
"MEU" => 'G',
"MGG" => 'R',
"MGN" => 'Q',
"MGY" => 'G',
"MHL" => 'L',
"MHO" => 'M',
"MHS" => 'H',
"MIS" => 'S',
"MK8" => 'L',
"ML3" => 'K',
"MLE" => 'L',
"MLL" => 'L',
"MLY" => 'K',
"MLZ" => 'K',
"MME" => 'M',
"MMO" => 'R',
"MND" => 'N',
"MNL" => 'L',
"MNV" => 'V',
"MP8" => 'P',
"MPQ" => 'G',
"MSA" => 'G',
"MSE" => 'M',
"MSL" => 'M',
"MSO" => 'M',
"MT2" => 'M',
"MTY" => 'Y',
"MVA" => 'V',
"N10" => 'S',
"N7P" => 'P',
"N80" => 'P',
"N8P" => 'P',
"NA8" => 'A',
"NAL" => 'A',
"NAM" => 'A',
"NB8" => 'N',
"NBQ" => 'Y',
"NC1" => 'S',
"NCB" => 'A',
"NDF" => 'F',
"NEM" => 'H',
"NEP" => 'H',
"NFA" => 'F',
"NHL" => 'E',
"NIY" => 'Y',
"NLE" => 'L',
"NLN" => 'L',
"NLO" => 'L',
"NLP" => 'L',
"NLQ" => 'Q',
"NLW" => 'L',
"NMC" => 'G',
"NMM" => 'R',
"NNH" => 'R',
"NPH" => 'C',
"NPI" => 'A',
"NTR" => 'Y',
"NTY" => 'Y',
"NVA" => 'V',
"NYS" => 'C',
"NZH" => 'H',
"OAR" => 'R',
"OAS" => 'S',
"OBS" => 'K',
"OCS" => 'C',
"OCY" => 'C',
"OHI" => 'H',
"OHS" => 'D',
"OLT" => 'T',
"OLZ" => 'S',
"OMT" => 'M',
"ONH" => 'A',
"OPR" => 'R',
"ORN" => 'A',
"ORQ" => 'R',
"OSE" => 'S',
"OTH" => 'T',
"OXX" => 'D',
"P1L" => 'C',
"P2Y" => 'P',
"PAQ" => 'Y',
"PAS" => 'D',
"PAT" => 'W',
"PAU" => 'A',
"PBB" => 'C',
"PBF" => 'F',
"PCA" => 'E',
"PCC" => 'P',
"PCS" => 'F',
"PEC" => 'C',
"PF5" => 'F',
"PFF" => 'F',
"PG1" => 'S',
"PG9" => 'G',
"PGY" => 'G',
"PH6" => 'P',
"PHA" => 'F',
"PHD" => 'D',
"PHE" => 'F',
"PHI" => 'F',
"PHL" => 'F',
"PHM" => 'F',
"PLE" => 'L',
"PM3" => 'F',
"POM" => 'P',
"PPN" => 'F',
"PR3" => 'C',
"PR9" => 'P',
"PRO" => 'P',
"PRS" => 'P',
"PSA" => 'F',
"PSH" => 'H',
"PTH" => 'Y',
"PTM" => 'Y',
"PTR" => 'Y',
"PVH" => 'H',
"PYA" => 'A',
"PYX" => 'C',
"QCS" => 'C',
"QMM" => 'Q',
"QPA" => 'C',
"QPH" => 'F',
"R1A" => 'C',
"R4K" => 'W',
"RE0" => 'W',
"RE3" => 'W',
"RPI" => 'R',
"RVX" => 'S',
"RZ4" => 'S',
"S1H" => 'S',
"S2C" => 'C',
"S2D" => 'A',
"S2P" => 'A',
"SAC" => 'S',
"SAH" => 'C',
"SAR" => 'G',
"SBL" => 'S',
"SCH" => 'C',
"SCS" => 'C',
"SCY" => 'C',
"SDP" => 'S',
"SE7" => 'A',
"SEB" => 'S',
"SEG" => 'A',
"SEL" => 'S',
"SEM" => 'S',
"SEN" => 'S',
"SEP" => 'S',
"SER" => 'S',
"SET" => 'S',
"SGB" => 'S',
"SHC" => 'C',
"SHP" => 'G',
"SHR" => 'K',
"SIB" => 'C',
"SLR" => 'P',
"SLZ" => 'K',
"SMC" => 'C',
"SME" => 'M',
"SMF" => 'F',
"SNC" => 'C',
"SNN" => 'N',
"SOC" => 'C',
"SOY" => 'S',
"SRZ" => 'S',
"STY" => 'Y',
"SUN" => 'S',
"SVA" => 'S',
"SVV" => 'S',
"SVW" => 'S',
"SVX" => 'S',
"SVY" => 'S',
"SVZ" => 'S',
"SYS" => 'C',
"T11" => 'F',
"TAV" => 'D',
"TBG" => 'V',
"TBM" => 'T',
"TCQ" => 'Y',
"TCR" => 'W',
"TDD" => 'L',
"TFQ" => 'F',
"TH6" => 'T',
"THC" => 'T',
"THR" => 'T',
"THZ" => 'R',
"TIH" => 'A',
"TMB" => 'T',
"TMD" => 'T',
"TNB" => 'C',
"TNR" => 'S',
"TOQ" => 'W',
"TPL" => 'W',
"TPO" => 'T',
"TPQ" => 'Y',
"TQI" => 'W',
"TQQ" => 'W',
"TRF" => 'W',
"TRG" => 'K',
"TRN" => 'W',
"TRO" => 'W',
"TRP" => 'W',
"TRQ" => 'W',
"TRW" => 'W',
"TRX" => 'W',
"TRY" => 'W',
"TTQ" => 'W',
"TTS" => 'Y',
"TXY" => 'Y',
"TY1" => 'Y',
"TY2" => 'Y',
"TY3" => 'Y',
"TY5" => 'Y',
"TYB" => 'Y',
"TYI" => 'Y',
"TYJ" => 'Y',
"TYN" => 'Y',
"TYO" => 'Y',
"TYQ" => 'Y',
"TYR" => 'Y',
"TYS" => 'Y',
"TYT" => 'Y',
"TYW" => 'Y',
"TYY" => 'Y',
"UMA" => 'A',
"VAD" => 'V',
"VAF" => 'V',
"VAL" => 'V',
"VB1" => 'K',
"WLU" => 'L',
"WPA" => 'F',
"WRP" => 'W',
"WVL" => 'V',
"X2W" => 'E',
"XCN" => 'C',
"XDT" => 'T',
"XPR" => 'P',
"XSN" => 'N',
"XX1" => 'K',
"YCM" => 'C',
"YOF" => 'Y',
"YTH" => 'T',
"Z01" => 'A',
"ZAL" => 'A',
"ZCL" => 'F',
"ZU0" => 'T',
"ZZJ" => 'A'
)


text:
# Used by fptas below
struct ScaleParams
    t::Vector{Int}
    ft::Vector{Float64}
    Ū::Int
    infty::Int

    function ScaleParams(mkt::VariedCostsMarket, ε::Float64)
        @assert 0 < ε < 1
        P = max(0,
            ceil(Int, log2(Int(mkt.m)^2 / (ε * sum(mkt.ft))))
        )
        t = mkt.t .* 2^P
        ft = mkt.f .* t
        Ū = ceil(Int, sum(ft))
        infty = sum(mkt.g) + 1

        return new(t, ft, Ū, infty)
    end
end

# Used by fptas below
@inbounds function G_recursor!(
    G_dict::Dict{Tuple{Int,Int},Int},
    j::Int,
    v::Int,
    mkt::VariedCostsMarket,
    sp::ScaleParams
)::Int
    return get(G_dict, (j, v)) do
        if v ≤ 0
            return 0
        elseif j == 0 || sp.t[j] < v || v ≥ sp.Ū
            return sp.infty
        else
            jmo = j - 1
            if mkt.f[j] < 1
                # Clamping prevents over/underflow: for any v<0 or v≥Ū the function
                # is trivially defined, so recording any more extreme number is meaningless
                v_minus_Δ = floor(Int, clamp((v - sp.ft[j]) / mkt.omf[j], -1, sp.Ū))

                push!(G_dict, (j, v) => min(
                    G_recursor!(G_dict, jmo, v, mkt, sp),
                    mkt.g[j] + G_recursor!(G_dict, jmo, v_minus_Δ, mkt, sp)
                ))
                return G_dict[(j, v)]
            else
                push!(G_dict, (j, v) => min(G_recursor!(G_dict, jmo, v, mkt, sp), mkt.g[j]))
                return G_dict[(j, v)]
            end
        end
    end
end

# Function I stole from Base.Sort that does midpoint using bitshift.
# Performance-optimized but safe only if `lo <= hi`.
midpoint(lo::T, hi::T) where {T<:Integer} = lo + ((hi - lo) >>> 0x01)

"""
    optimalportfolio_fptas(mkt, ε; verbose=false) -> X, v

Use the fully polynomial-time approximation scheme to produce a `1-ε`-optimal portfolio `X`,
with valuation `v`, for the [`VariedCostsMarket`](@ref) defined by `mkt`.
"""
function optimalportfolio_fptas(
    mkt::VariedCostsMarket,
    ε::Float64;
    verbose::Bool=false
)::Tuple{Vector{Int},Float64}
    sp = ScaleParams(mkt, ε)

    G_dict = Dict{Tuple{Int,Int},Int}()
    sizehint!(G_dict, mkt.m * mkt.m ÷ 2)

    # Binary search for max{ w :  G_recursor!(G_dict, mkt.m, w, mkt, sp ≤ H }
    # Should be able to use something like this with a proper by kw in searchsortedlast
    # vs = 0:sp.Ū
    # searchsortedlast(vs, by = ...)

    v = 0
    v_UB = sp.Ū
    @inbounds while v + 1 < v_UB
        # mid = (v + v_UB) ÷ 2
        mid = midpoint(v, v_UB)

        if G_recursor!(G_dict, mkt.m, mid, mkt, sp) > mkt.H
            v_UB = mid
        else
            v = mid
        end
    end

    X = Int[]
    @inbounds for j in reverse(1:mkt.m)
        # G_recursor!(G_dict, j, v, mkt, sp) < sp.infty &&
        if G_recursor!(G_dict, j, v, mkt, sp) < G_recursor!(G_dict, j - 1, v, mkt, sp)
            push!(X, j)
            v = floor(Int, clamp((v - sp.ft[j]) / mkt.omf[j], -1, sp.Ū))
        end
    end

    if verbose
        G_table = Array{Union{Missing,Int}}(missing, mkt.m, v_UB)
        for (j, v) in keys(G_dict)
            if 0 < j ≤ mkt.m && 0 < v ≤ v_UB
                G_table[j, v] = G_dict[(j, v)]
            end
        end

        issorted(mkt.perm) || @warn "t not sorted; table rows will differ from input"
        display(G_table)
    end

    return mkt.perm[X], valuation_nopermute_sorted(sort(X), mkt)
end

max_stars_repo_name: maxkapur/OptimalApplication

max_issues_repo_name: maxkapur/OptimalApplication

max_forks_repo_name: maxkapur/OptimalApplication


text:
#
# Julia. It supports UTF-8. (╯°□°）╯︵ ┻┻
#

using Plots
using StaticArrays
import Base.show

firstrun = true

struct FArrayDims
    fstride::Cptrdiff_t
    lower_bound::Cptrdiff_t
    upper_bound::Cptrdiff_t
end

struct FDataType
    len::Csize_t
    ver::Cint
    rank::Cchar
    type::Cchar
    attribute::Cshort
end

struct FArray1D
    base_addr::Ptr{Cvoid}
    offset::Csize_t
    dtype::FDataType
    span::Cptrdiff_t
    dim::FArrayDims

    FArray1D(vector::Array{Float64,1}) = new(
                pointer(vector),
                0,
                FDataType(8, 0, 1, 3, 0),
                length(vector)*8,
                FArrayDims(1, 1, length(vector))
        )
end

struct FArray2D
    base_addr::Ptr{Cvoid}
    offset::Csize_t
    dtype::FDataType
    span::Cptrdiff_t
    dim1::FArrayDims
    dim2::FArrayDims

    FArray2D(arr::Array{Float64,2}) = new(
                pointer(arr),
                0,
                FDataType(8, 0, 2, 3, 0),
                size(arr)[0]*size(arr)[1]*8,
                FArrayDims(1, 1, size(arr)[0]),
                FArrayDims(size(arr)[0], 1, size(arr)[1])
        )
end

function show(io::IO, farr::FArray1D)
    data = unsafe_wrap(Array{Float64,1}, Ptr{Float64}(farr.base_addr),
        (farr.dim.upper_bound - farr.dim.lower_bound + 1,))
    print(io, data)
end

mutable struct FieldCan
    Ath::Float64
    Aph::Float64
    hth::Float64
    hph::Float64
    Bmod::Float64

    dAth::SVector{3, Float64}
    dAph::SVector{3, Float64}
    dhth::SVector{3, Float64}
    dhph::SVector{3, Float64}
    dBmod::SVector{3, Float64}

    # second derivatives: drdr, drdhth, drdph, dthdth, dthdph, dphdph
    d2Ath::SVector{6, Float64}
    d2Aph::SVector{6, Float64}
    d2hth::SVector{6, Float64}
    d2hph::SVector{6, Float64}
    d2Bmod::SVector{6, Float64}

    H::Float64
    pth::Float64
    vpar::Float64

    dvpar::SVector{4, Float64}
    dH::SVector{4, Float64}
    dpth::SVector{4, Float64}

    # order of second derivatives:
    # d2dr2, d2drdth, d2drph, d2dth2, d2dthdph, d2dph2,
    # d2dpphdr, d2dpphdth, d2dpphdph, d2dpph2
    d2vpar::SVector{10, Float64}
    d2H::SVector{10, Float64}
    d2pth::SVector{10, Float64}

    mu::Float64
    ro0::Float64

    FieldCan() = new()
  end

mutable struct SymplecticIntegrator
    atol::Float64
    rtol::Float64

    z::SVector{4, Float64}
    pthold::Float64

    # Buffer for Lagrange polynomial interpolation
    kbuf::Int32
    kt::Int32
    k::Int32
    bufind::SVector{4, Int32}
    zbuf::SMatrix{4, 16, Float64, 64}
    coef::SMatrix{1, 4, Float64, 4}

    # Timestep and variables from z0
    ntau::Int32
    dt::Float64
    pabs::Float64

    # Integrator mode
    mode::Int32  # 1 = euler1, 2 = euler2, 3 = verlet

    SymplecticIntegrator() = new()
end

mutable struct MultistepIntegrator
  s::Int32
  alpha::SVector{64, Float64}
  beta::SVector{64, Float64}
  stages::SVector{42496, Int8}  # Static 64*664 array of SymplecticIntegrators
end

mutable struct NeoOrb
  fper::Float64
  dtau::Float64
  dtaumax::Float64
  v0::Float64
  n_e::Int32
  n_d::Int32
  integmode::Int32
  relerr::Float64

  NeoOrb() = new()
end

mutable struct CutDetector
  npl_half::Int32
  alam_prev::Float64
  par_inv::Float64
  iper::Int32
  itip::Int32
  kper::Int32

  orb_sten::SMatrix{6, 6, Float64, 36}
  coef::SMatrix{1, 6, Float64, 6}
  ipoi::SVector{6, Int32}

  norb::NeoOrb

#   fper::Float64
#   dtau::Float64
#   dtaumax::Float64
#   v0::Float64
#   n_e::Int32
#   n_d::Int32
#   integmode::Int32
#   relerr::Float64
#   firstrun::Bool

  CutDetector() = new()
end

norb = NeoOrb()
if !@isdefined(firstrun) | firstrun == true
      global firstrun, norb
      ccall((:__neo_orb_MOD_init_field, "lib/libneo_orb"),
      Cvoid, (Ref{NeoOrb}, Ref{Int32}, Ref{Int32}, Ref{Int32}, Ref{Int32}),
      norb, 5, 5, 3, 1)
      firstrun = false
end

Z_charge = 2
m_mass = 4
c=2.9979e10
e_charge=4.8032e-10
e_mass=9.1094e-28
p_mass=1.6726e-24
ev=1.6022e-12

bmod_ref = 5.0e4
Z_charge = 2
m_mass = 4
E_kin = 3.5e6
trace_time = 1.0e-2

v0 = sqrt(2.0*E_kin*ev/(m_mass*p_mass))    # alpha particle velocity, cm/s
rlarm=v0*m_mass*p_mass*c/(Z_charge*e_charge*bmod_ref) # reference gyroradius
tau=trace_time*v0

npoiper2 = 64

p_rmajor = cglobal((:__new_vmec_stuff_mod_MOD_rmajor, "lib/libneo_orb"), Float64)
rmajor = unsafe_load(p_rmajor)
println(rmajor)
rbig = rmajor*1.0e2
dtaumax = 2.0*pi*rbig/npoiper2
dtau = dtaumax
ccall((:__neo_orb_MOD_init_params, "lib/libneo_orb"), Cvoid,
      (Ref{NeoOrb}, Ref{Int32}, Ref{Int32}, Ref{Float64}, Ref{Float64},
       Ref{Float64}, Ref{Float64}),
      norb, Z_charge, m_mass, E_kin, dtau, dtaumax, 1e-8)

s = [0.5]
th = [0.0]
ph = [0.314]
lam = [0.22]

z = [0.5, 0.0, 0.314, 1.0, 0.22]
za = FArray1D(z)

ierr = 0

f = FieldCan()

## Multistage method
fac = 6.0

nstage = 2
stages = Array{SymplecticIntegrator}(undef, 2)

z[:] = [0.5, 0.0, 0.314, 1.0, 0.2]

# Explicit Euler
stages[1] = SymplecticIntegrator()
ccall((:__orbit_symplectic_MOD_orbit_sympl_init, "lib/libneo_orb"), Cvoid,
    (Ref{SymplecticIntegrator}, Ref{FieldCan}, Ref{FArray1D}, Ref{Float64},
     Ref{Float64}, Ref{Float64}, Ref{Int32}),
     stages[1], f, za, fac*norb.dtau/2.0, fac*norb.dtaumax/2.0, 1e-12, 2)

# ImplicitEuler
stages[2] = SymplecticIntegrator()
ccall((:__orbit_symplectic_MOD_orbit_sympl_init, "lib/libneo_orb"), Cvoid,
    (Ref{SymplecticIntegrator}, Ref{FieldCan}, Ref{FArray1D}, Ref{Float64},
     Ref{Float64}, Ref{Float64}, Ref{Int32}),
     stages[2], f, za, fac*norb.dtau/2.0, fac*norb.dtaumax/2.0, 1e-12, 2)

# Integrate 10000 timesteps

ntimstep = 10000
zpl = zeros(ntimstep, 5)

@time begin
    for i = 1:ntimstep
        ccall((:__neo_orb_MOD_timestep_sympl_z, "lib/libneo_orb"), Cvoid,
              (Ref{NeoOrb}, Ref{SymplecticIntegrator}, Ref{FieldCan}, Ref{FArray1D}, Ref{Int32}),
              norb, stages[1], f, za, ierr)
        stages[2].z = stages[1].z
        ccall((:__neo_orb_MOD_timestep_sympl_z, "lib/libneo_orb"), Cvoid,
              (Ref{NeoOrb}, Ref{SymplecticIntegrator}, Ref{FieldCan}, Ref{FArray1D}, Ref{Int32}),
              norb, stages[2], f, za, ierr)
        stages[1].z = stages[2].z
        zpl[i,:] = z
    end
end

plot(zpl[:,1].*cos.(zpl[:,2]), zpl[:,1].*sin.(zpl[:,2]))

## Cut detection
# println(za)
# println(si)

# ncut = 1000
# cutter = CutDetector()
# var_tip = zeros(6)
# var_tipa = FArray1D(var_tip)
# cut_type = 0

#println(cutter)

# ccall((:__cut_detector_MOD_init, "lib/libneo_orb"), Cvoid,
#        (Ref{CutDetector}, Ref{NeoOrb}, Ref{FArray1D}),
#        cutter, norb, za)

# println(cutter)

# @time begin
#         for i = 1:ncut
#                 ccall((:__cut_detector_MOD_trace_to_cut, "lib/libneo_orb"), Cvoid,
#                       (Ref{CutDetector}, Ref{FArray1D}, Ref{FArray1D}, Ref{Int32}, Ref{Int32}),
#                       cutter, za, var_tipa, cut_type, ierr)
#         end
# end

#%%

max_stars_repo_name: landreman/SIMPLE

max_issues_repo_name: landreman/SIMPLE

max_forks_repo_name: landreman/SIMPLE


text:
function getPOP2VerticalResolution()

# This z resolution is get from
    return - [
    0, 1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000, 9000, 10000,
    11000, 12000, 13000, 14000, 15000, 16000, 17019.681640625,
    18076.12890625, 19182.125, 20349.931640625, 21592.34375, 22923.3125,
    24358.453125, 25915.580078125, 27615.259765625, 29481.470703125,
    31542.373046875, 33831.2265625, 36387.47265625, 39258.046875,
    42498.88671875, 46176.65625, 50370.6875, 55174.91015625, 60699.66796875,
    67072.859375, 74439.8046875, 82960.6953125, 92804.3515625,
    104136.8203125, 117104.015625, 131809.359375, 148290.078125,
    166499.203125, 186301.4375, 207487.390625, 229803.90625, 252990.40625,
    276809.84375, 301067.0625, 325613.84375, 350344.875, 375189.1875,
    400101.15625, 425052.46875, 450026.0625, 475012, 500004.6875,
    525000.9375, 549999.0625
] / 100.0;

end

function getHOOMResolution(
    res
)

    local zs

    res = Symbol(res)

    if res == :POP2

        zs = getPOP2VerticalResolution()[1:41]

    elseif res == :Standard

        zs = collect(Float64, 0:-25:-500)[1:end-1]
        zs = vcat(zs, collect(Float64, -500:-50:-1000))
        if length(zs) != 31
            println(zs)
            throw(ErrorException("Not expected zs. Please check."))
        end

    else

        throw(ErrorException("Unknown resolution `" * string(res) * "`"))

    end

    return zs

end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
module WAVI

#Useful packages
using LinearAlgebra, SparseArrays, LinearMaps, Parameters,
      IterativeSolvers, Interpolations, BenchmarkTools, Reexport,
      NetCDF, JLD2, Setfield, MAT, ImageFiltering

#Import functions so they can be modified in this module.
import LinearAlgebra: ldiv!
import SparseArrays: spdiagm, spdiagm_internal, dimlub
import Setfield: @set

#This module will export these functions and types, allowing basic use of the model.
export
    #Structures
    Model, Params, TimesteppingParams, Grid, SolverParams, InitialConditions, OutputParams, Simulation,

    #Simulation controls
    update_state!, timestep!, run_simulation!,

    #Melt rates
    PlumeEmulator, BinfileMeltRate, UniformMeltRate, MISMIPMeltRateOne, PICO, QuadraticMeltRate,

    #Post-processing controls
    volume_above_floatation, height_above_floatation

#Reexport Modules useful for users of the WAVI module
@reexport using JLD2
@reexport using Setfield

#Abstract types
abstract type AbstractGrid{T <: Real, N <: Integer} end
abstract type AbstractMeltRate end
abstract type AbstractModel{T <: Real, N <: Integer, M <: AbstractMeltRate} end
abstract type AbstractPreconditioner{T <: Real, N <: Integer} end
#abstract type AbstractSimulation{T,N,R,A,W} end

#Type alias, just for abreviation
const KronType{T,N} = LinearMaps.KroneckerMap{T,Tuple{LinearMaps.WrappedMap{T,SparseMatrixCSC{T,N}},
                        LinearMaps.WrappedMap{T,SparseMatrixCSC{T,N}}}} where {T <: Real, N <: Integer}

#Concrete types

##################################################################################
#include all of the code
include("OutputParams/OutputParams.jl")
include("Grid.jl")
include("Params.jl")
include("SolverParams.jl")
include("TimesteppingParams.jl")
include("Clock.jl")
include("InitialConditions.jl")
include("Wavelets/Wavelets.jl")
include("Fields/Fields.jl")
include("Models/Model.jl")
include("MeltRate/MeltRate.jl")
include("Simulations/Simulation.jl")
include("utilities.jl")

"""
    spdiagm(m::Integer, n::Integer, kv::Pair{<:Integer,<:AbstractVector}...)

Method from Julia V1.4.2 to create non-square sparse matrix from diagonals. Included for backward compatibility.
"""
spdiagm(m::Integer, n::Integer, kv::Pair{<:Integer,<:AbstractVector}...) = _spdiagm((Int(m),Int(n)), kv...)
function _spdiagm(size, kv::Pair{<:Integer,<:AbstractVector}...)
    I, J, V = spdiagm_internal(kv...)
    mmax, nmax = dimlub(I), dimlub(J)
    mnmax = max(mmax, nmax)
    m, n = something(size, (mnmax,mnmax))
    (m ≥ mmax && n ≥ nmax) || throw(DimensionMismatch("invalid size=$size"))
    return sparse(I, J, V, m, n)
end
end


text:
# This file is part of Dieter: Utility functions

"Given a Dict of keys => filenames, check at files exist"
function check_files_exist(fileDict::Dict{String,String})
    for (k,v) in fileDict
        if !(isfile(fileDict[k]) || isdir(fileDict[k]))
            @error "File \"$(fileDict[k])\" not found."
            return false
        end
    end

    return true
end

"Parse a given data file and return the desired data in a `DataFrame`."
function parse_file(file::String; dataname::String="")
    filetype=split(lowercase(file), '.')[end]
    if filetype == "csv"
        df = CSV.read(file, DataFrame)
    elseif filetype == "sql" ## TODO
        # We expect dataname to be the path to a SQLite database file
        if !(isfile(dataname))
            @error "The given dataname \"$dataname\" is not a valid file."
        end
        try
            db = SQLite.DB(dataname)
        catch
            @error "Opening of SQLite database file \"$dataname\" unsuccessful."
        end
        db = SQLite.DB(dataname)
        queryString = read(file, String)
        df = DBInterface.execute(db,queryString) |> DataFrame
    else
        error("Unrecognised filetype (as parsed from the file extension).")
    end

    return df
end

tuple2_filter(func,A,B) = filter(x->func(x[1],x[2]), [(a,b) for a in A for b in B])
# tuple3_filter(func,A,B,C) = filter(x->func(x[1],x[2],x[3]), [(a,b,c) for a in A for b in B for c in C])
# tupleN_filter(func,X...) = filter(x->func(x...), "CartesianProduct"(X...))

"""
Using a DataFrame, `create_relation` returns a function that can be used to determine
a 2-valued relation, that is, a set-to-set correspondence.
DataFrame column-names (Symbols) are given that define the set tuples,
and a third Indicator column-name is used to indicate set membership.
The return type is a function that can be used with two "sets" A and B in
`tuple2_filter(func,A,B) = filter(x->func(x[1],x[2]), [(a,b) for a in A for b in B])`
"""
function create_relation(df::DataFrame,First::Symbol,Second::Symbol,
                Indicator::Symbol) # ; criterion=testfunc_returns_bool

    if !(String(First) in names(df))
        @error "The symbol $First must be a valid name of the DataFrame."
    elseif !(String(Second) in names(df))
        @error "The symbol $Second must be a valid name of the DataFrame."
    elseif !(String(Indicator) in names(df))
        @error "The symbol $Indicator must be a valid name of the DataFrame."
    end

    # FirstType = eltype(df[!,:($First)])
    # SecondType   = eltype(df[!,:($Second)])

    func(x,y) =
        begin
            r = @subset(df, :($First) .== x, :($Second) .== y)
            return sum(r[:,:($Indicator)]) > 0  ## This could be some other criteron for a relation existing:
            #  TODO allow user to pass the criterion function as an argument
        end

    return func
end

"Update a data Dict, merging the data if present, adding data if not present"
function update_dict!(dict::Dict{Symbol,Dict}, key::Symbol, val::Dict)
    # if !isa(val,Dict)
        # @warn "This function expects the value $val to be of type Dict"
    # end
    if haskey(dict, key)
        try
            merge!(dict[key], val) # Error if dict[key] & val are not Dict types
            # merge!(dict, Dict(key => val))
        catch e
            @warn "Merge of Dict failed for the key $(key)."
            throw(e)
        end
    else
        push!(dict, key => val)
    end
end

"""
Return a dictionary of `Dict`s from a `DataFrame` with column-names as keys,
and each `Dict` has a key-value pair of id_cols => data-value.
Column symbols listed in `skip_cols` are skipped.

# Examples

```julia
julia> df = DataFrame(x=1:4, y=11:14, z=2:5)
4×3 DataFrame
│ Row │ x     │ y     │ z     │
│     │ Int64 │ Int64 │ Int64 │
├─────┼───────┼───────┼───────┤
│ 1   │ 1     │ 11    │ 2     │
│ 2   │ 2     │ 12    │ 3     │
│ 3   │ 3     │ 13    │ 4     │
│ 4   │ 4     │ 14    │ 5     │

julia> map_idcol(df, 1, skip_cols=[:z])
Dict{Symbol,Dict} with 1 entry:
  :y => Dict(4=>14,2=>12,3=>13,1=>11)
```
"""
function map_idcol(df::DataFrame, id_cols::Int; skip_cols=Symbol[])

    push!(skip_cols, propertynames(df)[id_cols])

    if id_cols > 1
        ids = collect(zip([col for col in eachcol(df[:,1:id_cols])]...))
    else
        ids = df[!, 1]
    end

    non_id = df[:, id_cols+1:end]

    dict = Dict{Symbol, Dict}()

    for col in pairs(eachcol(non_id))
        if !(col[1] in skip_cols)
            vals = Dict(collect(zip(ids,col[2])))
            dict[col[1]] = vals
        end
    end
    return dict
end

map_idcol(df::DataFrame; skip_cols=Symbol[]) = map_idcol(df, 1, skip_cols=skip_cols)
map_idcol(df::DataFrame, id_col::Symbol; skip_cols=Symbol[]) = map_idcol(df, [id_col], skip_cols=skip_cols)

function map_idcol(df::DataFrame, id_cols::Array{Symbol,1}; skip_cols=Symbol[])

    for sym in skip_cols
        if sym in id_cols
            @warn "Skipped columm $sym is a member of the given id key columns.
                    It will be ignored."
        elseif !(sym in propertynames(df))
            @warn "Skipped columm $sym is a not column name in the DataFrame.
                    It will be ignored."
        end
    end

    push!(skip_cols, id_cols...)

    if length(id_cols) == 1 ## The Dict id will be a singleton
        Keys = df[!,:($(id_cols[1]))]   # Returns a Vector
    else # Make the Dict keys into Tuples
        Keys = [Tuple(row[id_cols]) for row in eachrow(df)]
    end

    non_id = setdiff(propertynames(df),skip_cols)  # Remove skipped columms

    dict = Dict{Symbol, Dict}()

    for col in pairs(eachcol(df))
        if col[1] in non_id
            dict[col[1]] = Dict(collect(zip(Keys,col[2])))
        end
    end

    return dict
end

"""
Create a `Dict` from a `DataFrame` with the `DataFrame` column-names as keys and columns as values.
"""
map_dfheader_to_col(df::DataFrame) = Dict(String(col[1]) => disallowmissing(col[2]) for col in pairs(eachcol(df)))

get_variable_names(m::JuMP.Model) = [v.name for (k,v) in m.varData]  # varData requires instantiation for JuMP Model

"""
Returns a `DataFrame` with the values of the variables from the JuMP container `var`.
The column names of the `DataFrame` can be specified for the indexing columns in `dim_names`,
and the name of the data value column by a Symbol `value_col` e.g. :Value
"""
function convert_jump_container_to_df(var::JuMP.Containers.DenseAxisArray;
    dim_names::Vector{Symbol}=Vector{Symbol}(),
    value_col::Symbol=:Value)

    if isempty(var)
        return DataFrame()
    end

    if length(dim_names) == 0
        dim_names = [Symbol("dim$i") for i in 1:length(var.axes)]
    end

    if length(dim_names) != length(var.axes)
        throw(ArgumentError("Length of given name list does not fit the number of variable dimensions"))
    end

    tup_dim = (dim_names...,)

    # With a product over all axis sets of size M, form an Mx1 Array of all indices to the JuMP container `var`
    ind = reshape([collect(k[i] for i in 1:length(dim_names)) for k in Base.Iterators.product(var.axes...)],:,1)

    var_val  = value.(var)

    df = DataFrame([merge(NamedTuple{tup_dim}(ind[i]), NamedTuple{(value_col,)}(var_val[(ind[i]...,)...])) for i in 1:length(ind)])

    return df
end

function convert_jump_container_to_df(var::JuMP.Containers.SparseAxisArray;
    dim_names::Vector{Symbol}=Vector{Symbol}(),
    value_col::Symbol=:Value)

    if isempty(var)
        return DataFrame()
    end

    var_data = (value.(var)).data

    num_dim = length(collect(keys(var_data))[1])  # Roundabout way of finding the dimension of the var array

    if length(dim_names) == 0
        dim_names = [Symbol("dim$i") for i in 1:num_dim]
    end

    if length(dim_names) != num_dim
        throw(ArgumentError("Length of given name list does not fit the number of variable dimensions"))
    end

    tup_dim = (dim_names..., value_col)

    df = DataFrame([NamedTuple{tup_dim}(([k[i] for i in 1:length(dim_names)]..., v)) for (k,v) in var_data])

    return df
end

function convert_jump_container_to_df(var::Array{VariableRef};
    dim_names::Vector{Symbol}=Vector{Symbol}(),
    value_col::Symbol=:Value)

    if isempty(var)
        return DataFrame()
    end

    num_dim = size(var,2)

    if length(dim_names) == 0
        dim_names = [Symbol("dim$i") for i in 1:num_dim]
    end

    if length(dim_names) != num_dim
        throw(ArgumentError("Length of given name list does not fit the number of variable dimensions"))
    end

    var_val = value.(var)

    df = DataFrame([var_val[:,i] for i in 1:num_dim], dim_names)

    return df
end

function convert_jump_container_to_df(var::VariableRef;
    dim_names::Vector{Symbol}=Vector{Symbol}(),
    value_col::Symbol=:Value)

    return convert_jump_container_to_df([var],dim_names=dim_names,value_col=value_col)
end

"""
Function to convert an Array{Any,2} with a header row of column name strings to a DataFrame
 - Thanks: https://stackoverflow.com/questions/25894634/dataframe-from-array-with-header
"""
table2df(mat) = DataFrame(mat[2:end,:], Symbol.(mat[1,:]))

"""
    tableToDict(tab::Array{Any,2}; keycols::Array{Int64,1}=[1])

Create a dictionary from an array in table form, where the table `tab` has:
 - key-index or indices given in the `keycol` array (integer column position), and
 - the column names for data appear in the first row.
"""
function tableToDict(tab::Array{<:Any,2}; keycols::Array{Int64,1}=[1])
    DataDict = Dict{Union{String,Tuple{String,Vararg{String}}},Any}() # initialise the dictionary that the function returns

    KeysName    = tab[1,keycols]      # the column label on the key column
    ColumnNames = tab[1,:]            # all the column names
    KeyArray    = tab[2:end,keycols]  # the keys for each row
    TableData   = tab[2:end,:]        # the data in the input array

    NumberRows    = size(TableData)[1]
    NumberColumns = size(TableData)[2]
    # Index columns with the key column removed:
    DataColumnIndices = setdiff(1:NumberColumns,keycols)

    if length(keycols) == 1 # # Make the DataDict key a singleton
        Keys = KeyArray
    else # Make the DataDict keys into Tuples
        Keys = [Tuple(KeyArray[i,:]) for i in 1:NumberRows]
    end

    for (k, key_value) in enumerate(Keys)
        DataDict[key_value] = Dict(ColumnNames[col] => TableData[k,col]
                                        for col in DataColumnIndices)
    end
    return DataDict
end

function SQLqueryToDict(sqlquery::SQLite.Query; keycols::Array{Int64,1}=[1])
    df = DataFrame(sqlquery)
    A = convert(Matrix, df)
    # or:
    # A = Tables.matrix(tb)s
    # Convert DataFrame to a table with column headers
    tableArray = [reshape(string.(names(df)),1,size(A)[2]); A]

    DataDict = tableToDict(tableArray; keycols=keycols)

    return DataDict
end

"""
`split_df_tuple` is a function that transforms a DataFrame with a Tuple{N,M...} type column.
It outputs an augmented DataFrame with the tuples split into two columns of type N, M... resp.

The function splits the column `InputCol` of type `Tuple` into individual columns named with `OutputCols` vector of Symbols.
Note that the `Symbol`s in `OutputCols` need to be of the name length and order to correspond to the splitting `Tuple`s.

 ### Arguments

 * `df` : a `DataFrame`
 * `InputCol` : a `Symbol` for a column of type `Tuple{Types...}`
 * `OutputCols`: a `Vector` of `Symbol`s for new column names,

 ### Returns

 * `::DataFrame`

 ### Examples

 ```
 julia> using DataFramesMeta, DataFrames

 julia> df = DataFrame(A_N = [("a",1),("b",2)], V = [8,9]);

 julia> split_df_tuple(df, :A_N, [:A, :N])
2×4 DataFrame
│ Row │ A_N      │ V     │ A      │ N     │
│     │ Tuple…   │ Int64 │ String │ Int64 │
├─────┼──────────┼───────┼────────┼───────┤
│ 1   │ ("a", 1) │ 8     │ a      │ 1     │
│ 2   │ ("b", 2) │ 9     │ b      │ 2     │
 ```

"""
function split_df_tuple(df::DataFrame, InputCol::Symbol,OutputCols::Vector{Symbol})
      return split_df_tuple!(copy(df),InputCol,OutputCols)
end

function split_df_tuple!(df::DataFrame, InputCol::Symbol,OutputCols::Vector{Symbol})
      for (i, newcol) in enumerate(OutputCols)
            df[!,newcol] = [tup[i] for tup in df[!,InputCol]]
      end
      return df
end

"Check if a `Dict`ionary has all missing values; if so, return `true`"
function isDictAllMissing(dict::Dict{T,S} where {T, S <: Union{Missing, Any}})
    return filter(x -> !ismissing(x.second), dict) |> isempty
end

# Filter a Dictionary for a matching string in the key.
dkeymatch(d::Dict,s::Regex) = filter(x -> (!ismissing(x.second) && occursin(s,x.first)),d)
# Filter a Dictionary for a matching string in the `pos`-th position in a Tuple key.
dkeymatch(d::Dict,s::Regex,pos::Int) = filter(x -> (!ismissing(x.second) && occursin(s,x.first[pos])),d)
# Filter a Dictionary for a matching string in the value.
dvalmatch(d::Dict,s::Regex) = filter(x -> (!ismissing(x.second) && occursin(s,x.second)),d)