text:
"""
Created on Thu Apr 09 2020

@author: Yoann Pradat

α and β divergences.
"""

"""
β_divergence(X, Y, β[, ϵ])

Compute the the β-divergence of X and Y, d_β(X||Y).

# Arguments
- `X::Matrix{<:Real}`: array
- `Y::Matrix{<:Real}`: array
- `β::Real`: real, parameter of the β-divergence.
    If β == 2, this is half the Frobenius *squared* norm.
    If β == 1, this is the Kullback-Leibler divergence.
    If β == 0, this is the Itakura-Saito divergence.
    Else, this is the general β-divergence.
"""
function β_divergence(X::Matrix{T}, Y::Matrix{T}, β::Real; ϵ::Float64=eps(T)) where T <: Real
    if β == 2
        #### Frobenius norm
        res = norm(X - Y, 2) ^ 2 / 2.
    else
        # vec() flattens the array column-wise
        Y_vec = vec(Y)
        X_vec = vec(X)

        if β == 1
            # drop 0 coordinates in X that contribute 0 to the divergence as (x -> xlog(x) continuous in 0)
            indices = X_vec .> ϵ
            Y_vec = Y_vec[indices]
            X_vec = X_vec[indices]

            #### Kullback-Leibler divergence
            sum_X_log_X_Y = sum(X_vec .* log.(X_vec ./ Y_vec))
            sum_X         = sum(X_vec)
            sum_Y         = sum(Y_vec)
            res           = abs(sum_X_log_X_Y - sum_X + sum_Y)

        elseif β == 0
            sum(X_vec .== 0) == 0 || throw(ArgumentError("for β=0, β-divergence is undefined if X has null entries"))
            sum(Y_vec .== 0) == 0 || throw(ArgumentError("for β=0, β-divergence is undefined if Y has null entries"))

            #### Itakura-Saito divergence
            X_Y         = X_vec ./ Y_vec
            sum_log_X_Y = sum(log.(X_Y))
            sum_X_Y_m1  = sum(X_Y .- 1)
            res         = abs(-sum_log_X_Y + sum_X_Y_m1)

        else
            if β > 0 && β < 1
                sum(Y_vec .== 0) == 0 || throw(ArgumentError("for 0<β<1, β-divergence is undefined if Y has null entries"))
                # drop 0 coordinates in X that contribute 0 to the divergence as (x -> xlog(x) continuous in 0)
                indices = X_vec .> ϵ
                Y_vec = Y_vec[indices]
                X_vec = X_vec[indices]
            elseif β < 0
                sum(X_vec .== 0) == 0 || throw(ArgumentError("for β<0, β-divergence is undefined if X has null entries"))
                sum(Y_vec .== 0) == 0 || throw(ArgumentError("for β<0, β-divergence is undefined if Y has null entries"))
            end

            #### General β-divergence for β not in (0, 1, 2)
            sum_X_β      = sum(X_vec .^ β)
            sum_Y_β      = sum(Y_vec .^ β)
            sum_X_Y_β_m1 = sum(X_vec .* Y_vec .^ (β-1))
            res          = abs((sum_X_β + (β-1) * sum_Y_β - β * sum_X_Y_β_m1)/(β*(β-1)))
        end
    end
    res
end

"""
α_divergence(X, Y, α[, ϵ])

Compute the the α-divergence of X and Y, d_α(X||Y).

# Arguments
- `X::Matrix{<:Real}`: array
- `Y::Matrix{<:Real}`: array
- `α::Real`: float, parameter of the α-divergence.
    If α == 2, this is Pearsoni's Chi-squared distance
    If α == 1, this is (at the limit) the Kullback-Leibler divergence
    If α == 0, this is (at the limit) the dual Kullback-Leibler divergence
    If α == 0.5, this is the squared Hellinger distance
    If α == -1, this is the Neyman's Chi-squared distance
    Else, this is the general α-divergence.
"""
function α_divergence(X::Matrix{T}, Y::Matrix{T}, α::Real; ϵ::Float64=eps(T)) where T <: Real
    # vec() flattens the array column-wise
    Y_vec = vec(Y)
    X_vec = vec(X)

    if α == 1
        # avoid division by zero
        Y_vec[Y_vec .== 0] .= ϵ

        #### Kullback-Leibler divergence
        sum_X_log_X_Y = sum(X_vec .* log.(X_vec ./ Y_vec))
        sum_X         = sum(X_vec)
        sum_Y         = sum(Y_vec)
        res           = abs(sum_X_log_X_Y - sum_X + sum_Y)
    elseif α == 0
        # avoid division by zero
        X_vec[X_vec .== 0] .= ϵ

        #### dual Kullback-Leibler divergence
        sum_Y_log_Y_X = sum(Y_vec .* log.(Y_vec ./ X_vec))
        sum_Y         = sum(Y_vec)
        sum_X         = sum(X_vec)
        res           = abs(sum_Y_log_Y_X - sum_Y + sum_X)
    else
        if α < 1
            # avoid division by zero
            Y_vec[Y_vec .== 0] .= ϵ
        end

        if α < 0
            # avoid division by zero
            X_vec[X_vec .== 0] .= ϵ
        end

        #### general α-divergence for α not in {0,1}
        sum_X_α_Y_α_m1 = sum(X_vec .^ α .* Y_vec .^ (1-α))
        sum_X          = sum(X_vec)
        sum_Y          = sum(Y_vec)
        res            = abs((sum_X_α_Y_α_m1 - α * sum_X + (α-1) * sum_Y)/(α*(α-1)))
    end
    res
end


text:
# ------------------------------------------------------------------------------------------
# # Getting to know Julia
#
#
# This notebook is meant to offer a crash course in Julia syntax to show you that Julia is
# lightweight and easy to use -- like your favorite high-level language!
#
# We'll talk about
# - Strings
# - Data structures
# - Loops
# - Conditionals
# - Functions
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## Strings
# ------------------------------------------------------------------------------------------

string1 = "How many cats "

string2 = "is too many cats?"

string(string1, string2)

😺 = 10
println("I don't know but $😺 are too few!")

# ------------------------------------------------------------------------------------------
# Note: Julia allows us to write super generic code, and 😺 is an example of this.
#
# This allows us to write code like
# ------------------------------------------------------------------------------------------

😺 = 1
😀 = 0
😞 = -1

😺 + 😞 == 😀

# ------------------------------------------------------------------------------------------
# ## Data structures
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ### Tuples
#
# We can create a tuple by enclosing an ordered collection of elements in `( )`.
#
# Syntax:

# ```julia
# (item1, item2, ...)```
# ------------------------------------------------------------------------------------------

myfavoriteanimals = ("penguins", "cats", "sugargliders")

myfavoriteanimals[1]

# ------------------------------------------------------------------------------------------
# ### Dictionaries
#
# If we have sets of data related to one another, we may choose to store that data in a
# dictionary. To do this, we use the `Dict()` function.
#
# Syntax:
# ```julia
# Dict(key1 => value1, key2 => value2, ...)```
#
# A good example of a dictionary is a contacts list, where we associate names with phone
# numbers.
# ------------------------------------------------------------------------------------------

myphonebook = Dict("Jenny" => "867-5309", "Ghostbusters" => "555-2368")

myphonebook["Jenny"]

# ------------------------------------------------------------------------------------------
# ### Arrays
#
# Unlike tuples, arrays are mutable. Unlike dictionaries, arrays contain ordered sequences
# of elements.

# We can create an array by enclosing this sequence of elements in `[ ]`.
#
# Syntax:

# ```julia
# [item1, item2, ...]```
#
#
# For example, we might create an array to keep track of my friends
# ------------------------------------------------------------------------------------------

myfriends = ["Ted", "Robyn", "Barney", "Lily", "Marshall"]

fibonacci = [1, 1, 2, 3, 5, 8, 13]

mixture = [1, 1, 2, 3, "Ted", "Robyn"]

# ------------------------------------------------------------------------------------------
# We can also create arrays of other data structures, or multi-dimensional arrays.
# ------------------------------------------------------------------------------------------

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

rand(4, 3)

# ------------------------------------------------------------------------------------------
# ## Loops
#
# ### `for` loops
#
# The syntax for a `for` loop is
#
# ```julia
# for *var* in *loop iterable*
#     *loop body*
# end
# ```
# ------------------------------------------------------------------------------------------

for n in 1:10
    println(n)
end

# ------------------------------------------------------------------------------------------
# ### `while` loops
#
# The syntax for a `while` is
#
# ```julia
# while *condition*
#     *loop body*
# end
# ```
# ------------------------------------------------------------------------------------------

n = 0
while n < 10
    n += 1
    println(n)
end

# ------------------------------------------------------------------------------------------
# ## Conditionals
#
# #### with `if`
#
# In Julia, the syntax
#
# ```julia
# if *condition 1*
#     *option 1*
# elseif *condition 2*
#     *option 2*
# else
#     *option 3*
# end
# ```
#
# allows us to conditionally evaluate one of our options.
# ------------------------------------------------------------------------------------------

x, y = # Enter two numbers here!
if x > y
    x
else
    y
end

# ------------------------------------------------------------------------------------------
# #### with ternary operators
#
# For this last block, we could instead use the ternary operator with the syntax
#
# ```julia
# a ? b : c
# ```
#
# which equates to
#
# ```julia
# if a
#     b
# else
#     c
# end
# ```
# ------------------------------------------------------------------------------------------

(x > y) ? x : y

# ------------------------------------------------------------------------------------------
# ## Functions
#
# Topics:
# 1. How to declare a function
# 2. Duck-typing in Julia
# 3. Mutating vs. non-mutating functions
# 4. Some higher order functions
#
# ### How to declare a function
#
# #### First way: with `function` and `end` keywords
# ------------------------------------------------------------------------------------------

function f(x)
    x^2
end

# ------------------------------------------------------------------------------------------
# #### Second way: with `=`
# ------------------------------------------------------------------------------------------

f2(x) = x^2

# ------------------------------------------------------------------------------------------
# Third way: as an anonymous function
# ------------------------------------------------------------------------------------------

f3 = x -> x^2

# ------------------------------------------------------------------------------------------
# #### Calling these functions
# ------------------------------------------------------------------------------------------

f(42)

f2(42)

f3(42)

# ------------------------------------------------------------------------------------------
# ### Duck-typing in Julia
# *"If it quacks like a duck, it's a duck."*

# Julia functions will just work on whatever inputs make sense.

# For example, `f` will work on a matrix.
# ------------------------------------------------------------------------------------------

A = rand(3, 3)
A

f(A)

# ------------------------------------------------------------------------------------------
# On the other hand, `f` will not work on a vector. Unlike `A^2`, which is well-defined, the
# meaning of `v^2` for a vector, `v`, is ambiguous.
# ------------------------------------------------------------------------------------------

v = rand(3)

f(v)

# ------------------------------------------------------------------------------------------
# ### Mutating vs. non-mutating functions
#
# By convention, functions followed by `!` alter their contents and functions lacking `!` do
# not.
#
# For example, let's look at the difference between `sort` and `sort!`.
# ------------------------------------------------------------------------------------------

v = [3, 5, 2]

sort(v)

v

# ------------------------------------------------------------------------------------------
# `sort(v)` returns a sorted array that contains the same elements as `v`, but `v` is left
# unchanged.

#
# On the other hand, when we run `sort!(v)`, the contents of v are sorted within the array
# `v`.
# ------------------------------------------------------------------------------------------

sort!(v)

v

# ------------------------------------------------------------------------------------------
# ### Some higher order functions
#
# #### map
#
# `map` is a "higher-order" function in Julia that *takes a function* as one of its input
# arguments.
# `map` then applies that function to every element of the data structure you pass it. For
# example, executing
#
# ```julia
# map(f, [1, 2, 3])
# ```
# will give you an output array where the function `f` has been applied to all elements of
# `[1, 2, 3]`
# ```julia
# [f(1), f(2), f(3)]
# ```
# ------------------------------------------------------------------------------------------

map(f, [1, 2, 3])

# ------------------------------------------------------------------------------------------
# Here we've squared all the elements of the vector `[1, 2, 3]`, rather than squaring the
# vector `[1, 2, 3]`.
#
# To do this, we could have passed to `map` an anonymous function rather than a named
# function, such as
# ------------------------------------------------------------------------------------------

x -> x^3

# ------------------------------------------------------------------------------------------
# via
# ------------------------------------------------------------------------------------------

map(x -> x^3, [1, 2, 3])

# ------------------------------------------------------------------------------------------
# and now we've cubed all the elements of `[1, 2, 3]`!
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ### broadcast
#
# `broadcast` is another higher-order function like `map`. `broadcast` is a generalization
# of `map`, so it can do every thing `map` can do and more. The syntax for calling
# `broadcast` is the same as for calling `map`
# ------------------------------------------------------------------------------------------

broadcast(f, [1, 2, 3])

# ------------------------------------------------------------------------------------------
# and again, we've applied `f` (squared) to all the elements of `[1, 2, 3]` - this time by
# "broadcasting" `f`!
#
# Some syntactic sugar for calling `broadcast` is to place a `.` between the name of the
# function you want to `broadcast` and its input arguments. For example,
#
# ```julia
# broadcast(f, [1, 2, 3])
# ```
# is the same as
# ```julia
# f.([1, 2, 3])
# ```
# ------------------------------------------------------------------------------------------

f.([1, 2, 3])

# ------------------------------------------------------------------------------------------
# Notice again how different this is from calling
# ```julia
# f([1, 2, 3])
# ```
# We can square every element of a vector, but we can't square a vector!
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# To drive home the point, let's look at the difference between
#
# ```julia
# f(A)
# ```
# and
# ```julia
# f.(A)
# ```
# for a matrix `A`:
# ------------------------------------------------------------------------------------------

A = [i + 3*j for j in 0:2, i in 1:3]

f(A)

# ------------------------------------------------------------------------------------------
# As before we see that for a matrix, `A`,
# ```
# f(A) = A^2 = A * A
# ```
#
# On the other hand,
# ------------------------------------------------------------------------------------------

B = f.(A)

# ------------------------------------------------------------------------------------------
# contains the squares of all the entries of `A`.
#
# This dot syntax for broadcasting allows us to write relatively complex compound
# elementwise expressions in a way that looks natural/closer to mathematical notation. For
# example, we can write
# ------------------------------------------------------------------------------------------

C = A .+ 2 .* f.(A) ./ A

# ------------------------------------------------------------------------------------------
# instead of
# ------------------------------------------------------------------------------------------

broadcast(x -> x + 2 * f(x) / x, A)

# ------------------------------------------------------------------------------------------
# and this will still compile down to code that runs as efficiently as `C`!
# ------------------------------------------------------------------------------------------

ext: jl

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
using Base.Filesystem
using Glob

@show pwd()

@show readdir(pwd())
@show @__FILE__
@show @__DIR__
@show @__LINE__

@show abspath("./fileio.jl")
@show abspath("../")

@show expanduser("~/work")
@show joinpath(expanduser("~/work"),"juliaLang","fileIO","fileio.jl")
jpath=joinpath(expanduser("~/work"),"juliaLang","fileIO","fileio.jl")
@show splitdir(jpath)
@show splitext(jpath)
@show basename(jpath)
@show dirname(jpath)
@show glob("*.jl")

for (root, dirs, files) in walkdir(".")
    println("Directories in $root")
    for dir in dirs
        println(joinpath(root, dir)) # path to directories
    end
    println("Files in $root")
    for file in files
        println(joinpath(root, file)) # path to files
    end
end

if !isdir("./myfolder/subfolder")
    mkpath("./myfolder/subfolder")
end

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
###
# Helpers for libarchive callbacks

"""
    ssize_t archive_read_callback(struct archive*,
                                  void *client_data, const void **buffer)

Returns pointer and size of next block of data from archive.
"""
to_read_callback{T}(func::Function, ::Type{T}) =
    cfunction(func, Cssize_t, Tuple{Ptr{Void},Ref{T},Ptr{Ptr{Void}}})

"""
    int64_t archive_skip_callback(struct archive*,
                                  void *client_data, int64_t request)

Skips at most request bytes from archive and returns the skipped amount.
This may skip fewer bytes than requested; it may even skip zero bytes.
If you do skip fewer bytes than requested, libarchive will invoke your
read callback and discard data as necessary to make up the full skip.
"""
to_skip_callback{T}(func::Function, ::Type{T}) =
    cfunction(func, Int64, Tuple{Ptr{Void},Ref{T},Int64})

"""
    int64_t archive_seek_callback(struct archive*,
                                  void *client_data, int64_t offset, int whence)

Seeks to specified location in the file and returns the position.
Whence values are SEEK_SET, SEEK_CUR, SEEK_END from stdio.h.
Return ARCHIVE_FATAL if the seek fails for any reason.
"""
to_seek_callback{T}(func::Function, ::Type{T}) =
    cfunction(func, Int64, Tuple{Ptr{Void},Ref{T},Int64,Cint})

"""
    ssize_t archive_write_callback(struct archive*,
                                   void *client_data,
                                   const void *buffer, size_t length);

Returns size actually written, zero on EOF, -1 on error.
"""
to_write_callback{T}(func::Function, ::Type{T}) =
    cfunction(func, Cssize_t, Tuple{Ptr{Void},Ref{T},Ptr{Void},Csize_t})

"""
    int archive_open_callback(struct archive*, void *client_data)
"""
to_open_callback{T}(func::Function, ::Type{T}) =
    cfunction(func, Cint, Tuple{Ptr{Void},Ref{T}})

"""
    int archive_close_callback(struct archive*, void *client_data)
"""
to_close_callback{T}(func::Function, ::Type{T}) =
    cfunction(func, Cint, Tuple{Ptr{Void},Ref{T}})

# """
#     int archive_switch_callback(struct archive*, void *client_data1,
#                                 void *client_data2)
#
# Switches from one client data object to the next/prev client data object.
# This is useful for reading from different data blocks such as a set of files
# that make up one large file.
# """
# to_switch_callback{T1,T2}(func::Function, ::Type{T1}=Void, ::Type{T2}=Void) =
#     cfunction(func, Cint, Tuple{Ptr{Void},Ptr{T1},Ptr{T2}})


text:
"""
Recursively replace all occurences of a field with dot notation:
```julia
type T
    a
    b
end
@with t::T a = b + 5
# becomes:
t.a = t.b + 5
```
"""
macro with(o, body::Expr)
    _with(o, body)
end

# --------------------------------------------------------------------------

# recurse
function replace_syms(expr::Expr, names, o)
    for i in eachindex(expr.args)
        expr.args[i] = replace_syms(expr.args[i], names, o)
    end
    expr
end

# replace?
function replace_syms(s::Symbol, names, o)
    if s in names
        :($(o).$(s))
    else
        s
    end
end

# identity
replace_syms(x, names, obj) = x

# if we pass in something like: type X; a; end; @with x f(x::X) = 10a
# then we look in the function definition to get the type of x
function _with(o::Symbol, ex::Expr)
    @assert ex.head == :function
    for farg in ex.args[1].args[2:end]
        # farg is an argument expression
        if farg == o
            error("You need to attach a type specification to know how to apply the @with macro, like: type X; a; end; @with x f(x::X) = 10a")
        elseif isa(farg, Expr) && farg.head == :(::) && farg.args[1] == o
            # found the right expression... apply it
            return _with(farg, ex)
        end
    end
    error("Couldn't match the symbol $o to a function argument in $(ex.args[1])")
end

function _with(oex::Expr, ex::Expr)
    # o is the symbol (x) and objtype is the type (X)
    @assert oex.head == :(::)
    @assert length(oex.args) == 2
    o, objtype = oex.args

    # if it's parameterized (x::X{T}) then just get the X-part
    if isa(objtype, Expr)
        @assert objtype.head == :curly
        objtype = objtype.args[1]
    end

    # get the fieldnames of that type, and replace them all
    names = fieldnames(eval(current_module(), objtype))
    esc(replace_syms(ex, names, o))
end

max_stars_repo_name: JuliaML/StochasticOptimization

max_issues_repo_name: JuliaML/StochasticOptimization

max_forks_repo_name: JuliaML/StochasticOptimization


text:
using ADCME
using LinearAlgebra
using AdFem
using PyPlot
using SparseArrays

function Q_exact(x, y)
    # (2*x^2*y*(x - 1)*(2*y - 1) + 1.0*(y^2 + 1)^2*(x*y*(x - 1) + x*y*(y - 1) + x*(x - 1)*(y - 1) + y*(x - 1)*(y - 1)) - (y^2 + 1)*(2*x*(x - 1)*(x + (x^2 + 1)*(y^2 + 1)) + 2*y*(x + (x^2 + 1)*(y^2 + 1))*(y - 1) + y*(2*x - 1)*(y - 1)*(2*x*(y^2 + 1) + 1)))/(y^2 + 1)^2
    # (2*x^2*y*(x - 1)*(2*y - 1) + x*(y^2 + 1)^2*(-y^2*(2*x - 1)*(y - 1)^2 + 1.0*y*(x - 1) + 1.0*(x - 1)*(y - 1)) - (y^2 + 1)*(2*x*(x - 1)*(x + (x^2 + 1)*(y^2 + 1)) + 2*y*(x + (x^2 + 1)*(y^2 + 1))*(y - 1) + y*(2*x - 1)*(y - 1)*(2*x*(y^2 + 1) + 1)))/(y^2 + 1)^2
    (2*x^2*y*(x - 1)*(2*y - 1)*(x + y + 1) - x*y^2*(2*x - 1)*(y - 1)^2*(y^2 + 1)^2*(x + y + 1) + x*(x - 1)*(2*y - 1)*(y^2 + 1)^2 - (y^2 + 1)*(x + y + 1)*(2*x*(x - 1)*(x + (x^2 + 1)*(y^2 + 1)) + 2*y*(x + (x^2 + 1)*(y^2 + 1))*(y - 1) + y*(2*x - 1)*(y - 1)*(2*x*(y^2 + 1) + 1)))/((y^2 + 1)^2*(x + y + 1))
end

function T_exact(x, y)
    x * (1-x) * y * (1-y)
end

function u_exact(x, y)
    x * y * (1-y)
end

function v_exact(x, y)
    1 / (1 + x + y)
end

function k_exact(x, y)
    1 + x^2 + x / (1 + y^2)
end

# function k_nn(xy) # xy shape N (=m*n) x 2
#     out = fc(xy, [20,20,20,1])^2 + 0.1 # N x 1
#     squeeze(out)
# end

#---------------------------------------------
# grid setup
m = 20
n = 20
h = 1/n

#---------------------------------------------
# discretized governing equation
#  J * u^{n+1} - J * u^n = - v ⋅ grad u + K * u^{n+1} + F^{n+1}
#---------------------------------------------
xy = fem_nodes(m, n, h)
x, y = xy[:,1], xy[:,2]
T0 = @. T_exact(x, y)
u = @. u_exact(x,y)
v = @. v_exact(x,y)
Q = @. Q_exact(x,y)
k = @. k_exact(x, y)
# k0 = @. k_nn(xy)

# ---------------------------------------------

bd = bcnode("all", m, n, h)

ugauss = fem_to_gauss_points(u, m, n, h)
vgauss = fem_to_gauss_points(v, m, n, h)
kgauss = fem_to_gauss_points(k, m, n, h)
Qgauss = fem_to_gauss_points(Q, m, n, h)

Advection = constant(compute_fem_advection_matrix1(constant(ugauss), constant(vgauss), m, n, h))
Laplace = constant(compute_fem_laplace_matrix1(kgauss, m, n, h))
A = Advection + Laplace
A, _ = fem_impose_Dirichlet_boundary_condition1(A, bd, m, n, h)
b = constant(compute_fem_source_term1(constant(Qgauss), m, n, h))
b = scatter_update(b, bd, zeros(length(bd)))
sol = A\b

# ---------------------------------------------------
# create a session and run
sess = Session(); init(sess)
T_computed = run(sess, sol)

#---------------------------------------------
# visualize numerical solution and exact solution
figure(figsize=(10,4))
subplot(121)
visualize_scalar_on_fem_points(T_computed, m, n, h)
subplot(122)
visualize_scalar_on_fem_points(T0, m, n, h)
savefig("forward_solution.png")
close("all")


text:
function twoSum(l::Vector{T}, target::T) where T
    len = length(l)
    for i in 1:(len - 1)
        for j in (i+1):len
            l[i] + l[j] == target && return [i,j]
        end
    end
    return []
end

max_stars_repo_name: ccqpein/Arithmetic-Exercises

max_issues_repo_name: ccqpein/Arithmetic-Exercises

max_forks_repo_name: ccqpein/Arithmetic-Exercises


text:
###
### Hash
###
###
### MurmurHash3 function of BioSequence.
###
### The hash function defined here cares about the starting position of a
### character sequence in the underlying data. That means, even if the starting
### positions of two sequences (`s1` and `s2`) are different in their `data`
### field, their hash values are identical if `s1 == s1` is true.  The
### implementation is based on the 128bit MurmurHash3 function, which was written
### by Austin Appleby, and the source code is distributed under the public domain:
### https://github.com/aappleby/smhasher/blob/master/src/MurmurHash3.cpp
###
### This file is a part of BioJulia.
### License is MIT: https://github.com/BioJulia/BioSequences.jl/blob/master/LICENSE.md

const c1 = 0x87c37b91114253d5
const c2 = 0x4cf5ad432745937f

@inline function rotl64(x::UInt64, r)
    return (x << (r & 63)) | (x >>> (-r & 63))
end

@inline function fmix64(k::UInt64)
    k = k ⊻ k >> 33
    k *= 0xff51afd7ed558ccd
    k = k ⊻ k >> 33
    k *= 0xc4ceb9fe1a85ec53
    k = k ⊻ k >> 33
    return k
end

@inline function murmur1(h1, k1)
    k1 *= c1
    k1 = rotl64(k1, 31)
    k1 *= c2
    h1 = h1 ⊻ k1
    return (h1, k1)
end

@inline function murmur2(h1, h2, k2)
    k2 *= c2
    k2 = rotl64(k2, 33)
    k2 *= c1
    h2 = h1 ⊻ k2
    return (h2, k2)
end

@inline function murmur(h1, h2, k1, k2)
    h1, k1 = murmur1(h1, k1)
    h1 = rotl64(h1, 27)
    h1 += h2
    h1 = h1 * 5 + 0x52dce729

    h2, k2 = murmur2(h1, h2, k2)
    h2 = rotl64(h2, 31)
    h2 += h1
    h2 = h2 * 5 + 0x38495ab5

    return (h1, h2, k1, k2)
end

function finalize(h1, h2, len)
    h1 = h1 ⊻ len
    h2 = h2 ⊻ len
    h1 += h2
    h2 += h1
    h1 = fmix64(h1)
    h2 = fmix64(h2)
    h1 += h2
    h2 += h1

    return h1
end

function tail(::Type{<:LongSequence}, data, next, stop, h1, h2)
    k1 = k2 = zero(UInt64)

    # Use this to mask any noncoding bits in last chunk
    mask = bitmask(offset(stop - bits_per_symbol(stop)) + bits_per_symbol(stop))
    @inbounds if next < stop
        k1 = data[index(next)]
        next += 64
    end
    @inbounds if next < stop
        k2 = data[index(next)] & mask
    else
        k1 &= mask
    end

    h1, k1 = murmur1(h1, k1)
    h2, k2 = murmur2(h1, h2, k2)
    return (h1, h2)
end

@inline function body(::Type{<:LongSequence}, next, stop, data, h1, h2)
    @inbounds while stop - next ≥ 128
        j = index(next)
        k1 = data[j]
        k2 = data[j+1]
        h1, h2, k1, k2 = murmur(h1, h2, k1, k2)
        next += 128
    end
    return (h1, h2, next)
end

# This version of the body loop code must take a nonzero offset into account
function body(::Type{<:LongSubSeq}, next, stop, data, h1, h2)
	off = offset(next)
	next < stop || return (h1, h2, next)
	# No offset, we can use the LongSequence one (more efficient)
	iszero(off) && return body(LongSequence, next, stop, data, h1, h2)

	@inbounds while stop - next ≥ 128
		k1 = data[index(next)]
		k2 = data[index(next) + 1]
		k3 = data[index(next) + 2]

		k1 = (k1 >>> off) | (k2 << ((64 - off) & 63))
		k2 = (k2 >>> off) | (k3 << ((64 - off) & 63))

		h1, h2, k1, k2 = murmur(h1, h2, k1, k2)
		next += 128
	end
	return h1, h2, next
end

function tail(::Type{<:LongSubSeq}, data, next, stop, h1, h2)
    next < stop || return (h1, h2)

    # Load in first up to 3 data elements where the last 128 bits may be stored
	firstindex = next
    k1 = data[index(next)]
    k2 = k3 = zero(UInt64)
    off = offset(next)
    next += 64 - off
    if next < stop
        k2 = @inbounds data[index(next)]
        next += 64
    end
    if next < stop
        k3 = @inbounds data[index(next)]
    end

    # Bitshift them to offset zero, only use k1 and k2
    if !iszero(off)
        k1 = (k1 >>> off) | (k2 << ((64 - off) & 63))
        k2 = (k2 >>> off) | (k3 << ((64 - off) & 63))
    end

	# Mask any noncoding bits
	mask = bitmask((stop-firstindex) & 63)
	if stop - firstindex > 64
		k2 &= mask
	else
		k1 &= mask
	end

    h1, k1 = murmur1(h1, k1)
    h2, k2 = murmur2(h1, h2, k2)

    return (h1, h2)
end

function Base.hash(seq::SeqOrView, seed::UInt64)
    # Mix sequence length so that dna"A" and dna"AA"
    # return the different hash values.
    h1::UInt64 = h2::UInt64 = hash(length(seq), seed)
    next = bitindex(seq, 1)
    stop = bitindex(seq, lastindex(seq) + 1)
    data = encoded_data(seq)

    h1, h2, next = body(typeof(seq), next, stop, data, h1, h2)
    h1, h2 = tail(typeof(seq), data, next, stop, h1, h2)
    h1 = finalize(h1, h2, length(seq))

    return h1
end


text:
using Documenter
using GraphQLClient
using StructTypes
using JSON3

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

makedocs(
    modules=[GraphQLClient],
    authors="Malcolm Miller",
    sitename="GraphQLClient.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://deloittedigitalapac.github.io/GraphQLClient.jl/stable",
    ),
    pages=[
        "Home" => "index.md",
        "Manual" => [
            "client.md",
            "operations.md",
            "struct_types_usage.md",
            "type_introspection.md",
            "low_level_execution.md",
            "limitations.md",
        ],
        "Library" => [
            "public.md",
            "private.md",
        ],
        "Contributing" => "contributing.md",
    ],
)

deploydocs(;
    repo="github.com/DeloitteDigitalAPAC/GraphQLClient.jl",
    push_preview=true,
    devbranch="main",
)


text:
#utilities for QUtilities
#some of these are candidates for Finometrics
const CRSP_MISSINGS = [
  "A", "B", "C", "D", "E", "S", "T", "P",
  -55.0, -66.0, -77.0, -88.0, -99.0,
  -55, -66, -77, -88, -99,
  "-55.0", "-66.0", "-77.0", "-88.0", "-99.0"]

const CRSP_PARSED_MISSINGS = (-55.0, -66.0, -77.0, -88.0, -99.0)
const CRSP_PARSED_MISSINGS = (-55.0, -66.0, -77.0, -88.0, -99.0)

const USED_COMPA_VALUE_FIELDS = [:at]
const USED_COMPQ_VALUE_FIELDS = [:atq, :ltq, :ibq, :txditcq, :pstkrq,
  :seqq, :ceqq, :pstkq]
const USED_CRSP_FIELDS = [:date, :ret, :mktcap]

#final columns for the decompressed file
const RETAINED_COLS_CRSP = [:permno, :date, :ret, :lret, :mktcap, :Lmktcap]

const RETAINED_COLS_COMP = [:begindate, :enddate,
  :gvkey, :adateq, :Nadateq,
  :fyear, :fqtr, :fyearqtr, :fdate, :fdateq,
    :cyearqtr, :ddate, :ddateq, :dyearq,
  :ia, :iaq,  :roeq]

const RETAINED_COLS_CCM = [:gvkey, :lpermno, :linkeffdate, :linkenddate]

#not sure if we really need much of these
const RETAINED_COLS_UNIV = [:rid, :permno, :date, :adateq, :ddateq, :fyear,
  :fyearqtr, :ret, :lret, :mktcap, :Lmktcap, :ia, :iaq, :roeq]

const RETAINED_VWCRSP_COLS = [:date, :vwindd]
const RETAINED_RFCRSP_COLS = [:date, :t30ind]

#these are the fields from which to build the factors
const Fia = :iaq
const Fmktcap = :Lmktcap
const Froe = :roeq

const USE_LOGS = false

const WEIGHT_ON = :Lmktcap
const RETURN_FIELD = USE_LOGS ? :lret : :ret

#NOTE: Below are portfolio construction parameters
const FACTOR_FIELDS = [Fmktcap, Fia, Froe]
const FACTOR_N_GROUPS = Dict(Fmktcap=>2, Fia=>3, Froe=>3)
const FACTOR_BREAKPOINTS = Dict(Fmktcap=>(0.5,1.0),
  Fia=>(0.3,0.7,1.0), Froe=>(0.3,0.7,1.0))
const FACTOR_GROUP_WEIGHTS = Dict(Fmktcap=>(1/9,-1/9),
  Fia=>(1/6, 0.0, -1/6), Froe=>(-1/6,0.0,1/6))
const FACTOR_VALUE_TYPES = Dict(Fia=>MFloat64, Froe=>MFloat64, Fmktcap=>MInt)

#NOTE: The below adjusts for logged vs unlogged returns
#const RETAINED_FACT_COLS = (( RETURN_FIELD == :lret) ?
#  [RETURN_FIELD; FACTOR_FIELDS; :lrfr, :lmkt] :
#  [RETURN_FIELD; FACTOR_FIELDS; :rfr, :mkt]

#converts gzip to lz4
function gz2lz4(file::String; filetype::String = "csv", location::String = "",
  gzextension::String = "gz", lz4extension::String = "lz4",
  deleteold::Bool = false)

  infile::String = "$location\\$file.$filetype.$gzextension"
  outfile::String = "$location\\$file.$filetype.$lz4extension"

  @assert isfile(infile)

  #Below doesnt work
  indata = open(infile)
  outdata = open(outfile, "w")
  s = GzipDecompressorStream(LZ4CompressorStream(outdata))
  write(s, indata)
  close(s)

  #run(``)

end

FILES_TO_RECOMPRESS = ["COMP-A", "COMP-Q", "CCM", "CRSP-M", "CRSP-D"]
function recompressqfactorfiles(files = FILES_TO_RECOMPRESS; datapath=DATA_PATH)::Nothing

  for f ∈ files
    gz2lz4(f, location=datapath)
  end

  return nothing
end

recompressqfactorfile(file; datapath=DATA_PATH)::Nothing = recompressqfactorfiles(
  [file]; datapath=datapath)

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
using Test
using DelimitedFiles, LinearAlgebra
import MLJBase
const Mlj = MLJBase
using StableRNGs
#rng = StableRNG(123)
using BetaML

#TESTRNG = FIXEDRNG # This could change...
TESTRNG = StableRNG(123)

println("*** Testing Decision trees/Random Forest algorithms...")

# ==================================
# NEW TEST
# ==================================
println("Testing basic classification of decision trees...")
# ---------------------
xtrain = [
    "Green"  3.0;
    "Yellow" 3.0;
    "Red"    1.0;
    "Red"    1.0;
    "Yellow" 3.0;
]

ytrain = ["Apple",  "Apple", "Grape", "Grape", "Lemon"]
myTree = buildTree(xtrain,ytrain,rng=copy(TESTRNG))

ŷtrain = predict(myTree, xtrain,rng=copy(TESTRNG))
@test accuracy(ŷtrain,ytrain,rng=copy(TESTRNG)) >= 0.8

xtest = [
    "Green"  3;
    "Yellow" 4;
    "Red"    2;
    "Red"    1;
    "Yellow" 3
]

ytest = ["Apple","Apple","Grape","Grape","Lemon"]
ŷtest = predict(myTree, xtest,rng=copy(TESTRNG))
@test accuracy(ŷtest,ytest,rng=copy(TESTRNG)) >= 0.8

#print(myTree)

# ==================================
# NEW TEST
println("Testing classification of the sepal database using decision trees...")
iris     = readdlm(joinpath(@__DIR__,"data","iris_shuffled.csv"),',',skipstart=1)
x = convert(Array{Float64,2}, iris[:,1:4])
y = convert(Array{String,1}, iris[:,5])

ntrain = Int64(round(size(x,1)*0.8))
xtrain = x[1:ntrain,:]
ytrain = y[1:ntrain]
xtest = x[ntrain+1:end,:]
ytest = y[ntrain+1:end]

myTree = buildTree(xtrain,ytrain, splittingCriterion=entropy,rng=copy(TESTRNG))
ŷtrain = predict(myTree, xtrain,rng=copy(TESTRNG))
@test accuracy(ŷtrain,ytrain,rng=copy(TESTRNG)) >= 0.98
ŷtest = predict(myTree, xtest,rng=copy(TESTRNG))
@test accuracy(ŷtest,ytest,rng=copy(TESTRNG))  >= 0.95

# ==================================
# NEW TEST
println("Testing decision trees regression...")

ϵtrain = [1.023,1.08,0.961,0.919,0.933,0.993,1.011,0.923,1.084,1.037,1.012]
ϵtest  = [1.056,0.902,0.998,0.977]
xtrain = [0.1 0.2; 0.3 0.5; 0.4 0.1; 0.5 0.4; 0.7 0.9; 0.2 0.1; 0.4 0.2; 0.3 0.3; 0.6 0.9; 0.3 0.4; 0.9 0.8]
ytrain = [(0.1*x[1]+0.2*x[2]+0.3)*ϵtrain[i] for (i,x) in enumerate(eachrow(xtrain))]
xtest  = [0.5 0.6; 0.14 0.2; 0.3 0.7; 20.0 40.0;]
ytest  = [(0.1*x[1]+0.2*x[2]+0.3)*ϵtest[i] for (i,x) in enumerate(eachrow(xtest))]

myTree = buildTree(xtrain,ytrain, minGain=0.001, minRecords=2, maxDepth=3,rng=copy(TESTRNG))
ŷtrain = predict(myTree, xtrain,rng=copy(TESTRNG))
ŷtest = predict(myTree, xtest,rng=copy(TESTRNG))
mreTrain = meanRelError(ŷtrain,ytrain)
@test mreTrain <= 0.06
mreTest  = meanRelError(ŷtest,ytest)
@test mreTest <= 0.3

# ==================================
# NEW TEST
println("Testing classification of the sepal database using random forests...")
iris     = readdlm(joinpath(@__DIR__,"data","iris_shuffled.csv"),',',skipstart=1)
x = convert(Array{Float64,2}, iris[:,1:4])
y = convert(Array{String,1}, iris[:,5])

ntrain = Int64(round(size(x,1)*0.8))
xtrain = x[1:ntrain,:]
ytrain = y[1:ntrain]
xtest = x[ntrain+1:end,:]
ytest = y[ntrain+1:end]

myForest = buildForest(xtrain,ytrain,β=0,maxDepth=20,oob=true,rng=copy(TESTRNG))

trees = myForest.trees
treesWeights = myForest.weights
oobError = myForest.oobError
ŷtrain = predict(myForest, xtrain,rng=copy(TESTRNG))
@test accuracy(ŷtrain,ytrain,rng=copy(TESTRNG)) >= 0.96
ŷtest = predict(myForest, xtest,rng=copy(TESTRNG))
@test accuracy(ŷtest,ytest,rng=copy(TESTRNG))  >= 0.96
updateTreesWeights!(myForest,xtrain,ytrain;β=1)
ŷtrain2 = predict(myForest, xtrain,rng=copy(TESTRNG))
@test accuracy(ŷtrain2,ytrain,rng=copy(TESTRNG)) >= 0.98
ŷtest2 = predict(myForest, xtest,rng=copy(TESTRNG))
@test accuracy(ŷtest2,ytest,rng=copy(TESTRNG))  >= 0.96
@test oobError <= 0.1

# ==================================
# NEW TEST
println("Testing random forest regression...")

ϵtrain = [1.023,1.08,0.961,0.919,0.933,0.993,1.011,0.923,1.084,1.037,1.012]
ϵtest  = [1.056,0.902,0.998,0.977]
xtrain = [0.1 0.2; 0.3 0.5; 0.4 0.1; 0.5 0.4; 0.7 0.9; 0.2 0.1; 0.4 0.2; 0.3 0.3; 0.6 0.9; 0.3 0.4; 0.9 0.8]
ytrain = [(0.1*x[1]+0.2*x[2]+0.3)*ϵtrain[i] for (i,x) in enumerate(eachrow(xtrain))]
xtest  = [0.5 0.6; 0.14 0.2; 0.3 0.7; 20.0 40.0;]
ytest  = [(0.1*x[1]+0.2*x[2]+0.3)*ϵtest[i] for (i,x) in enumerate(eachrow(xtest))]

myForest         = buildForest(xtrain,ytrain, minGain=0.001, minRecords=2, maxDepth=3,rng=copy(TESTRNG))
trees            = myForest.trees
treesWeights     = myForest.weights

ŷtrain           = predict(myForest, xtrain,rng=copy(TESTRNG))
ŷtest            = predict(myForest, xtest,rng=copy(TESTRNG))
mreTrain         = meanRelError(ŷtrain,ytrain)
@test mreTrain <= 0.08
mreTest  = meanRelError(ŷtest,ytest)
@test mreTest <= 0.4

updateTreesWeights!(myForest,xtrain,ytrain;β=50)
ŷtrain2 = predict(myForest, xtrain,rng=copy(TESTRNG))
ŷtest2 = predict(myForest, xtest,rng=copy(TESTRNG))
mreTrain = meanRelError(ŷtrain2,ytrain)
@test mreTrain <= 0.08
mreTest  = meanRelError(ŷtest2,ytest)
@test mreTest <= 0.4

# ==================================
# NEW TEST
println("Testing all possible combinations...")
xtrain = [1 "pippo" 1.5; 3 "topolino" 2.5; 1 "amanda" 5.2; 5 "zzz" 1.2; 7 "pippo" 2.2; 1 "zzz" 1.5; 3 "topolino" 2.1]
ytrain = [x[2][1] <= 'q' ? 5*x[1]-2*x[3] : -5*x[1]+2*x[3] for x in eachrow(xtrain)]
xtest = [2 "pippo" 3.4; 1 "amanda" 1.5; 4 "amanda" 0.5; 2 "topolino" 2.2; 7 "zzz" 3.2]
ytest = [x[2][1] <= 'q' ? 5*x[1]-2*x[3] : -5*x[1]+2*x[3] for x in eachrow(xtest)]
ytrainInt = Int64.(round.(ytrain))

myTree1 = buildTree(xtrain,ytrain,rng=copy(TESTRNG))
myForest = buildForest(xtrain,ytrain,oob=true,rng=copy(TESTRNG)) # TO.DO solved 20201130: If I add here β=1 I have no problem, but local testing gives a crazy error!!!
oobError = myForest.oobError
ŷtrain = predict(myForest,xtrain,rng=copy(TESTRNG))
ŷtest = predict(myForest,xtest,rng=copy(TESTRNG))
mreTrain = meanRelError(ŷtrain,ytrain)
mreTest  = meanRelError(ŷtest,ytest)

xtrain[3,3] = missing
xtest[3,2] = missing
myForest = buildForest(xtrain,ytrain,oob=true,β=1,rng=copy(TESTRNG))
oobError = myForest.oobError
ŷtrain = predict(myForest,xtrain,rng=copy(TESTRNG))
ŷtest = predict(myForest,xtest,rng=copy(TESTRNG))
mreTrain2 = meanRelError(ŷtrain,ytrain)
mreTest2  = meanRelError(ŷtest,ytest)
@test mreTest2 <= mreTest * 1.5

myTree2 = buildTree(xtrain,ytrainInt,rng=copy(TESTRNG))
myTree3 = buildTree(xtrain,ytrainInt, forceClassification=true,rng=copy(TESTRNG))
@test typeof(myTree1) <: Trees.DecisionNode && typeof(myTree2) <: Trees.DecisionNode && typeof(myTree3) <: Trees.DecisionNode

# ==================================
# NEW TEST
println("Testing trees with unsortable and missing X values...")

abstract type AType end
mutable struct SortableType<:AType
    x::Int64
    y::Int64
end

mutable struct UnsortableType<:AType
    x::Int64
    y::Int64
end
isless(x::SortableType,y::SortableType) = x.x < y.x

SortableVector = [SortableType(2,4),SortableType(1,5),SortableType(1,8),SortableType(12,5),
    SortableType(6,2),SortableType(2,2),SortableType(2,2),SortableType(2,4),
    SortableType(6,2),SortableType(1,5),missing,SortableType(2,4),
    SortableType(1,8),SortableType(12,5)]
UnSortableVector = [UnsortableType(2,5),UnsortableType(1,3),UnsortableType(1,8),UnsortableType(2,6),
    UnsortableType(6,3),UnsortableType(7,9),UnsortableType(2,5),UnsortableType(2,6),
    missing,UnsortableType(3,2),UnsortableType(6,3),UnsortableType(2,5),
    UnsortableType(7,9),UnsortableType(7,9)]

data = Union{Missing,Float64, String,AType}[
    0.9 0.6 "black" "monitor" 10.1
    0.3 missing "white" "paper sheet" 2.3
    4.0 2.2 missing "monitor"  12.5
    0.6 0.5 "white" "monitor" 12.5
    3.8 2.1 "gray" "car" 54.2
    0.3 0.2 "red" "paper sheet" 2.6
    0.1 0.1 "white" "paper sheet" 2.5
    0.3 0.2 "black" "monitor" 11.3
    0.1 0.2 "black" "monitor" 9.8
    0.31 0.2 "white" "paper sheet" 3.7
    3.2 1.9 "gray" "car" 64.3
    0.4 0.25 "white" "paper" 2.7
    0.9 0.4 "black" "monitor" 12.5
    4.1 2.1 "gray" "monitor" 13.2
]

X = hcat(data[:,[1,2,3]],UnSortableVector)
y = convert(Vector{String},  data[:,4])

((xtrain,xtest),(ytrain,ytest)) = Utils.partition([X,y],[0.7,0.3],shuffle=false,rng=copy(TESTRNG))

modelβ = buildForest(xtrain,ytrain,5,rng=copy(TESTRNG))
ŷtestβ = predict(modelβ,xtest,rng=copy(TESTRNG))
accβ   = accuracy(ŷtestβ,ytest,rng=copy(TESTRNG))
@test accβ >= 0.25

# ==================================
# NEW TEST
println("Testing MLJ interface for Trees models....")
X, y                           = Mlj.@load_boston
model_dtr                      = DecisionTreeRegressor(rng=copy(TESTRNG))
regressor_dtr                  = Mlj.machine(model_dtr, X, y)
(fitresult_dtr, cache, report) = Mlj.fit(model_dtr, 0, X, y)
yhat_dtr                       = Mlj.predict(model_dtr, fitresult_dtr, X)
@test meanRelError(yhat_dtr,y) < 0.02

model_rfr                      = RandomForestRegressor(rng=copy(TESTRNG))
regressor_rfr                  = Mlj.machine(model_rfr, X, y)
(fitresult_rfr, cache, report) = Mlj.fit(model_rfr, 0, X, y)
yhat_rfr                       = Mlj.predict(model_rfr, fitresult_rfr, X)
@test meanRelError(yhat_rfr,y) < 0.06

X, y                           = Mlj.@load_iris
model_dtc                      = DecisionTreeClassifier(rng=copy(TESTRNG))
regressor_dtc                  = Mlj.machine(model_dtc, X, y)
(fitresult_dtc, cache, report) = Mlj.fit(model_dtc, 0, X, y)
yhat_dtc                       = Mlj.predict(model_dtc, fitresult_dtc, X)
@test Mlj.mean(Mlj.LogLoss(tol=1e-4)(yhat_dtc, y)) < 0.0002

model_rfc                      = RandomForestClassifier(maxFeatures=3,rng=copy(TESTRNG))
regressor_rfc                  = Mlj.machine(model_rfc, X, y)
(fitresult_rfc, cache, report) = Mlj.fit(model_rfc, 0, X, y)
yhat_rfc                       = Mlj.predict(model_rfc, fitresult_rfc, X)
@test Mlj.mean(Mlj.LogLoss(tol=1e-4)(yhat_rfc, y)) < 0.04

# Other MLJ classifier models
#=
import MLJ
X, y                       = Mlj.@load_iris
MLJ.models(MLJ.matching(X,y))
Model                      = MLJ.@load XGBoostClassifier #  DecisionTreeClassifier    #    XGBoostClassifier
model                      = Model()
regressor                  = MLJ.machine(model, X, y)
(fitresult, cache, report) = MLJ.fit(model, 0, X, y)
yhat                       = MLJ.predict(model, fitresult, X)
MLJ.mean(MLJ.LogLoss(tol=1e-4)(yhat, y))
MLJ.evaluate!(regressor, measure=MLJ.LogLoss())
#XGBoostClassifier:
#- fit: https://github.com/alan-turing-institute/MLJModels.jl/blob/3687491b132be8493b6f7a322aedf66008caaab1/src/XGBoost.jl#L600
#- predict :
=#


text:
using Findpeaks
using Test

@testset "Single Gaussian peak" begin
    gaussian(x,μ,σ) = 1/sqrt(2*π)/σ*exp(-((x-μ)^2)/2/σ^2)
    x = collect(range(1,stop=1000,length=1000))
    μ = rand(x)
    data = gaussian.(x,μ,rand(x)/3)

    p = findpeaks(data)
    @test length(p) == 1
    @test p[1] == μ

    p = findpeaks(data,x)
    @test length(p) == 1
    @test p[1] == μ
end

@testset "Prominence" begin
    gaussian(x,μ,σ) = exp(-((x-μ)^2)/2/σ^2)

    x = range(0., stop = 1., length = 1000) |> collect

    positions = [1, 2, 3, 5, 7, 8]/10;
    amps = [3, 7, 5, 5, 4, 5];
    widths = [1, 3, 3, 4, 2, 3]/100;

    base = 4 * cos.(2π * x);
    signal = ( a * gaussian.(x, μ, σ) for (a, μ, σ) in zip(amps, positions, widths) )
    data = sum(signal) + base

    peaks = findpeaks(data, x, min_prom=4.)

    @test length(peaks) == 2

    expected_peak_1 = argmin(abs.(x .- positions[2])) # global highest peak
    expected_peak_2 = argmin(abs.(x .- positions[4])) # lowest peak

    # check peaks are around expected - may be shifted because of background
    @test abs(peaks[1] - expected_peak_1) < 20
    @test abs(peaks[2] - expected_peak_2) < 20
end

@testset "Threshold" begin
    y = [0., 1., 3., 1., 4., 5., 3., 0., 2., 5., 4., 0.]

    @test Set(findpeaks(y, threshold=3.1)) == Set([])
    @test Set(findpeaks(y, threshold=2.1)) == Set([])
    @test Set(findpeaks(y, threshold=1.1)) == Set([3])
    @test Set(findpeaks(y, threshold=0.1)) == Set([3, 6, 10])
    @test Set(findpeaks(y, threshold=0.0)) == Set([3, 6, 10])
end

@testset "Min. Distance" begin
    y = [0., 1., 3., 1., 4., 5., 3., 0., 2., 5., 4., 0.]

    @test Set(findpeaks(y, min_dist=4)) == Set([6])
    @test Set(findpeaks(y, min_dist=3)) == Set([6, 10])
    @test Set(findpeaks(y, min_dist=2)) == Set([3, 6, 10])
    @test Set(findpeaks(y, min_dist=0)) == Set([3, 6, 10])
    @test Set(findpeaks(y, min_dist=0)) == Set([3, 6, 10])
end

@testset "Min. Height" begin
    y = [0., 1., 3., 1., 4., 5., 3., 0., 2., 5., 4., 0.]

    @test Set(findpeaks(y, min_height=6.)) == Set([])
    @test Set(findpeaks(y, min_height=4.9)) == Set([6, 10])
    @test Set(findpeaks(y, min_height=2.9)) == Set([3, 6, 10])
end


text:
module Observables2

export Observable
export ObservingFunction
export observe!
export stop_observing!
export disable!
export is_disabled
export listeners
export notify!
export to_value
export on
export off
export onany
export connect! # obsid, async_latest, throttle
export n_ordinary_inputs
export n_observable_inputs

abstract type AbstractObservable{T} end

mutable struct Observable{V} <: AbstractObservable{V}
    val::V
    f::Any
    inputs::Vector{Any}
    listeners::Vector{Any}
    onlynew::Bool
end

struct ObservingFunction{F}
    f::F
    inputs::Vector{Any}
end

# creating an observable from a value only
# no inputs or listeners are set up
Observable(v; onlynew = false, type = nothing) =
    Observable{isnothing(type) ? typeof(v) : type}(v, nothing, [], Observable[], onlynew)

Observable{T}(v) where T = Observable(v; type = T)

# function Base.map(f, o::AbstractObservable, os...)
#     observe!(f, o, os...)
# end

# if we allow replacement of active observables with their values to inactivate them
# there needs to be a wrapper type to distinguish between an observable that is registered
# and an observable that is just a value
struct RegisteredObservable
    o::Observable
end

get_registered_value(r::RegisteredObservable) = r.o.val
get_registered_value(any) = any

# obs[!] = new_value always triggers a change message
function Base.setindex!(obs::Observable, value, ::typeof(!))
    obs.val = value
    for o in obs.listeners
        # ObservingFunctions don't get a new value set, just their function called
        if o isa ObservingFunction
            # call the stored function with the values of registered observables
            # or the inputs directly
            o.f((get_registered_value(input) for input in o.inputs)...)
        else
            # call the stored function with the values of registered observables
            # or the inputs directly and save the value, set it as the new observable value
            new_value = o.f((get_registered_value(input) for input in o.inputs)...)
            o[] = new_value
        end
    end
    obs
end

# obs[] = new_value only triggers a change message if onlynew is true
function Base.setindex!(obs::Observable, value)
    if !(obs.onlynew && obs.val == value)
        obs[!] = value
    end
end

Base.getindex(obs::Observable) = obs.val

function register!(with::Observable, o::Union{Observable, ObservingFunction})
    if o in with.listeners
        error("Observable already registered as observer.")
    end
    push!(with.listeners, o)
end

wrap_register(any) = any
wrap_register(o::Observable) = RegisteredObservable(o)

unwrap_register(any) = any
unwrap_register(ro::RegisteredObservable) = ro.o

# creating an observable from inputs and a function acting on those inputs

function observe!(f, inputs...; onlynew = false, type::Union{Type, Nothing} = nothing)

    # compute first value
    value = f((to_value(input) for input in inputs)...)

    # make a new observer that tracks who he observes but isn't yet registered with the other listeners
    param = if isnothing(type)
        # the parameter type is decided by the value
        typeof(value)
    else
        # the parameter type is manually given
        type
    end

    obs = Observable{param}(
        value,
        f,
        Any[wrap_register(i) for i in inputs],
        [],
        onlynew)

    # register the new observable with the inputs
    for o in inputs
        if o isa Observable
            register!(o, obs)
        end
    end

    obs
end

# these functions mimick the current observables api

# the difference for on / onany / off is that they return an ObservingFunction and not a function
# that is because a function doesn't know about observables that keep track of it
# so it becomes more difficult to unlink

function on(f, input)
    onany(f, input)
end

# here the argument order is different, the listener is removed from the input
# because only that direction is possible in Observables.jl
function off(input, listener)
    stop_observing!(listener, input)
end

function onany(f, inputs...)
    ofunc = ObservingFunction(f, Any[wrap_register(i) for i in inputs])

    for o in inputs
        if o isa Observable
            register!(o, ofunc)
        end
    end

    ofunc
end

Base.eltype(::AbstractObservable{T}) where {T} = T

# triggering an observable again with its current value
function notify!(o::Observable)
    o[!] = o[]
end

# extracting a value from an observable
to_value(o::Observable) = o.val
to_value(x) = x

listeners(o::Observable) = o.listeners

# function Base.copy(o::Observable{T}) where T
#     oc = Observable{T}(o.val)
#     on(o) do o
#         oc[] = o
#     end
#     oc
# end

function connect!(o1::Observable, o2::Observable)
    error("not implemented")
end

"""
    stop_observing!(observer::Observable, input::Observable)

Replace the `input` `Observable` in the input vector of the `observer` with its value
and remove the `observer` from the `listeners` of the `input` `Observable`.
"""
function stop_observing!(observer, input::Observable)
    # delete input in observer
    inputindex = findfirst(==(RegisteredObservable(input)), observer.inputs)

    if isnothing(inputindex)
        error("The input observable was not found as an input in the listener observable.")
    end

    observerindex = findfirst(==(observer), input.listeners)
    if isnothing(observerindex)
        error("The listener observable was not found as a listener in the input observable.")
    end

    observer.inputs[inputindex] = to_value(input)
    deleteat!(input.listeners, observerindex)
end

# if the input is not an observable
function stop_observing!(observer, input)
    # do nothing
end

"""
    stop_observing!(observer::Observable)

Replace all `Observable`s in the input vector of the `observer` with their values
and remove the `observer` from the `listeners` of each of these `Observable`s.
"""
function stop_observing!(observer)
    for input in observer.inputs
        stop_observing!(observer, unwrap_register(input))
    end
end

"""
    disable!(o::Observable; recursive = true)

Make an observable `o` stop observing all its inputs.
Set the stored function to `nothing`, so objects referenced in the
closure can possibly be released.

Then also make all listeners of `o` stop observing it.
If `recursive` is `true`, all listeners that don't have any `Observable` inputs
left after stopping are also disabled.

Returns the number of all disabled `Observable`s.
"""
function disable!(o::Observable; recursive = true)

    stop_observing!(o)
    o.f = nothing

    n_disabled = 1
    for observer in o.listeners
        stop_observing!(observer, o)
        # if the observer doesn't have any inputs left that are observables
        # we can disable it as well, because it will never be triggered again
        if recursive && n_observable_inputs(observer) == 0
            n_disabled += disable!(observer; recursive = true)
        end
    end
    n_disabled
end

"""
    is_disabled(o::Observable)

Checks that `o` has no observable inputs and no listeners.
Also, it should have no stored function.
"""
function is_disabled(o::Observable)
    n_observable_inputs(o) == 0 && isempty(listeners(o)) && isnothing(o.f)
end

n_ordinary_inputs(o) = sum((!isa).(o.inputs, RegisteredObservable))
n_observable_inputs(o) = sum(isa.(o.inputs, RegisteredObservable))

function Base.show(io::IO, o::Observable{V}) where V

    n_observable = n_observable_inputs(o)
    n_ordinary = n_ordinary_inputs(o)
    n_listeners = length(o.listeners)

    str = "Observable{$V} with $n_observable observable, $n_ordinary ordinary inputs, and $n_listeners listeners."
    println(io, str)
    print(io, "Value: $(o.val)")
    nothing
end

end # module


text:
Dims{N} = NTuple{N,Integer}

include("impl/pool.jl")
include("impl/conv.jl")

# Convolutions

function cdims(x::NTuple{N}, w::NTuple{N}, pad, stride) where N
  ntuple(Val{N}) do i
    if i < N-1
      1 + div(x[i] - w[i] + 2*pad[i], stride[i])
    elseif i == N-1
      w[N]
    else # i == N
      x[N]
    end
  end
end

# Interface

head(x) = reverse(Base.tail(reverse(x)))
padtuple(x::Tuple,p::Integer) = map(_->p, head(head(x)))
padtuple(x::Tuple,p::Tuple) = p
padtuple(x::AbstractArray,p) = padtuple(size(x),p)

function conv(x::A, w::A; pad = 0, stride = 1) where A<:AbstractArray
  pad_, stride_ = padtuple(x, pad), padtuple(x, stride)
  conv!(similar(x, cdims(size(x), size(w), pad_, stride_)),
        x, w, pad = pad_, stride = stride_)
end

∇conv_data(dy::A, x::A, w::A; pad = 0, stride = 1) where A<:AbstractArray =
  ∇conv_data!(zeros(x), dy, x, w; pad = pad, stride = stride)

∇conv_filter(dy::A, x::A, w::A; pad = 0, stride = 1) where A<:AbstractArray =
  ∇conv_filter!(zeros(w), dy, x, w; pad = pad, stride = stride)

# N-D dispatch

function conv!(y::AbstractArray{T,3}, x::AbstractArray{T,3}, w::AbstractArray{T,3};
               pad = 0, stride = 1) where T
    args = map(x -> reshape(x, size(x,1),1,size(x,2),size(x,3)), (y, x, w))
    conv!(args..., pad = (pad...,0), stride = (stride...,1))
    return y
end

function ∇conv_filter!(dw::AbstractArray{T,3}, dy::AbstractArray{T,3},
                       x::AbstractArray{T,3}, w::AbstractArray{T,3};
                       pad = 0, stride = 1) where T
    args = map(x -> reshape(x, size(x,1),1,size(x,2),size(x,3)), (dw, dy, x, w))
    ∇conv_filter!(args..., pad = (pad...,0), stride = (stride...,1))
    return dw
end

function ∇conv_data!(dx::AbstractArray{T,3}, dy::AbstractArray{T,3},
                     x::AbstractArray{T,3}, w::AbstractArray{T,3};
                     pad = 0, stride = 1) where T
    args = map(x -> reshape(x, size(x,1),1,size(x,2),size(x,3)), (dx, dy, x, w))
    ∇conv_data!(args..., pad = (pad...,0), stride = (stride...,1))
    return dx
end

conv!(y::AbstractArray{T,4}, x::AbstractArray{T,4}, w::AbstractArray{T,4};
      pad = 0, stride = 1) where T =
  conv2d!(y, x, w, padding = pad, stride = stride)

∇conv_filter!(dw::AbstractArray{T,4}, dy::AbstractArray{T,4}, x::AbstractArray{T,4}, w::AbstractArray{T,4};
              pad = 0, stride = 1) where T =
  conv2d_grad_w!(dw, x, w, dy, padding = pad, stride = stride)

∇conv_data!(dx::AbstractArray{T,4}, dy::AbstractArray{T,4}, x::AbstractArray{T,4}, w::AbstractArray{T,4};
            pad = 0, stride = 1) where T =
  conv2d_grad_x!(dx, x, w, dy, padding = pad, stride = stride)

conv!(y::AbstractArray{T,5}, x::AbstractArray{T,5}, w::AbstractArray{T,5};
      pad = 0, stride = 1) where T =
  conv3d!(y, x, w, padding = pad, stride = stride)

∇conv_filter!(dw::AbstractArray{T,5}, dy::AbstractArray{T,5}, x::AbstractArray{T,5}, w::AbstractArray{T,5};
              pad = 0, stride = 1) where T =
  conv3d_grad_w!(dw, x, w, dy, padding = pad, stride = stride)

∇conv_data!(dx::AbstractArray{T,5}, dy::AbstractArray{T,5}, x::AbstractArray{T,5}, w::AbstractArray{T,5};
            pad = 0, stride = 1) where T =
  conv3d_grad_x!(dx, x, w, dy, padding = pad, stride = stride)

# Pooling

function pdims(dims::Dims{N}, window, padding, stride) where N
  ntuple(Val{N}) do i
    if i < N-1
      1 + (dims[i] + 2*padding[i] - window[i])÷stride[i]
    else
      dims[i]
    end
  end
end

expand(::Type{Val{N}}, i::Integer) where N = ntuple(_ -> i, Val{N})
expand(::Type{Val{N}}, i::NTuple{N, Integer}) where N = i

# Interface

maxpool(x::AbstractArray, k; pad = map(_->0,k), stride = k) =
  maxpool!(similar(x, pdims(size(x), k, expand(Val{length(k)}, pad),
            expand(Val{length(k)}, stride))), x, k, pad = expand(Val{length(k)}, pad),
            stride = expand(Val{length(k)}, stride))

maxpool!(y::A, x::A, k; kw...) where A<:AbstractArray =
  maxpool_cpu!(y, x, k; kw...)

∇maxpool(dy::A, y::A, x::A, k; pad = map(_->0,k), stride = k) where A<:AbstractArray =
  ∇maxpool!(similar(x), dy, y, x, k, pad = expand(Val{length(k)}, pad),
            stride = expand(Val{length(k)}, stride))

∇maxpool!(dx::A, dy::A, y::A, x::A, k; kw...) where A<:AbstractArray =
  ∇maxpool_cpu!(dx, dy, y, x, k; kw...)

meanpool(x::AbstractArray, k; pad = map(_->0,k), stride = k) =
  meanpool!(similar(x, pdims(size(x), k, expand(Val{length(k)}, pad),
            expand(Val{length(k)}, stride))), x, k, pad = expand(Val{length(k)}, pad),
            stride = expand(Val{length(k)}, stride))

meanpool!(y::A, x::A, k; kw...) where A<:AbstractArray =
  meanpool_cpu!(y, x, k; kw...)

∇meanpool(dy::A, y::A, x::A, k; pad = map(_->0,k), stride = k) where A<:AbstractArray =
  ∇meanpool!(similar(x), dy, y, x, k, pad = expand(Val{length(k)}, pad),
            stride = expand(Val{length(k)}, stride))

∇meanpool!(dx::A, dy::A, y::A, x::A, k; kw...) where A<:AbstractArray =
  ∇meanpool_cpu!(dx, dy, y, x, k; kw...)

# N-D dispatch
# We use a separate function to avoid ambiguity issues
# (more specific array types vs. more specific dimensions)

maxpool_cpu!(y::A, x::A, k::Dims{2}; pad = (0,0), stride = k) where A<:AbstractArray{<:Real,4} =
  maxpool2d!(y, x, window = k, padding = pad, stride = stride)

∇maxpool_cpu!(dx::AbstractArray{<:Real,4}, dy::AbstractArray{<:Real,4}, y::AbstractArray{<:Real,4}, x::AbstractArray{<:Real,4},
              k::Dims{2}; pad = (0,0), stride = k) =
  maxpool2d_grad!(dx, dy, y, x,
                  window = k, padding = pad, stride = stride)

maxpool_cpu!(y::AbstractArray{<:Real,5}, x::AbstractArray{<:Real,5}, k::Dims{3}; pad = (0,0), stride = k) =
  maxpool3d!(y, x, window = k, padding = pad, stride = stride)

∇maxpool_cpu!(dx::AbstractArray{<:Real,5}, dy::AbstractArray{<:Real,5}, y::AbstractArray{<:Real,5}, x::AbstractArray{<:Real,5},
              k::Dims{3}; pad = (0,0), stride = k) =
  maxpool3d_grad!(dx, dy, y, x,
                  window = k, padding = pad, stride = stride)

meanpool_cpu!(y::AbstractArray{<:Real,4}, x::AbstractArray{<:Real,4}, k::Dims{2}; pad = (0,0), stride = k) =
  meanpool2d!(y, x, window = k, padding = pad, stride = stride)

∇meanpool_cpu!(dx::AbstractArray{<:Real,4}, dy::AbstractArray{<:Real,4}, y::AbstractArray{<:Real,4}, x::AbstractArray{<:Real,4},
              k::Dims{2}; pad = (0,0), stride = k) =
  meanpool2d_grad!(dx, dy, y, x,
                   window = k, padding = pad, stride = stride)

meanpool_cpu!(y::AbstractArray{<:Real,5}, x::AbstractArray{<:Real,5}, k::Dims{3}; pad = (0,0), stride = k) =
  meanpool3d!(y, x, window = k, padding = pad, stride = stride)

∇meanpool_cpu!(dx::AbstractArray{<:Real,5}, dy::AbstractArray{<:Real,5}, y::AbstractArray{<:Real,5}, x::AbstractArray{<:Real,5},
              k::Dims{3}; pad = (0,0), stride = k) =
  meanpool3d_grad!(dx, dy, y, x,
                   window = k, padding = pad, stride = stride)

# Deprecated 0.3

export conv2d, maxpool2d, avgpool2d

@deprecate conv2d(x, w; kw...) NNlib.conv(x, w; kw...)
@deprecate maxpool2d(x::AbstractArray{<:Real,4}, k::Integer) maxpool(x, (k,k))
@deprecate meanpool2d(x::AbstractArray{<:Real,4}, k::Integer) meanpool(x, (k,k))


text:
using Distributed
addprocs(4)
nprocs()

using MendelIHT
using SnpArrays
using DataFrames
using Distributions
using StatsFuns: logistic
using Random
using LinearAlgebra
using DelimitedFiles
using GLM
using RCall
R"library(glmnet)"

function iht_lasso_poisson(n :: Int64, p :: Int64, k :: Int64)
    glm = "poisson"

    #construct snpmatrix, covariate files, and true model b
    x, maf = simulate_random_snparray(n, p, "poisson.bed")
    xbm = SnpBitMatrix{Float64}(x, model=ADDITIVE_MODEL, center=true, scale=true);
    z = ones(n, 1) # non-genetic covariates, just the intercept
    true_b = zeros(p)
    true_b[1:k] = rand(Normal(0, 0.3), k)
    shuffle!(true_b)
    correct_position = findall(x -> x != 0, true_b)
    x_float = [convert(Matrix{Float64}, x, center=true, scale=true) z] #Float64 version of x

    # Simulate poisson data
    y_temp = xbm * true_b
    λ = exp.(y_temp) #inverse log link
    y = [rand(Poisson(x)) for x in λ]
    y = Float64.(y)

    #specify path and folds
    num_folds = 3
    folds = rand(1:num_folds, size(x, 1));

    #run glmnet via Rcall
    @rput x_float y folds num_folds #make variables visible to R
    R"lasso_cv_result = cv.glmnet(x_float, y, nfolds = num_folds, foldid = folds, family='poisson')"
    R"lasso_beta_tmp = glmnet(x_float, y, lambda=lasso_cv_result$lambda.min, family='poisson')$beta"
    R"lasso_beta = as.vector(lasso_beta_tmp)"
    @rget lasso_cv_result lasso_beta #pull result from R to Julia
    lasso_k_est = length(findall(!iszero, lasso_beta))

    #find non-zero entries returned by best lasso model as largest k estimate
    path = collect(1:50);

    #run IHT's cross validation routine
    mses = cv_iht_distributed(x, z, y, 1, path, folds, num_folds, glm, use_maf=false, debias=true)
    iht_k_est = argmin(mses)
    iht_result = L0_poisson_reg(x, xbm, z, y, 1, iht_k_est, glm = "poisson", debias=true, convg=false, show_info=false, true_beta=true_b, scale=false, init=false)
    iht_beta = iht_result.beta

    #show lasso and IHT's reconstruction result
    compare_model = DataFrame(
        true_β  = true_b[correct_position],
        IHT_β   = iht_beta[correct_position],
        lasso_β = lasso_beta[correct_position])
    @show compare_model

    #compute precision/recall for IHT and lasso
    iht_tp = length(findall(!iszero, iht_beta[correct_position]))
    iht_fp = iht_k_est - iht_tp
    iht_fn = k - iht_tp
    iht_precision = iht_tp / (iht_tp + iht_fp)
    iht_recall = iht_tp / (iht_tp + iht_fn)

    lasso_tp = length(findall(!iszero, lasso_beta[correct_position]))
    lasso_fp = lasso_k_est - lasso_tp
    lasso_fn = k - lasso_tp
    lasso_precision = lasso_tp / (lasso_tp + lasso_fp)
    lasso_recall = lasso_tp / (lasso_tp + lasso_fn)

    return iht_precision, iht_recall, lasso_precision, lasso_recall
end

function run()
    Random.seed!(2019)

    total_runs = 30
    iht_precision = zeros(total_runs)
    iht_recall = zeros(total_runs)
    lasso_precision = zeros(total_runs)
    lasso_recall = zeros(total_runs)
    for i in 1:total_runs
        println("current run = $i")
        n = 1000
        p = 10000
        k = 10
        ihtp, ihtr, lassop, lassor = iht_lasso_poisson(n, p, k)
        iht_precision[i] = ihtp
        iht_recall[i] = ihtr
        lasso_precision[i] = lassop
        lasso_recall[i] = lassor
    end

    writedlm("./poisson_results/iht_precision", iht_precision)
    writedlm("./poisson_results/iht_recall", iht_recall)
    writedlm("./poisson_results/lasso_precision", lasso_precision)
    writedlm("./poisson_results/lasso_precision", lasso_precision)
end

run()


text:
@doc """
The number of threads used by FFTW for CPU FFTs (default is the environment
variable `FFTW_NUM_THREADS`, or if that is not specified its
`Sys.CPU_THREADS÷2`). This must be set before creating any `FlatField` objects.
"""
FFTW_NUM_THREADS = nothing
@init global FFTW_NUM_THREADS = parse(Int,get(ENV,"FFTW_NUM_THREADS","$(Sys.CPU_THREADS÷2)"))

@doc """
Time-limit for FFT planning on CPU (default: 5 seconds). This must be set before
creating any `FlatField` objects.
"""
FFTW_TIMELIMIT = 5

# a set of wrapper FFT functions which use a @memoize'd plan
_fft_arr_type(arr) = basetype(typeof(parent(arr)))
m_rfft(arr::AbstractArray{T,N}, dims) where {T,N} = m_plan_rfft(_fft_arr_type(arr){T,N}, dims, size(arr)...) * arr
m_irfft(arr::AbstractArray{Complex{T},N}, d, dims) where {T,N} = m_plan_rfft(_fft_arr_type(arr){T,N}, dims, d, size(arr)[2:end]...) \ arr
m_rfft!(dst, arr::AbstractArray{T,N}, dims) where {T,N} = mul!(dst, m_plan_rfft(_fft_arr_type(arr){T,N}, dims, size(arr)...), arr)
m_irfft!(dst, arr::AbstractArray{Complex{T},N}, dims) where {T,N} = ldiv!(dst, m_plan_rfft(_fft_arr_type(arr){T,N}, dims, size(dst)...), copy_if_fftw(arr))
@memoize function m_plan_rfft(::Type{A}, dims::Dims, sz...) where {T, N, A<:AbstractArray{T,N}, Dims}
    FFTW.set_num_threads(FFTW_NUM_THREADS)
    plan_rfft(
        A(undef, sz...), dims; (A <: Array ? (timelimit=FFTW_TIMELIMIT,) : ())...
    ) #:: FFTW.rFFTWPlan{T,-1,false,N,Dims} # type assert to help inference in @memoized function
end
# FFTW (but not MKL) destroys the input array for inplace inverse real
# FFTs, so we need a copy. see https://github.com/JuliaMath/FFTW.jl/issues/158
copy_if_fftw(x) = FFTW.fftw_provider==:fftw ? copy(x) : x

"""
Convert a matrix A which is the output of a real FFT to a real vector, keeping
only unqiue real/imaginary entries of A
"""
function rfft2vec(A::AbstractMatrix)
    m,n = size(A)
    ireal,iimag = fftsyms(Val(m),Val(n))
    [real(A[ireal]); imag(A[iimag])]
end

"""
Convert a vector produced by rfft2vec back into a complex matrix.
"""
function vec2rfft(v::AbstractVector{<:Real})
    n = round(Int,sqrt(length(v)))
    m = n÷2+1
    nreal = (n^2)÷2 + (iseven(n) ? 2 : 1)
    ireal,iimag,inegks = fftsyms(Val(m),Val(n))
    A = fill(NaN+im*NaN,m,n)
    A[ireal] = v[1:nreal]
    A[iimag] += im*v[nreal+1:end]
    for i=1:m, j=1:n
        if isnan(A[i,j])
            A[i,j] = A[inegks[i,j]...]'
        end
    end
    A
end

"""
Convert an M×N matrix (with M=N÷2+1) which is the output a real FFT to a full
N×N one via symmetries.
"""
unfold(Tls::AbstractArray{<:Any,3}, Ny) = mapslices(X -> unfold(X, Ny), Array(Tls), dims=(1,2))
unfold(Tl::AbstractMatrix{<:Any}, Ny) = unfold(Array(Tl), Ny)
function unfold(Tl::Matrix{<:Any}, Ny::Int)
    m,n = size(Tl)
    @assert m==Ny÷2+1
    m2 = iseven(Ny) ? 2m : 2m+1
    n2 = iseven(n) ?  n+2 : n+3
    Tlu = similar(Tl,Ny,n)
    Tlu[1:m,1:n] = Tl
    @inbounds for i=m+1:Ny
        Tlu[i,1] = Tl[m2-i, 1]'
        @simd for j=2:n
            Tlu[i,j] = Tl[m2-i, n2-j]'
        end
    end
    Tlu
end

"""
Arguments `m` and `n` refer to the sizes of an `m`×`n` matrix (call it `A`) that is the output of a
real FFT (thus `m=n÷2+1`)

Returns a tuple of (ireal, iimag, negks) where these are

* `ireal` — `m`×`n` mask corrsponding to unique real entries of `A`
* `iimag` — `m`×`n` mask corrsponding to unique imaginary entries of `A`
* `negks` — `m`×`n` matrix of giving the index into A where the negative k-vector
            is, s.t. `A[i,j] = A[negks[i,j]]'`
"""
@generated function fftsyms(::Val{m},::Val{n}) where {m,n}
    k = ifftshift(-n÷2:(n-1)÷2)
    ks = tuple.(k',k)[1:n÷2+1,:]
    wrapk(k) = mod(k+n÷2,n) - n÷2
    negk(k) = @. wrapk(-k)
    k_in_ks(k) = -n÷2<=k[1]<=(n-1)÷2 && (0<=k[2]<=(n-1)÷2 || k[2]==-n÷2)
    ireal = map(k->(negk(k)==k || !k_in_ks(negk(k)) || (k[1]>0 || k[2]>0)), ks)
    iimag = map(k->(negk(k)!=k && !k_in_ks(negk(k)) || (k[1]>0 || k[2]>0)), ks)
    indexof(k) = (mod(k[2],n)+1, mod(k[1],n)+1)
    inegks = indexof.(negk.(ks))
    inegks[.!k_in_ks.(negk.(ks))] .= Ref((0,0))
    ireal,iimag,inegks#,ks,negk.(ks)#,k_in_ks.(negk.(ks)),map(k->k_in_ks(negk(k)),ks)
end

"""
    rfft_degeneracy_fac(n)

Returns an Array which is 2 if the complex conjugate of the
corresponding entry in the half-plane real FFT appears in the
full-plane FFT, and is 1 othewise. `n` is the length of the first
dimension of the full-plane FFT. The following identity holds:

    sum(abs2.(fft(x)) = sum(rfft_degeneracy_fac(size(x,1)) .* abs2.(rfft(x))
"""
function rfft_degeneracy_fac(n)
    if iseven(n)
        [1; fill(2,n÷2-1); 1]
    else
        [1; fill(2,n÷2)]
    end
end


text:
module SceneTrees

export SceneNode

mutable struct SceneNode
    object::Nullable{Vector{UInt8}}
    transform::Nullable{Vector{UInt8}}
    children::Dict{String, SceneNode}
end

SceneNode() = SceneNode(nothing, nothing, Dict{String, SceneNode}())

Base.getindex(s::SceneNode, name::AbstractString) = get!(SceneNode, s.children, name)
function Base.delete!(s::SceneNode, name::AbstractString)
    if haskey(s.children, name)
        delete!(s.children, name)
    end
end

function Base.getindex(s::SceneNode, path::AbstractVector{<:AbstractString})
    if length(path) == 0
        return s
    else
        return s[first(path)][path[2:end]]
    end
end

function Base.delete!(s::SceneNode, path::AbstractVector{<:AbstractString})
    parent = s[path[1:end-1]]
    child = path[end]
    delete!(parent, child)
end

function Base.foreach(f::Function, s::SceneNode)
    f(s)
    for (k, v) in s.children
        foreach(f, v)
    end
end

end


text:
####
# This file shows how the ∞-dimensional QL decomposition can be used
# for spectral theory of Toeplitz operators. We investigate the
# QL decomposition on each of the examples in Trefethen & Embree.
#

using InfiniteLinearAlgebra, BandedMatrices, PyPlot

###
# Basic routines for plotting
###

function ℓ11(A,λ; kwds...)
    try
        abs(ql(A-λ*I; kwds...).L[1,1])
    catch DomainError
        -1.0
    end
end

function branch(k)
    function(λ)
        j = sortperm(λ)[end-k+1]
        λ[j], j
    end
end

qlplot(A::AbstractMatrix; kwds...) = qlplot(BandedMatrix(A); kwds...)
function qlplot(A::BandedMatrix; branch=findmax, x=range(-4,4; length=200), y=range(-4,4;length=200), kwds...)
    z = ℓ11.(Ref(A), x' .+ y.*im; branch=branch)
    contourf(x,y,z; kwds...)
end

toepcoeffs(A::BandedMatrix) = InfiniteLinearAlgebra.rightasymptotics(A.data).args[1]
toepcoeffs(A::Adjoint) = reverse(toepcoeffs(A'))

function symbolplot(A::BandedMatrix; kwds...)
    l,u = bandwidths(A)
    a = toepcoeffs(A)
    θ = range(0,2π; length=1000)
    i = map(t -> dot(exp.(im.*(u:-1:-l).*t),a),θ)
    plot(real.(i), imag.(i); kwds...)
end

###
# Trefethen & Embree example
###

A = BandedMatrix(1 => Fill(2im,∞), 2 => Fill(-1,∞), 3 => Fill(2,∞), -2 => Fill(-4,∞), -3 => Fill(-2im,∞))
clf(); qlplot(A; x=range(-10,7; length=100), y=range(-7,8;length=100)); symbolplot(A; color=:black); title("Trefethen & Embree")
clf(); qlplot(transpose(A); x=range(-10,7; length=200), y=range(-7,8;length=200)); symbolplot(A; color=:black); title("Trefethen & Embree, transpose")

###
# limaçon
###

A = BandedMatrix(-2 => Fill(1.0,∞), -1 => Fill(1.0,∞), 1 => Fill(eps(),∞))
qlplot(A; x=range(-2,3; length=100), y=range(-2.5,2.5;length=100)); symbolplot(A; color=:black); title("Limacon")
clf(); qlplot(A; branch=findsecond, x=range(-2,3; length=100), y=range(-2.5,2.5;length=100)); symbolplot(A; color=:black); title("Limacon, second branch")
clf(); qlplot(transpose(A); x=range(-2,3; length=100), y=range(-2.5,2.5;length=100)); symbolplot(A; color=:black); title("Limacon, transpose")

###
# bull-head
###

A = BandedMatrix(-3 => Fill(7/10,∞), -2 => Fill(1,∞), 1 => Fill(2im,∞))
clf(); qlplot(A; x=range(-10,7; length=100), y=range(-7,8;length=100)); symbolplot(A; color=:black); title("Bull-head")
clf(); qlplot(transpose(A); x=range(-10,7; length=100), y=range(-7,8;length=100)); symbolplot(A; color=:black); title("Bull-head, transpose")

###
# Grcar
###

A = BandedMatrix(-3 => Fill(1,∞), -2 => Fill(1,∞), -1 => Fill(1,∞), 0 => Fill(1,∞), 1 => Fill(-1,∞))
clf(); qlplot(A; x=range(-4,5; length=100), y=range(-6,5;length=100)); symbolplot(A; color=:black); title("Grcar")
clf(); qlplot(A; branch=branch(2), x=range(-4,5; length=100), y=range(-6,5;length=100)); symbolplot(A; color=:black); title("Grcar, branch 2")
clf(); qlplot(A; branch=branch(3), x=range(-4,5; length=100), y=range(-6,5;length=100)); symbolplot(A; color=:black); title("Grcar, branch 3")
clf(); qlplot(transpose(A); x=range(-4,5; length=100), y=range(-6,5;length=100)); symbolplot(A; color=:black); title("Grcar, transpose")

###
# Triangle
###

A = BandedMatrix(-2 => Fill(1/4,∞), 1 => Fill(1,∞))
clf(); qlplot(A; x=range(-2,2; length=100), y=range(-2,2;length=100)); symbolplot(A; color=:black); title("Triangle")
clf(); qlplot(transpose(A); x=range(-2,2; length=100), y=range(-2,2;length=100)); symbolplot(A; color=:black); title("Triangle, transpose")
clf(); qlplot(transpose(A); branch=branch(2), x=range(-2,2; length=100), y=range(-2,2;length=100)); symbolplot(A; color=:black); title("Triangle, transpose, branch 2")

###
# Whale
###

A = BandedMatrix(-4 => Fill(im,∞), -3 => Fill(4,∞), -2 => Fill(3+im,∞), -1 => Fill(10,∞),
                    1 => Fill(1,∞), 2 => Fill(im,∞), 3 => Fill(-(3+2im),∞), 4=>Fill(-1,∞))

clf(); qlplot(A; x=range(-15,20; length=100), y=range(-20,20;length=100)); symbolplot(A; color=:black); title("Whale")
clf(); qlplot(transpose(A); x=range(-15,20; length=100), y=range(-20,20;length=100)); symbolplot(A; color=:black); title("Whale, transpose")

###
# Butterfly
###

A = BandedMatrix(-2 => Fill(1,∞), -1 => Fill(-im,∞), 1 => Fill(im,∞), 2 => Fill(-1,∞))
clf(); qlplot(A; x=range(-3,3; length=100), y=range(-3,3;length=100)); symbolplot(A; color=:black); title("Butterfly")
clf(); qlplot(A; branch=branch(2), x=range(-3,3; length=100), y=range(-3,3;length=100)); symbolplot(A; color=:black); title("Butterfly, branch 2")
clf(); qlplot(transpose(A); x=range(-3,3; length=100), y=range(-3,3;length=100)); symbolplot(A; color=:black); title("Butterfly")
clf(); qlplot(transpose(A); branch=branch(2), x=range(-3,3; length=100), y=range(-3,3;length=100)); symbolplot(A; color=:black); title("Butterfly")


text:
@blueprint struct Citation <: MustacheBlueprint
    template::String = TemplateFile(file="CITATION.bib")
    readme::Bool = false
end

blueprint_view(p::Citation, ctx::Context) = Dict(
    "AUTHORS" => join(ctx.kwargs.authors, ", "),
    "MONTH" => month(today()),
    "PKG" => ctx.kwargs.name,
    "URL" => "https://$(ctx.kwargs.host)/$(ctx.kwargs.user)/$(ctx.kwargs.name).jl",
    "YEAR" => year(today()),
)


text:
# Contains the configuration module, holding global settings for the package

"Configuration options"
module config
# Data directories
"The directory in which .jld data files are stored"
const datadir = normpath(joinpath(dirname(@__FILE__), "..", "data"))
"The directory in which raw data are stored in text format"
const rawdata = joinpath(datadir, "raw")
"The directory in which test data are stored"
const testdata = normpath(joinpath(dirname(@__FILE__), "..", "test", "data"))

# Gridded EOS parameters
"Grid resolution for making gridded equations of state"
const grid_resolution = 256
"Minimum pressure used in constructing the full EOS [Pa]"
const Pmin = 10^(4.5)
"Maximum pressure used in constructing the full EOS [Pa]"
const Pmax = 1e14
"Minimum temperature used in constructing the full EOS [K]"
const Tmin = 275
"Maximum temperature used in constructing the full EOS [K]"
const Tmax = 25000

end # module config


text:
# adapted from http://www.gis.usu.edu/~chrisg/python/2009/lectures/ospy_hw1a.py

import GDALUtils; const GU = GDALUtils

GU.registerdrivers() do
    open("tmp/hw1a.txt", "w") do file
        GU.read("ospy/data1/sites.shp") do datasource
            layer = GU.fetchlayer(datasource, 0)
            for feature in layer
                id = GU.fetchfield(feature, 0)
                cover = GU.fetchfield(feature, 1)
                geom = GU.fetchgeomfield(feature, 0)
                (x,y) = GU.getpoint(geom, 0)
                write(file, "$id $x $y $cover\n")
            end
        end
    end
end

GU.registerdrivers() do
    open("tmp/hw1a_v2.txt", "w") do file
        GU.read("ospy/data1/sites.shp") do datasource
            layer = GU.fetchlayer(datasource, 0)
            for feature in GU.DictIterator(layer)
                id = feature[:properties][0]
                cover = feature[:properties][1]
                (x,y) = GU.getpoint(feature[:geom][1], 0)
                write(file, "$id $x $y $cover\n")
            end
        end
    end
end


text:
import Base: firstindex, getindex, in, iterate, lastindex, length, show,
    summary, show_vector

using DataStructures

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

"""
    na(X)

The number of agents represented by (e.g., profile or allocation) `X`.
"""
function na end

"""
    agents(X)

Returns the set of agents associated with (e.g., profile or allocation)
`X`), as an iterable of `Int`s.
"""
agents(X) = 1:na(X)

"""
    ni(X)

The number of items represented by (e.g., profile or allocation) `X`.
"""
function ni end

"""
    items(X)

Returns the set of items associated with (e.g., profile or allocation) `X`,
as an iterable of `Int`s.
"""
items(X) = 1:ni(X)

# Used by Allocation and Additive
function show_agents_and_items(io::IO, X)
    n, m = na(X), ni(X)
    print(io,
        n, (n == 1 ? " agent" : " agents"),
        " and ",
        m, (m == 1 ? " item" : " items"))
end

## Allocations ###############################################################

"""
    struct Allocation <: Any

A mapping `A` from agents `i` to their assigned bundles `bundle(A, i)`. Agents
and items are represented as `Int`s, and bundles as `Set`s of `Int`s. The
`Allocation` also maintains an inverse mapping, from items `g` to their set of
owners, `owners(A, g)`. To keep these in sync, the bundles should be modified
using `give!` and `deny!`, rather than altering the bundle sets directly.
"""
struct Allocation
    # (Might make more sense to use `bundles` here, but to be consistent, we'd
    # then also need to rename `owners`...)
    bundle::Vector{Set{Int}}    # bundle[i]: Agent i's bundle
    owners::Vector{Set{Int}}    # owners[g]: Item g's owners
end

"""
    Allocation(n::Int, m::Int)

Construct an empty allocation with `n` agents and `m` items.
"""
Allocation(n::Int, m::Int) =
    Allocation([Set{Int}() for i = 1:n], [Set{Int}() for i = 1:m])

"""
    Allocation()

Construct an empty allocation with zero agents and items.
"""
Allocation() = Allocation(0, 0)

bundle(A) = A.bundle
owners(A) = A.owners

na(A::Allocation) = length(bundle(A))
ni(A::Allocation) = length(owners(A))

## Allocation printing

# Rather than reinventing printing code, we'll let built-in functions for
# AbstractDict and AbstractVector much of the heavy lifting, using a couple of
# one-off wrappers (AllocShowWrap and BundleShowWrap).

function show_bundle(io::IO, S)
    seq = sort(collect(S))
    # Supplying :typeinfo to "pretend" we've already printed the type info
    ctx = IOContext(io, :typeinfo => typeof(seq))
    show_vector(ctx, seq, "{", "}")
end

# Single-use -- assumes alloc will not be modified after construction, because
# this won't be reflected in pairs. This wrap is used to leverage the
# MIME"text/plain" printing of an AbstractDict; it is not used with single-line
# printing of Allocations (which is implemented more directly).
struct AllocShowWrap <: AbstractDict{Int,Set{Int}}
    alloc::Allocation
    pairs
end
AllocShowWrap(A) =
    AllocShowWrap(A, enumerate(BundleShowWrap.(A.bundle)))

summary(io::IO, w::AllocShowWrap) = summary(io, w.alloc)
show(io::IO, w::AllocShowWrap) = show(io, w.alloc)
length(w::AllocShowWrap) = na(w.alloc)
iterate(w::AllocShowWrap, args...) = iterate(w.pairs, args...)

struct BundleShowWrap
    bundle
end

show(io::IO, w::BundleShowWrap) = show_bundle(io, w.bundle)

function summary(io::IO, A::Allocation)
    print(io, "Allocation with ")
    show_agents_and_items(io, A)
    u = length([g for g in items(A) if !owned(A, g)])
    if u != 0
        print(io, ", ", u, " unallocated")
    end
end

function show(io::IO, A::Allocation)
    seq = BundleShowWrap.(A.bundle)
    # Supplying :typeinfo to "pretend" we've already printed the type info
    ctx = IOContext(io, :typeinfo => typeof(seq))
    show_vector(ctx, seq)
end

show(io::IO, m::MIME"text/plain", A::Allocation) = show(io, m, AllocShowWrap(A))

## Allocation accessors and manipulators

"""
    bundle(A, i)

The set of items allocated to agent `i` in the allocation `A`. The returned
`Set` should be treated as read-only.
"""
bundle(A, i) = bundle(A)[i]

"""
    owners(A, g)

The set of agents to which item `g` has been allocated in the allocation `A`.
The returned `Set` should be treated as read-only.
"""
owners(A, g) = owners(A)[g]

"""
    owner(A, g)

The agent to which item `g` has been allocated in the allocation `A`. Will
produce an error if `g` has been allocated to more than one agent.
"""
owner(A, g) = only(owners(A, g))

"""
    owned(A, g)

Whether or not the item `g` is owned by any agent in the allocation `A`.
"""
owned(A, g) = !isempty(owners(A, g))

"""
    give!(A, i, g::Int)

Give agent `i` the object `g` in the `Allocation` `A`.
"""
function give!(A, i, g::Int)
    push!(bundle(A, i), g)
    push!(owners(A, g), i)
    return A
end

"""
    give!(A, i, B)

Give agent `i` the bundle `B` in the `Allocation` `A`.
"""
function give!(A, i, B)
    union!(bundle(A, i), B)
    for g in B
        push!(owners(A, g), i)
    end
    return A
end

"""
    deny!(A, i, g)

Deny agent `i` the object `g`, which it has previously been given, in the
allocation `A`.
"""
function deny!(A, i, g)
    delete!(bundle(A, i), g)
    delete!(owners(A, g), i)
    return A
end

"""
    fill_even!(A)

Fill out the allocation by distributing the unallocated items evenly, by
repeatedly giving the next unallocated item to the agent with the fewest items
(ties broken arbitrarily).
"""
function fill_even!(A)
    n, m = na(A), ni(A)
    pq = PriorityQueue{Int, Int}()
    for i = 1:n
        enqueue!(pq, i, length(bundle(A, i)))
    end
    for g = 1:m
        owned(A, g) && continue
        i, u = peek(pq)
        give!(A, i, g)
        pq[i] = u + 1
    end
    return A
end

"""
    fill_random!(A)

Fill out the allocation by distributing the unallocated items randomly.
"""
function fill_random!(A)
    N, M = agents(A), items(A)
    for g in M
        owned(A, g) && continue
        i = rand(N)
        give!(A, i, g)
    end
end

## Valuation profiles ######################################################

"""
    abstract struct Profile <: Any

An abstract type representing an valuation profile. Which functions are used to
query it depends on the kind of valuation functions it represents. Additive
valuations act on individual objects, and simply sum those values over a bundle,
but profiles with quite different kinds of queries are possible for valuations
with other properties (see, e.g., [Fair Allocation of Indivisible Goods:
Improvements and
Generalizations](https://dl.acm.org/doi/10.1145/3219166.3219238) by Ghodsi et
al., 2018).
"""
abstract type Profile end

"""
    Allocation(V::Profile)

Construct an empty allocation with a number of agents and items equal to that of
the instance (i.e., profile) `V`.
"""
Allocation(V::Profile) = Allocation(na(V), ni(V))

"""
    value(V::Profile, i, S)
    value(V::Profile, i, g::Int)

The value agent `i` places on bundle `S`, according to the profile `V`. The
second form is a shortcut for `value(V, i, [g])`; the shortcut will generally be
more efficient. Note that the value of `S` may *not* in general be the sum of
the values of its items; that property is unique to `Additive` profiles.
"""
function value end

"""
    value(V::Profile, i, A::Allocation)

The value agent `i` receives in allocation `A`, under the profile `V`.
"""
value(V::Profile, i, A::Allocation) = bundle_value(V, i, A)

# Use for disambiguation definitions for subtypes:
bundle_value(V, i, A) = value(V, i, bundle(A, i))

"""
    value_1(V::Profile, i, S)

The value agent `i` places on bundle `S`, *up to one item*, that is, the
smallest value `i` can place on bundle `S` after removing (at most) one item,
according to the profile `V`.
"""
function value_1 end

"""
    value_x(V::Profile, i, S)

The value agent `i` places on bundle `S`, *up to any item*, that is, the
largest value `i` can place on bundle `S` after removing one item (or no
items, if the bundle is empty), according to the profile `V`.
"""
function value_x end

"""
    isintegral(V::Profile)

Test whether every value provided by `V` is an integer.
"""
function isintegral end

"""
    isnonnegative(V::Profile)

Test whether every value provided by `V` is nonnegative.
"""
function isnonnegative end

"""
    matrix(V::Profile)

Return a matrix `X` where `X[i, g]` is `value(V, i, g)`. May not be very useful
in general (especially if calculating single-item values isn't efficient to
begin with), but if such a matrix is available as part of the profile
implementation (as with `Additive`), it may be returned directly.
"""
matrix(V::Profile) = [value(V, i, g) for i in agents(V), g in items(V)]

"""
    struct Additive{T <: AbstractMatrix} <: Profile

An additive valuation profile, representing how each agent values all possible
bundles. Because of additivity, this is easily "lifted" from the values of
individual items, by addition, with an empty bundle given a value of zero. By
default, the profile is constructed from a real matrix `X`, supplied to the
default constructor, where `X[i, g]` is agent `i`'s value for item `g`.
"""
struct Additive{T <: AbstractMatrix{<:Real}} <: Profile
    values::T
end

"""
    Additive(n, m)

Create an additive profile for `n` agents and `m` items where all values are
set to zero.
"""
Additive(n, m) = Additive(zeros(n, m))

"""
    Profile(X::Matrix)

Alias for `Additive(X)`.
"""
Profile(X::Matrix) = Additive(X)

na(V::Additive) = size(V.values, 1)
ni(V::Additive) = size(V.values, 2)

function summary(io::IO, A::Additive{T}) where {T}
    print(io, "Additive{$T}")
end

function show(io::IO, ::MIME"text/plain", A::Additive)
    summary(io, A)
    print(io, " with ")
    show_agents_and_items(io, A)
    println(":")
    Base.print_matrix(io, A.values)
end

isintegral(V::Additive) = all(isinteger.(V.values))

isnonnegative(V::Additive) = all(V.values .>= zero(eltype(V.values)))

"""
    value(V::Additive, i, g::Int)

The value of item `g`, according to agent `i` under valuation profile `V`.
"""
value(V::Additive, i, g::Int) = V.values[i, g]

# Disambiguation:
value(V::Additive, i, A::Allocation) = bundle_value(V, i, A)

# The bundle value is "lifted" from item values by addition.
value(V::Additive, i, S) =
    isempty(S) ? zero(eltype(V.values)) : sum(value(V, i, g) for g in S)

# For the additive case, we can just subtract the highest item value.
value_1(V::Additive, i, S) =
    value(V, i, S) -
    (isempty(S) ? zero(eltype(V.values)) : maximum(value(V, i, g) for g in S))

# For the additive case, we can just subtract the lowest item value.
value_x(V::Additive, i, S) =
    value(V, i, S) -
    (isempty(S) ? zero(eltype(V.values)) : minimum(value(V, i, g) for g in S))

"""
    value!(V::Additive, i, g::Int, v)

Set the value of item `g`, according to agent `i`, to `v` in profile `V`.
"""
value!(V::Additive, i, g, v) = V.values[i, g] = v

"""
    matrix(V::Additive)

Return the underlying valuation matrix of `V`.
"""
matrix(V::Additive) = V.values

"""
    normalize(V::Additive)

Scale the values of `V` such that ``v_i(M) = n`` for all agents ``i``.
"""
normalize(V::Additive) =
    Additive(V.values .* [na(V) / value(V, i, items(V)) for i in agents(V)])

"""
    struct Submodular <: Profile

A submodular valuation profile, representing how each agent values all possible
bundles. The profile is constructed from a set of `n` submodular valuation
functions, one per agent, as well as the number of items, `m`. The profile
functions should, when supplied with a `Set` of items (a subset of `1:m`),
return the value of that set of items to the given agent (i.e., acting as
so-called *query oracles*).
"""
struct Submodular <: Profile
    oracles::Vector{Function}
    m::Int
end

na(V::Submodular) = length(V.oracles)
ni(V::Submodular) = V.m

value(V::Submodular, i, S) = V.oracles[i](Set(S))
value(V::Submodular, i, g::Int) = value(V, i, Set(g))

## Constraints ###############################################################

"""
    abstract type Constraint <: Any

Abstract supertype of various kinds of constraints. An instance of the
allocation problem is assumed to consist of a `Profile` object and at most one
`Constraint` object, embodying any and all constraints placed on feasible
solutions.
"""
abstract type Constraint end

"""
    mutable struct Category

One of the categories in a `Counts` constraint, from which each agent can hold
at most a given number of items. The category supports iteration (over its
members), and the threashold is available through the `threshold` accessor.
"""
mutable struct Category
    members::Set{Int}
    threshold::Int
end

iterate(c::Category, args...) = iterate(c.members, args...)
length(c::Category) = length(c.members)

"""
    threshold(c::Category)

The maximum number of items any agent can receive from the given category, as
part of a `Counts` constraint.
"""
threshold(c::Category) = c.threshold

"""
    mutable struct OrderedCategory

Used in place of `Category` when handling an ordered instance. The instance is
assumed to be such that items in the range `index:index + n_items - 1` belong to
the given category, i.e., the items of a category occupy a contiguous range of
integers.
"""
mutable struct OrderedCategory
    index::Int
    n_items::Int
    threshold::Int
end

in(g, c::OrderedCategory) = c.index <= g < c.index + c.n_items
lastindex(c::OrderedCategory) = length(c)
length(c::OrderedCategory) = c.n_items

getindex(c::OrderedCategory, i::Int) =
    getindex(c.index:c.index + c.n_items - 1, i)

getindex(c::OrderedCategory, v::UnitRange{Int64}) =
    getindex(c.index:c.index + c.n_items - 1, v)

iterate(c::OrderedCategory, args...) =
    iterate(c.index:c.index + c.n_items - 1, args...)

"""
    floor_n(c::OrderedCategory, n)

One `n`th of the number of items in the category, rounded down.
"""
floor_n(c::OrderedCategory, n) = floor(Int, c.n_items/n)

"""
    ceil_n(c::OrderedCategory, n)

One `n`th of the number of items in the category, rounded up.
"""
ceil_n(c::OrderedCategory, n) = ceil(Int, c.n_items/n)

"""
    required(c::OrderedCategory, n)

The number of items the next agent must take in order to keep the instance
valid, i.e., for there to be a maximum of `(n - 1) * threshold` remaining
items.
"""
required(c::OrderedCategory, n) = max(c.n_items - (n - 1) * c.threshold, 0)

"""
    struct Counts{T} <: Constraint

The *cardinality constraints* introduced by Biswas and Barman in their 2018
paper [Fair Division Under Cardinality
Constraints](https://www.ijcai.org/proceedings/2018/13). This is a form of
constraint consisting of several `Category` objects, available through
indexing or iteration. Any agent may hold at most a given number of items from
any given category.
"""
struct Counts{T} <: Constraint
    categories::Vector{T}
end

"""
    Counts(args::Pair...)

Create a `Counts` object where each pair `x => k` becomes a category with
members `Set(x)` and threshold `k`.
"""
Counts(args::Pair...) = Counts([Category(Set(p[1]), p[2]) for p in args])

getindex(C::Counts, i) = C.categories[i]
length(C::Counts) = length(C.categories)
iterate(C::Counts, args...) = iterate(C.categories, args...)

"""
    struct Conflicts{T <: AbstractGraph} <: Constraint

A kind of constraint -- or set of constraints -- that indicates that certain
items conflict, and thus cannot be allocated to the same agent. The
constraints are represented as a *conflict graph*
(`LightGraphs.AbstractGraph`), with items as nodes, and edges representing
conflicts. The `Conflicts` type is just a wrapper for dispatch purposes, with
the underlying graph available through the `graph` accessor.
"""
struct Conflicts{T <: AbstractGraph} <: Constraint
    graph::T
end

"""
    graph(C::Conflicts)

Return the conflict graph wrapped by a `Conflicts` object.
"""
graph(C::Conflicts) = C.graph

## Reductions ################################################################

"""
    mutable struct Reduction{S, T, U, V}

A reduction from one instance of a fair allocation problem to another. Contains
information about the profiles in the reduced instance, through an object of
type `S`. There must exist functions `agents(s::S)` and `items(s::S)` that
return iterators of, respectively, the agents and items in the reduced instance.
The reduction can also contain information about the constraints in the reduced
instance, through an object of type `T`.

In addition, the reduction contains two mappings (vectors), `λi` (of type `I`)
and `λg` (of type `G`). Both types should be indexable (for `i ∈ agents(s)` and
`g ∈ items(s)`, respectively). `λi[i]` and `λg[g]` should return the agent and
item identifier in the original instance of, respectively, agent `i` and item
`g` in the reduced instance.

The reduction also contains a function that can convert an allocation in the
reduced instance to one in the original instance.

The default constructor is `Reduction(V, C, λi, λg, transform::Function)`.
"""
mutable struct Reduction{S, T, I, G}
    V::S
    C::T
    λi::I
    λg::G
    transform::Function
end

"""
    Reduction(V, λi, λg, transform)

A simplified constructor for when there are no constraints.
"""
Reduction(V, λi, λg, transform) = Reduction(V, nothing, λi, λg, transform)

"""
    Reduction(V, C)

A simplified constructor for when either no changes have been performed or
changes only concern the profiles and/or constraints.
"""
Reduction(V, C) = Reduction(V, C, agents(V), items(V), identity)

"""
    Reduction(V)

A simplified constructor for when either no changes have been performed or
changes only concern the profiles.
"""
Reduction(V) = Reduction(V, nothing)

"""
    profile(R::Reduction)

Returns the valuation profile for the reduced instance.
"""
profile(R::Reduction) = R.V

"""
    constraint(R::Reduction)

Returns the constraint object for the reduced instance
"""
constraint(R::Reduction) = R.C

"""
    transform(R::Reduction, A::Allocation)

Converts the given allocation for the reduced instance to one for original
instance. The way the convertion occurs depends on the given reduction.
"""
transform(R::Reduction, A::Allocation) = R.transform(A)

"""
    agent(R::Reduction, i)

Converts the agent identifier `i` from the reduced instance to the agent
identifier of the same agent in the original instance.
"""
agent(R::Reduction, i) = R.λi[i]

"""
    item(R::Reduction, g)

Converts the item identifier `g` from the reduced instance to the item
identifier of the same item in the original instance.
"""
item(R::Reduction, g) = R.λg[g]

"""
    chain(R₁::Reduction, R₂::Reduction)

Assumes that R₂ is a reduction of the reduced instance of R₁. Combines the two
reductions, so that the original instance is the original instance of R₁ and the
reduced instance is the reduced instance of R₂ (essentially diagram-order
composition of the reductions).
"""
function chain(R₁::Reduction, R₂::Reduction)
    return Reduction(
            R₂.V, R₂.C,
            [agent(R₁, agent(R₂, i)) for i in agents(R₂.V)],
            [item(R₁, item(R₂, g)) for g in items(R₂.V)],
            (A) -> transform(R₁, transform(R₂, A))
        )
end