text:
using Plexifications
using Base.Test: @test, @test_throws

@test_throws ErrorException begin
    inv(TriPerplex(0, 0.0))
end

@test begin
    iszerodivisor(TriPerplex(0,1,1,0))
end

@test begin
    !iszerodivisor(TriPerplex(1,2,3,5))
end

@test begin
    a = Plexifications.randomBigFloat()
    isreal(TriPerplex(a))
end

@test begin
    a = 1
    !isreal(TriPerplex(a, a))
end

@test begin
    z = random(TriPerplex{Int})
    z == +(z)
end

@test begin
    a = Plexifications.randomBigFloat()
    real(TriPerplex(a)) == a
end

@test begin
    io = IOBuffer()
    show(io, TriPerplex(1, 2, 3, 4, 5, 6, 7, 8))
    l = String(take!(io))
    r = "[1: 1, R: 2, Q: 3, QR: 4, P: 5, PR: 6, PQ: 7, PQR: 8]"
    l == r
end

@test begin
    x = random(TriPerplex{BigFloat})
    y = random(TriPerplex{BigFloat})
    x != y
end

@test begin
    a = Plexifications.randomBigInt()
    z = random(TriPerplex{BigInt})
    l = a + z
    r = z + a
    l == r
end

@test begin
    a = Plexifications.randomBigInt()
    z = random(TriPerplex{BigInt})
    l = a - z
    r = -(z - a)
    l == r
end

@test begin
    x = random(TriPerplex{BigInt})
    y = random(TriPerplex{BigInt})
    l = x + y
    r = y + x
    l == r
end

@test begin
    x = random(BiPerplex{BigInt})
    y = random(TriPerplex{BigInt})
    l = x + y
    r = y + x
    l == r
end

@test begin
    x = random(TriPerplex{BigInt})
    y = random(TriPerplex{BigInt})
    l = x * y
    r = y * x
    l == r
end

@test begin
    x = random(BiPerplex{BigInt})
    y = random(TriPerplex{BigInt})
    l = x * y
    r = y * x
    l == r
end

@test begin
    x = random(TriPerplex{BigInt})
    y = random(TriPerplex{BigInt})
    l = x - y
    r = -(y - x)
    l == r
end

@test begin
    x = random(BiPerplex{BigInt})
    y = random(TriPerplex{BigInt})
    l = x - y
    r = -(y - x)
    l == r
end

@test begin
    x = random(TriPerplex{BigInt})
    y = random(TriPerplex{BigInt})
    z = random(TriPerplex{BigInt})
    l = (x + y) + z
    r = x + (y + z)
    l == r
end

@test begin
    x = random(TriPerplex{BigInt})
    y = random(TriPerplex{BigInt})
    z = random(TriPerplex{BigInt})
    l = (x * y) * z
    r = x * (y * z)
    l == r
end

@test begin
    x = random(TriPerplex{BigInt})
    barR(barR(x)) == x
end

@test begin
    x = random(TriPerplex{BigInt})
    barQ(barQ(x)) == x
end

@test begin
    x = random(TriPerplex{BigInt})
    barP(barP(x)) == x
end

@test begin
    x = random(TriPerplex{BigInt})
    barQ(barR(x)) == barR(barQ(x))
end

@test begin
    x = random(TriPerplex{BigInt})
    barP(barR(x)) == barR(barP(x))
end

@test begin
    x = random(TriPerplex{BigInt})
    barQ(barP(x)) == barP(barQ(x))
end

@test begin
    x = random(TriPerplex{Rational{BigInt}})
    inv(inv(x)) == x
end

@test begin
    x = random(TriPerplex{BigInt})
    y = random(TriPerplex{BigInt})
    l = barR(x * y)
    r = barR(x) * barR(y)
    l == r
end

@test begin
    x = random(TriPerplex{BigInt})
    y = random(TriPerplex{BigInt})
    l = barR(x + y)
    r = barR(x) + barR(y)
    l == r
end

@test begin
    x = random(TriPerplex{BigInt})
    y = random(TriPerplex{BigInt})
    l = barR(x - y)
    r = barR(x) - barR(y)
    l == r
end

@test begin
    x = random(TriPerplex{BigInt})
    y = random(TriPerplex{BigInt})
    l = barQ(x * y)
    r = barQ(x) * barQ(y)
    l == r
end

@test begin
    x = random(TriPerplex{BigInt})
    y = random(TriPerplex{BigInt})
    l = barQ(x + y)
    r = barQ(x) + barQ(y)
    l == r
end

@test begin
    x = random(TriPerplex{BigInt})
    y = random(TriPerplex{BigInt})
    l = barQ(x - y)
    r = barQ(x) - barQ(y)
    l == r
end

@test begin
    x = random(TriPerplex{BigInt})
    y = random(TriPerplex{BigInt})
    l = barP(x * y)
    r = barP(x) * barP(y)
    l == r
end

@test begin
    x = random(TriPerplex{BigInt})
    y = random(TriPerplex{BigInt})
    l = barP(x + y)
    r = barP(x) + barP(y)
    l == r
end

@test begin
    x = random(TriPerplex{BigInt})
    y = random(TriPerplex{BigInt})
    l = barP(x - y)
    r = barP(x) - barP(y)
    l == r
end

@test begin
    x = random(TriPerplex{Rational{BigInt}})
    y = random(TriPerplex{Rational{BigInt}})
    l = inv(x * y)
    r = inv(y) * inv(x)
    l == r
end

@test begin
    x = random(TriPerplex{Rational{BigInt}})
    y = random(TriPerplex{Rational{BigInt}})
    l = inv(x / y)
    r = y / x
    l == r
end

@test begin
    x = random(TriPerplex{Rational{BigInt}})
    y = random(TriPerplex{Rational{BigInt}})
    l = inv(x \ y)
    r = y \ x
    l == r
end

@test begin
    x = random(TriPerplex{BigInt})
    y = random(TriPerplex{BigInt})
    z = random(TriPerplex{BigInt})
    l = x * (y + z)
    r = (x * y) + (x * z)
    l == r
end

@test begin
    a = Plexifications.randomBigInt()
    x = random(TriPerplex{BigInt})
    y = random(TriPerplex{BigInt})
    l = a * (x + y)
    r = (a * x) + (a * y)
    l == r
end

@test begin
    a = Plexifications.randomBigInt()
    x = random(TriPerplex{BigInt})
    l = a * x
    r = x * a
    l == r
end

@test begin
    x = random(TriPerplex{BigInt})
    y = random(TriPerplex{BigInt})
    z = random(TriPerplex{BigInt})
    l = (x + y) * z
    r = (x * z) + (y * z)
    l == r
end

@test begin
    x = random(TriPerplex{BigInt})
    y = random(TriPerplex{BigInt})
    z = random(TriPerplex{BigInt})
    l = x * (y - z)
    r = (x * y) - (x * z)
    l == r
end

@test begin
    x = random(TriPerplex{BigInt})
    y = random(TriPerplex{BigInt})
    z = random(TriPerplex{BigInt})
    l = (x - y) * z
    r = (x * z) - (y * z)
    l == r
end

@test begin
    x = random(TriPerplex{BigInt})
    n = norm(x)
    n > big(0) || n < big(0) || n == big(0)
end

@test begin
    x = random(TriPerplex{BigInt})
    y = random(TriPerplex{BigInt})
    l = norm(x * y)
    r = norm(x) * norm(y)
    l == r
end

@test_throws ErrorException begin
    inv(TriPerplex(0, 0.0))
end

@test_throws ErrorException begin
    random(TriPerplex{Int}) / 0
end

@test_throws ErrorException begin
    0 \ random(TriPerplex{Int})
end

@test begin
    a = Plexifications.randomBigRational()
    z = random(TriPerplex{Rational{BigInt}})
    l = z / a
    r = a \ z
    l == r
end

@test begin
    a = Plexifications.randomBigRational()
    z = random(TriPerplex{Rational{BigInt}})
    l = a / z
    r = z \ a
    l == r
end

@test begin
    z = random(TriPerplex{BigInt})
    z == z + zero(TriPerplex{BigInt})
end

@test begin
    z = random(TriPerplex{BigInt})
    z == z * one(TriPerplex{BigInt})
end

@test begin
    z = random(TriPerplex{BigInt})
    z == z + zero(z)
end

@test begin
    z = random(TriPerplex{BigInt})
    z == z * one(z)
end

@test begin
    x = random(TriPerplex{BigInt})
    y = random(TriPerplex{BigInt})
    l = commutator(x, y)
    r = zero(TriPerplex{BigInt})
    l == r
end

@test begin
    x = random(TriPerplex{BigInt})
    y = random(TriPerplex{BigInt})
    l = alternatorL(x, y)
    r = zero(TriPerplex{BigInt})
    l == r
end

@test begin
    x = random(TriPerplex{BigInt})
    y = random(TriPerplex{BigInt})
    l = alternatorR(x, y)
    r = zero(TriPerplex{BigInt})
    l == r
end

@test begin
    x = random(TriPerplex{BigInt})
    y = random(TriPerplex{BigInt})
    l = flexator(x, y)
    r = zero(TriPerplex{BigInt})
    l == r
end

@test begin
    x = random(TriPerplex{BigInt})
    y = random(TriPerplex{BigInt})
    z = random(TriPerplex{BigInt})
    l = associator(x, y, z)
    r = zero(TriPerplex{BigInt})
    l == r
end

@test begin
    x = random(TriPerplex{BigInt})
    y = random(TriPerplex{BigInt})
    z = random(TriPerplex{BigInt})
    l = jacobiator(x, y, z)
    r = 3 * x * y * z
    l == r
end

@test begin
    l = TriPerplex(Prplex(1,2))
    r = TriPerplex(1,2)
    l == r
end

@test begin
    l = TriPerplex(BiPerplex(1,2))
    r = TriPerplex(1,2,0)
    l == r
end

@test begin
    l = TriPerplex(1, 2, 3.0)
    r = TriPerplex(1.0, 2.0, 3.0)
    l == r
end

@test begin
    l = TriPerplex(1, 2, 3, 4.0)
    r = TriPerplex(1.0, 2.0, 3.0, 4.0)
    l == r
end

@test begin
    l = TriPerplex(1, 2, 3, 4, 5.0)
    r = TriPerplex(1.0, 2.0, 3.0, 4.0, 5.0)
    l == r
end

@test begin
    l = TriPerplex(1, 2, 3, 4, 5, 6.0)
    r = TriPerplex(1.0, 2.0, 3.0, 4.0, 5.0, 6.0)
    l == r
end

@test begin
    l = TriPerplex(1, 2, 3, 4, 5, 6, 7.0)
    r = TriPerplex(1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0)
    l == r
end

@test begin
    l = TriPerplex(1, 2, 3, 4, 5, 6, 7, 8.0)
    r = TriPerplex(1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0)
    l == r
end


text:
"""
```
t = find_threshold(Moments(),  histogram, edges)
```

The following rule determines the threshold:  if one assigns all observations
below the threshold to a value z₀ and all observations above the threshold to a
value z₁, then the first three moments of the original histogram must match the
moments of this specially constructed bilevel histogram.

# Output

Returns a real number `t` that specifies the threshold.

# Arguments

The function arguments are described in more detail below.

##  `histogram`

An `AbstractArray` storing the frequency distribution.

##  `edges`

An `AbstractRange` specifying how the intervals for the frequency distribution
are divided.

# Example

Compute the threshold for the "camerman" image in the `TestImages` package.

```julia
using TestImages, ImageContrastAdjustment, HistogramThresholding

img = testimage("cameraman")
edges, counts = build_histogram(img,256)
#=
  The `counts` array stores at index 0 the frequencies that were below the
  first bin edge. Since we are seeking a threshold over the interval
  partitioned by `edges` we need to discard the first bin in `counts`
  so that the dimensions of `edges` and `counts` match.
=#
t = find_threshold(Moments(),counts[1:end], edges)
```

# Reference

[1] W.-H. Tsai, “Moment-preserving thresolding: A new approach,” Computer Vision, Graphics, and Image Processing, vol. 29, no. 3, pp. 377–393, Mar. 1985. [doi:10.1016/0734-189x(85)90133-1](https://doi.org/10.1016/0734-189x%2885%2990133-1)
"""
function find_threshold(algorithm::Moments, histogram::AbstractArray, edges::AbstractRange)
    n = sum(histogram)
    m₀ = 1.0
    m₁, m₂, m₃ =  compute_moments(histogram, edges, n)
    c₀, c₁ = compute_auxiliary_values(m₀, m₁, m₂, m₃)
    z₀, z₁ = compute_gray_levels(c₀, c₁)
    p₀, p₁ = compute_pixel_fractions(z₀, z₁, m₁)
    threshold_percentile = p₀
    target_count = ceil(Int,threshold_percentile * n)
    final_bin = determine_threshold(histogram, target_count)
    return edges[final_bin]
end

function determine_threshold(histogram, target_count)
    #= Finding the bin for which the bins below and including it contain p₀
     pixels. =#
    current_count = 0
    temp_bin = 0
    final_bin = 0
    holding_pixel_value = 0
    while current_count <= target_count
        temp_bin += 1
        current_count += histogram[temp_bin]
        if abs(target_count - current_count) <= abs(target_count - holding_pixel_value)
            holding_pixel_value = current_count
            final_bin = temp_bin
        end
    end
    final_bin
end

function compute_moments(histogram, edges, n)
    m₁ = 0.0
    m₂ = 0.0
    m₃ = 0.0
    # (2) Compute the moments from a histogram.
    for j in eachindex(histogram)
        m₁ += (histogram[j] / n) * edges[j]
        m₂ += (histogram[j] / n) * (edges[j])^2
        m₃ += (histogram[j] / n) * (edges[j])^3
    end
    m₁, m₂, m₃
end

function compute_auxiliary_values(m₀, m₁, m₂, m₃)
    # Appendix A.1.i calculates the auxiliary values c.
    cₐ = det([m₀ m₁;
              m₁ m₂])
    c₀ = (1/cₐ) * det([-m₂ m₁;
                       -m₃ m₂])
    c₁ = (1/cₐ) * det([m₀ -m₂;
                       m₁ -m₃])
    c₀, c₁
end

function compute_gray_levels(c₀, c₁)
    #= Appendix A.1.ii calculates grey values z which are the bi-level intensities
     associated with the hypothetical "unblurred" image. =#
    z₀ = (1/2) * ( -c₁ - sqrt(c₁^2 - 4*c₀))
    z₁ = (1/2) * ( -c₁ + sqrt(c₁^2 - 4*c₀))
    z₀, z₁
end

function compute_pixel_fractions(z₀, z₁, m₁)
    #= Appendix A.1.iii calculates the fraction of pixels below and above the
     threshold. =#
    pₐ = det([1 1;
              z₀ z₁])
    p₀ = (1/pₐ) * det([1 1;
                       m₁ z₁])
    p₁ = 1 - p₀
    p₀, p₁
end

function preserved_moments(p₀, z₀, p₁, z₁)
    # (3) Preserving the moments of the resulting image.
    m₁₀ = p₀ * (z₀)^0 + p₁ * (z₁)^0
    m₁₁ = p₀ * (z₀)^1 + p₁ * (z₁)^1
    m₁₂ = p₀ * (z₀)^2 + p₁ * (z₁)^2
    m₁₃ = p₀ * (z₀)^3 + p₁ * (z₁)^3
    m₁₀, m₁₁, m₁₂, m₁₃
end

ext: jl


text:
module Templates

    export DenseBNReLU, MLPEncoder, MLPDecoder, MLPAutoEncoder

    include("mlp.jl")

end

max_stars_repo_name: davnn/OutlierDetectionNetworks

max_issues_repo_name: davnn/OutlierDetectionNetworks

max_forks_repo_name: davnn/OutlierDetectionNetworks


text:
using SafeTestsets

@safetestset "basic symbolics" begin include("basic_symbolic_tests.jl") end
@safetestset "moment conversion" begin include("moment_convert_tests.jl") end
@safetestset "reaction systems" begin include("reaction_system_tests.jl") end
@safetestset "moment eqs" begin include("moment_equations_tests.jl") end
@safetestset "moment closure" begin include("closure_tests.jl") end
@safetestset "moment closure (bernoulli)" begin include("bernoulli_closure_tests.jl") end
@safetestset "conditional closure" begin include("conditional_closure_tests.jl") end
@safetestset "sample moments" begin include("sample_moment_tests.jl") end
@safetestset "linear mapping approximation" begin include("lma_tests.jl") end
@safetestset "moment eqs for SDEs" begin include("moment_equations_SDE_tests.jl") end
#@safetestset "polynomialization" begin include("polynomialization_tests.jl") end
@safetestset "latexify" begin include("latexify_tests.jl") end


text:
module LuaCutscenesLuaTalker

using ..Ahorn, Maple

@mapdef Entity "luaCutscenes/luaTalker" LuaTalker(x::Integer, y::Integer, width::Integer=Maple.defaultTriggerWidth, height::Integer=Maple.defaultTriggerHeight, onlyOnce::Bool=false, filename::String="", arguments::String="", unskippable::Bool=false)

const placements = Ahorn.PlacementDict(
    "Lua Talker (Lua Cutscenes)" => Ahorn.EntityPlacement(
        LuaTalker,
        "rectangle",
        Dict{String, Any}(),
        function(entity)
            entity.data["nodes"] = [(Int(entity.data["x"]) + Int(entity.data["width"]) + 8, Int(entity.data["y"]))]
        end
    )
)

Ahorn.nodeLimits(entity::LuaTalker) = 1, 1

Ahorn.minimumSize(entity::LuaTalker) = 8, 8
Ahorn.resizable(entity::LuaTalker) = true, true

const hoverTexture = "objects/LuaCutscenes/hover_idle"

function Ahorn.selection(entity::LuaTalker)
    x, y = Ahorn.position(entity)
    nx, ny = Int.(get(entity, "nodes", [(0, 0)])[1])

    width = Int(get(entity.data, "width", 8))
    height = Int(get(entity.data, "height", 8))

    return [
        Ahorn.Rectangle(x, y, width, height),
        Ahorn.getSpriteRectangle(hoverTexture, nx, ny, jx=0.5, jy=0.5)
    ]
end

function Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::LuaTalker)
    x, y = Ahorn.position(entity)
    nx, ny = Int.(get(entity, "nodes", [(0, 0)])[1])

    width = Int(get(entity.data, "width", 32))
    height = Int(get(entity.data, "height", 32))

    Ahorn.drawRectangle(ctx, x, y, width, height, (0.0, 1.0, 1.0, 0.4), (0.0, 1.0, 1.0, 1.0))
    Ahorn.drawSprite(ctx, hoverTexture, nx, ny, jx=0.5, jy=0.5)
end

end

max_stars_repo_name: rhelmot/LuaCutscenes

max_issues_repo_name: rhelmot/LuaCutscenes

max_forks_repo_name: rhelmot/LuaCutscenes


text:
"""
File: Helpers.jl
Authors: Reese Pathak, Stephen Boyd

This file defines a set of Helper functions to be used across the
MMA function modules. The implementations take advantage of advanced
Julia features (short-hand functions, macros, expressions, symbols, etc.). If
you're interested in learning about the metaprogramming features available in
Julia, you're encouraged to read more about it at the link below:
http://docs.julialang.org/en/release-0.5/manual/metaprogramming/
"""

module Utils

export @mma_function, demean

demean(x) = x - avg(x)

macro mma_function(ex)
    # get argument vector
    argv = ex.args[1]
    # get x parameter expression (is not evaluated)
    x_param = argv.args[2]
    if length(argv.args) == 2 # single argument funtion, e.g., avg, rms, etc.
      Expr(:function, argv, quote
          @assert !isempty($x_param) "The vector you provided is empty!"
          $(ex.args[2])
      end) |> esc
    else # two-argument functions, e.g., corrcoef or angle
      y_param = argv.args[3]
      Expr(:function, argv, quote
          @assert !isempty($x_param) && !isempty($y_param) "A vector you provided is empty!"
          @assert length($x_param) == length($y_param) "The vector dimensions did not agree!"
          $(ex.args[2])
      end) |> esc
    end
end

end # end module

max_stars_repo_name: reesepathak/mma-functions

max_issues_repo_name: reesepathak/mma-functions

max_forks_repo_name: reesepathak/mma-functions


text:
# BFS demo using Julia stdlib primitives

using sparse, BenchmarkTools

dim = 100

A = sparse_square_matrix(dim, Bool)

bench = @benchmark bfs_primitive(A,1)
display(bench)


text:
using QuadGK
gammaquad(t::Float64) = first(quadgk(x -> x ^ (t - 1) * exp(-x), zero(t), Inf, reltol = 100eps(t)))
@show gammaquad(1.0)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
function create_random_walk!(positions, T)
    for i in 1:length(positions)
        x = 0
        for _ in 1:T
            if rand() < 0.5
                x += 1
            else
                x -= 1
            end
        end
        positions[i] = x
    end
end

function create_random_walk_with_arrays!(positions, T)
    positions .= 0 # Reset to zero
    for _ in 1:T
        positions .+= rand([-1, 1],
         size(positions)...)
    end
end

function walk1D(T)
    x = 0
    for _ in 1:T
        x += rand() < 0.5 ? 1 : -1
    end
    return x
end

function run_walk_in_parallel(N, T)
    positions = zeros(Int, N)
    Threads.@threads for i in 1:N
        positions[i] = walk1D(T)
    end
    nothing
end

max_stars_repo_name: JamieMair/intro-to-julia-talk

max_issues_repo_name: JamieMair/intro-to-julia-talk

max_forks_repo_name: JamieMair/intro-to-julia-talk


text:
export ComplexNormal, Gaussian, Circular, Elliptic, MarchenkoPastur,
        stieltjes
        # re-export
export mean, var, pdf,rand

stieltjes(d::Semicircle) = z->(-z+sqrt(z^2-d.r^2))/2

"""
```julia
MarchenkoPastur <:ContinuousUnivariateDistribution

MarchenkoPastur(λ,σ)
pdf(d::MarchenkoPastur,x::Real)
```
- `λ` : default `0.5`
- `σ` : default `1`

# Examples

Generate a MP rv with ρ = 0.5, σ = 1
```julia
rand(MarchenkoPastur())

1.2656480923753979
```

Compute the desity for the MP distribution with λ=1.6 at the point 0
```julia
pdf(MarchenkoPastur(1.6),0)

0.375
```

Generate 100 MP rvs with ρ=0.1 and σ = 2
```julia
rand(MarchenkoPastur(0.1,2),100)

100-element Vector{Float64}:
 5.31001367107032
 2.3707745658317116
 2.5378523986772343
 2.585256212138476
 ⋮
 3.66799182506567
 6.517865226831382
 3.7628212250212423
 4.582520400697299
```
"""
struct MarchenkoPastur <:ContinuousUnivariateDistribution
    λ::Float64
    σ::Float64
    function MarchenkoPastur(λ=0.5,σ=1) new(λ,σ) end
end

function Distributions.pdf(d::MarchenkoPastur,x::Real)

    if d.λ > 1
        if x < 0 return 0 end
        if x == 0 return 1-1/d.λ end
    end

    λ₋,λ₊ = d.σ^2 * (1 ± sqrt(d.λ)).^2
    if x<λ₋ || x> λ₊
        return 0
    else
        return 1/(2π*d.σ^2)*(sqrt((λ₊-x)*(x-λ₋)))/(d.λ*x)
    end
end

function Base.rand(rng::AbstractRNG, d::MarchenkoPastur)
    if d.λ > 1
        if rand() <= 1-1/d.λ
            return 0
        end
    end

    λ₋,λ₊ = d.σ^2 * ±(1,sqrt(d.λ)).^2
    while true
        x = rand(Uniform(λ₋,λ₊))
        if rand() <= pdf(d,x)
            return x
        end
    end
end

"""
```julia
ComplexNormal(μ,σ)
```
- `μ` : mean, `μ = 0` by default
- `σ` : standard deviation, `σ  =1 ` by default

```julia
# Examples

# Generates 10 iid standard Complex Gaussian r.v.s
rand(ComplexNormal(), 10)

# Generates complex normal with mean 1+1im, variance 4
rand(ComplexNormal(1+1im,2))
```
"""
struct ComplexNormal<:ContinuousUnivariateDistribution
    μ::ComplexF64
    σ::Float64
    function ComplexNormal(μ=0,σ=1)
        @assert( σ >=0 , "σ must be non-negative")
        new(μ,σ)
    end
end

Base.eltype(::Type{ComplexNormal}) = ComplexF64
rand(rng::AbstractRNG, d::ComplexNormal) = randn(rng,ComplexF64)*d.σ+d.μ;

mean(d::ComplexNormal)=d.μ
var(d::ComplexNormal)=d.σ

"""
```julia
Gaussian(beta,μ,σ)
```
- `beta` : 1 (default) for Real Gaussian, 2 for Complex Gaussian
- `μ` : mean, `μ = 0` by default
-  `σ` : standard deviation, `σ  =1 ` by default

```julia
# Generates 10 iid standard normal r.v.s
rand(Gaussian(), 10)

# Generates a complex normal with mean 1+1im, variance 4
rand(Gaussian(2,1+1im,2))
```
"""
struct Gaussian <: ContinuousUnivariateDistribution
    beta::Int8
    μ::Union{ComplexF64,Float64}
    σ::Float64
    function Gaussian(beta=1::Int,μ=0::Int,σ=1::Int)
        if σ >=0
            (beta==1 ? Normal(μ,σ) : ComplexNormal(μ,σ))
        else
            error("σ must be non-negative")
        end
    end
end

mean(d::Gaussian)=d.μ
var(d::Gaussian)=d.σ
"""
```julia
Circular(c=0,R=1)
```
- The uniform distribution on the complex disk with center `c` and radius `R`

```julia
rand(Circular(1,10),200)
```
> Generates 200 uniform r.v.s on the complex disk of center 1 radius 10.
"""
struct  Circular <:ContinuousUnivariateDistribution
    c::ComplexF64
    R::Float64
    function Circular(c=0,R=1) new(c,R) end
end

Base.eltype(::Type{Circular}) = ComplexF64
function rand(rng::AbstractRNG, d::Circular)
    while true
        x_1 = rand(rng,Uniform(-1,1))
        x_2 = rand(rng,Uniform(-1,1))
        if x_1^2+x_2^2<=1
            return x_1*d.R+x_2*d.R*im + d.c
        end
    end
end

mean(d::Circular)=d.c

"""
```julia
Elliptic <:ContinuousUnivariateDistribution

Elliptic(ρ;c,R)
```
- the uniform rv on an ellipse of width `2(1+ρ)R`, height `2(1-ρ)R` centered at `c`
"""
struct  Elliptic <:ContinuousUnivariateDistribution
    ρ::Float64
    c::ComplexF64
    R::Float64

    function Elliptic(ρ=0.5;c=0,R=1)
        @assert(abs(ρ)<1, "ρ must be less than 1 in absolute value")
        new(ρ,c,R)
    end
end

Base.eltype(::Type{Elliptic}) = ComplexF64

function rand(rng::AbstractRNG, d::Elliptic)
    while true
        x_1 = rand(rng,Uniform(-(1+d.ρ),(1+d.ρ)))
        x_2 = rand(rng,Uniform(-(1-d.ρ),(1-d.ρ)))
        if x_1^2/(1+d.ρ)^2+x_2^2/(1-d.ρ)^2<=1
            return x_1*d.R+x_2*d.R*im + d.c
        end
    end
end

mean(d::Elliptic)=d.c


text:
# julia make_environment.jl

using Pkg

rm("./Project.toml"; force=true)
rm("./Manifest.toml"; force=true)

Pkg.activate(".")

Pkg.add([
    PackageSpec(;name="DataStructures",   version="0.18"),
    PackageSpec(;name="DisplayAs",        version="0.1"),
    PackageSpec(;name="Match",            version="1"),
    PackageSpec(;name="NativeSVG",        url="https://github.com/MarkNahabedian/NativeSVG.jl"),
    PackageSpec(;name="PanelCutting",
                url="https://github.com/MarkNahabedian/PanelCutting.jl"),
    PackageSpec(;name="Plots",            version="1.24.3"),  # later version might be broken.
    PackageSpec(;name="Revise",           version="3"),
    PackageSpec(;name="Unitful",          #=version="1"=#),
    # PackageSpec(;name="UnitfulCurrency", version="0.2"),
    PackageSpec(;name="UnitfulUS",        #=version="0.2"=#),
    PackageSpec(;name="VectorLogging",    url="https://github.com/MarkNahabedian/VectorLogging.jl")
    ])


text:
using Distributions

# A covariate expander is a function of one and more covariates and can describe
# its function in a string.
type CovariateExpander
	name::Function # Map the name of a covar to a string describing the expanded covar
	func::Function # Map covar value(s) to expanded values
end

expand(ce::CovariateExpander, data) = ce.func(data)
new_name(ce::CovariateExpander, name) = ce.name(name)

# A DataSet is a float array with named covariates. Each covariate
# corresponds to a row of the data set.
type DataSet
	names::Vector{String}
	data::Array{Float64, 2}
end

num_covars(d::DataSet) = size(d.data, 1)
length(d::DataSet) = size(d.data, 2)

function random_data(covars = 2, n = 30, offset = 10.0)
	names = [string("x", i) for i in 1:covars]
	DataSet(names, offset .+ randn(covars, n))
end

function expand!(d::DataSet, ce::CovariateExpander, indexToExpand)
	d.data = cat(1, d.data, expand(ce, d.data[indexToExpand,:]))
	push!(d.names, new_name(ce, d.names[indexToExpand]))
end
expand!(d::DataSet, ces::Vector{CovariateExpander}, indexToExpand) = map((ce) -> expand!(d, ce, indexToExpand), ces)

LogCE = CovariateExpander((n) -> "log($(n))", (d) -> log(d))
SqrtCE = CovariateExpander((n) -> "sqrt($(n))", (d) -> sqrt(d))
NLogNCE = CovariateExpander((n) -> "$(n)*log($(n))", (d) -> d .* log(d))
Pow2CE = CovariateExpander((n) -> "$(n)^2", (d) -> d.^2)
Pow3CE = CovariateExpander((n) -> "$(n)^3", (d) -> d.^3)

InterceptName = "(Intercept)"
InterceptCE = CovariateExpander((n) -> InterceptName, (d) -> ones(d))

StandardCES = [LogCE, SqrtCE, NLogNCE, Pow2CE, Pow3CE]

function rand_in_range(rows, columns, minRange = 0.0, maxRange = 100.0)
	minRange .+ (maxRange - minRange) .* rand(rows, columns)
end

# Create a random data set sampled from a uniform range.
function random_data(minRange = 0.0, maxRange = 100.0, numCovars = 2, n = 30;
	expandCovars = true
	)

	names = [string("x", i) for i in 1:numCovars]
	d = DataSet([InterceptName, names],
		cat(1, ones(Float64, 1, n), rand_in_range(numCovars, n, minRange, maxRange)))
	if expandCovars
		[expand!(d, StandardCES, i+1) for i in 1:numCovars]
	end
	d

end

# A regression problem is a dataset, a linear function of some subset of that data set and
# an error level (signal-to-noise ratio).
type RegressionProblem
	d::DataSet
	coefficients::Array{Float64,2}
	snr::Float64 # Mean of signal (y) over the standard deviation of the error
end

function random_regression_problem(numCovars = 2, numActiveCovars = 1;
	n = 30,
	minRange = 0.0,
	maxRange = 100.0,
	expandCovars = true,
	snr = 20.0
	)

	d = random_data(minRange, maxRange, numCovars, n; expandCovars = expandCovars)

	coefs = zeros(Float64, 1, num_covars(d))
	[coefs[i] = 10.0*rand() for i in shuffle(collect(1:num_covars(d)))[1:numActiveCovars]]
	p = RegressionProblem(d, coefs, snr)
	y, errors = sample_y(p)
	return p, y, errors

end

# Get the model encoded in a regression problem as a string.
function model(rp::RegressionProblem, coefs = rp.coefficients)
	mstr = ["y"]
	sep() = last(mstr) == "y" ? " = " : " + "
	for i in 1:size(coefs, 2)
		n = rp.d.names[i]
		c = coefs[i]
		if c != 0.0
			if n == InterceptName
				mstr = [mstr, sep(), @sprintf("%.4f", c)]
			else
				mstr = [mstr, sep(), @sprintf("%.4f", c), " * ", n]
			end
		end
	end
	join(mstr, "")
end

import Base.length
length(r::RegressionProblem) = size(r.d.data, 2)

y_actual(r::RegressionProblem) = r.coefficients * r.d.data

function sample_y(r::RegressionProblem)
	y = y_actual(r)
	mu = mean(y)
	sd = mu / r.snr
	errors = rand(Normal(0.0, sd), length(r))'
	return( y .+ errors, errors )
end

# Analyze a regression model proposed for a certain problem. Will calculate
# RSS, AIC, False positives (FP) and False negatives (FN) for the model.
function analyze_regression(rp::RegressionProblem, model::Array{Float64,2})
	report = [
		"M:             $(num_covars(rp.d))",
		"N:             $(length(rp))",
		"Active covars: $(num_nonzero(rp.coefficients))",
		"SNR (target):  $(rp.snr)"
	]
	println(join(report))
end


text:
"""
    mixing_rule_asymetric(op, op_asym, x, p, A, A_asym)

returns an efficient implementation of:
` sum(A[i,j] * x[i] * x[j] * op(p[i],p[j]) * op_asym(x[i],x[j],A_asym[i,j])) for i = 1:n , j = 1:n)`
where `op(p[i],p[j]) == op(p[j],p[i])` , op_asym doesn't follow this symmetry.

"""
function mixing_rule_asymetric(op, op_asym, x, p, A, A_asym)
    N = length(x)
    checkbounds(A, N, N)
    checkbounds(A_asym, N, N)
    @boundscheck checkbounds(p, N)
    @inbounds begin
        res1 = zero(eltype(x))
        for i = 1:N
            x[i] != 0 && begin
                res1 += p[i] * x[i]^2
                for j = 1:i - 1
                    res1 += 2*x[i]*x[j]*op(p[i], p[j])*A[i, j]*op_asym(x[i], x[j], A_asym[i, j])
                end
            end
        end
    end
    return res1
end

_gerg_asymetric_mix_rule(xi, xj, b) = b * (xi + xj) / (xi * b^2 + xj)

ith_index(pv::PackedVofV,i) = @inbounds begin (pv.p[i]):(pv.p[i+1]-1) end

@inline function ith_index(pol_i,exp_i,i)
    @inbounds begin
        kall_1::Int = pol_i[i]
        kall_end::Int = pol_i[i+1]#-1
        kexp_1::Int = exp_i[i]
        kexp_end::Int = exp_i[i+1]#-1
        lall = kall_end - kall_1
        lexp = kexp_end - kexp_1
        divider = lall - lexp
        kmid = kall_1+divider
        k1 = kall_1:(kmid - 1)
        k2 = kmid:(kall_end -1)
        kexp = kexp_1:(kexp_end-1)
        return k1,k2,kexp
    end
end

function _T_scale(model::MultiFluidModel,z=SA[1.],Σz = sum(z))
    Tc = model.properties.Tc.values
    #isone(length(z)) && return only(Tc)
    return mixing_rule_asymetric(
        (a,b)->sqrt(a * b),
        _gerg_asymetric_mix_rule,
        z,
        Tc,
        model.ideal.gamma_T.values,
        model.ideal.beta_T.values,
    )/(Σz*Σz)
end

T_scale(model::MultiFluidModel,z=SA[1.]) = _T_scale(model,z)

p_scale(model::MultiFluidModel,z=SA[1.]) = dot(z,model.properties.pc.values)/sum(z)

function T_scales(model::MultiFluidModel,z=SA[1.])
    return model.properties.Tc.values
end

function _v_scale(model::MultiFluidModel,z=SA[1.],Σz = sum(z))
    vc = model.properties.vc.values
    #isone(length(z)) && return only(vc)
    res = mixing_rule_asymetric(
        (a,b) -> ((cbrt(a) + cbrt(b))*0.5)^3,
        _gerg_asymetric_mix_rule,
        z,
        vc,
        model.ideal.gamma_v.values,
        model.ideal.beta_v.values,
    )
    return res/(Σz*Σz)
end

function lb_volume(model::MultiFluidModel,z=SA[1.])
    return dot(z,model.properties.lb_v.values)
end

function _delta(model::MultiFluidModel, rho, T, z=SA[1.],Σz = sum(z))
    vcmix = _v_scale(model,z,Σz)
    return rho * vcmix
end

function _tau(model::MultiFluidModel, rho, T, z=SA[1.],Σz = sum(z))
    Tcmix  = _T_scale(model,z,Σz)
    return Tcmix / T
end

function _f0(model::MultiFluidModel, ρ, T, z=SA[1.], Σz = sum(z))
    RR = 8.314472 / 8.314510
    #common_type = promote_type(typeof(ρ), typeof(T), eltype(x))
    _0 = zero(ρ + T + first(z))
    lnΣz = log(Σz)
    res = _0
    vc = model.properties.vc.values
    Tc = model.properties.Tc.values
    ϑ = model.ideal.theta.values
    n0 = model.ideal.n0.values
    Tinv = one(T)/T
    @inbounds for i ∈ @comps
        δi = ρ * vc[i]
        τi = Tc[i] * Tinv
        n1,n2,n3,n4,n5,n6,n7 = n0[i]
        ϑ1,ϑ2,ϑ3,ϑ4 = ϑ[i]
        ai = n1 + n2*τi + n3*log(τi)
        iszero(n4) || (ai += n4*log(abs(sinh(ϑ1*τi))))
        iszero(n5) || (ai -= n5*log(cosh(ϑ2*τi)))
        iszero(n6) || (ai += n6*log(abs(sinh(ϑ3*τi))))
        iszero(n7) || (ai -= n7*log(cosh(ϑ4*τi)))
        res += z[i]*(RR*ai + log(δi) + log(z[i]) - lnΣz)
    end
    return res
end

function _fr1(model::MultiFluidModel,δ,τ,z)
    _0 = zero(promote_type(typeof(δ), typeof(τ), eltype(z)))
    res = _0
    n  = model.single.n.values
    c = model.single.c.values
    k_all = n.p
    k_exp = c.p
    nᵢ = n.v
    dᵢ = model.single.d.values.v
    tᵢ = model.single.t.values.v
    cᵢ = c.v
    @inbounds for i ∈ @comps
        ai = _0
        k1,k2,kexp = ith_index(k_all,k_exp,i)
        for k ∈ k1
            ai += nᵢ[k]*(δ^dᵢ[k])*(τ^tᵢ[k])
        end

        for (k,k_) ∈ zip(k2,kexp)
            ai += nᵢ[k]*(δ^dᵢ[k])*(τ^tᵢ[k])*exp(-δ^cᵢ[k_])
        end
        res += z[i]*ai
    end
    return res
end
function _fr2(model::MultiFluidModel,δ,τ,z)
    _0 = zero(promote_type(typeof(δ), typeof(τ), eltype(z)))
    isone(length(z)) && return _0
    F = model.pair.Fij.values
    iszero(nnz(F)) && return _0
    res = _0
    Fij = nonzeros(F)
    n = model.pair.nij.values.storage
    η = model.pair.eta_ij.values.storage
    rows = rowvals(F)
    k_all = n.p
    k_exp = η.p
    nᵢⱼ = n.v
    tᵢⱼ = model.pair.tij.values.storage.v
    dᵢⱼ = model.pair.dij.values.storage.v
    ηᵢⱼ = η.v
    εᵢⱼ = model.pair.epsilon_ij.values.storage.v
    βᵢⱼ = model.pair.beta_ij.values.storage.v
    γᵢⱼ = model.pair.gamma_ij.values.storage.v
    @inbounds for j ∈ @comps
        for ii ∈ nzrange(F, j)
            i = rows[ii]
            Fᵢⱼ= Fij[ii]
            aij = _0
            k1,k2,kexp = ith_index(k_all,k_exp,ii)
            for k ∈ k1
                aij += nᵢⱼ[k]*(δ^(dᵢⱼ[k]))*(τ^(tᵢⱼ[k]))
            end
            for (k,k_) ∈ zip(k2,kexp)
                aij += nᵢⱼ[k]*(δ^(dᵢⱼ[k]))*(τ^(tᵢⱼ[k]))*
                exp(-ηᵢⱼ[k_]*(δ - εᵢⱼ[k_])^2 - βᵢⱼ[k_]*(δ -γᵢⱼ[k_]))
            end
           res +=z[i]*z[j]*Fᵢⱼ*aij
        end
     end
    return res
end

function a_ideal(model::MultiFluidModel, V, T, z=SA[1.0])
    Σz = sum(z)
    ρ = Σz*1.0e-3/V
    _f0(model, ρ, T, z,Σz)/Σz
end

function a_res(model::MultiFluidModel, V, T, z=SA[1.0])
    Σz = sum(z)
    invn = 1/Σz
    invn2 = invn*invn
    ρ = Σz*1.0e-3/V
    δ = _delta(model, ρ, T, z,Σz)
    τ = _tau(model, ρ, T, z,Σz)
    return (_fr1(model,δ,τ,z)*invn + _fr2(model,δ,τ,z)*invn2)
end

function eos(model::MultiFluidModel, V, T, z=SA[1.0])
    Σz = sum(z)
    invn = 1/Σz
    ρ = Σz*1.0e-3/V
    δ = _delta(model, ρ, T, z,Σz)
    τ = _tau(model, ρ, T, z,Σz)
    return  R̄*T*(_f0(model,ρ,T,z,Σz) + _fr1(model,δ,τ,z)+_fr2(model,δ,τ,z)*invn)
end

function eos_res(model::MultiFluidModel, V, T, z=SA[1.0])
    Σz = sum(z)
    invn = 1/Σz
    ρ = Σz*1.0e-3/V
    δ = _delta(model, ρ, T, z,Σz)
    τ = _tau(model, ρ, T, z,Σz)
    return  R̄*T*(_fr1(model,δ,τ,z)+_fr2(model,δ,τ,z)*invn)
end

function x0_sat_pure(model::MultiFluidModel,T)
    Ts = T_scale(model)
    vs = _v_scale(model)*0.001 #remember, vc constants in L/mol
    h = vs*5000.0
    T0 = 369.89*T/Ts
    vl = (1.0/_propaneref_rholsat(T0))*h
    vv = (1.0/_propaneref_rhovsat(T0))*h
    return [log10(vl),log10(vv)]
end

function x0_volume_liquid(model::MultiFluidModel,T,z)
    return 1.01*lb_volume(model,z)
end

molecular_weight(model::MultiFluidModel,z=SA[1.0]) = comp_molecular_weight(mw(model),z)
mw(model::MultiFluidModel) = model.properties.Mw.values


text:
### A Pluto.jl notebook ###
# v0.12.7

using Markdown
using InteractiveUtils

# ╔═╡ 6fd616f6-1ea1-11eb-3814-8bfb4a096c49
using Pkg, DrWatson

# ╔═╡ 6ffe2628-1ea1-11eb-24ea-57f985146a72
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using Distributions, Random, StatsBase, DataFrames, Plots
	Random.seed!(1)
end;

# ╔═╡ ed174bc4-1ea0-11eb-1e2f-a32874cec549
md"## Listing 3.18"

# ╔═╡ 70179b62-1ea1-11eb-13ce-7fb4d54bad4b
begin
	L, K, n  = 500, [450, 400, 250, 100, 50], 30
	hyperDists = [Hypergeometric(k,L-k,n) for k in K]
	xGrid = 0:1:n
	pmfs = [ pdf.(dist, xGrid) for dist in hyperDists ]
	labels = "Successes = " .* string.(K)
end

# ╔═╡ bfe8d6f0-2163-11eb-22e7-9b978bc1d407
begin
	bar( xGrid, pmfs,
		alpha=0.8, c=[:orange :purple :green :red :blue ],
		label=hcat(labels...), ylims=(0,0.25),
		xlabel="x", ylabel="Probability", legend=:top)
end

# ╔═╡ 70182cd8-1ea1-11eb-094d-8d0d49cf15f3
md"## End of listing 3.18"

# ╔═╡ Cell order:
# ╟─ed174bc4-1ea0-11eb-1e2f-a32874cec549
# ╠═6fd616f6-1ea1-11eb-3814-8bfb4a096c49
# ╠═6ffe2628-1ea1-11eb-24ea-57f985146a72
# ╠═70179b62-1ea1-11eb-13ce-7fb4d54bad4b
# ╠═bfe8d6f0-2163-11eb-22e7-9b978bc1d407
# ╟─70182cd8-1ea1-11eb-094d-8d0d49cf15f3


text:
export children, has_children, is_left_child, is_right_child

import AbstractTrees: children, has_children

function children(jlt::AbstractJLBoostTree)
    jlt.children
end

function has_children(jlt::AbstractJLBoostTree)
    length(children(jlt)) > 0
end

function is_left_child(jlt::AbstractJLBoostTree)
    jlt == jlt.parent.children[1]
end

is_right_child(jlt::AbstractJLBoostTree) = !is_left_child(jlt)


text:
using Statistics

function fastPooling!(
    cLayer::OneD,
    Ai::AbstractArray{T,3},
) where {OneD<:Union{MaxPool1D,AveragePool1D},T}

    Aip = padding(cLayer, Ai)
    n_Hi, c, m = size(Aip)
    s_H = cLayer.s
    f_H = cLayer.f
    n_H = (n_Hi - f_H) ÷ s_H + 1

    if cLayer isa MaxPoolLayer
        cLayer.A = reshape(
            maximum(reshape(Aip, f_H, n_H, c, m), dims = 1),
            n_H,
            c,
            m,
        )
    else
        cLayer.A =
            reshape(mean(reshape(Aip, f_H, n_H, c, m), dims = 1), n_H, c, m)
    end #if cLayer isa MaxPoolLayer

    return nothing
end #function fastPooling!(cLayer::OneD

function fastPooling!(
    cLayer::TwoD,
    Ai::AbstractArray{T,4},
) where {TwoD<:Union{MaxPool2D,AveragePool2D},T}

    Aip = padding(cLayer, Ai)
    n_Hi, n_Wi, c, m = size(Aip)
    s_H, s_W = cLayer.s
    f_H, f_W = cLayer.f
    n_H = (n_Hi - f_H) ÷ s_H + 1
    n_W = (n_Wi - f_W) ÷ s_W + 1

    if cLayer isa MaxPoolLayer
        cLayer.A = reshape(
            maximum(
                permutedims(
                    reshape(Aip, f_H, n_H, f_W, n_W, c, m),
                    [1, 3, 2, 4, 5, 6],
                ),
                dims = 1:2,
            ),
            n_H,
            n_W,
            c,
            m,
        )
    else
        cLayer.A = reshape(
            mean(
                permutedims(
                    reshape(Aip, f_H, n_H, f_W, n_W, c, m),
                    [1, 3, 2, 4, 5, 6],
                ),
                dims = 1:2,
            ),
            n_H,
            n_W,
            c,
            m,
        )
    end #if cLayer isa MaxPoolLayer

    return nothing
end #function fastPooling!(cLayer::TwoD,

function fastPooling!(
    cLayer::ThreeD,
    Ai::AbstractArray{T,5},
) where {ThreeD<:Union{MaxPool3D,AveragePool3D},T}

    Aip = padding(cLayer, Ai)
    n_Hi, n_Wi, n_Di, c, m = size(Aip)
    s_H, s_W, s_D = cLayer.s
    f_H, f_W, f_D = cLayer.f

    n_H = (n_Hi - f_H) ÷ s_H + 1
    n_W = (n_Wi - f_W) ÷ s_W + 1
    n_D = (n_Di - f_D) ÷ s_D + 1

    if cLayer isa MaxPoolLayer
        cLayer.A = reshape(
            maximum(
                permutedims(
                    reshape(Aip, f_H, n_H, f_W, n_W, f_D, n_D, c, m),
                    [1, 3, 5, 2, 4, 6, 7, 8],
                ),
                dims = 1:3,
            ),
            n_H,
            n_W,
            n_D,
            c,
            m,
        )
    else
        cLayer.A = reshape(
            mean(
                permutedims(
                    reshape(Aip, f_H, n_H, f_W, n_W, f_D, n_D, c, m),
                    [1, 3, 5, 2, 4, 6, 7, 8],
                ),
                dims = 1:3,
            ),
            n_H,
            n_W,
            n_D,
            c,
            m,
        )
    end #if cLayer isa MaxPoolLayer

    return nothing
end #function fastPooling!(cLayer::ThreeD,

export fastPooling!

### dfastPooling

function dfastPooling!(
    cLayer::OneD,
    Ai::AbstractArray{T,3},
    dAi::AbstractArray{T,3},
    Ao::AbstractArray{T,3},
    dAo::AbstractArray{T,3},
) where {OneD<:Union{MaxPool1D,AveragePool1D},T}

    padS = paddingSize(cLayer, Ai)
    Aip = padding(cLayer, Ai)
    n_Hi, c, m = size(Aip)
    s_H = cLayer.s
    f_H = cLayer.f
    n_H = (n_Hi - f_H) ÷ s_H + 1

    Aipr = reshape(Aip, f_H, n_H, c, m)
    if cLayer isa MaxPoolLayer
        mask = Aipr .== reshape(Ao, 1, n_H, c, m)
    else
        mask = similar(Aipr) .= 1 / prod((cLayer.f))
    end #if cLayer isa MaxPoolLayer
    dAi .= reshape(reshape(dAo, 1, n_H, c, m) .* mask, n_Hi, c, m)[1+padS[1]:end-padS[2], :, :]

    mask = nothing
    return nothing
end #function dfastPooling!(cLayer::OneD

function dfastPooling!(
    cLayer::TwoD,
    Ai::AbstractArray{T,4},
    dAi::AbstractArray{T,4},
    Ao::AbstractArray{T,4},
    dAo::AbstractArray{T,4},
) where {TwoD<:Union{MaxPool2D,AveragePool2D},T}

    padS = paddingSize(cLayer, Ai)
    Aip = padding(cLayer, Ai)
    n_Hi, n_Wi, c, m = size(Aip)
    s_H, s_W = cLayer.s
    f_H, f_W = cLayer.f
    n_H = (n_Hi - f_H) ÷ s_H + 1
    n_W = (n_Wi - f_W) ÷ s_W + 1

    Aipr =
        permutedims(reshape(Aip, f_H, n_H, f_W, n_W, c, m), [1, 3, 2, 4, 5, 6])

    if cLayer isa MaxPoolLayer
        mask = Aipr .== reshape(Ao, 1, 1, n_H, n_W, c, m)
    else
        mask = similar(Aipr) .= 1 / prod((cLayer.f))
    end
    dAi .= reshape(
        permutedims(
            reshape(dAo, 1, 1, n_H, n_W, c, m) .* mask,
            [1, 3, 2, 4, 5, 6],
        ),
        n_Hi,
        n_Wi,
        c,
        m,
    )[1+padS[1]:end-padS[2], 1+padS[3]:end-padS[4], :, :]

    mask = nothing
    return nothing
end #function dfastPooling!(cLayer::TwoD,

function dfastPooling!(
    cLayer::ThreeD,
    Ai::AbstractArray{T,5},
    dAi::AbstractArray{T,5},
    Ao::AbstractArray{T,5},
    dAo::AbstractArray{T,5},
) where {ThreeD<:Union{MaxPool3D,AveragePool3D},T}

    padS = paddingSize(cLayer, Ai)
    Aip = padding(cLayer, Ai)
    n_Hi, n_Wi, n_Di, c, m = size(Aip)
    s_H, s_W, s_D = cLayer.s
    f_H, f_W, f_D = cLayer.f

    n_H = (n_Hi - f_H) ÷ s_H + 1
    n_W = (n_Wi - f_W) ÷ s_W + 1
    n_D = (n_Di - f_D) ÷ s_D + 1

    Aipr = permutedims(
        reshape(Aip, f_H, n_H, f_W, n_W, f_D, n_D, c, m),
        [1, 3, 5, 2, 4, 6, 7, 8],
    )

    if cLayer isa MaxPoolLayer
        mask = Aipr .== reshape(Ao, 1, 1, 1, n_H, n_W, n_D, c, m)
    else
        mask = similar(Aipr) .= 1 / prod((cLayer.f))
    end

    dAi .= reshape(
        permutedims(
            reshape(dAo, 1, 1, 1, n_H, n_W, n_D, c, m) .* mask,
            [1, 4, 2, 5, 3, 6, 7, 8],
        ),
        n_Hi,
        n_Wi,
        n_Di,
        c,
        m,
    )[1+padS[1]:end-padS[2],
        1+padS[3]:end-padS[4],
        1+padS[5]:end-padS[6],
        :,
        :,
    ]

    mask = nothing
    return nothing
end #function dfastPooling!(cLayer::ThreeD,

export dfastPooling!


text:
using MLStyle

function foldast(leaf, branch; kwargs...)
    function go(ast)
        MLStyle.@match ast begin
            Expr(head, args...) => branch(head, map(go, args); kwargs...)
            x                   => leaf(x; kwargs...)
        end
    end

    return go
end

function matchcall(expr)
    @match expr begin
        Expr(:call, f, Expr(:parameters, kwargs...), args...) => (f, args, kwargs)
        Expr(:call, f, args...) => (f,args, [])
        _ => @error "matchcall called on $expr, not a function call"
    end
end

function callify(expr; call=:call)
    leaf(x) = x
    branch(head, newargs) = begin
        expr = Expr(head, newargs...)

        # If it's not a function call, just return it as-is
        head == :call || return expr

        (f, args, kwargs) = matchcall(expr)

        isempty(kwargs) && return Expr(:call, call, f, args...)

        return Expr(:call, call, Expr(:parameters, kwargs...), f, args...)
    end

    foldast(leaf, branch)(expr)
end

macro call(expr)
    callify(expr)
end

# julia> callify(:(f(g(x,y))))
# :(call(f, call(g, x, y)))

# julia> callify(:(f(x; a=3)))
# :(call(f, x; a = 3))

# julia> callify(:(a+b))
# :(call(+, a, b))

# julia> callify(:(call(f,3)))
# :(call(f, 3))

function call(f, args...; kwargs...)
    hasmethod(f, typeof.(args), keys(kwargs)) && return f(args...; kwargs...)

    term(f, args...)
end

# f(x) = x+1

# @call f(2)

# using SymbolicUtils

# @syms x::Vector{Float64} i::Int

# @call getindex(x,i)


text:
@info "Updating packages ..."
using Pkg
Pkg.instantiate()
Pkg.update()
Pkg.precompile()


text:
# Auxilary functions
# Arc functions ────────────────────────────────────────────────────────────────────────────────
"""
    cₐ(a::Arc, x=a.x)

Returns arc cost for arc `a` for arc flow `x`
"""
function cₐ(a::Arc, x=sum(a.xʳ))
    tₒ= a.tₒ
    α = a.α
    β = a.β
    V = a.V
    ϕ = a.ϕ

    t = tₒ * (1 + α * (x/V) ^ β)
    t′= ϕ == 0 || β == 0 ? 0.0 : tₒ * α * β * (x ^ (β - 1))/(V ^ β)

    c = t + ϕ * x * t′
    return c
end

"""
    cₐ′(a::Arc, x=a.x)

Returns first derivative of arc cost wrt arc flow for arc `a` at arc flow `x`
"""
function cₐ′(a::Arc, x=sum(a.xʳ))
    tₒ= a.tₒ
    α = a.α
    β = a.β
    V = a.V
    ϕ = a.ϕ

    t′= β == 0 ? 0.0 : tₒ * α * β * (x ^ (β - 1))/(V ^ β)
    t″ = ϕ == 0 || β == 0 || β == 1 ?  0.0 : tₒ * α * β * (β - 1) * (x ^ (β - 2))/(V ^ β)

    c′ = t′ + ϕ * x * t″
    return c′
end

# Segment functions ────────────────────────────────────────────────────────────────────────────────
"""
    cₑ(e::Vector{Arc})

Returns segment cost for segment `e`
"""
function cₑ(e::Vector{Arc})
    c = 0.0
    for a in e c += a.c end
    return c
end

"""
    cₑ′(e::Vector{Arc})

Returns first derivative of segment cost for segment `e`
"""
function cₑ′(e::Vector{Arc})
    c′ = 0.0
    for a in e c′ += a.c′ end
    return c′
end

"""
    fₑ(e::Vector{Arc}, o::Origin)

Return minimum flow on segment `e` from origin `o`
"""
function fₑ(e::Vector{Arc}, o::Origin)
    f = Inf
    k = o.k
    for a in e if a.xʳ[k] < f f = a.xʳ[k] end end
    return f
end

# Djikstra functions ────────────────────────────────────────────────────────────────────────────────
# NOTE: Node values must match node location (index) in the set of nodes N

"""
    djk(G::Graph, o::Origin)

Djikstra's label setting algorithm.
Returns predecessor arc label (index in vector `A`) `Lᵖ` for least cost paths from origin `o` on graph `G`.
"""
function djk(G::Graph, o::Origin)
    N = G.N
    A = G.A
    K = G.K
    r = o.n

    Lᵖ = Vector{Int64}(undef, length(N))     # Predecessor label
    Lᶜ = [Inf for _ in N]                    # Cost label
    Lᵖ[r.v] = -1
    Lᶜ[r.v] = 0.0

    X = copy(N)                              # Set of open nodes

    index = findfirst(x -> (x == r), X)::Int64
    t = X[index]
    deleteat!(X, index)
    while !isempty(X)
        i = t.v
        for j in t.H
            k = K[i,j]
            a = A[k]
            l = Lᶜ[i] + a.c
            if l < Lᶜ[j] && N[j] ∈ X Lᵖ[j], Lᶜ[j] = k, l
            end
        end
        index = argmin([Lᶜ[n.v] for n in X])
        t = X[index]
        deleteat!(X, index)
    end

    return Lᵖ
end

"""
    tree(G::Graph, o::Origin)

Returns tree for graph `G` rooted at origin `o`
"""
function tree(G::Graph, Lᵖ)
    N = G.N
    A = G.A

    T = Vector{Node}[[] for _ in N]

    for k in Lᵖ
        if k == -1 continue end
        a = A[k]
        t, h = a.t, a.h
        i = t.v
        push!(T[i], h)
    end

    return T
end

"""
    path(G::Graph, Lᵖ, r::Node, s::Node)

Returns path between origin node `r` and destination node `s` on graph `G` using predecessor label `Lᵖ` (for origin node `r`)
"""
function path(G::Graph, Lᵖ, r::Node, s::Node)
    A = G.A

    p = Arc[]

    h = s
    j = h.v
    while j ≠ r.v
        k = Lᵖ[j]
        a = A[k]
        push!(p, a)
        j = a.t.v
    end

    reverse!(p)
    return p
end
path(G::Graph, o::Origin, s::Node) = path(G, o.Lᵖ, o.n, s)  # Shorter version

# ────────────────────────────────────────────────────────────────────────────────

"""
    ispotential(a::Arc, o::Origin, G::Graph)

Identfies if arc `a` on graph `G` is a potential arc wrt flow from origin `o`
"""
function ispotential(a::Arc, o::Origin, G::Graph)
    k = o.k
    xʳₐ = a.xʳ[k]
    cₜₕ = a.c

    t = a.t
    pᵣₜ = path(G, o, t)
    uʳₜ = 0.0
    for a in pᵣₜ uʳₜ += a.c end

    h = a.h
    pᵣₕ = path(G, o, h)
    uʳₕ = 0.0
    for a in pᵣₕ uʳₕ += a.c end

    πʳₐ = uʳₜ + cₜₕ - uʳₕ

    bool = xʳₐ > 1e-12 && πʳₐ > 1e-16

    return bool
end

"""
    𝝳(p::PAS, λ)

Evaluates amount of flow `δ` to shift on pas `p`.
If `δ` is less than the threshold limit of `λ` then `δ` is assumed to be zero.
"""
function 𝝳(p::PAS, λ)
    e₁, e₂, o = p.e₁, p.e₂, p.o

    f₁, f₂ = fₑ(e₁, o), fₑ(e₂, o)
    c₁, c₂ = cₑ(e₁), cₑ(e₂)
    c₁′, c₂′ = cₑ′(e₁), cₑ′(e₂)

    Δ = (c₂ - c₁)/(c₁′ + c₂′)

    if abs(c₂ - c₁) < λ δ = 0.0 end
    if isnan(Δ) δ = 0.0
    elseif Δ ≥ 0 δ = min(Δ, f₂)
    else δ = max(Δ, -f₁) end

    return δ
end

"""
    shift(p::PAS, δ, ϕ)

Shifts flow `δ` on pas `p`.
Argument `ϕ` determines if the arc costs need to be updated for UE or SO assignment.
"""
function shift(p::PAS, δ)
    e₁, e₂, o = p.e₁, p.e₂, p.o
    k = o.k

    for a in e₁
        a.xʳ[k] += δ
        a.c = cₐ(a)
        a.c′= cₐ′(a)
    end

    for a in e₂
        a.xʳ[k] -= δ
        a.c = cₐ(a)
        a.c′= cₐ′(a)
    end

    return
end

"""
    MCS(a::Arc, o::Origin, G::Graph)

Develops pas for arc `a` wrt origin `o` using Maximum Cost Search method
"""
function MCS(a::Arc, o::Origin, G::Graph)
    depth, maxdepth = 1, 2

    i, j = a.t.v, a.h.v
    r, L₁ = o.n, o.Lᵖ
    N, A, K = G.N, G.A, G.K

    pᵣₕ = path(G, L₁, r, a.h)

    s = 1
    p = PAS(Arc[], Arc[], o)

    while depth ≤ maxdepth
        # Intialize
        l = zeros(Int64, length(N))
        for a in pᵣₕ l[a.t.v] = -1 end
        l[i] = 1
        l[j] = 1

        L₂ = Vector{Int64}(undef, length(N))
        L₂[j] = K[i,j]

        # Iterate
        t, h = a.t, a.h
        while true
            h = t

            # Maximum Cost Search
            f = 0.0
            for n in h.T
                k = K[n,h.v]
                x = A[k].xʳ[o.k]
                c = A[k].c
                if x > 1e-12 && c > f
                    f = c
                    t = N[n]
                    L₂[h.v] = k
                end
            end

            # PAS found
            if l[t.v] == -1
                e₁ = path(G, L₁, t, a.h)
                e₂ = path(G, L₂, t, a.h)

                s = l[t.v]
                p = PAS(e₁, e₂, o)

                δ = 𝝳(p, 0.0)
                shift(p, δ)
                bool = ispotential(a, o, G)
                if !bool
                    return s, p
                else
                    depth += 1
                    break
                end
            # Cycle found
            elseif l[t.v] == 1
                if depth < maxdepth
                    pₕₜ = path(G, L₂, h, t)
                    if h != t push!(pₕₜ, A[K[t.v, h.v]]) end
                    δ = Inf
                    k = o.k
                    for a in pₕₜ if a.xʳ[k] ≤ δ δ = a.xʳ[k] end end
                    for a in pₕₜ
                        a.xʳ[k] -= δ
                        a.c = cₐ(a)
                        a.c′= cₐ′(a)
                    end
                end
                depth += 1
                break
            # Continue
            else l[t.v] = 1
            end
        end
    end

    return s, p
end


text:
using Polynomials

d = 3
u = randn(d)
v = randn(d)
r = 3
s = 0.1

pmfg = PMFGaussian(r, s, d)

e  = (dot(u, v) - r)
lb = 0.5(e / s)^2
l  = loglik(pmfg, [u;v])

@test abs(l - lb) <= 1e-10

gub = -e * v / s^2
gvb = -e * u / s^2
g   = gradloglik(pmfg, [u; v])

@test norm(gub - g[1:d]) <= 1e-10 &&
      norm(gvb - g[d+1:end]) <= 1e-10

### test pmf_base corner cases
rex = 0.1
r   = -0.5 + sqrt(0.75)im
p   = poly([1.2 + 1e-9im, -1.7, r, conj(r)]) + rex
@test isapprox(PDSampler.pmf_base(rex, p), 1.2)

# the tests have crude accuracy (because the benchmark is drawn from visual
# inspections with Grapher). This doesn't matter much.

# CASE A # p = (s+1)(s^2+1)
p = Poly([1.0, 1.0]) * Poly([1.0, 0.0, 1.0])

@test abs(PDSampler.pmf_caseA(2.0, p) - 0.9782) <= 1e-3

# CASE B # p = (s-1)(s^2+1)
p = Poly([-1.0, 1.0]) * Poly([1.0, 0.0, 1.0])

@test abs(PDSampler.pmf_caseB(2.0, p, 1.0) - 2.0165) <= 1e-3

# CASE C # p (s+1)(s-1)(s-2)
p = Poly([1.0, 1.0]) * Poly([-1.0, 1.0]) * Poly([-2.0, 1.0])

@test abs(PDSampler.pmf_caseC(0.5014, p, 1.0, 2.0) - 0.2761) <= 1e-3
@test abs(PDSampler.pmf_caseC(1.4090, p, 1.0, 2.0) - 2.3968) <= 1e-3

# CASE D # p = (s-1)(s-2)(s-3)
p = Poly([-1.0, 1.0]) * Poly([-2.0, 1.0]) * Poly([-3.0, 1.0])

@test abs(PDSampler.pmf_caseD(0.1398, p, 1.0, 2.0, 3.0) - 1.4978) <= 1e-3
@test abs(PDSampler.pmf_caseD(0.8405, p, 1.0, 2.0, 3.0) - 3.5928) <= 1e-3


text:
using FactCheck
using Base.Test
using Primes

import OpenCL.cl

import CLFFT
const clfft = CLFFT

#macro throws_pred(ex) FactCheck.throws_pred(ex) end

const TOLERANCE = 1e-3

function allclose{T}(x::AbstractArray{T}, y::AbstractArray{T}; rtol=1e-5, atol=1e-8)
    @assert length(x) == length(y)
    @inbounds begin
        for i in length(x)
            xx, yy = x[i], y[i]
            if !(isapprox(xx, yy; rtol=rtol, atol=atol))
                return false
            end
        end
    end
    return true
end

#Note: these functions have been adapted from the clFFT tests (buffer.h)
function floats_about_equal{T<:clfft.clfftNumber}(a::T, b::T)
    tol  = convert(T, TOLERANCE)
    comp = 1e-5
    if T <: Real
        if abs(a) < comp && abs(b) < comp
            return true
        end
        return abs(a - b) > abs(a * tol) ? false : true
    else
        return (floats_about_equal(real(a), real(b)) &&
                floats_about_equal(imag(b), imag(b)))
    end
end

function allclose_clfft{T<:clfft.clfftNumber}(x::AbstractArray{T}, y::AbstractArray{T})
    @assert length(x) == length(y)
    @inbounds begin
        for i in length(x)
            xx, yy = x[i], y[i]
            if !(floats_about_equal(xx, yy))
                return false
            end
        end
    end
    return true
end

facts("Version") do
    @fact isa(CLFFT.version(), VersionNumber) --> true
    v = clfft.version()
    @fact v.major >= 2 --> true
    @fact v.minor >= 1 --> true
    @fact v.patch >= 0 --> true
end

## NOTE: the following 2 tests FAIL, testing/fixing is highly welcome
##
## facts("Plan") do
##     context("Constructor") do
##         ctx = cl.create_some_context()
##         @fact clfft.Plan(Complex64, ctx, (10, 10)) --> not(nothing)
##         # Plan's throw error on non-muliple 2,3,or 5 dims
##         for x in [2,3,5]
##             @fact (clfft.Plan(Complex64, ctx, (x^3,))) --> not(nothing)
##             @fact_throws (clfft.Plan(Complex64, ctx, (17^3,)))
##             for y in [2,3,5]
##                 @fact (clfft.Plan(Complex64, ctx, (x^3, y^3))) --> not(nothing)
##                 @fact_throws (clfft.Plan(Complex64, ctx, (17^3, y^3)))
##                 for z in [2,3,5]
##                 @fact (clfft.Plan(Complex64, ctx, (x^3, y^3, z^3))) --> not(nothing)
##                 @fact_throws (clfft.Plan(Complex64, ctx, (x^3, 17^3, z^3)))
##                 end
##             end
##         end
##         # FFT only for 1,2 or 3 dim
##         @fact_throws (clfft.Plan(Complex64, ctx, (2^2, 2^2, 2^2, 2^2)))
##         @fact_throws (clfft.Plan(Complex64, ctx, ()))
##         Base.gc()
##     end
## end

## facts("Example FFT Single") do
##     for N in [2^8,]# 3^7, 5^6]
##         @show N
##         X = rand(Complex64, N)
##         fftw_X = fft(X)
##         for device in cl.devices()
##             ctx = cl.Context(device)
##             queue = cl.CmdQueue(ctx)
##             bufX = cl.Buffer(Complex64, ctx, :copy, hostbuf=X)
##             p = clfft.Plan(Complex64, ctx, size(X))
##             clfft.set_layout!(p, :interleaved, :interleaved)
##             clfft.set_result!(p, :inplace)
##             clfft.bake!(p, queue)

##             @fact clfft.context(p) --> ctx
##             @fact clfft.precision(p) --> :single
##             @fact clfft.layout(p) --> (:interleaved, :interleaved)
##             @fact clfft.result(p) --> :inplace
##             @fact clfft.dim(p) --> 1
##             @fact length(clfft.lengths(p)) --> 1
##             @fact clfft.lengths(p)[1] --> length(X)
##             @fact clfft.transpose_result(p) --> false

##             @fact clfft.scaling_factor(p, :forward) --> Float32(1.0)
##             @fact clfft.scaling_factor(p, :backward) --> Float32(1.0 / length(X))
##             @fact clfft.batchsize(p) --> 1

##             clfft.enqueue_transform(p, :forward, [queue], bufX, nothing)
##             # read is blocking (waits on pending event for result)
##             R = cl.read(queue, bufX)
##             @fact allclose(R, fftw_X; rtol=1e-2, atol=1e-3) --> true
##             @fact allclose_clfft(R, fftw_X) --> true
##         end
##     end
## end

facts("Example FFT Double") do
    for N in [2^7,]# 3^6, 5^5]
        @show N
        X = rand(Complex128, N)
        fftw_X = fft(X)
        for device in cl.devices()
            ctx = cl.Context(device)
            queue = cl.CmdQueue(ctx)
            try
                bufX = cl.Buffer(Complex128, ctx, :copy, hostbuf=X)
                p = clfft.Plan(Complex128, ctx, size(X))
                clfft.set_layout!(p, :interleaved, :interleaved)
                clfft.set_result!(p, :inplace)
                clfft.bake!(p, queue)

                @fact clfft.context(p) --> ctx
                @fact clfft.precision(p) --> :double
                @fact clfft.layout(p) --> (:interleaved, :interleaved)
                @fact clfft.result(p) --> :inplace
                @fact clfft.dim(p) --> 1
                @fact length(clfft.lengths(p)) --> 1
                @fact clfft.lengths(p)[1] --> length(X)
                @fact clfft.transpose_result(p) --> false

                @fact clfft.scaling_factor(p, :forward) --> Float32(1.0)
                @fact clfft.scaling_factor(p, :backward) --> Float32(1.0 / length(X))
                @fact clfft.batchsize(p) --> 1

                clfft.enqueue_transform(p, :forward, [queue], bufX, nothing)
                R = cl.read(queue, bufX)
                @fact allclose(R, fftw_X; rtol=1e-2, atol=1e-3) --> true
                @fact allclose_clfft(R, fftw_X) --> true
            catch err
                if isa(err, clfft.CLFFTError)
                    if err.desc == :CLFFT_DEVICE_NO_DOUBLE
                        info("OpenCL.Device $device\ndoes not support double precision")
                    end
                else
                    throw(err)
                end
            end
        end
    end
end

facts("2D FFT Inplace") do
    transform_sizes = [2^6,]#3^4, 5^3]
    for N in transform_sizes
        for M in transform_sizes
            @show (N, M)
            X = rand(Complex64, (N, M))
            fftw_X = fft(X)
            for device in cl.devices()
                ctx = cl.Context(device)
                queue = cl.CmdQueue(ctx)
                p = clfft.Plan(Complex64, ctx, X)
                clfft.bake!(p, queue)
                bufX = cl.Buffer(Complex64, ctx, :copy, hostbuf=X)
                clfft.enqueue_transform(p, :forward, [queue], bufX, nothing)
                R = reshape(cl.read(queue, bufX), size(X))
                @fact allclose(R, fftw_X; rtol=1e-2, atol=1e-3) --> true
                @fact allclose_clfft(R, fftw_X) --> true
                Base.gc()
            end
        end
    end
end

facts("3D FFT Inplace") do
    const N = 64
    X = rand(Complex64, (N, N, N))
    fftw_X = fft(X)
    for device in cl.devices()
        ctx = cl.Context(device)
        queue = cl.CmdQueue(ctx)
        p = clfft.Plan(Complex64, ctx, X)
        clfft.bake!(p, queue)
        bufX = cl.Buffer(Complex64, ctx, :copy, hostbuf=X)
        clfft.enqueue_transform(p, :forward, [queue], bufX, nothing)
        R = reshape(cl.read(queue, bufX), size(X))
        @fact allclose(R, fftw_X; rtol=1e-2, atol=1e-3) --> true
        @fact allclose_clfft(R, fftw_X) --> true
        Base.gc()
    end
end

exitstatus()


text:
# list of industry codes used by Yahoo Finance API
Agricultural_Chemicals=112
Aluminum=132
Chemicals_Major_Diversified=110
Copper=131
Gold=134
Independent_Oil_Gas=121
Industrial_Metals_Minerals=133
Major_Integrated_Oil_Gas=120
Nonmetallic_Mineral_Mining=136
Oil_Gas_Drilling_Exploration=123
Oil_Gas_Equipment_Services=124
Oil_Gas_Pipelines=125
Oil_Gas_Refining_Marketing=122
Silver=135
Specialty_Chemicals=113
Steel_Iron=130
Synthetics=111
Conglomerates=210
Appliances=310
Auto_Manufacturers_Major=330
Auto_Parts=333
Beverages_Brewers=346
Beverages_Soft_Drinks=348
Beverages_Wineries_Distillers=347
Business_Equipment=313
Cigarettes=350
Cleaning_Products=326
Confectioners=345
Dairy_Products=344
Electronic_Equipment=314
Farm_Products=341
Food_Major_Diversified=340
Home_Furnishings_Fixtures=311
Housewares_Accessories=312
Meat_Products=343
Office_Supplies=327
Packaging_Containers=325
Paper_Paper_Products=324
Personal_Products=323
Photographic_Equipment_Supplies=318
Processed_Packaged_Goods=342
Recreational_Goods_Other=317
Recreational_Vehicles=332
Rubber_Plastics=322
Sporting_Goods=316
Textile_Apparel_Clothing=320
Textile_Apparel_Footwear_Accessories=321
Tobacco_Products_Other=351
Toys_Games=315
Trucks_Other_Vehicles=331
Accident_Health_Insurance=431
Asset_Management=422
Closed_End_Fund_Debt=425
Closed_End_Fund_Equity=426
Closed_End_Fund_Foreign=427
Credit_Services=424
Diversified_Investments=423
Foreign_Money_Center_Banks=417
Foreign_Regional_Banks=418
Insurance_Brokers=434
Investment_Brokerage_National=420
Investment_Brokerage_Regional=421
Life_Insurance=430
Money_Center_Banks=410
Mortgage_Investment=447
Property_Casualty_Insurance=432
Property_Management=448
REIT_Diversified=440
REIT_Healthcare_Facilities=442
REIT_Hotel_Motel=443
REIT_Industrial=444
REIT_Office=441
REIT_Residential=445
REIT_Retail=446
Real_Estate_Development=449
Regional_Mid_Atlantic_Banks=412
Regional_Midwest_Banks=414
Regional_Northeast_Banks=411
Regional_Pacific_Banks=416
Regional_Southeast_Banks=413
Regional_Southwest_Banks=415
Savings_Loans=419
Surety_Title_Insurance=433
Biotechnology=515
Diagnostic_Substances=516
Drug_Delivery=513
Drug_Manufacturers_Major=510
Drug_Manufacturers_Other=511
Drug_Related_Products=514
Drugs_Generic=512
Health_Care_Plans=522
Home_Health_Care=526
Hospitals=524
Long_Term_Care_Facilities=523
Medical_Appliances_Equipment=521
Medical_Instruments_Supplies=520
Medical_Laboratories_Research=525
Medical_Practitioners=527
Specialized_Health_Services=528
Aerospace_Defense_Major_Diversified=610
Aerospace_Defense_Products_Services=611
Cement=633
Diversified_Machinery=622
Farm_Construction_Machinery=620
General_Building_Materials=634
General_Contractors=636
Heavy_Construction=635
Industrial_Electrical_Equipment=627
Industrial_Equipment_Components=621
Lumber_Wood_Production=632
Machine_Tools_Accessories=624
Manufactured_Housing=631
Metal_Fabrication=626
Pollution_Treatment_Controls=623
Residential_Construction=630
Small_Tools_Accessories=625
Textile_Industrial=628
Waste_Management=637
Advertising_Agencies=720
Air_Delivery_Freight_Services=773
Air_Services_Other=772
Apparel_Stores=730
Auto_Dealerships=744
Auto_Parts_Stores=738
Auto_Parts_Wholesale=750
Basic_Materials_Wholesale=758
Broadcasting_Radio=724
Broadcasting_TV=723
Building_Materials_Wholesale=751
Business_Services=760
CATV_Systems=725
Catalog_Mail_Order_Houses=739
Computers_Wholesale=755
Consumer_Services=763
Department_Stores=731
Discount_Variety_Stores=732
Drug_Stores=733
Drugs_Wholesale=756
Education_Training_Services=766
Electronics_Stores=735
Electronics_Wholesale=753
Entertainment_Diversified=722
Food_Wholesale=757
Gaming_Activities=714
General_Entertainment=716
Grocery_Stores=734
Home_Furnishing_Stores=737
Home_Improvement_Stores=736
Industrial_Equipment_Wholesale=752
Jewelry_Stores=742
Lodging=710
Major_Airlines=770
Management_Services=769
Marketing_Services=721
Medical_Equipment_Wholesale=754
Movie_Production_Theaters=726
Music_Video_Stores=743
Personal_Services=762
Publishing_Books=729
Publishing_Newspapers=727
Publishing_Periodicals=728
Railroads=776
Regional_Airlines=771
Rental_Leasing_Services=761
Research_Services=768
Resorts_Casinos=711
Restaurants=712
Security_Protection_Services=765
Shipping=775
Specialty_Eateries=713
Specialty_Retail_Other=745
Sporting_Activities=715
Sporting_Goods_Stores=740
Staffing_Outsourcing_Services=764
Technical_Services=767
Toy_Hobby_Stores=741
Trucking=774
Wholesale_Other=759
Application_Software=821
Business_Software_Services=826
Communication_Equipment=841
Computer_Based_Systems=812
Computer_Peripherals=815
Data_Storage_Devices=813
Diversified_Communication_Services=846
Diversified_Computer_Systems=810
Diversified_Electronics=836
Healthcare_Information_Services=825
Information_Delivery_Services=827
Information_Technology_Services=824
Internet_Information_Providers=851
Internet_Service_Providers=850
Internet_Software_Services=852
Long_Distance_Carriers=843
Multimedia_Graphics_Software=820
Networking_Communication_Devices=814
Personal_Computers=811
Printed_Circuit_Boards=835
Processing_Systems_Products=842
Scientific_Technical_Instruments=837
Security_Software_Services=823
Semiconductor_Broad_Line=830
Semiconductor_Integrated_Circuits=833
Semiconductor_Specialized=832
Semiconductor_Equipment_Materials=834
Semiconductor_Memory_Chips=831
Technical_System_Software=822
Telecom_Services_Domestic=844
Telecom_Services_Foreign=845
Wireless_Communications=840
Diversified_Utilities=913
Electric_Utilities=911
Foreign_Utilities=910
Gas_Utilities=912
Water_Utilities=914


text:
# psf-gauss.jl

using SPECTrecon: psf_gauss
using Test: @test, @testset, @test_throws, @inferred

@testset "psf" begin
    psf = @inferred psf_gauss()
    @test psf isa Array{Float32,3}
end


text:
export generate_vortex, generate_patch

# Routine to generate the truth data
function generate_vortex(source, t0, tf, config::VortexConfig, path::String)
    Δt = config.Δt
    source = deepcopy(source)

    x0 = lagrange_to_state(source, config)
    freestream = PotentialFlow.Freestream(config.U)
    ϵinfl = RecipeInflation([config.ϵX; config.ϵΓ])
    cachevels = allocate_velocity(source)

    n0 = ceil(Int64, t0/Δt) + 1
    J = ceil(Int64, (tf-t0)/Δt)
    Acycle = n0:n0+J-1

    Ny = length(config.ss)
    Nx = 3*config.Nv

    xt = zeros(size(x0, 1), J+1)
    x = deepcopy(x0)
    yt = zeros(Ny, J+1)
    tt = zeros(J+1)

    tt[1] = t0
    xt[:,1] .= deepcopy(x0)
    ϵinfl(view(xt,:,1), config)

    yt[:,1] .= pressure(config.ss, source, freestream, t0)
    yt[:,1] .+= config.ϵY*randn(Ny)

    for i=1:length(Acycle)
        reset_velocity!(cachevels, source)
        self_induce_velocity!(cachevels, source, t0 + (i-1)*config.Δt)
        induce_velocity!(cachevels, source, freestream, t0 + (i-1)*config.Δt)

        # Advect the system
        advect!(source, source, cachevels, config.Δt)
        xt[:,i+1] .= deepcopy(lagrange_to_state(source, config))

        # Perturb the state
        ϵinfl(view(xt,:,i+1), config)
        source = deepcopy(state_to_lagrange(xt[:,i+1], config))
        # Compute pressure field
        yt[:,i+1] .= pressure(config.ss, source, freestream, t0+i*config.Δt)
        yt[:,i+1] .+= config.ϵY*randn(Ny)
        tt[i+1] = t0+i*config.Δt
    end

    # Generate the observation at the same time
    data = SyntheticData(tt, config.Δt, x0, xt, yt)

    save(path*"data_test.jld", "tt", tt, "Δt", config.Δt, "x0", x0, "xt", xt, "yt", yt)

    return data
end


text:
struct GoEnv <: GameEnv
  N::Int
  action_space::Int
  planes::Int
  max_action_space::Int
  EMPTY_BOARD::Array{Int8, 2}
  NEIGHBORS::Dict{NTuple{2, Int}, Array{NTuple{2, Int}, 1}}
  DIAGONALS::Dict{NTuple{2, Int}, Array{NTuple{2, Int}, 1}}

  function GoEnv(board_size::Int = 19)
    N = board_size
    action_space = N ^ 2 + 1 # plus 1 for pass move
    ALL_COORDS = [(i, j) for i = 1:N for j = 1:N]
    EMPTY_BOARD = zeros(Int8, N, N)
    planes = 1 * 2 * 8 + 1

    check_bounds(c) = 1 <= c[1] <= N && 1 <= c[2] <= N

    NEIGHBORS = Dict((x, y) => filter(k->check_bounds(k),
                              [(x+1, y), (x-1, y), (x, y+1), (x, y-1)]) for (x, y) in ALL_COORDS)
    DIAGONALS = Dict((x, y) => filter(k->check_bounds(k),
                              [(x+1, y+1), (x+1, y-1), (x-1, y+1), (x-1, y-1)]) for (x, y) in ALL_COORDS)
    new(N, action_space, planes, 361, EMPTY_BOARD, NEIGHBORS, DIAGONALS)
  end
end

include("board.jl")
include("coords.jl")


text:
module TapeTests

using ReverseDiff, Test
using ReverseDiff: SpecialInstruction, ScalarInstruction, NULL_TAPE

include(joinpath(dirname(@__FILE__), "utils.jl"))

for Instr in (SpecialInstruction, ScalarInstruction)
    x, y, k = rand(3), rand(2, 1), rand()
    z = rand()
    c = rand(1)
    instr = Instr(+, (x, y, k), z, c)
    @test isa(instr, Instr{typeof(+)})
    @test instr.func === +
    @test instr.input[1] !== x
    @test instr.input[2] !== y
    @test instr.input[3] === k
    @test instr.input[1] == x
    @test instr.input[2] == y
    @test instr.output === z
    @test instr.cache === c

    tp = InstructionTape()
    ReverseDiff.record!(tp, Instr, +, (x, y, k), z, c)
    @test tp[1] == instr
    @test tp[1].func === +
    @test tp[1].input[1] !== x
    @test tp[1].input[2] !== y
    @test tp[1].input[3] === k
    @test tp[1].input[1] == x
    @test tp[1].input[2] == y
    @test tp[1].output === z
    @test tp[1].cache === c

    ReverseDiff.record!(NULL_TAPE, Instr, +, (x, y, k), z, c)
    @test isempty(NULL_TAPE)
end

end # module


text:
module HartreeFockModule
export computeEnergyHartreeFock, evaluateSCF, evaluateHartreeFock
using ..BaseModule
using ..BasisModule
using ..BasisSetModule
using ..GeometryModule
using ..SpecialMatricesModule
using ..IntegralsModule
using ..ShellModule
using ..LibInt2Module
using ..InitialGuessModule

function computeEnergyHartreeFock(
  density::Matrix,
  kinetic::Matrix,
  nuclearAttraction::Matrix,
  coulomb::Matrix,
  exchange::Matrix,
  interatomicRepulsion::Float64;
  info::Bool=true)

  P = density
  info && println("=================================")
  eT = 2*trace(kinetic*P)
  info && println("Kinetic Energy: $eT")
  eV = 2*trace(nuclearAttraction*P)
  info && println("Nuclear Attraction Energy: $eV")
  info && println("Sum 1-Electron Energy: $(eT+eV)")
  info && println("---------------------------------")
  # G = 2J-K
  eJ = 2*trace(coulomb*P)
  info && println("Coulomb Energy: $eJ")
  eK = -trace(exchange*P)
  info && println("Exchange Energy: $eK")
  info && println("Total Electronic Energy: $(eT+eV+eJ+eK)")
  info && println("Total Energy: $(eT+eV+eJ+eK+interatomicRepulsion)")
  info && println("=================================")

  return eT+eV+eJ+eK
end

function computeEnergyHartreeFock(
  basis::GaussianBasis,
  geometry::Geometry,
  density::Matrix;
  info::Bool=true)

  P = density
  T = computeMatrixKinetic(basis)
  V = computeMatrixNuclearAttraction(basis,geometry)
  J = computeMatrixCoulomb(basis,density)
  K = computeMatrixExchange(basis,density)
  Vnn = computeEnergyInteratomicRepulsion(geometry)
  computeEnergyHartreeFock(P,T,V,J,K,Vnn;info=info)
end

function computeEnergyHartreeFock(
  density::Matrix,
  matrixKinetic::Matrix,
  matrixNuclearAttraction::Matrix,
  electronRepulsionIntegrals::Array{Float64,4},
  energyInteratomicRepulsion::Float64;
  info::Bool=true,
  # inceptions
  electronRepulsionIntegralsExchange::Array{Float64,4}=electronRepulsionIntegrals)

  J = computeMatrixCoulomb(electronRepulsionIntegrals,density)
  K = computeMatrixExchange(electronRepulsionIntegralsExchange,density)
  computeEnergyHartreeFock(density, matrixKinetic, matrixNuclearAttraction, J, K, energyInteratomicRepulsion;info=info)
end

function evaluateSCFStep(
  fock::Matrix,
  overlap::Matrix,
  electronNumber::Integer)

  eig = eigfact(Symmetric(fock),Symmetric(overlap)) # because Symmetric is faster and allows for sorted eigenvalues
  moEnergies = eig.values
  Cocc = eig.vectors[:,1:electronNumber]
  P = Cocc * Cocc.'

  return (moEnergies,P)
end

computeMatrixCoefficients(fock::Matrix, overlap::Matrix) = eigvecs(Symmetric(fock),Symmetric(overlap))

"""
    evaluateSCF(initialGuessDensity::Matrix, overlap::Matrix, kinetic::Matrix, nuclearAttraction::Matrix,
                coulomb::Function, exchange::Function,
                interatomicRepulsion::Float64, electronNumber::Integer)

All evaluateSCF calls take the following optional named arguments:
 - `info::Bool=true`,
 - `detailedinfo::Bool=info`,
 - `energyConvergenceCriterion::Float64 = 1e-8`

This is the fundamental routine to evaluate the Self-consistent field cycle of
Hartree-Fock theory and returns the converged (orbital energies, total energy,
density). This function call requires each element in its most primitive form
and is intended for use in unusual contexts or for performance reasons. For
more usable versions see below.
"""
function evaluateSCF(
  initialGuessDensity::Matrix,
  overlap::Matrix,
  kinetic::Matrix,
  nuclearAttraction::Matrix,
  coulomb::Function,
  exchange::Function,
  interatomicRepulsion::Float64,
  electronNumber::Integer;
  info::Bool=true,
  detailedinfo::Bool=info,
  energyConvergenceCriterion::Float64 = 1e-8)

  P::Matrix = initialGuessDensity
  totalEnergy = computeEnergyHartreeFock(P,kinetic,nuclearAttraction,coulomb(P),exchange(P),interatomicRepulsion;info=false)
  energies = Vector{Float64}()

  energyConverged = false
  i = 0
  J = coulomb(P)
  K = exchange(P)
  while (!energyConverged)
    i+=1
    F = kinetic+nuclearAttraction+2J-K
    energies,P = evaluateSCFStep(F,overlap,electronNumber)
    J = coulomb(P) # By requiring J,K to contain the new density, the variational principle guarantees no energy below the exact energy is obtained
    K = exchange(P) # think: E2 = P II P = tr(P*G[P]) != tr(P*G[P'])
    oldEnergy = totalEnergy
    totalEnergy = computeEnergyHartreeFock(P,kinetic,nuclearAttraction,J,K,interatomicRepulsion,info=detailedinfo)
    info && println("E[$i]: $totalEnergy")
    if (abs(totalEnergy-oldEnergy) < energyConvergenceCriterion)
      info && println("Converged!")
      energyConverged = true
    end
  end
  return energies, totalEnergy, P
end

"""
    evaluateSCF(basis, geometry, initialGuess, electronNumber)

allows to evaluate the Self-consistent field cycle of Hartree-Fock theory in the most
accessible way. Essentially, this function call combines its arguments in such a way
that it can then call upon the more fundamental `evaluateSCF` that takes only primitive
inputs. To this end, the arguments are very flexible:
- `basis`: One of `Basis`, `BasisSet`, Vector{Shell}, Vector{LibInt2Shell}
- `geo`: One of `Geometry`, `String` (the name of a file containing the geometry)
- `initialGuess`: Either a `Matrix` that is used as initial density matrix, or an `InitialGuess` (ZeroGuess, SADGuess)
- `electronNumber`: number of occupied closed shell orbitals (as the other forms are not yet implemented) or `Ionization` of the system

This evaluateSCF method additionally introduces the following optional named arguments:
- `computeTensorElectronRepulsionIntegralsCoulomb`: f: basis/shells/... -> ERIs (used for Coulomb part (and Exchange part if not overridden))
- `computeTensorElectronRepulsionIntegralsExchange`: f: basis/shells/... -> ERIs (used for Exchange part)
"""
function evaluateSCF(
  basis::Union{GaussianBasis,Vector{Shell},Vector{LibInt2Shell},BasisSet},
  geometry::Union{Geometry,String},
  initialGuessDensity::Union{Matrix,InitialGuess},
  electronNumber::Union{Integer,Ionization}=Ionization(0);
  detailedinfo::Bool=true,
  info::Bool=true,
  energyConvergenceCriterion::Float64 = 1e-8,
  computeTensorElectronRepulsionIntegralsCoulomb = computeTensorElectronRepulsionIntegrals,
  computeTensorElectronRepulsionIntegralsExchange = computeTensorElectronRepulsionIntegrals)

  if isa(geometry,String)
    geometry = Geometry(geometry)
  end

  input = basis
  if isa(input,Vector{LibInt2Shell}) && (LibInt2Shell != Shell)
    # keep working with LibInt2Shell
  elseif isa(input,BasisSet) && (LibInt2Shell != Shell)
    basis = computeBasisShellsLibInt2(input,geometry)
  elseif isa(input,BasisSet) && (LibInt2Shell == Shell)
    basis = computeBasis(input,geometry)
  elseif isa(input,GaussianBasis)
    # At the moment we simply use implicit Shell->Basis conversion which avoids ERI recomputation in every step and is therefore faster for the moment
    # later we will want to switch over to working with Shells directly as shell-structure allows for sparsity
    # evaluateSCF(Pinit, S, T, V, P->computeMatrixCoulomb(shells,P), P->computeMatrixExchange(shells,P), Vnn, electronNumber; energyConvergenceCriterion=energyConvergenceCriterion, info=info, detailedinfo=detailedinfo)
    basis = convert(GaussianBasis,input)
  end

  if isa(electronNumber,Ionization)
    electronNumber::Int64 = computeNumberElectrons(geometry,electronNumber.charge)/2
  end

  S = computeMatrixOverlap(basis)
  T = computeMatrixKinetic(basis)
  V = computeMatrixNuclearAttraction(basis,geometry)
  Pinit = InitialGuessModule.computeDensityMatrixGuess(initialGuessDensity,size(S)[1])
  Vnn = computeEnergyInteratomicRepulsion(geometry)

  if (LibInt2Shell != Shell) && isa(basis,Vector{LibInt2Shell})
    Xarg = basis
    Carg = basis
  else
    ERIsCoulomb = computeTensorElectronRepulsionIntegralsCoulomb(basis)
    ERIsExchange = computeTensorElectronRepulsionIntegralsCoulomb==computeTensorElectronRepulsionIntegralsExchange ? ERIsCoulomb : computeTensorElectronRepulsionIntegralsExchange(basis)
    Xarg = ERIsExchange
    Carg = ERIsCoulomb
  end

  evaluateSCF(Pinit, S, T, V, P->computeMatrixCoulomb(Carg,P), P->computeMatrixExchange(Xarg,P), Vnn, electronNumber; energyConvergenceCriterion=energyConvergenceCriterion, info=info, detailedinfo=detailedinfo)
end

evaluateHartreeFock = evaluateSCF

end # module


text:
include("betterment.jl")
include("../metric/patchDifference.jl")

#=
Result of a probe, a patch match for a given target point and corpus point.

!!! Not mutable.  Created as a return value.
=#
struct ProbeResult
    betterment::Betterment
    patchDifference::PatchDifferenceType
end

# Outer constructor for initial/empty/null result
function ProbeResult()
    return ProbeResult(NotBetter, Inf)
end

max_stars_repo_name: bootchk/resynthesizerJulia

max_issues_repo_name: bootchk/resynthesizerJulia

max_forks_repo_name: bootchk/resynthesizerJulia


text:
using NNlib: conv, ∇conv_data, depthwiseconv, output_size

# pad dims of x with dims of y until ndims(x) == ndims(y)
_paddims(x::Tuple, y::Tuple) = (x..., y[(end - (length(y) - length(x) - 1)):end]...)

expand(N, i::Tuple) = i
expand(N, i::Integer) = ntuple(_ -> i, N)

"""
    SamePad()

Passed as an option to convolutional layers (and friends), this causes
the padding to be chosen such that the input and output sizes agree
(on the first `N` dimensions, the kernel or window) when `stride==1`.

See also [`Conv`](@ref), [`MaxPool`](@ref).
"""
struct SamePad end

calc_padding(lt, pad, k::NTuple{N,T}, dilation, stride) where {T,N}= expand(Val(2*N), pad)
function calc_padding(lt, ::SamePad, k::NTuple{N,T}, dilation, stride) where {N,T}
  #Ref: "A guide to convolution arithmetic for deep learning" https://arxiv.org/abs/1603.07285

  # Effective kernel size, including dilation
  k_eff = @. k + (k - 1) * (dilation - 1)
  # How much total padding needs to be applied?
  pad_amt = @. k_eff - 1
  # In case amount of padding is odd we need to apply different amounts to each side.
  return Tuple(mapfoldl(i -> [cld(i, 2), fld(i,2)], vcat, pad_amt))
end

"""
    Conv(filter, in => out, σ = identity;
         stride = 1, pad = 0, dilation = 1, groups = 1, [bias, weight, init])

Standard convolutional layer. `filter` is a tuple of integers
specifying the size of the convolutional kernel;
`in` and `out` specify the number of input and output channels.

Image data should be stored in WHCN order (width, height, channels, batch).
In other words, a 100×100 RGB image would be a `100×100×3×1` array,
and a batch of 50 would be a `100×100×3×50` array.
This has `N = 2` spatial dimensions, and needs a kernel size like `(5,5)`,
a 2-tuple of integers.

To take convolutions along `N` feature dimensions, this layer expects as input an array
with `ndims(x) == N+2`, where `size(x, N+1) == in` is the number of input channels,
and `size(x, ndims(x))` is (as always) the number of observations in a batch.
Then:
* `filter` should be a tuple of `N` integers.
* Keywords `stride` and `dilation` should each be either single integer,
  or a tuple with `N` integers.
* Keyword `pad` can be:
  - a single integer for equal padding all around,
  - a tuple of `N` integers, to apply the same padding at begin/end of each spatial dimension,
  - a tuple of `2*N` integers, for asymmetric padding, or
  - the singleton `SamePad()`, to calculate padding such that
    `size(output,d) == size(x,d) / stride` (possibly rounded) for each spatial dimension.
* Keyword `groups` is expected to be an `Int`. It specifies the number of groups
  to divide a convolution into.

Keywords to control initialization of the layer:
* `init` - Function used to generate initial weights. Defaults to `glorot_uniform`.
* `weight` - Initial weights of the layer. Typically an array, and can be used to override
  other configurations. By default, these are generated using [`convfilter`](@ref).
* `bias` - Initial bias is zero by default, this can be disabled entirely by setting it to
  [`Flux.Zeros()`](@ref) or equivalently `false`, or another vector provided as
  `bias = randn(Float32, out)`.

See also [`ConvTranspose`](@ref), [`DepthwiseConv`](@ref), [`CrossCor`](@ref).

# Examples
```jldoctest
julia> xs = rand(Float32, 100, 100, 3, 50); # a batch of images

julia> layer = Conv((5,5), 3 => 7, relu; bias = false)
Conv((5, 5), 3 => 7, relu, bias=false)  # 525 parameters

julia> layer(xs) |> size
(96, 96, 7, 50)

julia> Conv((5,5), 3 => 7; stride = 2)(xs) |> size
(48, 48, 7, 50)

julia> Conv((5,5), 3 => 7; stride = 2, pad = SamePad())(xs) |> size
(50, 50, 7, 50)

julia> Conv((1,1), 3 => 7; pad = (20,10,0,0))(xs) |> size
(130, 100, 7, 50)

julia> Conv((5,5), 3 => 7; stride = 2, dilation = 4)(xs) |> size
(42, 42, 7, 50)
```
"""
struct Conv{N,M,F,A,V}
  σ::F
  weight::A
  bias::V
  stride::NTuple{N,Int}
  pad::NTuple{M,Int}
  dilation::NTuple{N,Int}
  groups::Int
end

"""
    Conv(weight::AbstractArray, [bias, activation; stride, pad, dilation])

Constructs a convolutional layer with the given weight and bias.
Accepts the same keywords (and has the same defaults) as the `Conv((4,4), 3 => 7, relu)`
method.

# Examples
```jldoctest
julia> weight = rand(3, 4, 5);

julia> bias = zeros(5);

julia> c1 = Conv(weight, bias, sigmoid)  # expects 1 spatial dimension
Conv((3,), 4 => 5, σ)  # 65 parameters

julia> c1(randn(100, 4, 64)) |> size
(98, 5, 64)

julia> params(c1) |> length
2
```
"""
function Conv(w::AbstractArray{T,N}, b = true, σ = identity;
              stride = 1, pad = 0, dilation = 1, groups = 1) where {T,N}
  stride = expand(Val(N-2), stride)
  dilation = expand(Val(N-2), dilation)
  pad = calc_padding(Conv, pad, size(w)[1:N-2], dilation, stride)
  bias = create_bias(w, b, size(w, N))
  return Conv(σ, w, bias, stride, pad, dilation, groups)
end

function Conv(k::NTuple{N,Integer}, ch::Pair{<:Integer,<:Integer}, σ = identity;
            init = glorot_uniform, stride = 1, pad = 0, dilation = 1, groups = 1,
            weight = convfilter(k, (ch[1] ÷ groups => ch[2]); init), bias = true) where N

  Conv(weight, bias, σ; stride, pad, dilation, groups)
end

"""
    convfilter(filter::Tuple, in => out[; init = glorot_uniform])

Constructs a standard convolutional weight matrix with given `filter` and
channels from `in` to `out`.

Accepts the keyword `init` (default: `glorot_uniform`) to control the sampling
distribution.

See also: [`depthwiseconvfilter`](@ref)
"""
convfilter(filter::NTuple{N,Integer}, ch::Pair{<:Integer,<:Integer};
          init = glorot_uniform) where N = init(filter..., ch...)

@functor Conv

function (c::Conv)(x::AbstractArray)
  σ, b = c.σ, reshape(c.bias, ntuple(_ -> 1, length(c.stride))..., :, 1)
  cdims = DenseConvDims(x, c.weight; stride = c.stride, padding = c.pad, dilation = c.dilation, groups = c.groups)
  σ.(conv(x, c.weight, cdims) .+ b)
end

function Base.show(io::IO, l::Conv)
  print(io, "Conv(", size(l.weight)[1:ndims(l.weight)-2])
  print(io, ", ", size(l.weight, ndims(l.weight)-1), " => ", size(l.weight, ndims(l.weight)))
  _print_conv_opt(io, l)
  print(io, ")")
end

function _print_conv_opt(io::IO, l)
  l.σ == identity || print(io, ", ", l.σ)
  all(==(0), l.pad) || print(io, ", pad=", _maybetuple_string(l.pad))
  all(==(1), l.stride) || print(io, ", stride=", _maybetuple_string(l.stride))
  all(==(1), l.dilation) || print(io, ", dilation=", _maybetuple_string(l.dilation))
  l.bias == Zeros() && print(io, ", bias=false")
end

"""
    ConvTranspose(filter, in => out, σ=identity; stride=1, pad=0, dilation=1, [bias, init])

Standard convolutional transpose layer. `filter` is a tuple of integers
specifying the size of the convolutional kernel, while
`in` and `out` specify the number of input and output channels.

Note that `pad=SamePad()` here tries to ensure `size(output,d) == size(x,d) * stride`.

Parameters are controlled by additional keywords, with defaults
`init=glorot_uniform` and `bias=true`.

See also [`Conv`](@ref) for more detailed description of keywords.

# Examples
```jldoctest
julia> xs = rand(Float32, 100, 100, 3, 50);  # a batch of 50 RGB images

julia> lay = ConvTranspose((5,5), 3 => 7, relu)
ConvTranspose((5, 5), 3 => 7, relu)  # 532 parameters

julia> lay(xs) |> size
(104, 104, 7, 50)

julia> ConvTranspose((5,5), 3 => 7, stride=2)(xs) |> size
(203, 203, 7, 50)

julia> ConvTranspose((5,5), 3 => 7, stride=3, pad=SamePad())(xs) |> size
(300, 300, 7, 50)
```
"""
struct ConvTranspose{N,M,F,A,V}
  σ::F
  weight::A
  bias::V
  stride::NTuple{N,Int}
  pad::NTuple{M,Int}
  dilation::NTuple{N,Int}
end

"""
    ConvTranspose(weight::AbstractArray, [bias, activation; stride, pad, dilation])

Constructs a layer with the given weight and bias arrays.
Accepts the same keywords as the `ConvTranspose((4,4), 3 => 7, relu)` method.
"""
function ConvTranspose(w::AbstractArray{T,N}, bias = true, σ = identity;
                      stride = 1, pad = 0, dilation = 1) where {T,N}
  stride = expand(Val(N-2), stride)
  dilation = expand(Val(N-2), dilation)
  pad = calc_padding(ConvTranspose, pad, size(w)[1:N-2], dilation, stride)
  b = create_bias(w, bias, size(w, N-1))
  return ConvTranspose(σ, w, b, stride, pad, dilation)
end

function ConvTranspose(k::NTuple{N,Integer}, ch::Pair{<:Integer,<:Integer}, σ = identity;
                      init = glorot_uniform, stride = 1, pad = 0, dilation = 1,
                      weight = convfilter(k, reverse(ch), init = init), bias = true) where N

  ConvTranspose(weight, bias, σ; stride, pad, dilation)
end

@functor ConvTranspose

function conv_transpose_dims(c::ConvTranspose, x::AbstractArray)
    # Calculate size of "input", from ∇conv_data()'s perspective...
    combined_pad = (c.pad[1:2:end] .+ c.pad[2:2:end])
    I = (size(x)[1:end-2] .- 1).*c.stride .+ 1 .+ (size(c.weight)[1:end-2] .- 1).*c.dilation .- combined_pad
    C_in = size(c.weight)[end-1]
    batch_size = size(x)[end]
    # Create DenseConvDims() that looks like the corresponding conv()
    return DenseConvDims((I..., C_in, batch_size), size(c.weight);
                        stride=c.stride,
                        padding=c.pad,
                        dilation=c.dilation,
    )
end

# TODO: Find proper fix for https://github.com/FluxML/Flux.jl/issues/900
@nograd conv_transpose_dims

function (c::ConvTranspose)(x::AbstractArray)
  σ, b = c.σ, reshape(c.bias, map(_->1, c.stride)..., :, 1)
  cdims = conv_transpose_dims(c, x)
  σ.(∇conv_data(x, c.weight, cdims) .+ b)
end

function Base.show(io::IO, l::ConvTranspose)
  print(io, "ConvTranspose(", size(l.weight)[1:ndims(l.weight)-2])
  print(io, ", ", size(l.weight, ndims(l.weight)), " => ", size(l.weight, ndims(l.weight)-1))
  _print_conv_opt(io, l)
  print(io, ")")
end

function calc_padding(::Type{ConvTranspose}, pad::SamePad, k::NTuple{N,T}, dilation, stride) where {N,T}
  calc_padding(Conv, pad, k .- stride .+ 1, dilation, stride)
end

"""
    DepthwiseConv(filter, in => out, σ=identity; stride=1, pad=0, dilation=1, [bias, init])

Depthwise convolutional layer. `filter` is a tuple of integers
specifying the size of the convolutional kernel, while
`in` and `out` specify the number of input and output channels.

Note that `out` must be an integer multiple of `in`.

Parameters are controlled by additional keywords, with defaults
`init=glorot_uniform` and `bias=true`.

See also [`Conv`](@ref) for more detailed description of keywords.

# Examples
```jldoctest
julia> xs = rand(Float32, 100, 100, 3, 50);  # a batch of 50 RGB images

julia> lay = DepthwiseConv((5,5), 3 => 6, relu; bias=false)
DepthwiseConv((5, 5), 3 => 6, relu, bias=false)  # 150 parameters

julia> lay(xs) |> size
(96, 96, 6, 50)

julia> DepthwiseConv((5,5), 3 => 9, stride=2, pad=2)(xs) |> size
(50, 50, 9, 50)
```
"""
struct DepthwiseConv{N,M,F,A,V}
  σ::F
  weight::A
  bias::V
  stride::NTuple{N,Int}
  pad::NTuple{M,Int}
  dilation::NTuple{N,Int}
end

"""
    DepthwiseConv(weight::AbstractArray, bias, [activation; stride, pad, dilation])

Constructs a layer with the given weight and bias arrays.
Accepts the same keywords as the `DepthwiseConv((4,4), 3 => 6, relu)` method.
"""
function DepthwiseConv(w::AbstractArray{T,N}, bias = true, σ = identity;
                      stride = 1, pad = 0, dilation = 1) where {T,N}
  stride = expand(Val(N-2), stride)
  dilation = expand(Val(N-2), dilation)
  pad = calc_padding(DepthwiseConv, pad, size(w)[1:N-2], dilation, stride)
  b = create_bias(w, bias, prod(size(w)[N-1:end]))
  return DepthwiseConv(σ, w, b, stride, pad, dilation)
end

function DepthwiseConv(k::NTuple{N,Integer}, ch::Pair{<:Integer,<:Integer}, σ = identity;
                init = glorot_uniform, stride = 1, pad = 0, dilation = 1,
                weight = depthwiseconvfilter(k, ch, init = init), bias = true) where N
  @assert ch[2] % ch[1] == 0 "Output channels must be integer multiple of input channels"
  return DepthwiseConv(weight, bias, σ; stride, pad, dilation)
end

@functor DepthwiseConv

"""
    depthwiseconvfilter(filter::Tuple, in => out)

Constructs a depthwise convolutional weight array defined by `filter` and channels
from `in` to `out`.

Accepts the keyword `init` (default: `glorot_uniform`) to control the sampling
distribution.

See also: [`convfilter`](@ref)
"""
depthwiseconvfilter(filter::NTuple{N,Integer}, ch::Pair{<:Integer,<:Integer};
                    init = glorot_uniform) where N = init(filter..., div(ch[2], ch[1]), ch[1])

function (c::DepthwiseConv)(x)
  σ, b = c.σ, reshape(c.bias, map(_->1, c.stride)..., :, 1)
  cdims = DepthwiseConvDims(x, c.weight; stride=c.stride, padding=c.pad, dilation=c.dilation)
  σ.(depthwiseconv(x, c.weight, cdims) .+ b)
end

function Base.show(io::IO, l::DepthwiseConv)
  print(io, "DepthwiseConv(", size(l.weight)[1:end-2])
  print(io, ", ", size(l.weight)[end], " => ", prod(size(l.weight)[end-1:end]))
  _print_conv_opt(io, l)
  print(io, ")")
end

"""
    CrossCor(filter, in => out, σ=identity; stride=1, pad=0, dilation=1, [bias, init])

Standard cross convolutional layer. `filter` is a tuple of integers
specifying the size of the convolutional kernel;
`in` and `out` specify the number of input and output channels.

Parameters are controlled by additional keywords, with defaults
`init=glorot_uniform` and `bias=true`.

See also [`Conv`](@ref) for more detailed description of keywords.

# Examples

```jldoctest
julia> xs = rand(Float32, 100, 100, 3, 50);  # a batch of 50 RGB images

julia> lay = CrossCor((5,5), 3 => 6, relu; bias=false)
CrossCor((5, 5), 3 => 6, relu, bias=false)  # 450 parameters

julia> lay(xs) |> size
(96, 96, 6, 50)

julia> CrossCor((5,5), 3 => 7, stride=3, pad=(2,0))(xs) |> size
(34, 32, 7, 50)
```
"""
struct CrossCor{N,M,F,A,V}
  σ::F
  weight::A
  bias::V
  stride::NTuple{N,Int}
  pad::NTuple{M,Int}
  dilation::NTuple{N,Int}
end

"""
    CrossCor(weight::AbstractArray, [bias, activation; stride, pad, dilation])

Constructs a layer with the given weight and bias arrays.
Accepts the same keywords as the `CrossCor((4,4), 3 => 7, relu)` method.
"""
function CrossCor(w::AbstractArray{T,N}, bias = true, σ = identity;
                  stride = 1, pad = 0, dilation = 1) where {T,N}
  stride = expand(Val(N-2), stride)
  dilation = expand(Val(N-2), dilation)
  pad = calc_padding(CrossCor, pad, size(w)[1:N-2], dilation, stride)
  b = create_bias(w, bias, size(w, N))
  return CrossCor(σ, w, b, stride, pad, dilation)
end

function CrossCor(k::NTuple{N,Integer}, ch::Pair{<:Integer,<:Integer}, σ = identity;
                  init = glorot_uniform, stride = 1, pad = 0, dilation = 1,
                  weight = convfilter(k, ch, init = init), bias = true) where N

  return CrossCor(weight, bias, σ; stride, pad, dilation)
end

@functor CrossCor

function crosscor(x, w, ddims::DenseConvDims)
  ddims = DenseConvDims(ddims, F=true)
  return conv(x, w, ddims)
end

function (c::CrossCor)(x::AbstractArray)
  σ, b = c.σ, reshape(c.bias, map(_->1, c.stride)..., :, 1)
  cdims = DenseConvDims(x, c.weight; stride=c.stride, padding=c.pad, dilation=c.dilation)
  σ.(crosscor(x, c.weight, cdims) .+ b)
end

function Base.show(io::IO, l::CrossCor)
  print(io, "CrossCor(", size(l.weight)[1:ndims(l.weight)-2])
  print(io, ", ", size(l.weight, ndims(l.weight)-1), " => ", size(l.weight, ndims(l.weight)))
  _print_conv_opt(io, l)
  print(io, ")")
end

"""
    AdaptiveMaxPool(out::NTuple)

Adaptive max pooling layer. Calculates the necessary window size
such that its output has `size(y)[1:N] == out`.

Expects as input an array with `ndims(x) == N+2`, i.e. channel and
batch dimensions, after the `N` feature dimensions, where `N = length(out)`.

See also [`MaxPool`](@ref), [`AdaptiveMeanPool`](@ref).

# Examples
```jldoctest
julia> xs = rand(Float32, 100, 100, 3, 50);  # batch of 50 RGB images

julia> AdaptiveMaxPool((25, 25))(xs) |> size
(25, 25, 3, 50)

julia> MaxPool((4,4))(xs) ≈ AdaptiveMaxPool((25, 25))(xs)
true
```
"""
struct AdaptiveMaxPool{S, O}
  out::NTuple{O, Int}
  AdaptiveMaxPool(out::NTuple{O, Int}) where O = new{O + 2, O}(out)
end

function (a::AdaptiveMaxPool{S})(x::AbstractArray{T, S}) where {S, T}
  insize = size(x)[1:end-2]
  outsize = a.out
  stride = insize .÷ outsize
  k = insize .- (outsize .- 1) .* stride
  pad = 0
  pdims = PoolDims(x, k; padding=pad, stride=stride)
  return maxpool(x, pdims)
end

function Base.show(io::IO, a::AdaptiveMaxPool)
  print(io, "AdaptiveMaxPool(", a.out, ")")
end

"""
    AdaptiveMeanPool(out::NTuple)

Adaptive mean pooling layer. Calculates the necessary window size
such that its output has `size(y)[1:N] == out`.

Expects as input an array with `ndims(x) == N+2`, i.e. channel and
batch dimensions, after the `N` feature dimensions, where `N = length(out)`.

See also [`MaxPool`](@ref), [`AdaptiveMaxPool`](@ref).

# Examples
```jldoctest
julia> xs = rand(Float32, 100, 100, 3, 50);  # batch of 50 RGB images

julia> AdaptiveMeanPool((25, 25))(xs) |> size
(25, 25, 3, 50)

julia> MeanPool((4,4))(xs) ≈ AdaptiveMeanPool((25, 25))(xs)
true
```
"""
struct AdaptiveMeanPool{S, O}
  out::NTuple{O, Int}
  AdaptiveMeanPool(out::NTuple{O, Int}) where O = new{O + 2, O}(out)
end

function (a::AdaptiveMeanPool{S})(x::AbstractArray{T, S}) where {S, T}
  insize = size(x)[1:end-2]
  outsize = a.out
  stride = insize .÷ outsize
  k = insize .- (outsize .- 1) .* stride
  pad = 0
  pdims = PoolDims(x, k; padding=pad, stride=stride)
  return meanpool(x, pdims)
end

function Base.show(io::IO, a::AdaptiveMeanPool)
  print(io, "AdaptiveMeanPool(", a.out, ")")
end

"""
    GlobalMaxPool()

Global max pooling layer.

Transforms (w,h,c,b)-shaped input into (1,1,c,b)-shaped output,
by performing max pooling on the complete (w,h)-shaped feature maps.

See also [`MaxPool`](@ref), [`GlobalMeanPool`](@ref).

```jldoctest
julia> xs = rand(Float32, 100, 100, 3, 50);

julia> m = Chain(Conv((3,3), 3 => 7), GlobalMaxPool());

julia> m(xs) |> size
(1, 1, 7, 50)

julia> GlobalMaxPool()(rand(3,5,7)) |> size  # preserves 2 dimensions
(1, 5, 7)
```
"""
struct GlobalMaxPool end

function (g::GlobalMaxPool)(x)
  # Input size
  x_size = size(x)
  # Kernel size
  k = x_size[1:end-2]
  # Pooling dimensions
  pdims = PoolDims(x, k)

  return maxpool(x, pdims)
end

function Base.show(io::IO, g::GlobalMaxPool)
  print(io, "GlobalMaxPool()")
end

"""
    GlobalMeanPool()

Global mean pooling layer.

Transforms (w,h,c,b)-shaped input into (1,1,c,b)-shaped output,
by performing mean pooling on the complete (w,h)-shaped feature maps.

```jldoctest
julia> xs = rand(Float32, 100, 100, 3, 50);

julia> m = Chain(Conv((3,3), 3 => 7), GlobalMeanPool());

julia> m(xs) |> size
(1, 1, 7, 50)
```
"""
struct GlobalMeanPool end

function (g::GlobalMeanPool)(x)
  # Input size
  x_size = size(x)
  # Kernel size
  k = x_size[1:end-2]
  # Pooling dimensions
  pdims = PoolDims(x, k)

  return meanpool(x, pdims)
end

function Base.show(io::IO, g::GlobalMeanPool)
  print(io, "GlobalMeanPool()")
end

"""
    MaxPool(window::NTuple; pad=0, stride=window)

Max pooling layer, which replaces all pixels in a block of
size `window` with one.

Expects as input an array with `ndims(x) == N+2`, i.e. channel and
batch dimensions, after the `N` feature dimensions, where `N = length(window)`.

By default the window size is also the stride in each dimension.
The keyword `pad` accepts the same options as for the `Conv` layer,
including `SamePad()`.

See also [`Conv`](@ref), [`MeanPool`](@ref), [`AdaptiveMaxPool`](@ref), [`GlobalMaxPool`](@ref).

# Examples

```jldoctest
julia> xs = rand(Float32, 100, 100, 3, 50);  # batch of 50 RGB images

julia> m = Chain(Conv((5, 5), 3 => 7, pad=SamePad()), MaxPool((5, 5), pad=SamePad()))
Chain(
  Conv((5, 5), 3 => 7, pad=2),          # 532 parameters
  MaxPool((5, 5), pad=2),
)

julia> m[1](xs) |> size
(100, 100, 7, 50)

julia> m(xs) |> size
(20, 20, 7, 50)

julia> lay = MaxPool((5,), pad=2, stride=(3,))  # one-dimensional window
MaxPool((5,), pad=2, stride=3)

julia> lay(rand(Float32, 100, 7, 50)) |> size
(34, 7, 50)
```
"""
struct MaxPool{N,M}
  k::NTuple{N,Int}
  pad::NTuple{M,Int}
  stride::NTuple{N,Int}
end

function MaxPool(k::NTuple{N,Integer}; pad = 0, stride = k) where N
  stride = expand(Val(N), stride)
  pad = calc_padding(MaxPool ,pad, k, 1, stride)
  return MaxPool(k, pad, stride)
end

function (m::MaxPool)(x)
  pdims = PoolDims(x, m.k; padding=m.pad, stride=m.stride)
  return maxpool(x, pdims)
end

function Base.show(io::IO, m::MaxPool)
  print(io, "MaxPool(", m.k)
  all(==(0), m.pad) || print(io, ", pad=", _maybetuple_string(m.pad))
  m.stride == m.k || print(io, ", stride=", _maybetuple_string(m.stride))
  print(io, ")")
end

_maybetuple_string(pad) = string(pad)
_maybetuple_string(pad::Tuple) = all(==(pad[1]), pad) ? string(pad[1])  : string(pad)

"""
    MeanPool(window::NTuple; pad=0, stride=window)

Mean pooling layer, averaging all pixels in a block of size `window`.

Expects as input an array with `ndims(x) == N+2`, i.e. channel and
batch dimensions, after the `N` feature dimensions, where `N = length(window)`.

By default the window size is also the stride in each dimension.
The keyword `pad` accepts the same options as for the `Conv` layer,
including `SamePad()`.

See also [`Conv`](@ref), [`MaxPool`](@ref), [`AdaptiveMeanPool`](@ref).

# Examples

```jldoctest
julia> xs = rand(Float32, 100, 100, 3, 50);

julia> m = Chain(Conv((5,5), 3 => 7), MeanPool((5,5), pad=SamePad()))
Chain(
  Conv((5, 5), 3 => 7),                 # 532 parameters
  MeanPool((5, 5), pad=2),
)

julia> m[1](xs) |> size
(96, 96, 7, 50)

julia> m(xs) |> size
(20, 20, 7, 50)
```
"""
struct MeanPool{N,M}
  k::NTuple{N,Int}
  pad::NTuple{M,Int}
  stride::NTuple{N,Int}
end

function MeanPool(k::NTuple{N,Integer}; pad = 0, stride = k) where N
  stride = expand(Val(N), stride)
  pad = calc_padding(MeanPool, pad, k, 1, stride)
  return MeanPool(k, pad, stride)
end

function (m::MeanPool)(x)
  pdims = PoolDims(x, m.k; padding=m.pad, stride=m.stride)
  return meanpool(x, pdims)
end

function Base.show(io::IO, m::MeanPool)
  print(io, "MeanPool(", m.k)
  all(==(0), m.pad) || print(io, ", pad=", _maybetuple_string(m.pad))
  m.stride == m.k || print(io, ", stride=", _maybetuple_string(m.stride))
  print(io, ")")
end