text:
using BenchmarkTools
using TimeZones
using Dates
using HTTP
using JSON

include("benchmark_utils.jl")

const CI_GITHUB_REPO = "nep-pack/julia-ci"
const SOURCE_GITHUB_REPO = "nep-pack/NonlinearEigenproblems.jl"
const NR_COMMITS_TO_SHOW = 16

abstract type BenchmarkFile end

struct GitHubBenchmarkFile <: BenchmarkFile
    build_nr
    url
end

struct LocalBenchmarkFile <: BenchmarkFile
    path
end

build_nr(file_name) = parse(Int, split(file_name, "-")[end-1])

GitHubBenchmarkFile(entry) =
    GitHubBenchmarkFile(build_nr(entry["name"]), entry["download_url"])

get_benchmark_json(file::GitHubBenchmarkFile) = (println("Loading $file"); String(HTTP.request("GET", file.url).body))
get_benchmark_json(file::LocalBenchmarkFile) = read(file.path, String)

import Base.show

show(io::IO, f::GitHubBenchmarkFile) = print(io, f.url)
show(io::IO, f::LocalBenchmarkFile) = print(io, f.path)

function get_github_files(branch, count)
    req = HTTP.request("GET", "https://api.github.com/repos/$CI_GITHUB_REPO/contents", [("User-Agent", "julia-ci")])
    json = JSON.parse(String(req.body))
    files = [GitHubBenchmarkFile(f) for f in json if f["type"] == "file" && startswith(f["name"], "benchmark-$branch-")]
    sort!(files, by = x -> x.build_nr)
    return files[max(1, end-count+1):end]
end

function create_report(report_file, branch, files)
    file_name = coalesce(report_file, tempname() * ".html")
    files = ismissing(branch) ? LocalBenchmarkFile.(files) : get_github_files(branch, NR_COMMITS_TO_SHOW)

    open(file_name, "w") do io
        write_report(io, files, ismissing(branch))
    end

    # if no report name was specified, open report in browser
    ismissing(report_file) && open_file(file_name)
end

function get_tooltip(run)
    runstr = run["time"]
    haskey(run, "runtime") && (runstr *= @sprintf("; runtime=%.0fs", run["runtime"]))
    haskey(run, "config") && (runstr *= "; $(run["config"])")
    haskey(run, "julia-version") && (runstr *= "; Julia $(run["julia-version"])")

    ghz, arch, cpus, mem = match(r"([\d\.]+)GHz \(([^\)]+)\), (\d+) threads, (\d+) GB", run["cpu"]).captures

    commit_id_and_branch, git_timestamp, commit_msg = split(run["git"], ", ", limit=3)
    length(commit_msg) > 100 && (commit_msg = commit_msg[1:97] * "...")

    return "Run: $runstr
" *
        "CPU: $cpus × $(ghz)GHz ($arch), $mem GB memory
" *
        "Git: $commit_id_and_branch, $git_timestamp
" *
        "$commit_msg"
end

function write_report(io, files, sort_benchmarks_by_date)
    print_report_header(io)

    # write table header
    println(io, "
")
    println(io, "
")
    println(io, "Test	")

    println(io, "	")
    println(io, "Baseline	")

    runs = map(f -> load_benchmark_from_string(get_benchmark_json(f)), files)
    sort_benchmarks_by_date && sort!(runs, by = r -> r["time"])
    bruns = map(r -> r["benchmark"][1], runs)
    println("Loaded $(length(runs)) runs")

    for run in runs
        println(io, "	")
        commit_id = split(run["git"], " ")[1]
        id = "$commit_id";
        println(io, "
$(get_tooltip(run))$id
	")
    end

    println(io, "
")

    root = get_testset_hierarchy(bruns[end])

    # write table itself
    iterate_hierarchy(root) do test
        println(io, "
")
        spaces = findfirst(!isequal(' '), test.indented_name)
        indented_name = " "^((spaces-1)*2) * test.indented_name[spaces:end]
        println(io, "$indented_name	")

        key = test.trial[1]

        min_time = mapreduce(run -> haskey(run, key) ? run[key].time : Inf, min, bruns)
        min_memory = mapreduce(run -> haskey(run, key) ? run[key].memory : Inf, min, bruns)

        println(io, "	")
        println(io, "$(prettytime(min_time))	")
        println(io, "$(prettymemory(min_memory))	")

        for run in bruns
            println(io, "	")

            if haskey(run, key)
                b = run[key]

                time_increase_pct = 100 * (b.time / min_time - 1)
                time_str = prettypercent(time_increase_pct)
                time_class = time_increase_pct >= 10 ? "bad" : time_increase_pct >= 5 ? "warn" : "ok"

                memory_increase_pct = 100 * (b.memory / min_memory - 1)
                memory_str = prettypercent(memory_increase_pct)
                memory_class = memory_increase_pct >= 5 ? "bad" : memory_increase_pct >= 1 ? "warn" : "ok"
            else
                time_str = prettypercent(NaN)
                memory_str = prettypercent(NaN)
                time_class = "missing"
                memory_class = "missing"
            end

            println(io, "$time_str	")
            println(io, "$memory_str	")
        end

        println(io, "
")
    end

    println(io, "
")

    print_report_footer(io)
end

function open_file(filename)
    if Sys.isapple()
        run(`open $filename`)
    elseif Sys.islinux() || Sys.isbsd()
        run(`xdg-open $filename`)
    elseif Sys.iswindows()
        run(`$(ENV["COMSPEC"]) /c start $filename`)
    else
        @warn "Can not open $filename on OS $(Sys.KERNEL) ($(Sys.MACHINE))"
    end
end

function print_report_header(io)
    println(io, """
<!DOCTYPE html>

<meta charset='utf-8'>
<title>Regression Analysis</title>

Regression Analysis

""")
end

function print_report_footer(io)
    println(io, """

""")
end

function parse_args(args)
    report_file = branch = missing
    while !isempty(args)
        if args[1] == "-o"
            report_file = args[2]
            args = args[3:end]
        elseif args[1] == "-b"
            branch = args[2]
            args = args[3:end]
        else break end
    end
    return (report_file, branch, args)
end

if length(ARGS) < 1
    println("Arguments required:\n\n$(basename(@__FILE__)) [-o html_output] [-b branch] [json_file ...]")
else
    create_report(parse_args(ARGS)...)
end


text:
using DataFrames
using CSV
using Statistics

filenames = [

    ### Put result table files here for reviewing the results, e.g.

    # LD 10
    # "Experiments/results/ld10_2020-02-19_003851.csv"

    # LD 3
    # "Experiments/results/ld3_2020-02-19_050959.csv"

    # CLD 10
    # "Experiments/results/cld10_2020-02-19_015442.csv"

    # CLD 3
    # "Experiments/results/cld3_2020-02-19_030803.csv"

    ]

data = vcat([CSV.read(fname) for fname in filenames]...)

r = by(data, :solver, N = :reward => length, mean = :reward => mean, std = :reward => r->std(r)/sqrt(length(r)))
sort!(r, :mean, rev=true)
display(filenames)
@show r

ext: jl


text:
__precompile__()

"""
    Phylo package

The `Phylo` package provides some simple phylogenetics types (e.g. NamedTree)
to interface to the `Diversity` package for measuring phylogenetic diversity.
It also provides an interface to `R` for copying trees to and from that
language and can read newick and nexus tree files (including `TreeSet`s that
contain multiple trees).

Finally it also provides a standard abstract interface to
phylogenetic trees, by defining `AbstractNode`, `AbstractBranch` and
`AbstractTree` supertypes, and methods to interface to them. It also
provides (through the `Phylo.API` submodule) methods to
(re)define to write your own phylogenetic type in a way that will
interact cleanly with other phylogenetic packages.
"""
module Phylo

import Base: Pair, Tuple, show, eltype, length, getindex
abstract type Rootedness end
struct Unrooted <: Rootedness end
abstract type Rooted <: Rootedness end
struct OneRoot <: Rooted end
struct ManyRoots <: Rooted end
export Unrooted, OneRoot, ManyRoots

abstract type TreeType end
struct OneTree <: TreeType end
struct ManyTrees <: TreeType end
export OneTree, ManyTrees

abstract type AbstractNode{RootType <: Rootedness, NodeLabel} end
abstract type AbstractBranch{RootType <: Rootedness, NodeLabel} end

abstract type AbstractTree{TT <: TreeType, RT <: Rootedness, NL,
                           N <: AbstractNode{RT, NL},
                           B <: AbstractBranch{RT, NL}} end
export AbstractTree

@enum TraversalOrder anyorder preorder inorder postorder breadthfirst
export TraversalOrder, anyorder, preorder, inorder, postorder, breadthfirst

"""
    Phylo.API submodule

The `Phylo.API` submodule should be `import`ed if you want to
create a new phylogeny, node or branch subtype. Otherwise it can be
ignored.
"""
module API
include("API.jl")
# AbstractTree methods
export _ntrees, _gettrees, _nroots, _getroots, _getroot
export _treenametype, _gettreenames, _gettree, _gettreename
export _createbranch!, _deletebranch!, _createnode!, _deletenode!
export _getnodenames, _hasnode, _getnode, _getnodes
export _getbranchnames, _getbranchname, _hasbranch, _getbranch, _getbranches
export _hasrootheight, _getrootheight, _setrootheight!, _clearrootheight!
export _getleafinfo, _setleafinfo!, _leafinfotype, _gettreeinfo
export _getnodedata, _setnodedata!, _nodedatatype
export _getbranchdata, _setbranchdata!, _branchdatatype
export _hasheight, _getheight, _setheight!
export _hasparent, _getparent, _getancestors
export _haschildren, _getchildren, _getdescendants
export _validate!, _traversal, _branchdims
export _getleafnames, _getleaves, _resetleaves!, _nleaves, _nnodes, _nbranches
export HoldsNodeData, MatchTreeNameType

# AbstractNode methods
export _isleaf, _isroot, _isinternal, _isunattached
export _indegree, _hasinboundspace, _outdegree, _hasoutboundspace, _degree
export _hasinbound, _getinbound, _addinbound!, _removeinbound!
export _getoutbounds, _addoutbound!, _removeoutbound!
export _getconnections, _addconnection!, _removeconnection!
export MatchNodeType, MatchNodeTypes, PreferNodeObjects, _prefernodeobjects

# AbstractBranch methods
export _src, _dst, _getlength
export MatchBranchType, PreferBranchObjects, _preferbranchobjects
export MatchBranchNodeType

# Label names
export _newnodelabel, _newbranchlabel

end

include("Interface.jl")
# AbstractTree methods
export ntrees, gettrees, nroots, getroots, getroot
export treenametype, gettreenames, gettreename #, getonetree #unimplemented
export roottype, nodetype, nodedatatype, nodenametype
export branchtype, branchdatatype, branchnametype
export createbranch!, deletebranch!
export createnode!, createnodes!, deletenode!
export getnodenames, getnodename, hasnode, getnode, getnodes
export getbranchnames, getbranchname, hasbranch, getbranch, getbranches
export hasrootheight, getrootheight, setrootheight!
export getparent, getancestors #, hasparent # unimplemented
export getchildren, getdescendants #, haschildren # unimplemented
export validate!, traversal, branchdims

@deprecate addnode! createnode!
@deprecate addnodes! createnodes!
@deprecate addbranch! createbranch!
@deprecate(branch!(tree, source, length = missing;
                   destination = _newnodelabel(tree),
                   branchname = _newbranchlabel(tree)),
           createbranch!(tree, source, createnode!(tree, destination),
                         length; name = branchname))

# AbstractTree / AbstractNode methods
export isleaf, isroot, isinternal, isunattached
export indegree, outdegree, hasinbound, getinbound, getoutbounds
export hasoutboundspace, hasinboundspace
export getleafnames, getleaves, nleaves, nnodes, nbranches
export getleafinfo, setleafinfo!, leafinfotype
export getnodedata, setnodedata!
export hasheight, getheight, setheight!

# AbstractTree / AbstractBranch methods
export src, dst, getlength
export hasrootheight, getrootheight, setrootheight! #, clearrootheight! #unimplemented
export getbranchdata, setbranchdata!
# export getrootdistance # unimplemented

@deprecate getnoderecord getnodedata
@deprecate setnoderecord! setnodedata!
@deprecate getbranchrecord getbranchdata
@deprecate setbranchrecord! setbranchdata!
@deprecate getbranchinfo getbranchdata
@deprecate setbranchinfo! setbranchdata!

include("Branch.jl")
export Branch

include("Node.jl")
export BinaryNode, Node

include("Tree.jl")
export BinaryTree, NamedBinaryTree, NamedTree
export PolytomousTree, NamedPolytomousTree

include("LinkTree.jl")
export LinkBranch, LinkNode, LinkTree
export RootedTree, ManyRootTree, UnrootedTree

include("routes.jl")
export branchhistory, branchfuture, branchroute
export nodehistory, nodefuture, noderoute

# Iterator methods expanded
include("Iterators.jl")
export nodeiter, nodefilter, nodenameiter, nodenamefilter,
    branchiter, branchfilter, branchnameiter, branchnamefilter
@deprecate treeiter gettrees
@deprecate treenameiter gettreenames

# A set of multiple trees
include("TreeSet.jl")
export TreeSet, gettreeinfo

# Random tree generator
include("rand.jl")
export Nonultrametric, Ultrametric
export BrownianTrait, DiscreteTrait, SymmetricDiscreteTrait

# Read Newick Tree
include("newick.jl")
export parsenewick, parsenexus

# Display methods expanded
include("show.jl")

# Method for trimming trees
include("trim.jl")
export droptips!, keeptips!

# Plot recipes
include("plot.jl")
export map_depthfirst

# Metrics from the tree
include("metrics.jl")
export mrca, nodeheights
export distance, distances, heighttoroot, heightstoroot

# Path into package
path(path...; dir::String = "test") = joinpath(@__DIR__, "..", dir, path...)

using Requires
function __init__()
    @require RCall="6f49c342-dc21-5d91-9882-a32aef131414" begin
        println("Creating Phylo RCall interface...")
        include("rcall.jl")
    end
end

end # module


text:
# ======================================
#    Tables
# ======================================

#"""
#AbstractTable{Index,Key} represents tables with a given set of indices (a
#FieldIndex) and a key (either one of the fields in FieldIndex, or else
#DefaultKey(), which represents the intrinsic row number of a table)
#"""
#abstract AbstractTable{Index,Key}
#@inline index{Index,Key}(::AbstractTable{Index,Key}) = Index
#@inline eltypes{Index,Key}(::AbstractTable{Index,Key}) = eltypes(Index)
#@inline Base.names{Index,Key}(::AbstractTable{Index,Key}) = names(Index)
#@inline key{Index,Key}(::AbstractTable{Index,Key}) = Key
#@inline Base.keytype{Index,Key}(::AbstractTable{Index,Key}) = eltype(Key)
#@inline keyname{Index,Key}(::AbstractTable{Index,Key}) = name(Key)

"""
A table stores columns of data accessible by their field names.
"""
immutable Table{Names, StorageTypes <: Tuple} #<: AbstractTable{Index,DefaultKey()}
    data::StorageTypes

    function Table(data_in::Tuple, check_sizes::Type{Val{true}} = Val{true})
        check_table(Val{Names}, StorageTypes)
        ls = map(length, data_in)
        for i in 2:length(data_in)
            if ls[i] != ls[1]
                error("Column inputs must be same length.")
            end
        end
        new(data_in)
    end

    function Table(data_in::Tuple, check_sizes::Type{Val{false}})
        check_table(Val{Names}, StorageTypes)
        new(data_in)
    end
end

@generated function check_table{Names, Types}(::Type{Val{Names}}, ::Type{Types})
    if !isa(Names, Tuple) || eltype(Names) != Symbol || length(Names) != length(unique(Names))
        str = "Table parameter 1 (Names) is expected to be a tuple of unique symbols, got $Names"
        return :(error($str))
    end
    if :Row ∈ Names
        return :( error("Field name cannot be :Row") )
    end
    N = length(Names)
    if length(Types.parameters) != N || reduce((a,b) -> a | !isa(b, DataType), false, Types.parameters)
        str = "Table parameter 2 (Types) is expected to be a Tuple{} of $N DataTypes, got $Types"
        return :(error($str))
    end
    for j = 1:N
        if eltype(Types.parameters[j]) == Types.parameters[j]
            warn("Column :$(Names[j]) storage type $(Types.parameters[j]) doesn't appear to be a storage container")
        end
    end

    return nothing
end

@compat @generated function (::Type{Table{Names}}){Names, CheckSizes}(data::Tuple, ::Type{Val{CheckSizes}} = Val{true})
    if !isa(Names, Tuple) || eltype(Names) != Symbol || length(Names) != length(unique(Names))
        str = "Table parameter 1 (Names) is expected to be a tuple of unique symbols, got $Names"
        return :(error($str))
    end

    if length(Names) == length(data.parameters)
        return quote
            $(Expr(:meta,:inline))
            Table{Names, $data}(data, Val{CheckSizes})
        end
    else
        return :(error("Can't construct Table with $(length(Names)) columns with input $data"))
    end
end

@compat @generated function (::Type{Table{Names,StorageTypes}}){Names, StorageTypes <: Tuple}()
    exprs = [:($(StorageTypes.parameters[j])()) for j = 1:length(Names)]
    return Expr(:call, Table{Names, StorageTypes}, Expr(:tuple, exprs...))
end

# TODO think about conversions...
@compat (::Type{Table{Names,Types_new}}){Names, Types, Types_new}(t::Table{Names,Types}) = convert(Table{Names,Types_new}, t)

@generated function Base.convert{Names, Names_new, Types, Types_new <: Tuple}(::Type{Table{Names_new,Types_new}}, t::Table{Names,Types})
    if !isa(Names_new, Tuple) || eltype(Names_new) != Symbol || length(Names_new) != length(Names) || length(Names_new) != length(unique(Names_new))
        str = "Cannot convert $(length(Names)) columns to new names $(Names_new)."
        return :(error($str))
    end
    if length(Types_new.parameters) != length(Types.parameters)
        str = "Cannot convert $(length(Types.pareters)) columns to $(length(Types_new.parameters)) new types."
        return :(error($str))
    end
    exprs = [:(convert(Types_new.parameters[$j], t.data[$j])) for j = 1:length(Names)]
    return Expr(:call, Table{Names_new,Types_new}, Expr(:tuple, exprs...), Val{false})
end

@compat Base.:(==){Names, Types1, Types2}(table1::Table{Names, Types1}, table2::Table{Names, Types2}) = (table1.data == table2.data)
@compat @generated function Base.:(==){Names1, Types1, Names2, Types2}(table1::Table{Names1,Types1}, table2::Table{Names2,Types2})
    try
        order = permutator(Names1, Names2)
        expr = :( table1.data[$(order[1])] == table2.data[1] )
        for j = 2:length(Names1)
            expr = Expr(:call, :(&), expr, :( table1.data[$(order[j])] == table2.data[$j] ))
        end
        return expr
    catch
        return false
    end
end

rename{Names, NewNames}(table::Table{Names}, ::Type{Val{NewNames}}) = Table{NewNames}(table.data)

@generated function rename{Names, OldName, NewName}(t::Table{Names}, ::Type{Val{OldName}}, ::Type{Val{NewName}})
    j = columnindex(Names, OldName)

    NewNames = [Names...]
    NewNames[j] = NewName
    NewNames = (NewNames...)

    return :(Table{$NewNames}(t.data, Val{false}))
end

@inline Base.names{Names}(::Table{Names}) = Names
@inline Base.names{Names, Types <: Tuple}(::Type{Table{Names,Types}}) = Names
@inline Base.names{Names}(::Type{Table{Names}}) = Names
@generated function eltypes{Names, Types}(::Union{Table{Names,Types}, Type{Table{Names,Types}}})
    elem_types = map(eltype, Types.parameters)
    quote
        $(Expr(:meta, :inline))
        $(Expr(:curly, :Tuple, elem_types...))
    end
end
@inline storagetypes{Names, Types <: Tuple}(::Table{Names,Types}) = Types
@inline storagetypes{Names, Types <: Tuple}(::Type{Table{Names,Types}}) = Types

@inline nrow(t::Table) = length(t.data[1])
@inline ncol{Names}(::Table{Names}) = length(Names)
@inline ncol{Names,Types}(::Type{Table{Names,Types}}) = length(Names)
@inline ncol{Names}(::Type{Table{Names}}) = length(Names)

#############################################################################
# Indexing columns
#############################################################################
@generated function Base.getindex{Names, GetName}(t::Table{Names}, ::Type{Val{GetName}})
    if isa(GetName, Symbol)
        if GetName == :Row
            return quote
                $(Expr(:meta, :inline))
                1:nrow(t)
            end
        end

        j = columnindex(Names, GetName)

        return quote
            $(Expr(:meta, :inline))
            t.data[$j]
        end
    elseif isa(GetName, Tuple)
        inds = columnindex(Names, GetName)
        exprs = [:(t.data[$(inds[j])]) for j = 1:length(inds)]
        expr = Expr(:call, Table{GetName}, Expr(:tuple, exprs...))

        return quote
            $(Expr(:meta, :inline))
            $expr
        end
    else
        str = "Can't get column(s) named $Name"
        return :(error($str))
    end
end

#############################################################################
# iterating over rows
#############################################################################

# TODO don't assume all containers are compatible with each other... (e.g. different implementations of dictionaries).
# Probably should define some kind of table key for this

# TODO very strong assumption about columns being of same length. Probably
# should check in start? And hope the user doesn't change the column sizes
# differently? Or do the safe thing and get the user to wrap it in @inbounds??

Base.start(t::Table) = 1
@generated function Base.next{Names}(t::Table{Names}, state)
    exprs = [:(Base.unsafe_getindex(t.data[$i],state)) for i = 1:length(Names)]
    return Expr(:tuple, Expr(:call, Row{Names}, Expr(:tuple, exprs...)), :(state + 1))
end
Base.done(t::Table, state) = state > nrow(t)

#############################################################################
# indexing rows
#############################################################################

@inline Base.length(t::Table) = nrow(t)
@inline Base.endof(t::Table) = nrow(t)
@inline Base.size(t::Table) = (nrow(t),)
@inline Base.size(t::Table, i::Integer) = i == 1 ? nrow(t) : error("Tables are one-dimensional storage containers of Rows. Consider using `ncol()`")
@inline Base.ndims(t::Table) = 1
@inline Base.isempty(t::Table) = isempty(t.data[1])

# head/tail
function head(t::Table, n = 5)
    if nrow(t) <= n
        return t
    else
        return t[1:n]
    end
end

function tail(t::Table, n = 5)
    if nrow(t) <= n
        return t
    else
        return t[end-n+1:end]
    end
end

@generated function Base.getindex{Names}(t::Table{Names}, i::Integer)
    exprs = [:(getindex(t.data[$c], i)) for c = 1:length(Names)]
    return Expr(:call, Row{Names}, Expr(:tuple, exprs...))
end
@generated function Base.getindex{Names}(t::Table{Names}, inds)
    exprs = [:(getindex(t.data[$c], inds)) for c = 1:length(Names)]
    return Expr(:call, Table{Names}, Expr(:tuple, exprs...))
end
Base.getindex{Names}(t::Table{Names}, ::Colon) = t

@generated function Base.setindex!{Names1,Names2}(t::Table{Names1}, v::Row{Names2}, i::Integer)
    if Names1 == Names2
        exprs = [:(setindex!(t.data[$c], v.data[$c], i)) for c = 1:length(Names1)]
        return Expr(:block, exprs...)
    else
        if length(Names1) != length(Names2)
            str = "Cannot assign $(length(v.parameters)) columns to $(length(Names)) columns"
            return :(error($str))
        end

        order = permutator(Names1, Names2)
        exprs = [:(setindex!(t.data[$(order[c])], v.data[$c], i)) for c = 1:length(Names1)]
        return Expr(:block, exprs...)
    end
end
@generated function Base.setindex!{Names}(t::Table{Names}, v::Tuple, i)
    if length(v.parameters) != length(Names)
        str = "Cannot assign a $(length(v.parameters))-tuple to $(length(Names)) columns"
        return :(error($str))
    end
    exprs = [:(setindex!(t.data[$c], v.data[$c], i)) for c = 1:length(Names)]
    return Expr(:block, exprs...)
end
@generated function Base.setindex!{Names1,Names2}(t::Table{Names1}, v::Table{Names2}, inds)
    if Names1 == Names2
        exprs = [:(setindex!(t.data[$c], v.data[$c], inds)) for c = 1:length(Names1)]
        return Expr(:block, exprs...)
    else
        if length(Names1) != length(Names2)
            str = "Cannot assign $(length(v.parameters)) columns to $(length(Names)) columns"
            return :(error($str))
        end

        order = permutator(Names1, Names2)
        exprs = [:(setindex!(t.data[$(order[c])], v.data[$c], inds)) for c = 1:length(Names1)]
        return Expr(:block, exprs...)
    end
end

@generated function Base.unsafe_getindex{Names}(t::Table{Names}, i::Integer)
    exprs = [:(Base.unsafe_getindex(t.data[$c], i)) for c = 1:length(Names)]
    return Expr(:call, Row{Names}, Expr(:tuple, exprs...))
end
@generated function Base.unsafe_getindex{Names}(t::Table{Names}, inds)
    exprs = [:(Base.unsafe_getindex(t.data[$c], inds)) for c = 1:length(Names)]
    return Expr(:call, Table{Names}, Expr(:tuple, exprs...))
end

@generated function Base.unsafe_setindex!{Names1,Names2}(t::Table{Names1}, v::Row{Names2}, i::Integer)
    if Names1 == Names2
        exprs = [:(Base.unsafe_setindex!(t.data[$c], v.data[$c], i)) for c = 1:length(Names1)]
        return Expr(:block, exprs...)
    else
        if length(Names1) != length(Names2)
            str = "Cannot assign $(length(v.parameters)) columns to $(length(Names)) columns"
            return :(error($str))
        end

        order = permutator(Names1, Names2)
        exprs = [:(Base.unsafe_setindex!(t.data[$(order[c])], v.data[$c], i)) for c = 1:length(Names1)]
        return Expr(:block, exprs...)
    end
end
@generated function Base.unsafe_setindex!{Names}(t::Table{Names}, v::Tuple, i)
    if length(v.parameters) != length(Names)
        str = "Cannot assign a $(length(v.parameters))-tuple to $(length(Names)) columns"
        return :(error($str))
    end
    exprs = [:(Base.unsafe_setindex!(t.data[$c], v.data[$c], i)) for c = 1:length(Names)]
    return Expr(:block, exprs...)
end
@generated function Base.unsafe_setindex!{Names1,Names2}(t::Table{Names1}, v::Table{Names2}, inds)
    if Names1 == Names2
        exprs = [:(Base.unsafe_setindex!(t.data[$c], v.data[$c], inds)) for c = 1:length(Names1)]
        return Expr(:block, exprs...)
    else
        if length(Names1) != length(Names2)
            str = "Cannot assign $(length(v.parameters)) columns to $(length(Names)) columns"
            return :(error($str))
        end

        order = permutator(Names1, Names2)
        exprs = [:(Base.unsafe_setindex!(t.data[$(order[c])], v.data[$c], inds)) for c = 1:length(Names1)]
        return Expr(:block, exprs...)
    end
end

#############################################################################
# Indexing with two variables
#############################################################################

@generated function Base.getindex{Names, GetName}(t::Table{Names}, rowinds, ::Type{Val{GetName}})
    if isa(GetName, Symbol)
        if GetName == :Row
            return quote
                $(Expr(:meta, :inline))
                (1:nrow(t))[rowinds]
            end
        end

        j = columnindex(Names, GetName)

        return quote
            $(Expr(:meta, :inline))
            t.data[$j][rowinds]
        end
    elseif isa(GetName, Tuple) && eltype(GetName) == Symbol && length(GetName) == length(unique(GetName))
        inds = columnindex(Names, GetName)
        exprs = [:(t.data[$(inds[j])][rowinds]) for j = 1:length(inds)]
        if rowinds <: Integer
            expr = Expr(:call, Row{GetName}, Expr(:tuple, exprs...))
        else
            expr = Expr(:call, Table{GetName}, Expr(:tuple, exprs...))
        end

        return quote
            $(Expr(:meta, :inline))
            $expr
        end
    else
        str = "Can't get column(s) named $Name"
        return :(error($str))
    end
end

@generated function Base.unsafe_getindex{Names, GetName}(t::Table{Names}, rowinds, ::Type{Val{GetName}})
    if isa(GetName, Symbol)
        if GetName == :Row
            return quote
                $(Expr(:meta, :inline))
                Base.unsafe_getindex(1:nrow(t), rowinds)
            end
        end

        j = columnindex(Names, GetName)

        return quote
            $(Expr(:meta, :inline))
            Base.unsafe_getindex(t.data[$j], rowinds)
        end
    elseif isa(GetName, Tuple) && eltype(GetName) == Symbol && length(GetName) == length(unique(GetName))
        inds = columnindex(Names, GetName)
        exprs = [:(Base.unsafe_getindex(t.data[$(inds[j])], rowinds)) for j = 1:length(inds)]
        if rowinds <: Integer
            expr = Expr(:call, Row{GetName}, Expr(:tuple, exprs...))
        else
            expr = Expr(:call, Table{GetName}, Expr(:tuple, exprs...))
        end

        return quote
            $(Expr(:meta, :inline))
            $expr
        end
    else
        str = "Can't get column(s) named $Name"
        return :(error($str))
    end
end

@generated function Base.setindex!{Names, GetName}(t::Table{Names}, value, rowinds, ::Type{Val{GetName}})
    if isa(GetName, Symbol)
        j = columnindex(Names, GetName)
        if value <: Union{Cell, Column}
            return quote
                $(Expr(:meta, :inline))
                t.data[$j][rowinds] = value.data
            end
        else
            return quote
                $(Expr(:meta, :inline))
                t.data[$j][rowinds] = value
            end
        end
    elseif isa(GetName, Tuple) && eltype(GetName) == Symbol && length(GetName) == length(unique(GetName))
        inds = columnindex(Names, GetName)
        if value <: Union{Row, Table}
            order = permutator(names(value), GetName)
            exprs = [:(t.data[$(inds[j])][rowinds] = value.data[$(order[j])]) for j = 1:length(inds)]
            return Expr(:block, Expr(:meta, :inline), exprs...)
        elseif value <: Tuple && length(value.parameters) == length(inds)
            exprs = [:(t.data[$(inds[j])][rowinds] = value[$j]) for j = 1:length(inds)]
            return Expr(:block, Expr(:meta, :inline), exprs...)
        else
            str = "Can't set columns $GetName with a $value"
            return :(error($str))
        end
    else
        str = "Can't set column(s) named $Name"
        return :(error($str))
    end
end

@generated function Base.unsafe_setindex!{Names, GetName}(t::Table{Names}, value, rowinds, ::Type{Val{GetName}})
    if isa(GetName, Symbol)
        j = columnindex(Names, GetName)
        if value <: Union{Cell, Column}
            return quote
                $(Expr(:meta, :inline))
                Base.unsafe_setindex!(t.data[$j], value.data, rowinds)
            end
        else
            return quote
                $(Expr(:meta, :inline))
                Base.unsafe_setindex!(t.data[$j], value, rowinds)
            end
        end
    elseif isa(GetName, Tuple) && eltype(GetName) == Symbol && length(GetName) == length(unique(GetName))
        inds = columnindex(Names, GetName)
        if value <: Union{Row, Table}
            order = permutator(names(value), GetName)
            exprs = [:(Base.unsafe_setindex!(t.data[$(inds[j])], value.data[$(order[j])], rowinds)) for j = 1:length(inds)]
            return Expr(:block, Expr(:meta, :inline), exprs...)
        elseif value <: Tuple && length(value.parameters) == length(inds)
            exprs = [:(Base.unsafe_setindex!(t.data[$(inds[j])], value[$j], rowinds)) for j = 1:length(inds)]
            return Expr(:block, Expr(:meta, :inline), exprs...)
        else
            str = "Can't set columns $GetName with a $value"
            return :(error($str))
        end
    else
        str = "Can't set column(s) named $Name"
        return :(error($str))
    end
end

Base.getindex{Names}(t::Table{Names}, inds, ::Colon) = t[inds]
Base.unsafe_getindex{Names}(t::Table{Names}, inds, ::Colon) = Base.unsafe_getindex(t, inds)
Base.setindex!{Names}(t::Table{Names}, value, inds, ::Colon) = Base.setindex!(t, value, inds)
Base.unsafe_setindex!{Names}(t::Table{Names}, value, inds, ::Colon) = Base.unsafe_setindex!(t, value, inds)

#############################################################################
# push, pop, etc
#############################################################################

@generated function Base.pop!{Names}(t::Table{Names})
    exprs = [:(pop!(t.data[$c])) for c = 1:length(Names)]
    return Expr(:call, Row{Names}, Expr(:tuple, exprs...))
end

@generated function Base.shift!{Names}(t::Table{Names})
    exprs = [:(shift!(t.data[$c])) for c = 1:length(Names)]
    return Expr(:call, Row{Names}, Expr(:tuple, exprs...))
end

@generated function Base.push!{Names}(t::Table{Names}, v::Tuple)
    if length(Names) != length(v.parameters)
        str = "Cannot assign a $(length(v.parameters))-tuple to $(length(Names)) columns"
        return :(error($str))
    end
    exprs = [:(push!(t.data[$c], v[$c])) for c = 1:length(Names)]
    return Expr(:block, exprs..., :t)
end
@generated function Base.push!{Names1,Names2}(t::Table{Names1}, r::Row{Names2})
    order = permutator(Names1, Names2)
    exprs = [:(push!(t.data[$(order[c])], r.data[$c])) for c = 1:length(Names1)]
    return Expr(:block, exprs..., :t)
end

@generated function Base.unshift!{Names}(t::Table{Names}, v::Tuple)
    if length(Names) != length(v.parameters)
        str = "Cannot assign a $(length(v.parameters))-tuple to $(length(Names)) columns"
        return :(error($str))
    end
    exprs = [:(unshift!(t.data[$c], v[$c])) for c = 1:length(Names)]
    return Expr(:block, exprs..., :t)
end
@generated function Base.unshift!{Names1,Names2}(t::Table{Names1}, r::Row{Names2})
    order = permutator(Names1, Names2)
    exprs = [:(unshift!(t.data[$(order[c])], r.data[$c])) for c = 1:length(Names1)]
    return Expr(:block, exprs..., :t)
end

@generated function Base.append!{Names}(t::Table{Names}, v::Tuple)
    if length(Names) != length(v.parameters)
        str = "Cannot assign a $(length(v.parameters))-tuple to $(length(Names)) columns"
        return :(error($str))
    end
    exprs = [:(append!(t.data[$c], v[$c])) for c = 1:length(Names)]
    return Expr(:block, exprs..., :t)
end
@generated function Base.append!{Names1,Names2}(t::Table{Names1}, t2::Table{Names2})
    order = permutator(Names1, Names2)
    exprs = [:(append!(t.data[$(order[c])], t2.data[$c])) for c = 1:length(Names1)]
    return Expr(:block, exprs..., :t)
end

@generated function Base.prepend!{Names}(t::Table{Names}, v::Tuple)
    if length(Names) != length(v.parameters)
        str = "Cannot assign a $(length(v.parameters))-tuple to $(length(Names)) columns"
        return :(error($str))
    end
    exprs = [:(prepend!(t.data[$c], v[$c])) for c = 1:length(Names)]
    return Expr(:block, exprs..., :t)
end
@generated function Base.prepend!{Names1,Names2}(t::Table{Names1}, t2::Table{Names2})
    order = permutator(Names1, Names2)
    exprs = [:(prepend!(t.data[$(order[c])], t2.data[$c])) for c = 1:length(Names1)]
    return Expr(:block, exprs..., :t)
end

@generated function Base.insert!{Names}(t::Table{Names}, i::Integer, v::Tuple)
    if length(Names) != length(v.parameters)
        str = "Cannot assign a $(length(v.parameters))-tuple to $(length(Names)) columns"
        return :(error($str))
    end
    exprs = [:(insert!(t.data[$c], i, v[$c])) for c = 1:length(Names)]
    return Expr(:block, exprs..., :t)
end
@generated function Base.insert!{Names1,Names2}(t::Table{Names1}, i::Integer, r::Row{Names2})
    order = permutator(Names1, Names2)
    exprs = [:(insert!(t.data[$(order[c])], i, r.data[$c])) for c = 1:length(Names1)]
    return Expr(:block, exprs..., :t)
end

@generated function Base.empty!{Names}(t::Table{Names})
    exprs = [:(empty!(t.data[$c])) for c = 1:length(Names)]
    return Expr(:block, exprs..., :t)
end

@generated function Base.deleteat!{Names}(t::Table{Names}, i)
    exprs = [:(deleteat!(t.data[$c], i)) for c = 1:length(Names)]
    return Expr(:block, exprs..., :t)
end

@generated function Base.splice!{Names}(t::Table{Names}, i::Integer)
    exprs = [:(splice!(t.data[$c], i)) for c = 1:length(Names)]
    return Expr(:call, Row{Names}, Expr(:tuple, exprs...))
end
@generated function Base.splice!{Names}(t::Table{Names}, i::Integer, v::Tuple)
    exprs = [:(splice!(t.data[$c], i, v[$c])) for c = 1:length(Names)]
    return Expr(:call, Row{Names}, Expr(:tuple, exprs...))
end
@generated function Base.splice!{Names1,Names2}(t::Table{Names1}, i::Integer, v::Union{Row{Names2},Table{Names2}})
    order = permutator(Names2, Names1)
    exprs = [:(splice!(t.data[$c], i, v.data[$(order[c])])) for c = 1:length(Names1)]
    return Expr(:call, Row{Names1}, Expr(:tuple, exprs...))
end
@generated function Base.splice!{Names}(t::Table{Names}, i)
    exprs = [:(splice!(t.data[$c], i)) for c = 1:length(Names)]
    return Expr(:call, Table{Names}, Expr(:tuple, exprs...))
end
@generated function Base.splice!{Names}(t::Table{Names}, i, v::Tuple)
    exprs = [:(splice!(t.data[$c], i, v[$c])) for c = 1:length(Names)]
    return Expr(:call, Table{Names}, Expr(:tuple, exprs...))
end
@generated function Base.splice!{Names1,Names2}(t::Table{Names1}, i, v::Union{Row{Names2},Table{Names2}})
    order = permutator(Names2, Names1)
    exprs = [:(splice!(t.data[$c], i, v.data[$(order[c])])) for c = 1:length(Names1)]
    return Expr(:call, Table{Names1}, Expr(:tuple, exprs...))
end

# Vertically concatenate rows and tables into tables
@generated function Base.vcat{Names}(t1::Union{Row{Names}, Table{Names}})
    exprs = [:(vcat(t1.data[$c])) for c = 1:length(Names)]
    return Expr(:call, Table{Names}, Expr(:tuple, exprs...))
end

# TODO still problem with mixed Row-Table vcat inherited from Base scalar-Vector vcat
@generated function Base.vcat{Names1, Names2}(t1::Union{Row{Names1}, Table{Names1}}, t2::Union{Row{Names2}, Table{Names2}})
    if Names1 == Names2
        exprs = [:(vcat(t1.data[$c], t2.data[$c])) for c = 1:length(Names1)]
        return Expr(:call, Table{Names1}, Expr(:tuple, exprs...))
    else
        if length(Names1) != length(Names2)
            str = "Cannot match $(length(v.parameters)) columns to $(length(Names)) columns"
        end

        order = permutator(Names1, Names2)
        exprs = [:(vcat(t1.data[$c], t2.data[$(order[c])])) for c = 1:length(Names1)]
        return Expr(:call, Table{Names1}, Expr(:tuple, exprs...))
    end
end

Base.vcat{Names1, Names2}(t1::Union{Row{Names1}, Table{Names1}}, t2::Union{Row{Names2}, Table{Names2}}, ts::Union{Row, Table}...) = vcat(vcat(t1, t2), ts...)

# Horizontally concatenate columns and tables into tables
@generated Base.hcat{Name}(c::Column{Name}) = :( Table{$((Name,))}((c.data,)) )
Base.hcat(t::Table) = t

@generated function Base.hcat(r1::Union{Column,Table}, r2::Union{Column,Table})
    names1 = (r1 <: Column ? (name(r1),) : names(r1))
    names2 = (r2 <: Column ? (name(r2),) : names(r2))

    if length(intersect(names1, names2)) != 0
        str = "Column names are not distinct. Got $names1 and $names2"
        return :(error($str))
    end

    newnames = (names1..., names2...)
    exprs1 = r1 <: Column ? [:(r1.data)] : [:(r1.data[$j]) for j = 1:length(names1)]
    exprs2 = r2 <: Column ? [:(r2.data)] : [:(r2.data[$j]) for j = 1:length(names2)]
    exprs = vcat(exprs1, exprs2)

    return Expr(:call, Table{newnames}, Expr(:tuple, exprs...))
end

Base.hcat(t1::Union{Column,Table}, t2::Union{Column,Table}, ts::Union{Column,Table}...) = hcat(hcat(t1, t2), ts...)

# copy
@generated function Base.copy{Names}(t::Table{Names})
    exprs = [:(copy(t.data[$j])) for j = 1:length(Names)]
    return Expr(:call, Table{Names}, Expr(:tuple, exprs...))
end

macro Table(exprs...)
    N = length(exprs)
    names = Vector{Any}(N)
    values = Vector{Any}(N)
    for i = 1:N
        expr = exprs[i]
        if expr.head != :(=) && expr.head != :(kw) # strange Julia bug, see issue 7669
            error("A Expecting expression like @Table(name1::Type1 = value1, name2::Type2 = value2) or @Table(name1 = value1, name2 = value2)")
        end
        if isa(expr.args[1],Symbol)
            names[i] = (expr.args[1])
            values[i] = esc(expr.args[2])
        elseif isa(expr.args[1],Expr)
            if expr.args[1].head != :(::) || length(expr.args[1].args) != 2
                error("A Expecting expression like @Table(name1::Type1 = value1, name2::Type2 = value2) or @Table(name1 = value1, name2 = value2)")
            end
            names[i] = (expr.args[1].args[1])
            values[i] = esc(Expr(:call, :convert, expr.args[1].args[2], expr.args[2]))
        else
            error("A Expecting expression like @Table(name1::Type1 = value1, name2::Type2 = value2) or @Table(name1 = value1, name2 = value2)")
        end
    end

    tabletype = TypedTables.Table{(names...)}
    return Expr(:call, tabletype, Expr(:tuple, values...))
end

#=

@generated Table{Index<:FieldIndex,StorageTypes<:Tuple}(::Index,data_in::StorageTypes,check_sizes::Union{Type{Val{true}},Type{Val{false}}} = Val{true}) = :(Table{$(Index()),$(eltypes(Index)),$StorageTypes}(data_in,check_sizes) )
@generated Table{Index<:FieldIndex}(::Index,check_sizes::Union{Type{Val{true}},Type{Val{false}}} = Val{true}) = :(Table{$(Index()),$(eltypes(Index)),$(makestoragetypes(eltypes(Index)))}(instantiate_tuple($(makestoragetypes(eltypes(Index)))),check_sizes) )
@generated Table{Index,ElTypes}(row::Row{Index,ElTypes},check_sizes::Union{Type{Val{true}},Type{Val{false}}} = Val{true}) = :(Table{Index,ElTypes,$(makestoragetypes(ElTypes))}(ntuple(i->[row.data[i]],$(length(Index))),check_sizes) )
@generated Table{F,ElType,StorageType}(col::Column{F,ElType,StorageType},check_sizes::Union{Type{Val{true}},Type{Val{false}}} = Val{true}) = :(Table{$(FieldIndex(F)),Tuple{ElType},Tuple{StorageType}}((col.data,),check_sizes) )

@generated function Base.call{Fields,StorageTypes<:Tuple}(Index::FieldIndex{Fields},x::StorageTypes)
    if StorageTypes == eltypes(Index)
        return :(Row{$(Index()),$StorageTypes}(x))
    end

    try
        check_table(Index(),eltypes(Index),StorageTypes)
        return :(Table{$(Index()),$(eltypes(Index)),$StorageTypes}(x))
    catch
        str = "Can't instantiate a Row or Table having index $(Index()) with data of types $StorageTypes."
        error(str)
    end
end
Base.call{Fields}(::FieldIndex{Fields},x...) = error("Must instantiate Row with a tuple of type $(eltypes(Index)) or a Table with a tuple of appropriate storage containers")

@generated function makestoragetypes{T<:Tuple}(::Type{T})
    eltypes = T.parameters
    storagetypes = Vector{DataType}(length(eltypes))
    for i = 1:length(eltypes)
        storagetypes[i] = makestoragetype(eltypes[i])
    end
    return :(Tuple{$(storagetypes...)} )
end

Base.convert{Index,DataTypes<:Tuple,StorageTypes<:Tuple}(::Type{StorageTypes},x::Table{Index,DataTypes,StorageTypes}) = x.data

=={Index,ElTypes,StorageType}(table1::Table{Index,ElTypes,StorageType},table2::Table{Index,ElTypes,StorageType}) = (table1.data == table2.data)
function =={Index1,ElTypes1,StorageType1,Index2,ElTypes2,StorageType2}(table1::Table{Index1,ElTypes1,StorageType1},table2::Table{Index2,ElTypes2,StorageType2})
    if samefields(Index1,Index2)
        idx = Index2[Index1]
        for i = 1:length(Index1)
            if table1.data[i] != table2.data[idx[i]]
                return false
            end
        end
        return true
    else
        return false
    end
end

# Conversion between Dense and normal Tables??
# Some more convient versions. (a) One that takes pairs of fields and storage.
#                              (b) A macro @table(x=[1,2,3],y=["a","b","c"]) -> Table(Field{:x,eltype([1,2,3)}()=>[1,2,3], Field{:y,eltype{["a","b","c"]}()=>["a","b","c"])

# Data from the index
Base.names{Index}(table::Table{Index}) = names(Index)
Base.names{Index,ElTypes,StorageTypes}(table::Type{Table{Index,ElTypes,StorageTypes}}) = names(Index)
eltypes{Index}(table::Table{Index}) = eltypes(Index)
eltypes{Index,ElTypes,StorageTypes}(table::Type{Table{Index,ElTypes,StorageTypes}}) = eltypes(Index)
storagetypes{Index,ElTypes,StorageTypes}(table::Table{Index,ElTypes,StorageTypes}) = StorageTypes
storagetypes{Index,ElTypes,StorageTypes}(table::Type{Table{Index,ElTypes,StorageTypes}}) = StorageTypes

@generated rename{Index,ElTypes,StorageTypes,NewIndex<:FieldIndex}(table::Table{Index,ElTypes,StorageTypes}, new_names::NewIndex) = :( Table{$(rename(Index,NewIndex())),$ElTypes,$StorageTypes}(table.data, Val{false}) )
@generated rename{Index,ElTypes,DataTypes,OldFields,NewFields}(table::Table{Index,ElTypes,DataTypes}, old_names::OldFields, ::NewFields) = :(Table($(rename(Index,OldFields(),NewFields())),table.data, Val{false}))

# Vector-like introspection
Base.length{Index,ElTypes,StorageTypes}(table::Table{Index,ElTypes,StorageTypes}) = length(table.data[1])
Base.length{Index}(table::Table{Index,Tuple{},Tuple{}}) = 0

@generated ncol{Index,ElTypes,StorageTypes}(table::Table{Index,ElTypes,StorageTypes}) = :($(length(Index)))
nrow{Index,ElTypes,StorageTypes}(table::Table{Index,ElTypes,StorageTypes}) = length(table.data[1])
Base.size{Index,ElTypes,StorageTypes}(table::Table{Index,ElTypes,StorageTypes}) = (length(table.data[1]),length(Index))
Base.size{Index,ElTypes,StorageTypes}(table::Table{Index,ElTypes,StorageTypes},i::Int) = i == 2 ? length(Index) : (i == 1 ? length(table.data[1]) : error("Tables are two-dimensional"))
@generated Base.eltype{Index,ElTypes,StorageTypes}(table::Table{Index,ElTypes,StorageTypes}) = :($(eltypes(Index)))
Base.isempty{Index,ElTypes,StorageTypes}(table::Table{Index,ElTypes,StorageTypes}) = isempty(table.data[1])
Base.endof{Index,ElTypes,StorageTypes}(table::Table{Index,ElTypes,StorageTypes}) = endof(table.data[1])
index{Index}(table::Table{Index}) = Index

samefields{Index1,Index2}(::Table{Index1},::Table{Index2}) = samefields(Index1,Index2)
samefields{Index1,Index2}(::Row{Index1},::Table{Index2}) = samefields(Index1,Index2)
samefields{Index1,Index2}(::Table{Index1},::Row{Index2}) = samefields(Index1,Index2)
samefields{Index1<:FieldIndex,Index2}(::Index1,::Table{Index2}) = samefields(Index1(),Index2)
samefields{Index1,Index2<:FieldIndex}(::Table{Index1},::Index2) = samefields(Index1,Index2())

# Iterators
Base.start{Index,ElTypes,StorageTypes}(table::Table{Index,ElTypes,StorageTypes}) = 1
Base.next{Index,ElTypes,StorageTypes}(table::Table{Index,ElTypes,StorageTypes},i) = (table[i],i+1)
Base.done{Index,ElTypes,StorageTypes}(table::Table{Index,ElTypes,StorageTypes},i) = (i-1 == length(table))

# get/set index (not fast... use generated functions?)
Base.getindex{Index,ElTypes,StorageTypes}(table::Table{Index,ElTypes,StorageTypes},idx::Int) = Row{Index,ElTypes}(ntuple(i->getindex(table.data[i],idx),length(Index)))
Base.getindex{Index,ElTypes,StorageTypes}(table::Table{Index,ElTypes,StorageTypes},idx) = Table{Index,ElTypes,StorageTypes}(ntuple(i->getindex(table.data[i],idx),length(Index)))

Base.setindex!{Index,ElTypes,StorageTypes}(table::Table{Index,ElTypes,StorageTypes},val::Row{Index},idx::Int) = (for i = 1:length(Index); setindex!(table.data[i],val.data[i],idx); end; val)
Base.setindex!{Index,ElTypes,StorageTypes}(table::Table{Index,ElTypes,StorageTypes},val::Table{Index,ElTypes, StorageTypes},idx) = (for i = 1:length(Index); setindex!(table.data[i],val.data[i],idx); end; val)

Base.unsafe_getindex{Index,ElTypes,StorageTypes}(table::Table{Index,ElTypes,StorageTypes},idx::Int) = Row{Index, ElTypes}(ntuple(i->unsafe_getindex(table.data[i],idx),length(Index)))
Base.unsafe_getindex{Index,ElTypes,StorageTypes}(table::Table{Index,ElTypes,StorageTypes},idx) = Table{Index,ElTypes,StorageTypes}(ntuple(i->Base.unsafe_getindex(table.data[i],idx),length(Index)))
Base.unsafe_setindex!{Index,ElTypes,StorageTypes}(table::Table{Index,ElTypes,StorageTypes},val::Row{Index},idx::Int) = for i = 1:length(Index); Base.unsafe_setindex!(table.data[i],val.data[i],idx); end
Base.unsafe_setindex!{Index,ElTypes,StorageTypes}(table::Table{Index,ElTypes,StorageTypes},val::Table{Index,StorageTypes},idx) = for i = 1:length(Index); Base.unsafe_setindex!(table.data[i],val.data[i],idx); end

# Push, append, pop
function Base.push!{Index,ElTypes,StorageTypes}(table::Table{Index,ElTypes,StorageTypes},row::Row{Index,ElTypes})
    for i = 1:length(Index)
        push!(table.data[i],row.data[i])
    end
    table
end
@inline Base.push!{Index,ElTypes,StorageTypes,Index2,ElTypes2}(table::Table{Index,ElTypes,StorageTypes},row::Row{Index2,ElTypes2}) = push!(table,row[Index])
@inline Base.push!{Index,ElTypes,StorageTypes}(table::Table{Index,ElTypes,StorageTypes},data_in::ElTypes) = push!(table, Row{Index,ElTypes}(data_in))
function Base.append!{Index,ElTypes,StorageTypes}(table::Table{Index,ElTypes,StorageTypes},table_in::Table{Index,ElTypes,StorageTypes})
    for i in 1:length(Index)
        append!(table.data[i],table_in.data[i])
    end
end
@inline Base.append!{Index,ElTypes,StorageTypes,Index2,ElTypes2,StorageTypes2}(table::Table{Index,ElTypes,StorageTypes},table_in::Table{Index2,ElTypes2,StorageTypes2}) = append!(table,table_in[Index])
function Base.append!{Index,ElTypes,StorageTypes}(table::Table{Index,ElTypes,StorageTypes},data_in::StorageTypes)
    ls = map(length,data_in)
    for i in 2:length(Index)
        if ls[i] != ls[1]
            error("Column inputs must be same length.")
        end
    end

    for i in 1:length(data_in)
        append!(table.data[i],data_in[i])
    end
end
function Base.prepend!{Index,ElTypes,StorageTypes}(table::Table{Index,ElTypes,StorageTypes},table_in::Table{Index,ElTypes,StorageTypes})
    for i in 1:length(Index)
        prepend!(table.data[i],table_in.data[i])
    end
end
@inline Base.prepend!{Index,ElTypes,StorageTypes,Index2,ElTypes2,StorageTypes2}(table::Table{Index,ElTypes,StorageTypes},table_in::Table{Index2,ElTypes2,StorageTypes2}) = prepend!(table,table_in[Index])
function Base.prepend!{Index,ElTypes,StorageTypes}(table::Table{Index,ElTypes,StorageTypes},data_in::StorageTypes)
    ls = map(length,data_in)
    for i in 2:length(Index)
        if ls[i] != ls[1]
            error("Column inputs must be same length.")
        end
    end

    for i in 1:length(data_in)
        prepend!(table.data[i],data_in[i])
    end
end
function Base.unshift!{Index,ElTypes,StorageTypes}(table::Table{Index,ElTypes,StorageTypes},row::Row{Index,ElTypes})
    for i = 1:length(Index)
        unshift!(table.data[i],row.data[i])
    end
    table
end
@inline Base.unshift!{Index,ElTypes,StorageTypes,Index2,ElTypes2}(table::Table{Index,ElTypes,StorageTypes},row::Row{Index2,ElTypes2}) = unshift!(table,row[Index])
@inline Base.unshift!{Index,ElTypes,StorageTypes}(table::Table{Index,ElTypes,StorageTypes},data_in::ElTypes) = unshift!(table, Row{Index,ElTypes}(data_in))
function Base.pop!{Index,ElTypes,StorageTypes}(table::Table{Index,ElTypes,StorageTypes})
    return Row{Index,eltypes(Index)}(ntuple(i->pop!(table.data[i]),length(Index)))
end
function Base.shift!{Index,ElTypes,StorageTypes}(table::Table{Index,ElTypes,StorageTypes})
    return Row{Index,eltypes(Index)}(ntuple(i->shift!(table.data[i]),length(Index)))
end
function Base.empty!{Index,ElTypes,StorageTypes}(table::Table{Index,ElTypes,StorageTypes})
    for i = 1:length(Index)
        empty!(table.data[i])
    end
end

# copies
Base.copy{Index,ElTypes,StorageTypes}(table::Table{Index,ElTypes,StorageTypes}) = Table{Index,ElTypes,StorageTypes}(copy(table.data))
Base.deepcopy{Index,ElTypes,StorageTypes}(table::Table{Index,ElTypes,StorageTypes}) = Table{Index,ElTypes,StorageTypes}(deepcopy(table.data))

# Some output
function Base.summary{Index,ElTypes,StorageTypes}(io::IO,table::Table{Index,ElTypes,StorageTypes})
    println(io,"$(ncol(table))-column × $(nrow(table))-row Table{$Index,$StorageTypes}")
end

head(t::Table, n = 5) = length(t) >= n ? t[1:n] : t
tail(t::Table, n = 5) = length(t) >= n ? t[end-n+1:end] : t

# Can create a subtable easily by selecting columns
@generated Base.getindex{Index,ElTypes,StorageTypes,F<:Field}(table::Table{Index,ElTypes,StorageTypes},::F) = :( table.data[$(Index[F()])] )
Base.getindex{Index,ElTypes,StorageTypes,F<:DefaultKey}(table::Table{Index,ElTypes,StorageTypes},::F) = 1:length(table.data[1])

@generated function Base.getindex{Index,ElTypes,StorageTypes,NewIndex<:FieldIndex}(table::Table{Index,ElTypes,StorageTypes},::NewIndex) # Need special possibility for DefaultKey, which may or many not exist in array
    if issubset(NewIndex(),Index+DefaultKey())
        str = "("
        for i = 1:length(NewIndex())
            if NewIndex()[i] == DefaultKey()
                str *= "TableKey{Index,ElTypes,StorageTypes}(Ref(table))"
            else
                str *= "table.data[$(Index[NewIndex()[i]])]"
            end
            if i < length(NewIndex()) || i == 1
                str *= ","
            end
        end
        str *= ")"

        return :(Table(NewIndex(),$(parse(str))))
    else
        str = "Cannot index $NewIndex from table with indices $Index"
        return :(error($str))
    end
end

@generated function Base.getindex{Index,DataTypes,StorageTypes,I}(table::Table{Index,DataTypes,StorageTypes},::Type{Val{I}})
    if isa(I, Int) # TODO: Figure out if this one can be made faster... seems to be some overhead in creating the new Column (similarly for Table below)
        return :( $(Expr(:meta,:inline)); $(Column{Index[I],DataTypes.parameters[I],StorageTypes.parameters[I]})(table.data[$I]) )
    elseif isa(I, Symbol)
        if I == name(DefaultKey)
            return :( 1:length(table) )
        else
            return :( table.data[$(Index[Val{I}])] )
        end
    elseif isa(I, Tuple)
        l_I = length(I)
        if isa(I, NTuple{l_I, Int})
            expr = Expr(:tuple, ntuple(i-> :(table.data[$(I[i])]), l_I)...)
            return :( Table($(Index[Val{I}]), $expr, Val{false}) )
        elseif isa(I, NTuple{l_I, Symbol})
            expr = Expr(:tuple, ntuple(i -> I[i] == name(DefaultKey) ? :( 1:length(table) ) : :(table.data[$(Index[Val{I[i]}])]), l_I)...)
            return :( Table($(Index[Val{Index[Val{I}]}]), $expr, Val{false}) )
        else
            str = "Can't index Table with fields $Fields with a Val{$I}"
            return :(error($str))
        end
    else # e.g. UnitRange{Int} and other methods of indexing a Tuple
        str = "Can't index Table with fields $Fields with a Val{$I}"
        return :(error($str))
    end
end

# Simultaneously indexing by single row and column
@generated Base.getindex{Index,ElTypes,StorageTypes,F<:Field}(table::Table{Index,ElTypes,StorageTypes},idx::Int,::F) = :( table.data[$(Index[F()])][idx] )
Base.getindex{Index,ElTypes,StorageTypes,F<:DefaultKey}(table::Table{Index,ElTypes,StorageTypes},idx::Int,::F) = idx

@generated function Base.getindex{Index,ElTypes,StorageTypes,NewIndex<:FieldIndex}(table::Table{Index,ElTypes,StorageTypes},idx::Int,::NewIndex) # Need special possibility for DefaultKey, which may or many not exist in array
    if issubset(NewIndex(),Index+DefaultKey())
        str = "("
        for i = 1:length(NewIndex())
            if NewIndex()[i] == DefaultKey()
                str *= "idx"
            else
                str *= "table.data[$(Index[NewIndex()[i]])][idx]"
            end
            if i < length(NewIndex()) || i == 1
                str *= ","
            end
        end
        str *= ")"

        return :(Row(NewIndex(),$(parse(str))))
    else
        str = "Cannot index $NewIndex from table with indices $Index"
        return :(error($str))
    end
end

@generated function Base.getindex{Index,DataTypes,StorageTypes,I}(table::Table{Index,DataTypes,StorageTypes},idx::Int,::Type{Val{I}})
    if isa(I, Int) # TODO: Figure out if this one can be made faster... seems to be some overhead in creating the new Column (similarly for Table below)
        return :( $(Expr(:meta,:inline)); $(Cell{Index[I],DataTypes.parameters[I]})(table.data[$I][idx]) )
    elseif isa(I, Symbol)
        if I == name(DefaultKey)
            return :( (1:length(table))[idx] )
        else
            return :( table.data[$(Index[Val{I}])][idx] )
        end
    elseif isa(I, Tuple)
        l_I = length(I)
        if isa(I, NTuple{l_I, Int})
            expr = Expr(:tuple, ntuple(i-> :(table.data[$(I[i])][idx]), l_I)...)
            return :( Row($(Index[Val{I}]), $expr) )
        elseif isa(I, NTuple{l_I, Symbol})
            expr = Expr(:tuple, ntuple(i -> :(table.data[$(Index[Val{I[i]}])][idx]), l_I)...)
            return :( Row($(Index[Val{Index[Val{I}]}]), $expr) )
        else
            str = "Can't index Table with fields $Fields with a Val{$I}"
            return :(error($str))
        end
    else # e.g. UnitRange{Int} and other methods of indexing a Tuple
        str = "Can't index Table with fields $Fields with a Val{$I}"
        return :(error($str))
    end
end

# Simultaneously indexing by multiple rows and column
@generated Base.getindex{Index,ElTypes,StorageTypes,F<:Field}(table::Table{Index,ElTypes,StorageTypes},idx,::F) = :( table.data[$(Index[F()])][idx] )
Base.getindex{Index,ElTypes,StorageTypes,F<:DefaultKey}(table::Table{Index,ElTypes,StorageTypes},idx,::F) = (1:length(table.data[1]))[idx]

@generated function Base.getindex{Index,ElTypes,StorageTypes,NewIndex<:FieldIndex}(table::Table{Index,ElTypes,StorageTypes},idx,::NewIndex) # Need special possibility for DefaultKey, which may or many not exist in array
    if issubset(NewIndex(),Index+DefaultKey())
        str = "("
        for i = 1:length(NewIndex())
            if NewIndex()[i] == DefaultKey()
                str *= "TableKey{Index,StorageTypes}(Ref(table))"
            else
                str *= "table.data[$(Index[NewIndex()[i]])][idx]"
            end
            if i < length(NewIndex()) || i == 1
                str *= ","
            end
        end
        str *= ")"

        return :(Table(NewIndex(),$(parse(str))))
    else
        str = "Cannot index $NewIndex from table with indices $Index"
        return :(error($str))
    end
end

@generated function Base.getindex{Index,DataTypes,StorageTypes,I}(table::Table{Index,DataTypes,StorageTypes},idx,::Type{Val{I}})
    if isa(I, Int) # TODO: Figure out if this one can be made faster... seems to be some overhead in creating the new Column (similarly for Table below)
        return :( $(Expr(:meta,:inline)); $(Column{Index[I],DataTypes.parameters[I],StorageTypes.parameters[I]})(table.data[$I][idx]) )
    elseif isa(I, Symbol)
        if I == name(DefaultKey)
            return :( (1:length(table))[idx] )
        else
            return :( table.data[$(Index[Val{I}])][idx] )
        end
    elseif isa(I, Tuple)
        l_I = length(I)
        if isa(I, NTuple{l_I, Int})
            expr = Expr(:tuple, ntuple(i-> :(table.data[$(I[i])][idx]), l_I)...)
            return :( Table($(Index[Val{I}]), $expr, Val{false}) )
        elseif isa(I, NTuple{l_I, Symbol})
            expr = Expr(:tuple, ntuple(i -> :(table.data[$(Index[Val{I[i]}])][idx]), l_I)...)
            return :( Table($(Index[Val{Index[Val{I}]}]), $expr, Val{false}) )
        else
            str = "Can't index Table with fields $Fields with a Val{$I}"
            return :(error($str))
        end
    else # e.g. UnitRange{Int} and other methods of indexing a Tuple
        str = "Can't index Table with fields $Fields with a Val{$I}"
        return :(error($str))
    end
end

Base.getindex{Index,ElTypes,StorageTypes}(table::Table{Index,ElTypes,StorageTypes},idx,::Colon) = table[idx]

function Base.filter!{Index,ElTypes,StorageTypes}(idx::Vector{Bool}, table::Table{Index,ElTypes,StorageTypes})
    for i = 1:ncol(table)
        tmp = table.data[i][idx]
        resize!(table.data[i],length(tmp))
        table.data[i][:] = tmp
    end
    table
end

# TODO - sub() returns table/row with sub as data elements.

# Concatenate rows and tables into tables
@inline Base.vcat(row::Row) = Table(row)
@inline Base.vcat(table::Table) = table

@generated Base.vcat{Index}(row1::Row{Index},row2::Union{Row{Index},Table{Index}}) = :( Table{Index,$(eltypes(Index)),$(makestoragetypes(eltypes(Index)))}(ntuple(i->vcat(row1.data[i],row2.data[i]),$(length(Index)))) )
@generated Base.vcat{Index}(table1::Table{Index},table2::Union{Row{Index},Table{Index}}) = :( Table{Index,$(eltypes(table1)),$(storagetypes(table1))}(ntuple(i->vcat(table1.data[i],table2.data[i]),length(Index))) )

# Otherwise, fields don't match...
function Base.vcat{Index1,Index2}(x1::Union{Row{Index1},Table{Index1}}, x2::Union{Row{Index2},Table{Index2}})
    if samefields(Index1, Index2)
        vcat(x1, x2[Index1])
    else
        error("Indices $Index1 and $Index2 don't match")
    end
end

# More than two inputs
Base.vcat(c1::Union{Row,Table},c2::Union{Row,Table},c3::Union{Row,Table}) = vcat(vcat(c1,c2),c3)
Base.vcat(c1::Union{Row,Table},c2::Union{Row,Table},c3::Union{Row,Table},cs...) = vcat(vcat(c1,c2),c3,cs...)

# Concatenate columns and tables into tables
# Generated functions appear to be needed for speed...
Base.hcat{F,ElType,StorageType}(col::Column{F,ElType,StorageType}) = Table(col)
Base.hcat{Index,ElTypes,StorageTypes}(table::Table{Index,ElTypes,StorageTypes}) = table

@generated Base.hcat{F1,ElType1,Storage1,F2,ElType2,Storage2}(col1::Column{F1,ElType1,Storage1},col2::Column{F2,ElType2,Storage2}) = :( Table{$(FieldIndex{(F1,F2)}()),$(Tuple{ElType1,ElType2}),$(Tuple{Storage1,Storage2})}((col1.data,col2.data)) )
@generated function Base.hcat{F1,ElType1,Storage1,Index2,ElTypes2,StorageTypes2}(col1::Column{F1,ElType1,Storage1},table2::Table{Index2,ElTypes2,StorageTypes2})
    Index = F1 + Index2
    ElTypes = eltypes(Index)
    StorageTypes = Tuple{Storage1 ,StorageTypes2.parameters...}

    :(Table{$Index,$ElTypes,$StorageTypes}((col1.data, table2.data...)) )
end
@generated function Base.hcat{Index1,ElTypes1,StorageTypes1,F2,ElType2,Storage2}(table1::Table{Index1,ElTypes1,StorageTypes1},col2::Column{F2,ElType2,Storage2})
    Index = Index1 + F2
    ElTypes = eltypes(Index)
    StorageTypes = Tuple{StorageTypes1.parameters..., Storage2}

    :(Table{$Index,$ElTypes,$StorageTypes}((table1.data..., col2.data)) )
end
@generated function Base.hcat{Index1,ElTypes1,StorageTypes1,Index2,ElTypes2,StorageTypes2}(table1::Table{Index1,ElTypes1,StorageTypes1},table2::Table{Index2,ElTypes2,StorageTypes2})
    Index = Index1 + Index2
    ElTypes = eltypes(Index)
    StorageTypes = Tuple{StorageTypes1.parameters..., StorageTypes2.parameters...}

    :(Table{$Index,$ElTypes,$StorageTypes}((table1.data..., table2.data...)) )
end

# Strangely I get dispatch problems for Column->Table but not Cell->Row. Solution is to make sure this is at least as spececific in its templating as those above...
@inline Base.hcat{F1,ElType1,Storage1,F2,ElType2,Storage2}(in1::Column{F1,ElType1,Storage1},in2::Column{F2,ElType2,Storage2},ins::Union{Column,Table}...) = hcat(hcat(in1,in2),ins...)
@inline Base.hcat{F1,ElType1,Storage1,Index2,ElTypes2,StorageTypes2}(in1::Column{F1,ElType1,Storage1},in2::Table{Index2,ElTypes2,StorageTypes2},ins::Union{Column,Table}...) = hcat(hcat(in1,in2),ins...)
@inline Base.hcat{Index1,ElTypes1,StorageTypes1,F2,ElType2,Storage2}(in1::Table{Index1,ElTypes1,StorageTypes1},in2::Column{F2,ElType2,Storage2},ins::Union{Column,Table}...) = hcat(hcat(in1,in2),ins...)
@inline Base.hcat{Index1,ElTypes1,StorageTypes1,Index2,ElTypes2,StorageTypes2}(in1::Table{Index1,ElTypes1,StorageTypes1},in2::Table{Index2,ElTypes2,StorageTypes2},ins::Union{Column,Table}...) = hcat(hcat(in1,in2),ins...)

macro table(exprs...)
    N = length(exprs)
    field = Vector{Any}(N)
    value = Vector{Any}(N)
    for i = 1:N
        expr = exprs[i]
        if expr.head != :(=) && expr.head != :(kw) # strange Julia bug, see issue 7669
            error("A Expecting expression like @table(name::Type = value) or @table(field = value)")
        end
        if isa(expr.args[1],Symbol)
            field[i] = expr.args[1]
        elseif isa(expr.args[1],Expr)
            if expr.args[1].head != :(::) || length(expr.args[1].args) != 2
                field[i] = expr.args[1]
            else
                field[i] = :(TypedTables.Field{$(Expr(:quote,expr.args[1].args[1])),$(expr.args[1].args[2])}())
            end
        else
            error("C Expecting expression like @table(name::Type = value) or @table(field = value)")
        end
        value[i] = expr.args[2]
    end

    fields = Expr(:tuple,field...)
    values = Expr(:tuple,value...)

    return :(TypedTables.Table(TypedTables.FieldIndex{$(esc(fields))}(),$(esc(values))))
end

"This is a fake 'storage container' for the field DefaultKey() in a Table"
immutable TableKey{Index,ElTypes,StorageTypes}
    parent::Ref{Table{Index,ElTypes,StorageTypes}}
end
Base.getindex{Index,ElTypes,StorageTypes,T}(k::TableKey{Index,ElTypes,StorageTypes},i::T) = getindex(1:length(k.parent.x),i)
Base.length{Index,ElTypes,StorageTypes}(k::TableKey{Index,ElTypes,StorageTypes}) = length(k.parent.x)
Base.endof{Index,ElTypes,StorageTypes}(k::TableKey{Index,ElTypes,StorageTypes}) = endof(k.parent.x)
Base.eltype(k::TableKey) = Int
Base.eltype{Index,ElTypes,StorageTypes}(k::Type{TableKey{Index,ElTypes,StorageTypes}}) = Int
Base.first{Index,ElTypes,StorageTypes}(::TableKey{Index,ElTypes,StorageTypes},i::Int) = 1
Base.next{Index,ElTypes,StorageTypes}(::TableKey{Index,ElTypes,StorageTypes},i::Int) = (i, i+1)
Base.done{Index,ElTypes,StorageTypes}(k::TableKey{Index,ElTypes,StorageTypes},i::Int) = (i-1) == length(k.parent.x)
Base.show(io::IO,k::TableKey) = show(io::IO,1:length(k.parent.x))
Base.copy(k::TableKey) = 1:length(k)

=#


text:
using PointsReconstruction
using Test
using Distances

@testset "summary characteristics" begin include("summary_characteristics.jl") end


text:
export build_conv

function _compute_kernel_size(receptive_field, points_per_unit, num_layers)
    receptive_points = receptive_field * points_per_unit
    return ceil_odd(1 + (receptive_points - 1) / num_layers)
end

_compute_padding(kernel_size) = Integer(floor(kernel_size / 2))

_glorot_normal(dims...) = randn(Float32, dims...) .* sqrt(2f0 / sum(dims[1:2]) / dims[3])

_init_conv_fixed_bias(k, ch) =
    (_glorot_normal(k..., ch...), fill(1f-3, ch[2]))
_init_depthwiseconv_fixed_bias(k, ch) =
    (_glorot_normal(k..., div(ch[2], ch[1]), ch[1]), fill(1f-3, ch[2]))

_init_conv_random_bias(k, ch) =
    (_glorot_normal(k..., ch...), 1f-3 .* randn(Float32, ch[2]))
_init_depthwiseconv_random_bias(k, ch) =
    (_glorot_normal(k..., div(ch[2], ch[1]), ch[1]), 1f-3 .* randn(Float32, ch[2]))

_expand_kernel(n, ::Val{1}) = (n, 1)
_expand_kernel(n, ::Val{2}) = (n, n)

_expand_padding(n, ::Val{1}) = (n, 0)
_expand_padding(n, ::Val{2}) = (n, n)

"""
    build_conv(
        receptive_field::Float32,
        num_layers::Integer,
        num_channels::Integer;
        points_per_unit::Float32=30f0,
        multiple::Integer=1,
        num_in_channels::Integer=2,
        num_out_channels::Integer=2,
        dimensionality::Integer=1,
        init_conv::Function=_init_conv_fixed_bias,
        init_depthwiseconv::Function=_init_depthwiseconv_fixed_bias,
        act=x -> leakyrelu(x, 0.1f0)
    )

Build a CNN with a specified receptive field size.

# Arguments
- `receptive_field::Float32`: Width of the receptive field.
- `num_layers::Integer`: Number of layers of the CNN, excluding an initial
    and final pointwise convolutional layer to change the number of channels
    appropriately.
- `num_channels::Integer`: Number of channels of the CNN.

# Keywords
- `points_per_unit::Float32=30f0`: Points per unit for the discretisation. See
     `UniformDiscretisation1D`.
- `multiple::Integer=1`: Multiple for the discretisation. See `UniformDiscretisation1D`.
- `num_in_channels::Integer=2`: Number of input channels.
- `num_out_channels::Integer=2`: Number of output channels.
- `dimensionality::Integer=1`: Dimensionality of the filters.
- `init_conv::Function=_init_conv_fixed_bias`: Initialiser for dense convolutions.
- `init_depthwiseconv::Function=_init_depthwiseconv_fixed_bias`: Initialiser for depthwise
    separable convolutions.
- `act=x -> leakyrelu(x, 0.1f0)`: Activation function to use.

# Returns
- `BatchedConv`: CNN.
"""
function build_conv(
    receptive_field::Float32,
    num_layers::Integer,
    num_channels::Integer;
    points_per_unit::Float32=30f0,
    multiple::Integer=1,
    num_in_channels::Integer=2,
    num_out_channels::Integer=2,
    dimensionality::Integer=1,
    init_conv::Function=_init_conv_fixed_bias,
    init_depthwiseconv::Function=_init_depthwiseconv_fixed_bias,
    act=x -> leakyrelu(x, 0.1f0)
)
    # Appropriate expand the kernels to the right dimensionality.
    kernel = _expand_kernel(
        _compute_kernel_size(receptive_field, points_per_unit, num_layers),
        Val(dimensionality)
    )
    padding = _expand_padding(_compute_padding(kernel[1]), Val(dimensionality))

    # Build layers of the conv net.
    layers = Any[Conv(init_conv((1, 1), num_in_channels=>num_channels)..., act)]
    for i = 1:num_layers
        push!(layers, DepthwiseConv(
            init_depthwiseconv(kernel, num_channels=>num_channels)...,
            pad=padding
        ))
        push!(layers, Conv(
            init_conv((1, 1), num_channels=>num_channels)...,
            act
        ))
    end
    push!(layers, Conv(init_conv((1, 1), num_channels=>num_out_channels)...))

    return BatchedConv(
        Chain(layers...),
        points_per_unit,
        multiple,
        dimensionality
    )
end

"""
    struct BatchedConv

A batched CNN bundled with contextual information.

# Fields
- `conv`: Batched CNN.
- `points_per_unit::Float32`: Points per unit for the discretisation. See
     `UniformDiscretisation1D`.
- `multiple::Integer`: Multiple for the discretisation. See `UniformDiscretisation1D`.
- `dimensionality::Integer`: Dimensionality of the filters.
"""
struct BatchedConv
    conv
    points_per_unit::Float32
    multiple::Integer
    dimensionality::Integer
end

@Flux.functor BatchedConv

(layer::BatchedConv)(x) =  layer(x, Val(layer.dimensionality))

function (layer::BatchedConv)(x, dimensionality::Val{1})
    x, back = to_rank(3, x)        # Compress batch dimensions.
    x = with_dummy(layer.conv, x)  # Dummy dimension required!
    return back(x)                 # Uncompress batch dimensions
end

function (layer::BatchedConv)(x, dimensionality::Val{2})
    x, back = to_rank(4, x)  # Compress batch dimensions.
    x = layer.conv(x)        # No dummy dimension required!
    return back(x)           # Uncompress batch dimensions
end


text:
import Downloads
import JSON
url="http://fanyi.youdao.com/translate?smartresult=dict&smartresult=rule"
data=Dict{String,String}(
    "from"=>"AUTO",
    "to"=>"AUTO",
    "smartresult"=>"dict",
    "client"=>"fanyideskweb",
    "salt"=>"15801391750396",
    "sign"=>"74bbb50b1bd6c62fbff24be5f3787e2f",
    "ts"=>"1580139175039",
    "bv"=>"e2a78ed30c66e16a857c5b6486a1d326",
    "doctype"=>"json",
    "version"=>"2.1",
    "keyfrom"=>"fanyi.web",
    "action"=>"FY_BY_CLICKBUTTION",
    "i"=>"Hello world!"
)
head=Dict{String,String}(
    "User-Agent"=>"Mozilla/5.0 (Windows NT 10.0; Win64; x64) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/79.0.3945.130 Safari/537.36"
)
while true
    i::String=readline()
    if i=="exit"
        break
    end
    data["i"]=i
    r=Downloads.request(url;method="POST",input=IOBuffer(JSON.json(data)),headers=head)
    dump(r)
end

max_stars_repo_name: Rratic/notes

max_issues_repo_name: Rratic/notes

max_forks_repo_name: Rratic/notes


text:
# Description: test script for robot planning algorithms for a robot of radious 1m
using Pkg
Pkg.activate(@__DIR__)
Pkg.instantiate()
using GRECO, Revise
using Plots, Random
using MeshCat, GeometryTypes, CoordinateTransformations
# use include("src/GRECO.jl") # use if module alejandro is modified

# Some functions provided by Dr. Hassan Poonawala for visualization (pathvisualization.jl)
function sampleh() # random homogenous matrix
  randvec=rand(1,3)*2*pi
  rotmat = [1 0.0 0.0; 0.0 cos(randvec[1]) -sin(randvec[1]);0.0 sin(randvec[1]) cos(randvec[1])] *[cos(randvec[1]) 0.0 -sin(randvec[1]); 0.0 1.0 0.0;sin(randvec[1]) 0.0 cos(randvec[1])]*[1 0.0 0.0; 0.0 cos(randvec[2]) -sin(randvec[2]);0.0 sin(randvec[2]) cos(randvec[2])]
  h = zeros(Float64,4,4)
  h[1:3,1:3] = rotmat;
  h[1:3,4] = randn(3,1)*100
  return h
end

function plot_cylinder(vis,c1,c2,radius,mat,name="")
    geom = Cylinder(Point3f0(c1),Point3f0(c2),convert(Float32,radius))
    setobject!(vis["cyl"][name],geom,MeshPhongMaterial(color=RGBA(1, 0, 0, 1.0)))
end

function plot_sphere(vis,c1,radius,mat,name="")
    geom = HyperSphere(Point3f0(c1), convert(Float32,radius))
    setobject!(vis["sph"][name],geom,mat)
end

function add_obs!(vis,obs) # visualize obstacles
    for i in 1:size(obs,1)

        if obs[i,5] == 0.0

          plot_sphere(vis,[obs[i,1],obs[i,2],obs[i,3]],obs[i,4],MeshPhongMaterial(color=RGBA(0, 0, 1, 1.0)),"sph$i")
        else
          plot_cylinder(vis,[obs[i,1],obs[i,2],0],[obs[i,1],obs[i,2],obs[i,3]],obs[i,4],MeshPhongMaterial(color=RGBA(0, 0, 1, 1.0)),"cyl_$i")
        end
    end
end

function draw_robot(vis,h,name="")  # visualize frame as quadrotor
    robot_obj = HyperRectangle(Vec(0.0,0.0,0.0),Vec(1.0,1.0,1.0))
    setobject!(vis["robot"][name],robot_obj,MeshPhongMaterial(color=RGBA(0, 0, 0, 1.0)));
    settransform!(vis["robot"][name], compose(Translation(h[1:3,4]),LinearMap(h[1:3,1:3])))
end

# Example selector: 1 to 3
SF=3;
if SF==1
    s=[0.0 0.0 0.0]';
    g=[20.0 20.0 5.0]';
    O=[10.0 10.0 5.0 2 0;
        5.0 5.0 8.0 3 0;
        14.0 7.0 10.0 2.0 0;
        6.0 16.0 9.0 3.0 0;
        10.0 10.0 3.0 2.0 1;
        18.0 18.0 12.0 1.0 1;
        17.0 19.0 12.0 1.0 1
        19.0 17.0 12.0 1.0 1]
elseif SF==2
    s=[0.0 0.0 0.0]';
    g=[0.0 0.0 23.0]';

    O=[5.0 5.0 10.0 1.0 1;
        5.0 0.0 10.0 1.0 1;
        5.0 -5.0 10.0 1.0 1;
        -5.0 5.0 10.0 1.0 1;
        -5.0 0.0 10.0 1.0 1;
        -5.0 -5.0 10.0 1.0 1;
        0.0 5.0 10.0 1.0 1;
        0.0 -5.0 10.0 1.0 1;
        0.0 0.0 15.0 5.0 0]
    elseif SF==3
        # Start node
        s=[0.0 0.0 0.0]';

        # End node
        g=[4.0 4.0 1.0]';

        # Obstacles: each row represents an obstacle
        O=[2 2 2 0.2 0;         # sphere
            2 2 1.5 0.5 1];      # cylinder
end

# Path with PRM
Frames,Cost=PRM(s,g,O);
traj_sol=Frames;

# Visualize everything: obs + solution path
vis = Visualizer()
open(vis)
# obstacles:
add_obs!(vis,O)
# visualize robot at each coordinate frame on path:
for i in 1:length(traj_sol)
  draw_robot(vis,traj_sol[i],"$i")
end

max_stars_repo_name: anniexen/GRECO

max_issues_repo_name: anniexen/GRECO

max_forks_repo_name: anniexen/GRECO


text:
using Fermi.ConfigurationInteraction: Determinant, excitation_level, αexcitation_level, βexcitation_level, αexclusive, βexclusive, phase, create, annihilate

"""
    Fermi.CoupledCluster.ecRCCSD

Fermi struct that holds information about ecRCCSD wavefunctions

# Fields

    CorrelationEnergy   CCSD correlation energy
    T1                  T1 amplitudes
    T2                  T2 amplitudes

_struct tree:_

**RCCSD** <: AbstractCCWavefunction <: AbstractCorrelatedWavefunction <: AbstractWavefunction
"""
struct ecRCCSD{T} <: AbstractCCWavefunction
    GuessEnergy::T
    CorrelationEnergy::T
    T1::_T where _T <: Fermi.AbstractTensor
    T2::_T where _T <: Fermi.AbstractTensor
end

# Most general function. It defines the precision and call a precision-specific function
#"""
#    Fermi.CoupledCluster.RCCSD()
#
#Compute a RCCSD wave function using Fermi.CurrentOptions data.
#"""
function ecRCCSD()
    @output "Selecting precision...\n"
    prec = select_precision(Fermi.CurrentOptions["precision"])
    ecRCCSD{prec}()
end

"""
    Fermi.CoupledCluster.RCCSD{T}()

Compute a RCCSD wave function for a given precision T (Float64 or Float32)
"""
function ecRCCSD{T}() where T <: AbstractFloat
    @output "Selecting Algorithm...\n"
    alg = select_algorithm(Fermi.CurrentOptions["cc_alg"])
    ecRCCSD{T}(alg)
end

# kernel (loop logic) from main RCCSD
include("../RCCSD/kernel.jl")

#implementations
include("CTF.jl")
include("ClusterDecomposition.jl")


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


text:
using FixedPointNumbers
using ColorTypes

struct TextureArray <: AbstractArray{RGBA,2}
    backing::Array{RGBA{N0f8},2}
end

Base.size(A::TextureArray)::Tuple{Int, Int} = (size(A.backing)[1]*2, floor(size(A.backing)[2]/2))
Base.getindex(A::TextureArray, I::Vararg{Int, N}) where N = A.backing[Int(ceil(I[1] / 2)), iseven(I[1]) ? I[2] + size(A)[2] : I[2]]
Base.setindex!(A::TextureArray, v::RGBA{N0f8}, I::Vararg{Int, 2}) = (A.backing[floor(I[1] / 2, iseven(I[1]) ? I[2] + size(A)[1] : I[2])] = v)
Base.IndexStyle(::Type{<:TextureArray}) = IndexCartesian()

gray(c) = RGBA{N0f8}(c, c, c, 1)

testtexture() = TextureArray(gray.(reshape(0:15, 4, 4) ./ 15))

max_stars_repo_name: npdmfixup/Fuji

max_issues_repo_name: npdmfixup/Fuji

max_forks_repo_name: npdmfixup/Fuji


text:
import Base: show

using Dates: AbstractTime, Nanosecond

"""
    Clock{T<:Number}

    Clock{T}(time, iteration)

Keeps track of the current `time` and `iteration` number. The `time::T` can be either a number of a `DateTime` object.
"""
mutable struct Clock{T}
         time :: T
    iteration :: Int
end

Clock(; time, iteration=0) = Clock(time, iteration)

Base.show(io::IO, c::Clock{FT}) where FT =
    println(io, "Clock{$FT}: time = ", prettytime(c.time), ", iteration = ", c.iteration)

function tick!(clock, Δt)
    clock.time += Δt
    clock.iteration += 1
    return nothing
end

function tick!(clock::Clock{<:AbstractTime}, Δt)
    clock.time += Nanosecond(round(Int, 1e9 * Δt))
    clock.iteration += 1
    return nothing
end

"Adapt `Clock` to work on the GPU via CUDAnative and CUDAdrv."
Adapt.adapt_structure(to, clock::Clock) = (time=clock.time, iteration=clock.iteration)


text:
import Base: push!, getindex, setindex!

using OrderedCollections: OrderedDict
using Oceananigans: AbstractOutputWriter, AbstractDiagnostic

#####
##### Utilities shared between diagnostics and output writers
#####

defaultname(::AbstractDiagnostic, nelems) = Symbol(:diag, nelems+1)
defaultname(::AbstractOutputWriter, nelems) = Symbol(:writer, nelems+1)

const DiagOrWriterDict = OrderedDict{S, <:Union{AbstractDiagnostic, AbstractOutputWriter}} where S

function push!(container::DiagOrWriterDict, elem)
    name = defaultname(elem, length(container))
    container[name] = elem
    return nothing
end

getindex(container::DiagOrWriterDict, inds::Integer...) = getindex(container.vals, inds...)
setindex!(container::DiagOrWriterDict, newvals, inds::Integer...) = setindex!(container.vals, newvals, inds...)

function push!(container::DiagOrWriterDict, elems...)
    for elem in elems
        push!(container, elem)
    end
    return nothing
end


text:
# Julia code
# Copyright (c) 2022 the SLICOTMath.jl developers
# Portions extracted from SLICOT-Reference distribution:
# Copyright (c) 2002-2020 NICONET e.V.

function run_mb03bz(datfile, io=stdout)
    NIN = 5
    NOUT = 6
    KMAX = 6
    NMAX = 50
    LDA1 = NMAX
    LDA2 = NMAX
    LDQ1 = NMAX
    LDQ2 = NMAX
    LDWORK = NMAX
    LZWORK = NMAX
    S = Array{BlasInt,1}(undef, KMAX)
    A = Array{ComplexF64,3}(undef, LDA1,LDA2,KMAX)
    Q = Array{ComplexF64,3}(undef, LDQ1,LDQ2,KMAX)
    ALPHA = Array{ComplexF64,1}(undef, NMAX)
    BETA = Array{ComplexF64,1}(undef, NMAX)
    SCAL = Array{BlasInt,1}(undef, NMAX)
    ZWORK = Array{ComplexF64,1}(undef, LZWORK)
    DWORK = Array{Float64,1}(undef, LDWORK)
    f = open(datfile,"r")
    readline(f)
    vs = split(readline(f))
    JOB =  vs[1][1]
    COMPQ =  vs[2][1]
    K = parse(BlasInt, vs[3])
    N = parse(BlasInt, vs[4])
    ILO = parse(BlasInt, vs[5])
    IHI = parse(BlasInt, vs[6])
if ( N<0 || N>NMAX )
else

    vs = String[]
    _isz = K
    while length(vs) < _isz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    S[1:_isz] .= parsex.(BlasInt, vs)

    vs = String[]
    _isz,_jsz,_ksz = (N,N,K)
    while length(vs) < _isz*_jsz*_ksz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for k in 1:_ksz
      for i in 1:_isz
        _i0 = (i-1)*_jsz + (k-1)*_jsz*_isz
        A[i,1:_jsz,k] .= parsex.(ComplexF64, vs[_i0+1:_i0+_jsz])
      end
    end
# interp call 1

    INFO = SLICOT.mb03bz!(JOB, COMPQ, K, N, ILO, IHI, S, A, Q, ALPHA, BETA, SCAL, LDWORK, LZWORK)
    @test INFO == 0
    INFO == 0 || return

if ( INFO!=0 )
else
if ( LSAME( JOB, 'S' ) )
# interp output 1
    println(io,"A:")
    _nc = N
    _nr = N
    _nk = K
    show(io,"text/plain",A[1:_nr,1:_nc,1:_nk])
    println(io,)

end # if
if ( !LSAME( COMPQ, 'N' ) )
# interp output 2
    println(io,"Q:")
    _nc = N
    _nr = N
    _nk = K
    show(io,"text/plain",Q[1:_nr,1:_nc,1:_nk])
    println(io,)

end # if
# interp output 3
    println(io,"ALPHA:")
    _nr = N
    show(io,"text/plain",ALPHA[1:_nr])
    println(io,)

# interp output 4
    println(io,"BETA:")
    _nr = N
    show(io,"text/plain",BETA[1:_nr])
    println(io,)

# interp output 5
    println(io,"SCAL:")
    _nr = N
    show(io,"text/plain",SCAL[1:_nr])
    println(io,)

end # if
end # if
    close(f)
end # run_X()


text:
using Test
using Dates

##################################
## Bruteforce version
##################################

function constructSet(line)
    path = []
    x = 0
    y = 0
    instructions = split(line, ",")
    for inst in instructions
        dir = inst[1]
        number = parse(Int, inst[2:end])

        if dir == 'R'
            positions = map(i -> (x, y + i), 1:number)
        elseif dir == 'L'
            positions = map(i -> (x, y - i), 1:number)
        elseif dir == 'U'
            positions = map(i -> (x + i , y), 1:number)
        elseif dir == 'D'
            positions = map(i -> (x - i , y), 1:number)
        end
        x = positions[end][1]
        y = positions[end][2]
        append!(path, positions)
    end
    return path
end

function run(input)
    input = split(input, "\n")
    set1 = constructSet(input[1])
    set2 = constructSet(input[2])
    inter = intersect(Set(set1), Set(set2))
    return  minimum(pos -> abs(pos[1]) + abs(pos[2]), inter)
end

##############################################################
### Another faster version
### We don't store all the points, but only the
### vertices (begin, end) and compute their intersections.
##############################################################

function getVertices(line)
    vertices = []
    (x, y) = (0, 0)
    instructions = split(line, ",")
    for inst in instructions
        dir = inst[1]
        steps = parse(Int, inst[2:end])
        if dir == 'R'
            append!(vertices, [((x, y), (x, y+steps))])
            y = y + steps
        elseif dir == 'L'
            append!(vertices,[((x, y-steps), (x, y))])
            y = y - steps
        elseif dir == 'U'
            append!(vertices, [((x, y), (x + steps, y))])
            x = x + steps
        elseif dir == 'D'
            append!(vertices, [((x - steps, y), (x, y))])
            x = x - steps
        end
    end
    return vertices
end

function getDir(ver)
    (x1, y1), (x2, y2) = ver
    if x1 == x2
        return "LR"
    elseif y1 == y2
        return "UD"
    end
end

@test getDir(((1, 2), (1, 3))) == "LR"
@test getDir(((1, 2), (3, 2))) == "UD"

function getIntersectionVertices(ver1, ver2)
    ((x1, y1), (x2, y2)) = ver1
    ((x3, y3), (x4, y4)) = ver2
    dir1 = getDir(ver1)
    dir2 = getDir(ver2)
    if dir1 == dir2 == "LR"
        if x1 == x3
            return [(x1, y) for y in max(y1, y3):min(y2, y4)]
        else
            return []
        end
    elseif dir1 == dir2 == "UD"
        if y1 == y3
            return [(x, y1) for x in max(x1, x3):min(x2, x4)]
        else
            return []
        end
    elseif dir1 == "UD" && dir2 == "LR"
        if x1 <= x3 <= x2 && y3 <= y1 <= y4
            return [(x3, y1)]
        else
            return []
        end
    elseif dir1 == "LR" && dir2 == "UD"
        # invert the two
        ((x1, y1), (x2, y2)) = ver2
        ((x3, y3), (x4, y4)) = ver1
        if x1 <= x3 <= x2 && y3 <= y1 <= y4
            return [(x3, y1)]
        else
            return []
        end
    end
end

@test getIntersectionVertices(
    ((1, 2), (1, 3)),
    ((1, 2), (1, 4))
) == [(1, 2), (1, 3)]

@test getIntersectionVertices(
    ((1, 2), (1, 7)),
    ((1, 4), (1, 20))
) == [(1, 4), (1, 5), (1, 6),  (1, 7)]

@test getIntersectionVertices(
    ((2, 1), (7, 1)),
    ((4, 1), (20, 1))
) == [(4, 1), (5, 1), (6, 1),  (7, 1)]

@test getIntersectionVertices(
    ((2, 5), (2, 10)),
    ((1, 6), (5, 6))
) == [(2, 6)]

@test getIntersectionVertices(
    ((1, 6), (5, 6)),
    ((2, 5), (2, 10))
) == [(2, 6)]

function getIntersections(vertices1, vertices2)
    intersections = []
    for ver1 in vertices1
        for ver2 in vertices2
            inter = getIntersectionVertices(ver1, ver2)
            append!(intersections, inter)
        end
    end
    return intersections
end

function run2(input)
    input = split(input, "\n")
    vertices1 = getVertices(input[1])
    vertices2 = getVertices(input[2])
    inter = getIntersections(vertices1, vertices2)
    inter = filter(x -> x ≠ (0, 0), inter)
    return minimum(pos -> abs(pos[1]) + abs(pos[2]), inter)
end

@test run2("R75,D30,R83,U83,L12,D49,R71,U7,L72\nU62,R66,U55,R34,D71,R55,D58,R83") == 159
@test run2("R98,U47,R26,D63,R33,U87,L62,D20,R33,U53,R51\nU98,R91,D20,R16,D67,R40,U7,R15,U6,R7") == 135

#########################################

function main()
    input = ARGS[1]
    # input = split(input, "\n")
    res, time, memory = @timed run2(input)
    println("_duration:$(time*1000)")
    println(res)
end

main()


text:
using CSV
using FilePaths;
using FilePathsBase: /
"""
    read_graph(dir::String, name::String, sparse_size_threshold = 10000)

    A utility to parse input csv files into a RedBlackGraph. The first argument is the directory that holds the csv files.
    The second is the basename for the csv files. By convention, there are two csv files <basename>.canonical.edges.csv and
    <basename>.canonical.vertices.csv.

    <basename>.canonical.edges.csv may have more than two columns, but the first two are significant for these utilities.
    Column 1 is the source vetex number, Column 2 is the destination vertex number.

    <basename>.canonical.vertices.csv may have more than 2 columns, but two are significant for these utilites. It must contain
    a Column labeled "vertex_number" and one labeled "color". Both columns contain integers.

    This utility will return either a Matrix{AInt64} or a SparseMatrixCSC{AInt64, Int64} if the number of vertices in the graph
    is greater than the value of sparse_size_threshold.

    (Note: https://github.com/rappdw/fs-crawler can be used to generate the required csv files.)
"""
function read_graph(dir::String, name::String, sparse_size_threshold = 10000)
    dir = Path(dir)
    vertex_rows = CSV.Rows(dir / "$name.canonical.vertices.csv"; select=["vertex_number"], reusebuffer=true, types=Dict("vertex_number"=>Int64))

    lastarg(_,x) = x
    n = foldl(lastarg, vertex_rows)[1] # got this from the discussion https://discourse.julialang.org/t/getting-the-last-element-of-an-iterator/49696/44
    use_spare = n > 10000
    graph = nothing
    r = red_one(AInt64)
    if use_spare
        # figure this one out later...
    else
        graph = zeros(AInt64, n, n)
    end
    vertex_rows = CSV.Rows(dir / "$name.canonical.vertices.csv"; select=["vertex_number" ,"color"], reusebuffer=true, types=Dict("vertex_number"=>Int64, "color"=>Int64))
    for vertex in vertex_rows
        # get the vertex number (first column) and color (second column)
        i = vertex[1]
        val = r                 # TODO: support AInt64(-1) and then just use that for value, e.g. graph[i, i] = AInt64(vertex[2])
        if vertex[2] == 1
            val = 1
        end
        graph[i, i] = val
    end
    edge_rows = CSV.Rows(dir / "$name.canonical.edges.csv"; types=Dict("Column1"=>Int64, "Column2"=>Int64), header=0)
    for edge in edge_rows
        # get the source (first column) and destination (second column)
        i = edge[1]
        j = edge[2]
        # set the relationship to 2 (father), 3 (mother) based on color of vertex
        val = 2
        if graph[j, j] == 1
            val = 3
        end
        graph[i, j] = val
    end
    return graph
end


text:
include(joinpath("conjugates", "fallbacks.jl"))
include(joinpath("conjugates", "beta_binom.jl"))
include(joinpath("conjugates", "dirichlet_multi.jl"))
include(joinpath("conjugates", "gamma_exp.jl"))

include(joinpath("conjugates", "normalgamma.jl"))
include(joinpath("conjugates", "normalinversegamma.jl"))
include(joinpath("conjugates", "normalwishart.jl"))
include(joinpath("conjugates", "normalinversewishart.jl"))
include(joinpath("conjugates", "normal.jl"))
include(joinpath("conjugates", "mvnormal.jl"))


text:
# Beta marginal distribution

"""
	$(SIGNATURES)

Constructor for `BetaMarginal`.
"""
function BetaMarginal(lb :: T1, ub :: T1, alpha :: T1, beta :: T1) where
    T1 <: AbstractFloat

    betaDistr = Distributions.Beta(alpha, beta);
    return BetaMarginal(lb, ub, betaDistr);
end

draw_test_endowments(m :: BetaMarginal{T1}, n, rng :: AbstractRNG) where T1 =
    lb(m) .+ (ub(m) .- lb(m)) .* rand(rng, m.betaDistr, n);

validate_draws(m :: BetaMarginal{T1}, x :: AbstractVector{T1})  where T1 =
    validate_bounded_draws(m, x);

marginal_quantile(m :: BetaMarginal{T1}, pct) where T1 =
    lb(m) .+ (ub(m) .- lb(m)) .* quantile.(m.betaDistr, pct);

# ----------


text:
# using LazyArrays
# const AbstractOuterProductArray{T,N} = BroadcastArray{T,N,typeof(*),I} where {T,N,I}

# OuterProductArray(arrays::Vararg{<:AbstractVector{T},N}) where {T,N} =
#                  BroadcastArray(*,ntuple(k->reshape(arrays[k],ntuple(k1-> k==k1 ? length(arrays[k]) : 1,Val(N)))  ,Val(N))...)

tensorproduct(v::AbstractVector...)  =
    OuterProductArray(v...)
tensorproduct(v::GridArrays.AbstractIntervalGrid...)  =
    ProductGrid(v...)
abstract type AbstractOuterProductArray{T,N} <: AbstractArray{T,N} end

components(opa::AbstractOuterProductArray) = opa.vectors
size(A::AbstractOuterProductArray) = map(length, components(A))
axes(A::AbstractOuterProductArray) = map(x->axes(x,1), components(A))

@propagate_inbounds function Base.getindex(F::AbstractOuterProductArray, k::Integer)
    @boundscheck checkbounds(F, k)
    unsafe_array_getindex(F, (k,))
end

@propagate_inbounds function Base.getindex(F::AbstractOuterProductArray{T, N}, kj::Vararg{<:Integer, N}) where {T, N}
    @boundscheck checkbounds(F, kj...)
    unsafe_array_getindex(F, kj)
end

unsafe_array_getindex(F::AbstractOuterProductArray{T,N}, kj::NTuple{N,Int}) where {T,N} =
    prod(map(getindex, components(F), kj))

struct OuterProductArray{T,N,V} <: AbstractOuterProductArray{T,N}
    vectors   ::  V
    OuterProductArray(v::AbstractVector...) =
        new{promote_type(map(eltype,v)...),length(v),typeof(v)}(v)
end

struct OPArrayStyle{N} <: Broadcast.AbstractArrayStyle{N} end
Broadcast.BroadcastStyle(::Type{<:OuterProductArray{T,N}}) where {T,N} = OPArrayStyle{N}()
Base.copyto!(dest::OuterProductArray, bc::Broadcast.Broadcasted{<:OPArrayStyle}) =
    OuterProductArray(map(argsi->broadcast(bc.f, argsi...), zip(map(components, bc.args)...))...)
Base.copy(a::OuterProductArray) = OuterProductArray(map(copy, components(a))...)
Base.similar(a::OuterProductArray) = OuterProductArray(map(similar, components(a))...)
Base.similar(::Type{OuterProductArray{T}}, shape::Tuple{Union{Integer, Base.OneTo},Vararg{Union{Integer,Base.OneTo},N} where N}) where T =
    OuterProductArray(map(x->Base.similar(Vector{T}, x), shape)...)
Base.similar(bc::Broadcast.Broadcasted{OPArrayStyle{N}},::Type{ElType}) where {N,ElType} =
    Base.similar(OuterProductArray{ElType}, axes(bc))
function Base.copyto!(dest::OuterProductArray, src::OuterProductArray)
    for (deste,srce) in zip(components(dest),components(src))
        copy!(deste,srce)
    end
end


text:
module layers

using Flux
using Flux: glorot_uniform, @functor#, destructure
using Zygote: @adjoint, @nograd
using LinearAlgebra, SparseArrays
using Statistics
using SimpleWeightedGraphs
using ChemistryFeaturization
#using DifferentialEquations, DiffEqSensitivity

# regularized norm fcn, cut out the dims part
function reg_norm(x::AbstractArray, ϵ=sqrt(eps(Float32)))
    μ′ = mean(x)
    σ′ = std(x, mean = μ′, corrected=false)
    return Float32.((x .- μ′) ./ (σ′ + ϵ))
end

"""
    AGNConv(in=>out)

Atomic graph convolutional layer. Almost identical to GCNConv from GeometricFlux but adapted to be most similar to Tian's original AGNN structure, so explicitly has self and convolutional weights separately.

# Fields
- `selfweight::Array{T,2}`: weights applied to features at a node
- `convweight::Array{T,2}`: convolutional weights
- `bias::Array{T,2}`: additive bias (second dimension is always 1 because only learnable per-feature, not per-node)
- `σ::F`: activation function (will be applied before `reg_norm` to outputs), defaults to softplus

# Arguments
- `in::Integer`: the dimension of input features.
- `out::Integer`: the dimension of output features.
- `σ=softplus`: activation function
- `initW=glorot_uniform`: initialization function for weights
- `initb=zeros`: initialization function for biases

"""
struct AGNConv{T,F}
    selfweight::Array{T,2}
    convweight::Array{T,2}
    bias::Array{T,2}
    σ::F
end

function AGNConv(ch::Pair{<:Integer,<:Integer}, σ=softplus; initW=glorot_uniform, initb=zeros, T::DataType=Float32)
    selfweight = T.(initW(ch[2], ch[1]))
    convweight = T.(initW(ch[2], ch[1]))
    b = T.(initb(ch[2], 1))
    AGNConv(selfweight, convweight, b, σ)
end

@functor AGNConv

"""
 Define action of layer on inputs: do a graph convolution, add this (weighted by convolutional weight) to the features themselves (weighted by self weight) and the per-feature bias (concatenated to match number of nodes in graph).

# Arguments
- input: AtomGraph object
"""
function (l::AGNConv)(ag::AtomGraph)
    lapl = ag.lapl
    X = ag.features
    out_mat = Float32.(reg_norm(l.σ.(l.convweight * X * lapl + l.selfweight * X + reduce(hcat,l.bias for i in 1:size(X, 2)))))
    AtomGraph(ag.graph, ag.elements, ag.lapl, out_mat, AtomFeat[], ag.id)
end

# fixes from Dhairya so backprop works
@adjoint function SparseMatrixCSC{T,N}(arr) where {T,N}
  SparseMatrixCSC{T,N}(arr), Δ -> (collect(Δ),)
end
@nograd LinearAlgebra.diagm

@adjoint function Broadcast.broadcasted(Float32, a::SparseMatrixCSC{T,N}) where {T,N}
  Float32.(a), Δ -> (nothing, T.(Δ), )
end
@nograd issymmetric

@adjoint function softplus(x::Real)
  y = softplus(x)
  return y, Δ -> (Δ * σ(x),)
end

"""
Custom pooling layer that outputs a fixed-length feature vector irrespective of input dimensions, for consistent handling of different-sized graphs feeding to fully-connected dense layers afterwards. Adapted from Flux's MeanPool.

It accepts a pooling width and will adjust stride and/or padding such that the output vector length is correct.
"""
struct AGNPool
    pool_func::Function
    dim::Int64
    str::Int64
    pad::Int64
    function AGNPool(pool_type::String, in_num_features::Int64, out_num_features::Int64, pool_width_frac::Float64)
        @assert in_num_features >= out_num_features "I don't think you actually want to pool to a LONGER vector, do you?"
        dim, str, pad = compute_pool_params(in_num_features, out_num_features, Float32(pool_width_frac))
        if pool_type=="max"
            pool_func = Flux.maxpool
        elseif pool_type=="mean"
            pool_func = Flux.meanpool
        end
        new(pool_func, dim, str, pad)
    end
end

pool_out_features(num_f::Int64, dim::Int64, stride::Int64, pad::Int64) = Int64(floor((num_f + 2 * pad - dim) / stride + 1))

"""
Helper function to work out dim, pad, and stride for desired number of output features, given a fixed pooling width.
"""
function compute_pool_params(num_f_in::Int64, num_f_out::Int64, dim_frac::Float32; start_dim=Int64(round(dim_frac*num_f_in)), start_str=Int64(floor(num_f_in/num_f_out)))
    # take starting guesses
    dim = start_dim
    str = start_str
    p_numer = str*(num_f_out-1) - (num_f_in - dim)
    if p_numer < 0
        p_numer == -1 ? dim = dim + 1 : str = str + 1
    end
    p_numer = str*(num_f_out-1) - (num_f_in - dim)
    if p_numer < 0
        error("problem, negative p!")
    end
    if p_numer % 2 == 0
        pad = Int64(p_numer/2)
    else
        dim = dim - 1
        pad = Int64((str*(num_f_out-1) - (num_f_in - dim))/2)
    end
    out_fea_len = pool_out_features(num_f_in, dim, str, pad)
    if !(out_fea_len==num_f_out)
        print("problem, output feature wrong length!")
    end
    # check if pad gets comparable to width...
    if pad >= 0.8*dim
        @warn "specified pooling width was hard to satisfy without nonsensically large padding relative to width, had to increase from desired width"
        dim, str, pad  = compute_pool_params(num_f_in, num_f_out, dim_frac, start_dim=Int64(round(1.2*start_dim)))
    end
    dim, str, pad
end

function (m::AGNPool)(ag::AtomGraph)
      # compute what pad and stride need to be...
      x = ag.features
      x = reshape(x, (size(x)..., 1, 1))
      # do mean pooling across feature direction, average across all nodes in graph
      # TODO: decide if this approach makes sense or if there's a smarter way
      pdims = PoolDims(x, (m.dim,1); padding=(m.pad,0), stride=(m.str,1))
      mean(m.pool_func(x, pdims), dims=2)[:,:,1,1]
end

# following commented out for now because it only runs suuuuper slowly but slows down precompilation a lot
"""
# DEQ-style model where we treat the convolution as a SteadyStateProblem
struct AGNConvDEQ{T,F}
    conv::AGNConv{T,F}
end

function AGNConvDEQ(ch::Pair{<:Integer,<:Integer}, σ=softplus; initW=glorot_uniform, initb=glorot_uniform, T::DataType=Float32, bias::Bool=true)
    conv = AGNConv(ch, σ; initW=initW, initb=initb, T=T)
    AGNConvDEQ(conv)
end

@functor AGNConvDEQ

# set up SteadyStateProblem where the derivative is the convolution operation
# (we want the "fixed point" of the convolution)
# need it in the form f(u,p,t) (but t doesn't matter)
# u is the features, p is the parameters of conv
# re(p) reconstructs the convolution with new parameters p
function (l::AGNConvDEQ)(gr::AtomGraph)
    p,re = Flux.destructure(l.conv)
    # do one convolution to get initial guess
    guess = l.conv(gr).features

    f = function (dfeat,feat,p,t)
        input = gr
        input.features = reshape(feat,size(guess))
        output = re(p)(input)
        dfeat .= vec(output.features) .- vec(input.features)
    end

    prob = SteadyStateProblem{true}(f, vec(guess), p)
    #return solve(prob, DynamicSS(Tsit5())).u
    alg = SSRootfind()
    #alg = SSRootfind(nlsolve = (f,u0,abstol) -> (res=SteadyStateDiffEq.NLsolve.nlsolve(f,u0,autodiff=:forward,method=:anderson,iterations=Int(1e6),ftol=abstol);res.zero))
    out_mat = reshape(solve(prob, alg, sensealg = SteadyStateAdjoint(autodiff = false, autojacvec = ZygoteVJP())).u,size(guess))
    return AtomGraph(gr.graph, gr.elements, out_mat, gr.featurization)
end
"""

end