text:
using Diagonalizations, LinearAlgebra, PosDefManifold, Test

##  Create data for testing the case k>1, m>1 ##
# `t` is the number of samples,
# `m` is the number of datasets,
# `k` is the number of observations,
# `n` is the number of variables,
# `noise` must be smaller than 1.0. The smaller the noise, the more data are correlated
# Output k vectors of m data data matrices
function getData(t, m, k, n, noise)
    # create m identical data matrices and rotate them by different
    # random orthogonal matrices V_1,...,V_m
    𝐕=[randU(n) for i=1:m] # random orthogonal matrices
    # variables common to all subjects with unique variance profile across k
    X=[(abs2.(randn(n))).*randn(n, t) for s=1:k]
    # each subject has this common part plus a random part
    𝐗=[[𝐕[i]*((1-noise)*X[s] + noise*randn(n, t)) for i=1:m] for s=1:k]
    return 𝐗, 𝐕
end

function getData(::Type{Complex{T}}, t, m, k, n, noise) where {T<:AbstractFloat}
    # create m identical data matrices and rotate them by different
    # random orthogonal matrices V_1,...,V_m
    𝐕=[randU(ComplexF64, n) for i=1:m] # random orthogonal matrices
    # variables common to all subjects with unique variance profile across k
    X=[(abs2.(randn(n))).*randn(ComplexF64, n, t) for s=1:k]
    # each subject has this common part plus a random part
    𝐗=[[𝐕[i]*((1-noise)*X[s] + noise*randn(ComplexF64, n, t)) for i=1:m] for s=1:k]
    return 𝐗, 𝐕
end

# REAL data
# do joint blind source separation of non-stationary data
t, m, n, k, noise = 200, 5, 4, 6, 0.1
Xset, Vset=getData(t, m, k, n, noise)
𝒞=Array{Matrix}(undef, k, m, m)
for s=1:k, i=1:m, j=1:m 𝒞[s, i, j]=(Xset[s][i]*Xset[s][j]')/t end

aX=majd(Xset; fullModel=true, algorithm=:OJoB)
# the spForm index of the estimated demixing matrices times the true
# mixing matrix must be low
@test mean(spForm(aX.F[i]'*Vset[i]) for i=1:m)<0.1

# test the same using NoJoB algorithm
aX=majd(Xset; fullModel=true, algorithm=:NoJoB)
@test mean(spForm(aX.F[i]'*Vset[i]) for i=1:m)<0.1

# plot the original cross-covariance matrices and the rotated
# cross-covariance matrices

# Get all products 𝐔[i]' * 𝒞[l, i, j] * 𝐔[j]
function _rotate_crossCov(𝐔, 𝒞, m, k)
    𝒮=Array{Matrix}(undef, k, m, m)
    @inbounds for l=1:k, i=1:m, j=1:m 𝒮[l, i, j]=𝐔[i]'*𝒞[l, i, j]*𝐔[j] end
    return 𝒮
end

# Put all `k` cross-covariances in a single matrix
# of dimension m*n x m*n for visualization
function 𝒞2Mat(𝒞::AbstractArray, m, k)
    n=size(𝒞[1, 1, 1], 1)
    C=Matrix{Float64}(undef, m*n, m*n)
    for i=1:m, j=1:m, x=1:n, y=1:n C[i*n-n+x, j*n-n+y]=𝒞[k, i, j][x, y] end
    return C
end

using Plots

Cset=[𝒞2Mat(𝒞, m, s) for s=1:k]
 Cmax=maximum(maximum(abs.(C)) for C ∈ Cset)
 h1 = heatmap(Cset[1], clim=(-Cmax, Cmax), yflip=true, c=:bluesreds, title="all cross-cov, k=1")
 h2 = heatmap(Cset[2], clim=(-Cmax, Cmax), yflip=true, c=:bluesreds, title="all cross-cov, k=2")
 h3 = heatmap(Cset[2], clim=(-Cmax, Cmax), yflip=true, c=:bluesreds, title="all cross-cov, k=3")
 h4 = heatmap(Cset[2], clim=(-Cmax, Cmax), yflip=true, c=:bluesreds, title="all cross-cov, k=4")
 h5 = heatmap(Cset[2], clim=(-Cmax, Cmax), yflip=true, c=:bluesreds, title="all cross-cov, k=5")
 h6 = heatmap(Cset[2], clim=(-Cmax, Cmax), yflip=true, c=:bluesreds, title="all cross-cov, k=6")
 📈=plot(h1, h2, h3, h4, h5, h6, size=(1200,550))
# savefig(📈, homedir()*"\\Documents\\Code\\julia\\Diagonalizations\\docs\\src\\assets\\FigmAJD1.png")

𝒮=_rotate_crossCov(aX.F, 𝒞, m, k)
 Sset=[𝒞2Mat(𝒮, m, s) for s=1:k]
 Smax=maximum(maximum(abs.(S)) for S ∈ Sset)
 h11 = heatmap(Sset[1], clim=(-Smax, Smax), yflip=true, c=:bluesreds, title="rotated cross-cov, k=1")
 h12 = heatmap(Sset[2], clim=(-Smax, Smax), yflip=true, c=:bluesreds, title="rotated cross-cov, k=2")
 h13 = heatmap(Sset[2], clim=(-Smax, Smax), yflip=true, c=:bluesreds, title="rotated cross-cov, k=3")
 h14 = heatmap(Sset[2], clim=(-Smax, Smax), yflip=true, c=:bluesreds, title="rotated cross-cov, k=4")
 h15 = heatmap(Sset[2], clim=(-Smax, Smax), yflip=true, c=:bluesreds, title="rotated cross-cov, k=5")
 h16 = heatmap(Sset[2], clim=(-Smax, Smax), yflip=true, c=:bluesreds, title="rotated cross-cov, k=6")
 📉=plot(h11, h12, h13, h14, h15, h16, size=(1200,550))
# savefig(📉, homedir()*"\\Documents\\Code\\julia\\Diagonalizations\\docs\\src\\assets\\FigmAJD2.png")

# COMPLEX data
# do joint blind source separation of non-stationary data
t, m, n, k, noise = 200, 5, 4, 6, 0.1
Xcset, Vcset=getData(ComplexF64, t, m, k, n, noise)
𝒞=Array{Matrix}(undef, k, m, m)
for s=1:k, i=1:m, j=1:m 𝒞[s, i, j]=(Xcset[s][i]*Xcset[s][j]')/t end

aXc=majd(Xcset; fullModel=true, algorithm=:OJoB)
# the spForm index of the estimated demixing matrices times the true
# mixing matrix must be low
@test mean(spForm(aXc.F[i]'*Vcset[i]) for i=1:m)<0.1

# test the same using NoJoB algorithm
aXc=majd(Xcset; fullModel=true, algorithm=:NoJoB)
@test mean(spForm(aXc.F[i]'*Vcset[i]) for i=1:m)<0.1


text:
export update_particles!

"""
    update_particles!(brownian; period_x, period_y)

Advances the particles stored under `brownian.particles` by a timestep
`brownian.dt`.  Apply periodic boundaries with periodicity `period_x`,
`period_y` if `period_x > 0` or `period_y > 0`.
"""
function update_particles!(brownian::Brownian; period_x::Float64 = -1.0, period_y::Float64 = -1.0)
    @use_threads brownian.multithreaded for particle in brownian.particles
        if !isnothing(particle.active_force)
            particle.θ = update_active_force(particle.active_force; particle_orientation = particle.θ, dt = brownian.dt)

            af_x, af_y = get_active_force(particle.active_force)
            particle.f_x += af_x
            particle.f_y += af_y
        end

        dtγ_trans = brownian.dt / particle.γ_trans
        amplitude = sqrt(2 * particle.D_trans * brownian.dt)
        particle.x = wrap_position(particle.x + dtγ_trans * particle.f_x + amplitude * randn(); period = period_x)
        particle.y = wrap_position(particle.y + dtγ_trans * particle.f_y + amplitude * randn(); period = period_y)

        particle.f_x = 0.0
        particle.f_y = 0.0

        if brownian.rotations
            dtγ_rot = brownian.dt / particle.γ_rot
            particle.θ += dtγ_rot * particle.t_θ + sqrt(2 * particle.D_rot * brownian.dt) * randn()

            particle.t_θ = 0.0
        end
    end
end

ext: jl


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ fe9f9510-6055-11eb-2064-fb0d68f76e64
using PlutoUI

# ╔═╡ b7178270-6055-11eb-1d9d-1558690fb171
include("hodgkin_huxley.jl")

# ╔═╡ 3630e4c0-6056-11eb-3ab8-9d5c27ca9099
@bind input_current Slider(0:0.01:10, default=2)

# ╔═╡ 7f15a132-6056-11eb-1f09-8d04fefe528b
generate_plots(input_current)

# ╔═╡ aee4424e-605a-11eb-1b67-09c0df2fe9a0
input_current

# ╔═╡ Cell order:
# ╠═b7178270-6055-11eb-1d9d-1558690fb171
# ╠═fe9f9510-6055-11eb-2064-fb0d68f76e64
# ╠═3630e4c0-6056-11eb-3ab8-9d5c27ca9099
# ╠═7f15a132-6056-11eb-1f09-8d04fefe528b
# ╠═aee4424e-605a-11eb-1b67-09c0df2fe9a0

max_stars_repo_name: coinslab/ComputationalCognitiveModeling

max_issues_repo_name: coinslab/ComputationalCognitiveModeling

max_forks_repo_name: coinslab/ComputationalCognitiveModeling


text:
module TestChecker
    using Calculus2
    using Base.Test

    # ---
    # Univariate
    # ---

    f(x::Real) = sin(x)
    g(x::Real) = cos(x)
    h(x::Real) = -sin(x)

    @test isgrad(g, f, 0.0)
    @test isgrad(g, f, 1.0)
    @test isgrad(g, f, 10.0)
    @test isgrad(g, f, 100.0)
    @test isgrad(g, f, 1000.0)

    @test ishess(h, f, 0.0)
    @test ishess(h, f, 1.0)
    @test ishess(h, f, 10.0)
    @test ishess(h, f, 100.0)
    @test ishess(h, f, 1000.0)

    # ---
    # Multivariate
    # ---

    f(x::Vector) = sin(x[1]) + cos(x[2])
    g(x::Vector) = [cos(x[1]), -sin(x[2])]
    function g!(gr::Vector, x::Vector)
    	gr[1] = cos(x[1])
    	gr[2] = -sin(x[2])
        return gr
    end
    h(x::Vector) = [
    	-sin(x[1]) 0.0;
    	0.0        -cos(x[2]);
    ]
    function h!(H::Matrix, x::Vector)
    	H[1, 1] = -sin(x[1])
    	H[1, 2] = 0.0
    	H[2, 1] = 0.0
    	H[2, 2] = -cos(x[2])
        return H
    end

    @test isgrad(g, f, [0.0, 0.0])
    @test isgrad(g, f, [1.0, 1.0])
    @test isgrad(g, f, [10.0, 10.0])
    @test isgrad(g, f, [100.0, 100.0])
    @test isgrad(g, f, [1000.0, 1000.0])
    @test isgrad!(g!, f, [0.0, 0.0])
    @test isgrad!(g!, f, [1.0, 1.0])
    @test isgrad!(g!, f, [10.0, 10.0])
    @test isgrad!(g!, f, [100.0, 100.0])
    @test isgrad!(g!, f, [1000.0, 1000.0])

    @test ishess(h, f, [0.0, 0.0])
    @test ishess(h, f, [1.0, 1.0])
    @test ishess(h, f, [10.0, 10.0])
    @test ishess(h, f, [100.0, 100.0])
    @test ishess(h, f, [1000.0, 1000.0])
    @test ishess!(h!, f, [0.0, 0.0])
    @test ishess!(h!, f, [1.0, 1.0])
    @test ishess!(h!, f, [10.0, 10.0])
    @test ishess!(h!, f, [100.0, 100.0])
    @test ishess!(h!, f, [1000.0, 1000.0])
end


text:
module Hungarian

using LinearAlgebra
using SparseArrays

# Zero markers used in hungarian algorithm
# 0 => NON   => Non-zero
# 1 => Z     => ordinary zero
# 2 => STAR  => starred zero
# 3 => PRIME => primed zero
const NON = Int8(0)
const Z = Int8(1)
const STAR = Int8(2)
const PRIME = Int8(3)

include("./Munkres.jl")

"""
    hungarian(costMat) -> (assignment, cost)

Find an optimal solution of the rectangular assignment problem represented by
the ``N x M`` matrix `costMat`. Return the optimal column indices and
corresponding minimal cost.

The `costMat[n,m]` denotes the `cost` to assign the `n`th worker to the `m`th job.
The zero element in the return value `assignment` means that these workers have
no assigned job.

Elements in the matrix can be set to `missing`. In this case, the corresponding
matching cannot be considered by the algorithm.

# Examples
```julia
julia> A = [ 24     1     8;
              5     7    14;
              6    13    20;
             12    19    21;
             18    25     2];

julia> assignment, cost = hungarian(A)
([2,1,0,0,3],8)

julia> assignment, cost = hungarian(A')
([2,1,5],8)

julia> using Missings

julia> costMat = [ missing  1   1
                      1     0   1
                      1     1   0 ]
3×3 Array{Union{Float64, Missings.Missing},2}:
missing  1.0  1.0
  1.0    0.0  1.0
  1.0    1.0  0.0

julia> hungarian(costMat)
([2, 1, 3], 2)
```

"""
function hungarian(costMat::AbstractMatrix)
    rowNum, colNum = size(costMat)
    # currently, the function `hungarian` automatically transposes `cost matrix` when there are more workers than jobs.
    costMatrix = rowNum ≤ colNum ? costMat : costMat'
    matching = munkres(costMatrix)
    assignment = zeros(Int, rowNum)
    rows = rowvals(matching)
    for c = 1:size(matching,2), i in nzrange(matching, c)
        r = rows[i]
        if matching[r,c] == STAR
            if rowNum ≤ colNum
                assignment[r] = c
            else
                assignment[c] = r
            end
        end
    end
    # calculate minimum cost
    cost = sum(costMat[i...] for i in zip(1:rowNum, assignment) if i[2] != 0)
    return assignment, cost
end

export hungarian

end # module


text:
using Test

include("instructs.jl")
include("functions.jl")
include("simple_hessian.jl")
include("taylor.jl")


text:
import FireflyModule.Firefly

function path_cost(firefly::Firefly, distance_matrix)::Float32
    """Calculates the total path cost from given distance matrix.
    """
    path = firefly.path
    path_len = length(path) - 1
    p_cost = distance_matrix[path[end].id, path[1].id]
    for i in range(1, length=path_len)
        p_cost +=  distance_matrix[path[i].id, path[i+1].id]
    end

    return p_cost
end

max_stars_repo_name: zahidkizmaz/firefly-algorithm-tsp

max_issues_repo_name: zahidkizmaz/firefly-algorithm-tsp

max_forks_repo_name: zahidkizmaz/firefly-algorithm-tsp


text:
"""
Transform the results of GenGalaxy to Euclidean coordinates.
"""
function converttoeuclidean(coordinate::Array{G,1}) where G <: AbstractFloat
    radius = coordinate[1]
    theta = acos(1.0 - 2.0 * rand())
    phi = 2.0 * pi * rand()

    velocity = sqrt(coordinate[2]^2 + coordinate[3]^2)
    velSign = sign(rand() - 0.5)
    velPhi = 2.0 * pi * rand()

    x = radius * sin(theta) * cos(phi)
    y = radius * sin(theta) * sin(phi)
    z = radius * cos(theta)

    vz2 = velSign * coordinate[2]
    vx2 = coordinate[3] * cos(velPhi)
    vy2 = coordinate[3] * sin(velPhi)

    vx = cos(theta) * cos(phi) * vx2 - sin(phi) * vy2 + sin(theta) * cos(phi) * vz2
    vy = cos(theta) * sin(phi) * vx2 + cos(phi) * vy2 + sin(theta) * sin(phi) * vz2
    vz = -sin(theta) * vx2 + cos(theta) * vz2
    return x, y, z, vx, vy, vz
end


text:
__precompile__(true)

module LeastSquaresOptim

##############################################################################
##
## Dependencies
##
##############################################################################

import Base: A_mul_B!, Ac_mul_B!, A_ldiv_B!, copy!, fill!, scale!, norm, axpy!, eltype, length, size, call
import Base.SparseArrays.CHOLMOD: VTypes, ITypes, Sparse, Factor, C_Sparse, SuiteSparse_long, transpose_, @cholmod_name, common, defaults, set_print_level, common_final_ll, analyze, factorize_p!, check_sparse
using ForwardDiff

##############################################################################
##
## Exported methods and types
##
##############################################################################

export optimize!,
LeastSquaresProblem,
LeastSquaresProblemAllocated,
LeastSquaresResult

##############################################################################
##
## Load files
##
##############################################################################

include("utils/lsmr.jl")
include("utils/utils.jl")

include("types.jl")
include("optimizer/levenberg_marquardt.jl")
include("optimizer/dogleg.jl")

include("solver/dense_qr.jl")
include("solver/dense_cholesky.jl")
include("solver/sparse_cholesky.jl")
include("solver/iterative_lsmr.jl")

end


text:
_is_incremental(::Type{Ipopt.Optimizer}) = true

function set_default_config!(ext::ExtensionType, d::GlobalOptimizer, m::Incremental{Ipopt.Optimizer}, local_solver::Bool)
    #=
    c_tol = _constraint_tol(d)*1E-3
    MOI.set(m, MOI.RawOptimizerAttribute("tol"), _absolute_tol(d)*1E-3)
    MOI.set(m, MOI.RawOptimizerAttribute("print_level"), 0)
    MOI.set(m, MOI.RawOptimizerAttribute("constr_viol_tol"), c_tol)
    if local_solver
        MOI.set(m, MOI.RawOptimizerAttribute("max_iter"), 20000)
        MOI.set(m, MOI.RawOptimizerAttribute("acceptable_tol"), 1E30)
        MOI.set(m, MOI.RawOptimizerAttribute("acceptable_iter"), 10000)
        MOI.set(m, MOI.RawOptimizerAttribute("acceptable_compl_inf_tol"), c_tol)
        MOI.set(m, MOI.RawOptimizerAttribute("acceptable_dual_inf_tol"), 1.0)
        MOI.set(m, MOI.RawOptimizerAttribute("acceptable_constr_viol_tol"), c_tol)
    else
        MOI.set(m, MOI.RawOptimizerAttribute("max_iter"), 1E5)
        MOI.set(m, MOI.RawOptimizerAttribute("acceptable_iter"), 1E5+1)
    end
    =#
    return
end

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
module POMDPLinter

using Logging

export
    implemented,
    @implemented,
    RequirementSet,
    check_requirements,
    show_requirements,
    get_requirements,
    requirements_info,
    @POMDP_require,
    @POMDP_requirements,
    @requirements_info,
    @get_requirements,
    @show_requirements,
    @warn_requirements,
    @req,
    @subreq

include("requirements_internals.jl")
include("requirements_printing.jl")
include("requirements_interface.jl")

end


text:
module SimplePadics

using Base: String, prec_arrow, get_preferences
using Nemo, SimplePolynomials

_PF = PadicField

export Padic, Dyadic

const _DEFAULT_PRECISION = 10
_PREC = _DEFAULT_PRECISION    # standard precision for p-adic numbers

export set_precision, get_precision

"""
    set_precision(prec::Integer)
Controls the number of digits of precision for `Padic` numbers.
Without an argument, set the precision to a default value (`_DEFAULT_PRECISION`).
"""
function set_precision(prec::Integer = _DEFAULT_PRECISION)
    @assert prec > 1 "Precision must be set to 2 or higher"
    global _PREC = prec
end

"""
    get_precision()
Return the number of digits of precision for `Padic` numbers.
"""
get_precision() = _PREC

"""
    Padic

Constructor for `p`-adic numbers:
* `Padic{P}(n::Integer)`
* `Padic{P}(r::Rational)`
"""
struct Padic{P} <: Number
    x::padic
    function Padic{P}(s::Integer) where {P}
        F = _PF(P, _PREC)
        new(F(s))
    end

    function Padic{P}(r::Rational) where {P}
        F = _PF(P, _PREC)
        a = F(numerator(r)) // F(denominator(r))
        new(a)
    end

    function Padic(s::padic)   # conver Nemo padic to our Padic
        F = s.parent
        digits = F.prec_max
        P = F.p
        new{P}(s)
    end

    function Padic(s::Padic)
        return s
    end
end

"""
    Dyadic

This is an abbrevation for `Padic{2}`.
"""
Dyadic = Padic{2}

function padic2str(a::padic)::String
    r = lift(QQ, a)
    t = numerator(r)
    b = BigInt(denominator(r))
    p = a.parent.p
    prefix = "…"

    digs = base(t, p)

    if b > 1   # we need to shift to
        ndigs = length(digs)
        shift = Int(round(log(p, b)))
        if shift < ndigs
            digs = digs[1:end-shift] * "." * digs[end-shift+1:end]
        elseif shift == ndigs
            digs = "0." * digs
        else
            pad = "0"^(shift - ndigs + 1)
            digs = pad * digs
            digs = digs[1:end-shift] * "." * digs[end-shift+1:end]
        end
    end

    suffix = valuation(a) >= 0 ? ".0" : ""

    return prefix * digs * suffix
end

function padic2str(a::Padic{P})::String where {P}
    s = padic2str(a.x)
    return s * "_{$P}"
end

import Base.Multimedia.display
display(a::Padic) = println(padic2str(a))

import Base.show
show(io::IO, a::Padic) = print(io, padic2str(a))

import Base: digits

"""
    digits(a::Padic)

Return the digits of `a`, starting from the right.
```
julia> a = Padic{7}(124//49)
…2.35_{7}

julia> digits(a)
3-element Vector{Int64}:
 5
 3
 2
```
"""
function digits(a::Padic{P})::Vector{Int} where {P}
    z = numerator(lift(QQ, a.x))
    return digits(z, base = P)
end

include("arithmetic.jl")
include("functions.jl")
include("hensel.jl")

end # module


text:
using Flux: testmode!

@testset "Dropout" begin
  x = [1.,2.,3.]
  @test x == testmode!(Dropout(0.1))(x)
  @test x == Dropout(0)(x)
  @test zeros(x) == Dropout(1)(x)

  x = rand(100)
  m = Dropout(0.9)
  y = m(x)
  @test count(a->a==0, y) > 50
  testmode!(m)
  y = m(x)
  @test count(a->a==0, y) == 0
  testmode!(m, false)
  y = m(x)
  @test count(a->a==0, y) > 50

  x = rand(100)
  m = Chain(Dense(100,100),
            Dropout(0.9))
  y = m(x)
  @test count(a->a == 0, y) > 50
  testmode!(m)
  y = m(x)
  @test count(a->a == 0, y) == 0
end

@testset "BatchNorm" begin
  let m = BatchNorm(2), x = param([1 2; 3 4; 5 6]')

    @test m.β.data == [0, 0]  # initβ(2)
    @test m.γ.data == [1, 1]  # initγ(2)
    # initial m.σ is 1
    # initial m.μ is 0
    @test m.active

    # @test m(x).data ≈ [-1 -1; 0 0; 1 1]'
    m(x)

    # julia> x
    #  2×3 Array{Float64,2}:
    #  1.0  3.0  5.0
    #  2.0  4.0  6.0
    #
    # μ of batch will be
    #  (1. + 3. + 5.) / 3 = 3
    #  (2. + 4. + 6.) / 3 = 4
    #
    # ∴ update rule with momentum:
    #  .1 * 3 + 0 = .3
    #  .1 * 4 + 0 = .4
    @test m.μ ≈ reshape([0.3, 0.4], 2, 1)

    # julia> .1 .* std(x, 2, corrected=false) .* (3 / 2).+ .9 .* [1., 1.]
    # 2×1 Array{Float64,2}:
    #  1.14495
    #  1.14495
    @test m.σ ≈ .1 .* std(x.data, 2, corrected=false) .* (3 / 2).+ .9 .* [1., 1.]

    testmode!(m)
    @test !m.active

    x′ = m(x).data
    @test x′[1] ≈ (1 - 0.3) / 1.1449489742783179
  end

  # with activation function
  let m = BatchNorm(2, σ), x = param([1 2; 3 4; 5 6]')
    @test m.active
    m(x)

    testmode!(m)
    @test !m.active

    x′ = m(x).data
    @test x′[1] ≈ σ((1 - 0.3) / 1.1449489742783179)
  end

  let m = BatchNorm(2), x = param(reshape(1:6, 3, 2, 1))
    y = reshape(permutedims(x, [2, 1, 3]), 2, :)
    y = permutedims(reshape(m(y), 2, 3, 1), [2, 1, 3])
    @test m(x) == y
  end

  let m = BatchNorm(2), x = param(reshape(1:12, 2, 3, 2, 1))
      y = reshape(permutedims(x, [3, 1, 2, 4]), 2, :)
    y = permutedims(reshape(m(y), 2, 2, 3, 1), [2, 3, 1, 4])
    @test m(x) == y
  end

  let m = BatchNorm(2), x = param(reshape(1:24, 2, 2, 3, 2, 1))
    y = reshape(permutedims(x, [4, 1, 2, 3, 5]), 2, :)
    y = permutedims(reshape(m(y), 2, 2, 2, 3, 1), [2, 3, 4, 1, 5])
    @test m(x) == y
  end
end


text:
using Combinatorics
using StatsBase

# Bounds: puzzle input
const LOWER = 387638
const UPPER = 919123

function int_from_digits(digits::AbstractVector)
    range = length(digits)-1:-1:0
    sum(@. digits * 10^range)
end

function countcombs(f)
    combs = Iterators.filter(f, with_replacement_combinations(1:9, 6))
    # Just iterate to avoid materializing a large array
    n = 0
    for c in combs
        n += 1
    end
    n
end

println("Part 1 answer: ", countcombs() do c
    issorted(c) && any(iszero, diff(c)) && LOWER ≤ int_from_digits(c) ≤ UPPER
end)

println("Part 2 answer: ", countcombs() do c
    issorted(c) && LOWER ≤ int_from_digits(c) ≤ UPPER && any(==(2), last(rle(c)))
end)

max_stars_repo_name: ararslan/advent-solutions

max_issues_repo_name: ararslan/advent-solutions

max_forks_repo_name: ararslan/advent-solutions


text:
# Copyright 2017-2018 Gandalf Software, Inc. (Scott Paul Jones)
# Licensed under MIT License, see LICENSE.md

# These use the "Holy Trait Trick", which was created by @timholy (Tim Holy),
# which was made possible by the type system in Julia from @jeff.bezanson (Jeff Bezanson).

# One is proud to stand on the shoulders of such people!

## Traits for string and character types ##

"""
    ValidatedStyle(A)
    ValidatedStyle(typeof(A))

`ValidatedStyle` specifies the whether a string or character type is always valid or not
When you define a new `AbstractString` or `AbstractChar` type, you can choose to implement it
as always validated, or validation state unknown.

    StrAPI.ValidatedStyle(::Type{<:MyCharType}) = StrAPI.AlwaysValid()

The default is `UnknownValidity`
"""
abstract type ValidatedStyle end
struct AlwaysValid     <: ValidatedStyle end
struct UnknownValidity <: ValidatedStyle end

# single or multiple codeunits per codepoint

"""
    EncodingStyle(A)
    EncodingStyle(typeof(A))

`EncodingStyle` specifies the whether a character set encoding uses one or multiple codeunits to
encode a single codepoint.
When you define a new `AbstractString` type, you can choose to implement it with either
single or multi-codeunit indexing.

    StrAPI.EncodingStyle(::Type{MyCharSetEncoding}) = StrAPI.MultiCodeUnitEncoding()

The default is `SingleCodeUnitEncoding()`
"""
abstract type EncodingStyle end

struct SingleCodeUnitEncoding <: EncodingStyle end
struct MultiCodeUnitEncoding  <: EncodingStyle end

# Type of character set

"""
    CharSetStyle(A)
    CharSetStyle(typeof(A))

`CharSetStyle` specifies the information about the character set used by the string or
characters.
When you define a new `AbstractString` or `AbstractChar` type,
you can choose to implement it with

    StrAPI.CharSetStyle(::Type{<:MyString}) = StrAPI.CharSetISOCompat()

The default is `CharSetUnicode()`
"""
abstract type CharSetStyle end

"""Codepoints are not in Unicode compatible order"""
struct CharSetOther         <: CharSetStyle end
"""Characters 0-0x7f same as ASCII"""
struct CharSetASCIICompat   <: CharSetStyle end
"""Characters 0-0x9f follows ISO 8859"""
struct CharSetISOCompat     <: CharSetStyle end
"""Characters 0-0xd7ff, 0xe000-0xffff follow Unicode BMP"""
struct CharSetBMP           <: CharSetStyle end
"""Full Unicode character set, no additions"""
struct CharSetUnicode       <: CharSetStyle end
"""Unicode character set, plus encodings of invalid characters"""
struct CharSetUnicodePlus   <: CharSetStyle end
"""8-bit Binary string, not text"""
struct CharSetBinary        <: CharSetStyle end
"""Raw bytes, words, or character string, unknown encoding/character set"""
struct CharSetUnknown       <: CharSetStyle end

# Comparison traits

"""
    CompareStyle(Union{A, typeof(A)}, Union{B, typeof(B)})

`CompareStyle` specifies how to compare two strings with character set encodings A and B

    StrAPI.CompareStyle(::Type{<:MyString}, ::Type{String}) = StrAPI.ByteCompare()

The default is `CodePointCompare`
"""
abstract type CompareStyle end

struct NoCompare        <: CompareStyle end # For equality checks, can't be equal
struct ByteCompare      <: CompareStyle end # Compare bytewise
struct ASCIICompare     <: CompareStyle end # Compare bytewise for ASCII subset, else codepoint
struct WordCompare      <: CompareStyle end # Compare first not equal word with <
struct UTF16Compare     <: CompareStyle end # Compare first not equal word adjusting > 0xd7ff
struct WidenCompare     <: CompareStyle end # Narrower can be simply widened for comparisons
struct CodePointCompare <: CompareStyle end # Compare CodePoints

promote_rule(::Type{NoCompare},    ::Type{<:CompareStyle})   = NoCompare
promote_rule(::Type{ByteCompare},  ::Type{CodePointCompare}) = ByteCompare
promote_rule(::Type{WordCompare},  ::Type{CodePointCompare}) = WordCompare
promote_rule(::Type{UTF16Compare}, ::Type{CodePointCompare}) = UTF16Compare
promote_rule(::Type{ASCIICompare}, ::Type{CodePointCompare}) = ASCIICompare
promote_rule(::Type{WidenCompare}, ::Type{CodePointCompare}) = WidenCompare

"""Determine if a string has multiple codeunit encoding"""
function is_multi end

is_multi(::Type{T}) where {T<:AbstractString} = EncodingStyle(T) === MultiCodeUnitEncoding()
is_multi(::T) where {T<:AbstractString} = is_multi(T)

@api develop AlwaysValid, UnknownValidity,
             SingleCodeUnitEncoding, MultiCodeUnitEncoding,
             CharSetOther, CharSetBinary, CharSetASCIICompat, CharSetISOCompat,
             CharSetBMP, CharSetUnicode, CharSetUnicodePlus, CharSetUnknown,
             NoCompare, CodePointCompare, ByteCompare, ASCIICompare, WordCompare,
             UTF16Compare, WidenCompare

@api develop! ValidatedStyle, EncodingStyle, CharSetStyle, CompareStyle, is_multi


text:
function visible_asteroids_at_coord(grid, coord)
    # Search outwards in "shells" from coord
    remaining = copy(grid)
    for r in 1:maximum(size(grid))
        shell = map(x -> coord + x, setdiff(CartesianIndices((-r:r, -r:r)), CartesianIndices((-(r-1):(r-1), -(r-1):(r-1)))))
        for i in filter(x -> x ∈ CartesianIndices(remaining) && remaining[x], shell)
            # We can see this asteroid. Find the minimum integer direction vector.
            delta = CartesianIndex(tuple((x -> x .÷ gcd(x))(collect((i - coord).I))...))
            # Delete all occluded locations by stepping along the direction vector
            occluded = i + delta
            while occluded ∈ CartesianIndices(grid)
                remaining[occluded] = false
                occluded += delta
            end
        end
    end
    # Remember to remove the asteroid we're sitting on
    remaining[coord] = false
    return remaining
end

# Part 1: read asteroid map as bool matrix and find out which one can detect most asteroids
grid = hcat(map(line -> '#' .== collect(line), open(f -> readlines(f), "day10_input.txt"))...)'
visible_counts = map(x -> grid[x] ? count(visible_asteroids_at_coord(grid, x)) : 0, CartesianIndices(grid))
println("Max detected asteroids: $(maximum(visible_counts))")

# Part 2: get coord of best station, then work out which is the 200th asteroid to be vaporized
# We actually have >200 detected, so we need less than one rotation. Sort detected by angle.
coord = argmax(visible_counts)
detected = visible_asteroids_at_coord(grid, coord)
item = sort(map(x -> (π - atan(reverse((x - coord).I)...), x), filter(x -> detected[x], CartesianIndices(detected))))[200]
println("200th to be vaporized: $item, checksum: $((item[2][2] - 1) * 100 + item[2][1] - 1))")

max_stars_repo_name: ericgreveson/adventofcode

max_issues_repo_name: ericgreveson/adventofcode

max_forks_repo_name: ericgreveson/adventofcode


text:
using IntArrays
using FactCheck

srand(12345)

const Ts = (UInt8, UInt16, UInt32, UInt64)

facts("conversion") do
    context("IntVector") do
        data = [0x00, 0x01, 0x02]
        ivec = IntArray{2}(data)
        @fact typeof(ivec) --> IntArray{2,UInt8,1}
        @fact typeof(ivec) --> IntVector{2,UInt8}

        ivec = convert(IntVector{2}, data)
        @fact typeof(ivec) --> IntArray{2,UInt8,1}
        @fact typeof(ivec) --> IntVector{2,UInt8}
    end
    context("IntMatrix") do
        data = [0x00 0x01; 0x02 0x03; 0x04 0x05]
        imat = IntArray{3}(data)
        @fact typeof(imat) --> IntArray{3,UInt8,2}
        @fact typeof(imat) --> IntMatrix{3,UInt8}

        imat = convert(IntMatrix{3}, data)
        @fact typeof(imat) --> IntArray{3,UInt8,2}
        @fact typeof(imat) --> IntMatrix{3,UInt8}
    end
    context("three-dimensional array") do
        data = rand(0x00:0x03, 2, 3, 4)
        iarr = IntArray{2}(data)
        @fact typeof(iarr) --> IntArray{2,UInt8,3}
    end
end

facts("construction") do
    context("IntVector") do
        ivec = IntArray{9,UInt16}(10)
        @fact typeof(ivec) --> IntArray{9,UInt16,1}
        @fact size(ivec) --> (10,)
        @fact length(ivec) --> 10
    end
    context("IntMatrix") do
        imat = IntMatrix{3,UInt8}(2, 3)
        @fact typeof(imat) --> IntArray{3,UInt8,2}
        @fact size(imat) --> (2, 3)
    end
    context("invalid width") do
        @fact_throws Exception IntArray{9,UInt8,1}(4)
        @fact_throws Exception IntArray{20,UInt16,1}(4)
        @fact_throws Exception IntArray{50,UInt32,1}(4)
    end
    context("mmap") do
        @fact typeof(IntArray{2,UInt8}(10, true)) --> IntArray{2,UInt8,1}
    end
end

facts("similar") do
    context("IntVector") do
        data = rand(0x00:0x02, 10)
        ivec = IntVector{2}(data)
        @fact similar(ivec) === ivec --> false
        @fact size(similar(ivec)) --> (10,)
        @fact typeof(similar(ivec)) --> IntVector{2,UInt8}
        @fact size(similar(ivec, UInt16)) --> (10,)
        @fact typeof(similar(ivec, UInt16)) --> IntVector{2,UInt16}
        @fact size(similar(ivec, UInt16, (20,))) --> (20,)
        @fact typeof(similar(ivec, UInt16, (20,))) --> IntVector{2,UInt16}
    end
end

facts("getindex") do
    context("IntVector") do
        data = [0x00, 0x01, 0x02, 0x03, 0x04]
        ivec = IntArray{3}(data)
        @fact_throws BoundsError ivec[0]
        @fact ivec[1] --> 0x00
        @fact ivec[2] --> 0x01
        @fact ivec[3] --> 0x02
        @fact ivec[4] --> 0x03
        @fact ivec[5] --> 0x04
        @fact_throws BoundsError ivec[6]
    end
    context("IntMatrix") do
        # 2x3
        data = [0x00 0x01 0x02; 0x03 0x04 0x05]
        imat = IntArray{4}(data)
        # linear indexing
        @fact_throws BoundsError imat[0]
        @fact imat[1] --> 0x00
        @fact imat[2] --> 0x03
        @fact imat[3] --> 0x01
        @fact imat[4] --> 0x04
        @fact imat[5] --> 0x02
        @fact imat[6] --> 0x05
        @fact_throws BoundsError imat[7]
        # tuples
        @fact_throws BoundsError imat[0,1]
        @fact imat[1,1] --> 0x00
        @fact imat[1,2] --> 0x01
        @fact imat[1,3] --> 0x02
        @fact imat[2,1] --> 0x03
        @fact imat[2,2] --> 0x04
        @fact imat[2,3] --> 0x05
        @fact_throws BoundsError imat[1,4]
        @fact_throws BoundsError imat[2,4]
        @fact_throws BoundsError imat[3,1]
    end
end

facts("setindex!") do
    context("IntVector") do
        data = [0x00, 0x01, 0x02]
        ivec = IntArray{2}(data)
        @fact_throws BoundsError ivec[0] = 0x00
        ivec[1] = 0x01
        @fact ivec[1] --> 0x01
        ivec[2] = 0x02
        @fact ivec[2] --> 0x02
        ivec[3] = 0x03
        @fact ivec[3] --> 0x03
        @fact_throws BoundsError ivec[4] = 0x01
    end
    context("IntMatrix") do
        data = [0x00 0x00 0x00; 0x00 0x00 0x00]
        imat = IntArray{2}(data)
        # linear
        @fact_throws BoundsError imat[0] = 0x00
        imat[2] = 0x01
        @fact imat[2] --> 0x01
        imat[2] = 0x00
        @fact imat[2] --> 0x00
        @fact_throws BoundsError imat[7] = 0x00
        # tuple
        imat[1,1] = 0x01
        @fact imat[1,1] --> 0x01
        imat[1,3] = 0x03
        @fact imat[1,3] --> 0x03
        @fact_throws BoundsError imat[1,4] = 0x00
        @fact_throws BoundsError imat[3,1] = 0x01
    end
end

facts("comparison") do
    a = IntVector{2}([0x00, 0x01, 0x02, 0x03])
    b = IntVector{2}([0x00, 0x01, 0x02, 0x03])
    c = IntVector{4}([0x00, 0x01, 0x02, 0x03])
    d = IntVector{2}([0x00, 0x02, 0x02, 0x03])
    e = IntMatrix{2}([0x00  0x02; 0x01  0x03])
    @fact a == b --> true
    @fact a == c --> true
    @fact a == d --> false
    @fact a == e --> false
end

facts("sizeof") do
    context("smaller bits") do
        n = 100
        data = rand(0x00:0x01, n)
        @fact sizeof(IntVector{1}(data)) --> less_than(sizeof(data))
        @fact sizeof(IntVector{2}(data)) --> less_than(sizeof(data))
        @fact sizeof(IntVector{3}(data)) --> less_than(sizeof(data))
        @fact sizeof(IntVector{4}(data)) --> less_than(sizeof(data))
    end
end

facts("copy") do
    context("same length") do
        data = rand(0x00:0x03, 10)
        ivec = IntVector{2,UInt8}(data)
        @fact copy(ivec) --> ivec
        @fact typeof(copy(ivec)) --> typeof(ivec)
        @fact size(copy(ivec)) --> size(ivec)
        ivec′ = IntVector{2,UInt8}(10)
        @fact copy!(ivec′, ivec) === ivec′ --> true
        @fact ivec′ --> ivec
    end
    context("larger") do
        data = rand(0x00:0x03, 10)
        ivec = IntVector{2,UInt8}(data)
        ivec′ = IntVector{2,UInt8}(20)
        copy!(ivec′, ivec)
        @fact ivec′[1:10] --> ivec
    end
    context("smaller") do
        data = rand(0x00:0x03, 10)
        ivec = IntVector{2,UInt8}(data)
        ivec′ = IntVector{2,UInt8}(5)
        @fact_throws BoundsError copy!(ivec′, ivec)
    end
end

facts("fill!") do
    context("IntVector") do
        n = 100
        data = rand(0x00:0x03, n)
        for x in 0x00:0x03
            ivec = IntVector{2}(data)
            @fact fill!(ivec, x) === ivec --> true
            @fact ivec --> ones(UInt8, n) * x
        end
    end
    context("IntMatrix") do
        m, n = 10, 11
        data = rand(0x00:0x03, (m, n))
        for x in 0x00:0x03
            #imat = IntMatrix{2}(m, n)
            #@fact fill!(imat, x) === imat --> true
            #@fact imat --> ones(UInt8, (m, n)) * x
        end
    end
end

facts("reverse") do
    ivec = IntVector{2,UInt8}()
    @fact reverse!(ivec) === ivec --> true
    @fact ivec --> isempty

    ivec = IntVector{2}([0x00])
    @fact reverse!(ivec) === ivec --> true
    @fact ivec --> [0x00]

    ivec = IntVector{2}([0x00, 0x01])
    @fact reverse!(ivec) === ivec --> true
    @fact ivec --> [0x01, 0x00]

    ivec = IntVector{2}([0x00, 0x01, 0x02])
    @fact reverse!(ivec) === ivec --> true
    @fact ivec --> [0x02, 0x01, 0x00]

    ivec = IntVector{2}([0x00, 0x01, 0x02, 0x03])
    @fact reverse!(ivec) === ivec --> true
    @fact ivec --> [0x03, 0x02, 0x01, 0x00]
end

facts("push!/pop!") do
    ivec = IntVector{4,UInt8}()
    @fact length(ivec) --> 0
    @fact push!(ivec, 3) === ivec --> true
    @fact length(ivec) --> 1
    @fact ivec[end] --> 3
    @fact pop!(ivec) --> 3
    @fact length(ivec) --> 0
end

facts("append!") do
    ivec = IntVector{4}([0x00])
    @fact append!(ivec, [0x01, 0x02]) === ivec --> true
    @fact ivec --> [0x00, 0x01, 0x02]
    append!(ivec, [3, 4])
    @fact ivec --> [0x00, 0x01, 0x02, 0x03, 0x04]
end

facts("radixsort") do
    data = UInt8[]
    ivec = IntVector{1}(data)
    @fact radixsort(ivec) --> issorted

    data = [0x00]
    ivec = IntVector{2}(data)
    @fact radixsort(ivec) --> issorted

    data = [0x01, 0x00]
    ivec = IntVector{2}(data)
    @fact radixsort(ivec) --> issorted

    n = 101
    data = rand(0x00:0x01, n)
    ivec = IntVector{1}(data)
    @fact radixsort(ivec) --> issorted
    @fact radixsort!(ivec) === ivec --> true
    @fact issorted(ivec) --> true
    data = rand(0x00:0x03, n)
    ivec = IntVector{2}(data)
    @fact radixsort(ivec) --> issorted
    @fact radixsort!(ivec) === ivec --> true
    @fact issorted(ivec) --> true
    data = rand(0x00:0x07, n)
    ivec = IntVector{3}(data)
    @fact radixsort(ivec) --> issorted
    @fact radixsort!(ivec) === ivec --> true
    @fact issorted(ivec) --> true
end

# thorough and time-consuming tests for each combination of width and element type

facts("conversion") do
    context("empty") do
        for T in Ts, w in 1:sizeof(T)*8
            data = Vector{T}(0)
            @fact typeof(IntArray{w}(data)) --> IntArray{w,T,1}
            @fact length(IntArray{w}(data)) --> 0

            data = Matrix{T}(0, 0)
            @fact typeof(IntArray{w}(data)) --> IntArray{w,T,2}
            @fact size(IntArray{w}(data)) --> (0, 0)
        end
    end
    context("small") do
        for T in Ts, w in 1:sizeof(T)*8
            data = T[0x00, 0x01]
            @fact typeof(IntArray{w}(data)) --> IntArray{w,T,1}
            @fact length(IntArray{w}(data)) --> 2

            data = T[0x00 0x01 0x00; 0x01 0x00 0x01]
            @fact typeof(IntArray{w}(data)) --> IntArray{w,T,2}
            @fact size(IntArray{w}(data)) --> (2, 3)
        end
    end
    context("large") do
        n = 1000
        for T in Ts, w in 1:sizeof(T)*8
            data = rand(T, n)
            @fact typeof(IntArray{w}(data)) --> IntArray{w,T,1}
            @fact length(IntArray{w}(data)) --> n
        end
    end
end

facts("random run") do
    context("IntVector") do
        n = 123
        for T in Ts, w in 1:sizeof(T)*8
            data = rand(T(0):T(2)^w-T(1), n)
            ivec = IntVector{w,T}(data)
            for i in 1:n
                @fact ivec[i] --> data[i]
            end

            # random update
            for _ in 1:100
                i = rand(1:n)
                x::T = rand(T) % w
                data[i] = x
                ivec[i] = x
                @fact ivec[i] --> data[i]
            end
            for i in 1:n
                @fact ivec[i] --> data[i]
            end
            @fact_throws BoundsError ivec[0]
            @fact_throws BoundsError ivec[n+1]

            # sort
            sort!(data)
            radixsort!(ivec)
            for i in 1:n
                @fact ivec[i] --> data[i]
            end

            # random push!/pop!
            for _ in 1:100
                if rand() < 0.5
                    x::T = rand(T) % w
                    push!(ivec, x)
                    push!(data, x)
                else
                    pop!(ivec)
                    pop!(data)
                end
                @fact length(ivec) --> length(data)
                @fact ivec[end] --> data[end]
            end
            while !isempty(ivec)
                x = pop!(ivec)
                y = pop!(data)
                @fact x --> y
            end
            @fact isempty(ivec) --> true
        end
    end
    context("IntMatrix") do
        m, n = 11, 28
        for T in Ts, w in 1:sizeof(T)*8
            data = rand(T(0):T(2)^w-T(1), (m, n))
            imat = IntMatrix{w,T}(data)
            for i in 1:m, j in 1:n
                @fact imat[i,j] --> data[i,j]
            end

            # random update
            for _ in 1:100
                i = rand(1:m)
                j = rand(1:n)
                x::T = rand(T) % w
                data[i,j] = x
                imat[i,j] = x
                @fact imat[i,j] --> data[i,j]
            end
            for i in 1:m, j in 1:n
                @fact imat[i,j] --> data[i,j]
            end
        end
    end
end


text:
module NamedAxesArrays

import Base: eltype, getindex, ndims, setindex!

immutable NamedAxesArray{T,N,A<:AbstractArray,axisnames<:(Symbol...)} <: AbstractArray{T,N}
    data::A
end
NamedAxesArray{T,N}(A::AbstractArray{T,N}, axisnames::NTuple{N,Symbol}) = NamedAxes{T,N,typeof(A),axisnames}(A)

   eltype{T,N,A,axisnames}(::Type{NamedAxesArray{T,N,A,axisnames}}) = T
    ndims{T,N,A,axisnames}(::Type{NamedAxesArray{T,N,A,axisnames}}) = N
axisnames{T,N,A,axisnames}(::Type{NamedAxesArray{T,N,A,axisnames}}) = axisnames
   eltype{T,N,A,axisnames}(A::NamedAxesArray{T,N,A,axisnames}) = T
    ndims{T,N,A,axisnames}(A::NamedAxesArray{T,N,A,axisnames}) = N
axisnames{T,N,A,axisnames}(A::NamedAxesArray{T,N,A,axisnames}) = axisnames

immutable Ax{axisname<:Symbol,T}
    I::T
end
Ax{T}(name::Symbol, I::T) = Ax{name,T}(I)

axisname{S,T}(i::Ax{S,T}) = S
axistype{S,T}(i::Ax{S,T}) = T
axisname{S,T}(::Type{Ax{S,T}}) = S
axistype{S,T}(::Type{Ax{S,T}}) = T

### Traditional position-based indexing ###
typealias RegularIndex{R} Union(R, AbstractVector{R})
# Efficient unsplatted variants
let params = Expr[], args = Expr[], indexes = Symbol[]
for i = 1:4
    Rsym = symbol("R$i")
    isym = symbol("i$i")
    push!(params, :($Rsym<:Real))
    push!(args, :($isym::RegularIndex{$Rsym}))
    push!(indexes, isym)
    @eval begin
        getindex{$(params...)}(A::NamedAxesArray, $(args...)) = A.data[$(indexes...)]
        setindex!{$(params...)}(A::NamedAxesArray, v, $(args...)) = A.data[$(indexes...)] = v
    end
end
# Splatted variants for higher dimensions
getindex(A::NamedAxesArray, I::Union(Real, AbstractVector)...) = A.data[I...]
setindex!(A::NamedAxesArray, v, I::Union(Real, AbstractVector)...) = A.data[I...] = v

### Indexing using the axis names ###
# For integer-valued dimensions, this has slicing semantics
stagedfunction getindex(A::NamedAxesArray, i1::Ax)
    getindex_gen(A, i1)
end
stagedfunction getindex(A::NamedAxesArray, i1::Ax, i2::Ax)
    getindex_gen(A, i1, i2)
end
stagedfunction getindex(A::NamedAxesArray, i1::Ax, i2::Ax, i3::Ax)
    getindex_gen(A, i1, i2, i3)
end
stagedfunction getindex(A::NamedAxesArray, i1::Ax, i2::Ax, i3::Ax, i4::Ax)
    getindex_gen(A, i1, i2, i3, i4)
end

function getindex_gen(A, indexes...)
    anames = [map(axisname, indexes)...]
    Anames = [axisnames(A)...]
    nd = ndims(A)
    ind = indexin(anames, Anames)
    for i = 1:length(ind)
        if ind[i] == 0
        return :(error("axis name $(anames[i]) is not in $(axisnames(A))"))
    end
    indexexprs = [:(1:size(A,d)) for d = 1:nd]
    indextypes = [Range{Int} for d = 1:nd]
    for i = 1:length(ind)
        indexexprs[ind[i]] = :(indexes[i].I)
        indextypes[ind[i]] = axistype(indexes[i])
    end
    sliceindexes = filter(d -> axistypes[d]<:Real, 1:nd)
    deleteat!(Anames, sliceindexes)
    AvT = viewapl1(A, indextypes...)  # viewapl1 throws an error when indexing with higher-dimensional array indexes
    :(NamedAxesArray{$(eltype(AvT)),$(ndims(AvT)), $AvT, $(tuple(Anames...))}(viewapl1(A, $(indexexprs...))))
end

end


text:
function run_rayleigh_benard_regression_test(arch)

    #####
    ##### Parameters
    #####
          α = 2                 # aspect ratio
          n = 1                 # resolution multiple
         Ra = 1e6               # Rayleigh number
    Nx = Ny = 8n * α            # horizontal resolution
    Lx = Ly = 1.0 * α           # horizontal extent
         Nz = 16n               # vertical resolution
         Lz = 1.0               # vertical extent
         Pr = 0.7               # Prandtl number
          a = 1e-1              # noise amplitude for initial condition
         Δb = 1.0               # buoyancy differential

    # Rayleigh and Prandtl determine transport coefficients
    ν = sqrt(Δb * Pr * Lz^3 / Ra)
    κ = ν / Pr

    #####
    ##### Model setup
    #####

    # Force salinity as a passive tracer (βS=0)
    c★(x, z) = exp(4z) * sin(2π/Lx * x)
    Fc(i, j, k, grid, time, U, C, params) = 1/10 * (c★(grid.xC[i], grid.zC[k]) - C.c[i, j, k])

    bbcs = HorizontallyPeriodicBCs(   top = BoundaryCondition(Value, 0.0),
                                   bottom = BoundaryCondition(Value, Δb))

    model = Model(
               architecture = arch,
                       grid = RegularCartesianGrid(size=(Nx, Ny, Nz), length=(Lx, Ly, Lz)),
                    closure = ConstantIsotropicDiffusivity(ν=ν, κ=κ),
                    tracers = (:b, :c),
                   buoyancy = BuoyancyTracer(),
        boundary_conditions = HorizontallyPeriodicSolutionBCs(b=bbcs),
                    forcing = ModelForcing(c=Fc)
    )

    # The type of the underlying data, not the offset array.
    ArrayType = typeof(model.velocities.u.data.parent)

    Δt = 0.01 * min(model.grid.Δx, model.grid.Δy, model.grid.Δz)^2 / ν

    spinup_steps = 1000
      test_steps = 100

    output_U(model) = datatuple(model.velocities)
    output_Φ(model) = datatuple(model.tracers)
    output_G(model) = datatuple(model.timestepper.Gⁿ)
    outputfields = Dict(:U=>output_U, :Φ=>output_Φ, :G=>output_G)

    prefix = "data_rayleigh_benard_regression"
    outputwriter = JLD2OutputWriter(model, outputfields; dir=joinpath(dirname(@__FILE__), "data"),
                                    prefix=prefix, frequency=test_steps, including=[])

    #####
    ##### Initial condition and spinup steps for creating regression test data
    #####

    #=
    @warn ("Generating new data for the Rayleigh-Benard regression test.
           New regression test data generation will fail unless the JLD2
           file $prefix.jld2 is manually deleted.")

    ξ(z) = a * rand() * z * (Lz + z) # noise, damped at the walls
    b₀(x, y, z) = (ξ(z) - z) / Lz
    set!(model, b=b₀)

    time_step!(model, spinup_steps-test_steps, Δt)
    push!(model.output_writers, outputwriter)
    time_step!(model, 2test_steps, Δt)
    =#

    #####
    ##### Regression test
    #####

    # Load initial state
    u₀, v₀, w₀ = get_output_tuple(outputwriter, spinup_steps, :U)
    b₀, c₀ = get_output_tuple(outputwriter, spinup_steps, :Φ)
    Gu₀, Gv₀, Gw₀, Gb₀, Gc₀ = get_output_tuple(outputwriter, spinup_steps, :G)

    model.velocities.u.data.parent .= ArrayType(u₀.parent)
    model.velocities.v.data.parent .= ArrayType(v₀.parent)
    model.velocities.w.data.parent .= ArrayType(w₀.parent)
    model.tracers.b.data.parent    .= ArrayType(b₀.parent)
    model.tracers.c.data.parent    .= ArrayType(c₀.parent)

    model.timestepper.Gⁿ.u.data.parent .= ArrayType(Gu₀.parent)
    model.timestepper.Gⁿ.v.data.parent .= ArrayType(Gv₀.parent)
    model.timestepper.Gⁿ.w.data.parent .= ArrayType(Gw₀.parent)
    model.timestepper.Gⁿ.b.data.parent .= ArrayType(Gb₀.parent)
    model.timestepper.Gⁿ.c.data.parent .= ArrayType(Gc₀.parent)

    model.clock.iteration = spinup_steps
    model.clock.time = spinup_steps * Δt
    length(model.output_writers) > 0 && pop!(model.output_writers)

    # Step the model forward and perform the regression test
    time_step!(model, test_steps, Δt; init_with_euler=false)

    u₁, v₁, w₁ = get_output_tuple(outputwriter, spinup_steps+test_steps, :U)
    b₁, c₁ = get_output_tuple(outputwriter, spinup_steps+test_steps, :Φ)

    field_names = ["u", "v", "w", "b", "c"]
    fields = [model.velocities.u.data.parent, model.velocities.v.data.parent,
              model.velocities.w.data.parent, model.tracers.b.data.parent,
              model.tracers.c.data.parent]
    fields_correct = [u₁.parent, v₁.parent, w₁.parent,
                      b₁.parent, c₁.parent]
    summarize_regression_test(field_names, fields, fields_correct)

    # Now test that the model state matches the regression output.
    @test all(Array(model.velocities.u.data.parent) .≈ u₁.parent)
    @test all(Array(model.velocities.v.data.parent) .≈ v₁.parent)
    @test all(Array(model.velocities.w.data.parent) .≈ w₁.parent)
    @test all(Array(model.tracers.b.data.parent)    .≈ b₁.parent)
    @test all(Array(model.tracers.c.data.parent)    .≈ c₁.parent)

    return nothing
end


text:
using Pkg
Pkg.activate(pwd())

# # Linear regression with sparse constraints
# ## Ridge regression

using LinearAlgebra
using Random
using Plots

n = 10000
m = 1000

Random.seed!(666)

#+

X = randn(n, m)
y = 10*X[:,1] + X[:,2] + randn(n)

# ### Exercise:
# Implement the methods for the `ridge_reg` function. Verify that the result in the same
# result.
#
# **Hints:**
# - The eigendecomposition can be found by `eigen(A)` or `eigen(A).values`.
# - The identity matrix is implemented by `I` in the `LinearAlgebra` package.
#
# ---
# ### Solution:

# ---

using BenchmarkTools

@btime ridge_reg(X, y, 10);
@btime ridge_reg(X, y, 10, Q, Q_inv, λ);

#+

μs = range(0, 1000; length=50)
ws = hcat(ridge_reg.(Ref(X), Ref(y), μs, Ref(Q), Ref(Q_inv), Ref(λ))...)

plot(μs, abs.(ws');
    label="",
    yscale=:log10,
    xlabel="mu",
    ylabel="weights: log scale",
)

# ## Lasso

S(x, η) = max(x-η, 0) - max(-x-η, 0)

function lasso(X, y, μ, Q, Q_inv, λ;
        max_iter = 100,
        ρ = 1e3,
        w = zeros(size(X,2)),
        u = zeros(size(X,2)),
        z = zeros(size(X,2)),
    )

    for i in 1:max_iter
        w = Q * ( (Diagonal(1 ./ (λ .+ ρ)) * ( Q_inv * (X'*y + ρ*(z-u)))))
        z = S.(w + u, μ / ρ)
        u = u + w - z
    end
    return w, u, z
end

#+

ws = zeros(size(X,2), length(μs))

for (i, μ) in enumerate(μs)
    global w, u, z
    w, u, z = i > 1 ? lasso(X, y, μ, Q, Q_inv, λ; w, u, z) : lasso(X, y, μ, Q, Q_inv, λ)
    ws[:,i] = w
end

plot(μs, abs.(ws');
    label="",
    yscale=:log10,
    xlabel="mu",
    ylabel="weights: log scale",
)


text:
####
#### Chebyshev on [-1,1]
####

@testset "Chebyshev" begin
    @test_throws ArgumentError Chebyshev(InteriorGrid(), 0)
    @test_throws ArgumentError Chebyshev(EndpointGrid(), 1)

    for grid_kind in (InteriorGrid(), EndpointGrid())
        for N in (grid_kind ≡ InteriorGrid() ? 1 : 2):10

            basis = Chebyshev(grid_kind, N)
            @test is_function_basis(basis)
            @test is_function_basis(typeof(basis))
            @test dimension(basis) == N

            for i in 1:100
                x = rand_in_domain(i, -1, 1)
                bx = @inferred basis_at(basis, x)

                @test length(bx) == N
                @test eltype(bx) == Float64
                @test collect(bx) ≈ [chebyshev_cos(x, i) for i in 1:N]

                θ = rand(N)
                @test linear_combination(basis, θ, x) ≈
                    sum(chebyshev_cos(x, i) * θ for (i,θ) in enumerate(θ))
                @test derivative(linear_combination(basis, θ), x) ≈
                    sum(chebyshev_cos_deriv(x, i) * θ for (i,θ) in enumerate(θ))
            end

            g = @inferred grid(basis)
            @test length(g) == N
            if grid_kind ≡ InteriorGrid()
                @test all(x -> isapprox(chebyshev_cos(x, N + 1), 0, atol = 1e-14), g)
            else
                @test all(x -> isapprox(chebyshev_cos_deriv(x, N), 0, atol = 1e-13),
                          g[2:(end-1)])
                @test g[1] == -1
                @test g[end] == 1
            end
        end
    end
end


text:
module IPFT

using Distributions
using Random
using Statistics
using StatsBase: sample, mean, weights
import StatsBase: entropy
using LinearAlgebra

# using Parameters
using CPUTime

using POMDPs
using POMDPSimulators
using POMDPModelTools
using MCTS
using MCTS: convert_estimator
using BasicPOMCP
using ParticleFilters
using QMDP

using Test
# using InteractiveUtils

include("information_measures/InformationMeasures.jl")

export
    AbstractMultiObjectiveCriterion,
    MultiObjective,
    MultiObjectiveUCB,
    q_value
include("util/action_selection.jl")

export
    InformationMCTSPlanner,
    InformationBelief,
    pdf,
    resample
include("util/util.jl")

export
    IPFTPlanner,
    IPFTSolver,
    IPFTStateNode,
    IPFTree,
    action_info,
    InformationRewardBeliefMDP
include("core/ipft_types.jl")
include("core/ipft.jl")
include("core/information_mdp.jl")

export
    InformationFilter,
    update,
    initialize_belief
include("util/information_filter.jl")

include("util/rollout.jl")

export
    IBMDPSolver
include("util/ibmdp_solver.jl")

end # module


text:
using Distributions
using Gadfly
using PMTK4Stats

"""
Plot CI and HPD for beta posterior
"""
function drawBetaHPD()
    a = 3
    b = 9
    α = 0.05
    β = Beta(a, b)

    CI = quantile(β, [α/2, 1-α/2])

    xs = linspace(0.001, 0.999, 40)
    ps = pdf(β, xs)

    icdf(p::Real) = quantile(β, p)
    HPD = hdiFromIcdf(icdf)

    ints = hcat(CI, HPD)
    names = ["CI"; "HPD"]

    function plot_ci(int, name)
        l, u = int
        pl = pdf(β, l)
        pu = pdf(β, u)
        vert_up_l = layer(x=[l, l], y=[0, pl], Geom.line, Theme(default_color=colorant"blue", line_width=3px))
        hor_l     = layer(x=[l, u], y=[pl, pu], Geom.line, Theme(default_color=colorant"blue", line_width=3px))
        vert_dn_l = layer(x=[u, u], y=[pu, 0], Geom.line, Theme(default_color=colorant"blue", line_width=3px))
        pdf_l     = layer(x=xs, y=ps, Geom.line, Theme(default_color=colorant"black", line_width=3px))

        pp = plot(vert_up_l, hor_l, vert_dn_l, pdf_l)
    end

    [plot_ci(slicedim(ints, 2, i), names[i]) for i=1:length(names)]
end


text:
import Base: length,show,push!
struct StateSpace{N,T}

    roots::AbstractVector{<:AbstractVector{<:AbstractVector}}
    parameter::AbstractRange
    targets::AbstractVector

end

"""
	data = StateSpace( dimension, parameter, targets; nRoots=2, eltype=Float64 )

Define state space with targets to be used in optimisation

# Positional Arguments
- `dimension` dimensionality of state space `u`
- `parameter` one dimensional bifurcation parameter grid `p`
- `targets` vector of target locations

# Keyword Arguments
- `nRoots` number of roots to continue solutions from
- `eltype` numeric type for vector elements

# Struct Fields
- `roots` vector of roots `F(u.p)=0` to continue solutions from
- `parameter` one dimensional bifurcation parameter grid `p`
- `targets` vector of target locations
"""
function StateSpace(dimension::Integer,parameter::AbstractRange,targets::AbstractVector; nRoots::Integer=2, eltype::DataType=Float64)
    roots, nTargets = fill( [zero(SizedVector{dimension,eltype})] ,nRoots ), length(targets)
    return StateSpace{dimension,eltype}( SizedVector{nRoots}(roots), StepRangeLen{eltype}(parameter), SVector{nTargets,eltype}(targets) )
end

"""
	branch = Branch{V,T}()

Initialises vector of named tuples that contain the following fields

# Fields
- `z` steady state solutions `(u,p)` along branch
- `λ` vector of eigenvalues
- `ds` arclength steps sizes between continuation points
- `bif` boolean telling us if point is a bifurcation
"""
BranchPoint{V,T} = NamedTuple{(:z,:λ,:ds,:bif),Tuple{BorderedArray{V,T},Vector{Complex{T}},T,Bool}}
Branch{V,T} = Vector{BranchPoint{V,T}}
dim(branch::Branch) = length(branch) > 0 ? length(first(branch).z)-1 : 0

function push!(branch::Branch,state::ContState)
    z,∂z = copy(state.z_old), [state.tau.u;state.tau.p]
    push!(branch, ( z=z, λ=state.eigvals, ds=norm(∂z)*abs(state.ds), bif=detectBifucation(state) ))
end

import Base: zero,isapprox,unique,+ ################################## methods for BorderedArray
import LinearAlgebra: norm

+(x::BorderedArray, y::BorderedArray) = BorderedArray(x.u+y.u,x.p+y.p)
norm(x::Branch, y::Branch) = mapreduce( (x,y)->norm(x.z.u-y.z.u)+norm(x.z.p-y.z.p),+,x,y)/min(length(x),length(y))

zero(::Type{BorderedArray{T,U}}; ϵ::Bool=true) where {T<:Union{Number,StaticArray},U<:Union{Number,StaticArray}} = BorderedArray(zero(T).+ϵ*eps(),zero(U).+ϵ*eps())
isapprox( x::BorderedArray, y::BorderedArray; kwargs... ) = isapprox( [x.u;x.p], [y.u;y.p] ; kwargs...)
function unique(X::AbstractVector{T}; kwargs...) where T<:BorderedArray
    z = T[]
    for x ∈ X
        if all( zi -> ~isapprox( x, zi; kwargs... ), z )
            push!(z,x)
        end
    end
    return z
end

function unique(X::AbstractVector{T}; kwargs...) where T<:Branch
    branches = T[]
    for branch ∈ X
        if all( branchᵢ -> ~isapprox( norm(branch,branchᵢ), 0; kwargs... ), branches )
            push!(branches,branch)
        end
    end
    return branches
end

function kernel(A::AbstractMatrix; nullity::Int=0)

    if iszero(nullity) # default to SVD
        return nullspace(A)

    else
        return qr(A').Q[:,end+1-nullity:end]
    end
end

#################################################### display methods
show(io::IO, branches::Vector{Branch{V,T}}) where {V,T} = print(io,
    "Vector{Branch}[dim=$(dim(first(branches))) bifurcations=$(length(unique([ s.z for branch ∈ branches for s ∈ branch if s.bif ]; atol=0.1))) branches=$(length(branches)), states=$(sum(branch->length(branch),branches))]")
show(io::IO, M::MIME"text/plain", branches::Vector{Branch{V,T}}) where {V,T} = show(io,branches)
show(io::IO, states::StateSpace{N,T}) where {N,T} = print(io,"StateSpace{$N,$T}(parameters=$(states.parameter),targets=$(states.targets))")


text:
"""
    BufferState(;busy, error, empty)

An object representing the state of the buffer.

# Fields

- `busy=false`: busy if true
- `error=false`: error if true
- `empty=true`: empty if true
"""
mutable struct BufferState
    busy::Bool
    error::Bool
    empty::Bool
end

BufferState(;busy=false, error=false, empty=true) = BufferState(busy, error, empty)

abstract type AbstractParms end

"""
    Parms(; kwargs...) -> Parms

ACT-R parameters with default values. Default values are overwritten with keyword arguments.

# Fields

- `d=0.5`: activation decay
- `τ=0.0`: retrieval threshold
- `s=0.2`: logistic scalar for activation noise.
- `γ=1.6`: maximum associative strength
- `blc=0.0`: base level constant
- `δ=0.0`: mismatch penalty
- `ter=0.0`: a constant for encoding and responding time
- `mmp_fun`: a mismatch penalty function. By default, `mmp_fun` subtracts `δ` from each non-matching slot value
- `sa_fun`: a function for spreading activation which requires arguments for actr and chunk
- `util_mmp_fun`: utility mismatch penalty function applied to each condition
- `lf=1.0:` latency factor parameter
- `u0=0.0`: initial utility value
- `σu=.2`: standard deviation of utility noise
- `δu=1.0`: mismatch penalty parameter for utility
- `bll=false`: base level learning on
- `mmp=false`: mismatch penalty on
- `sa=false`: spreading activatin on
- `noise=false`: noise on
- `mmp_utility=false`: mismatch penalty for procedural memory
- `utility_noise=false`: utility noise for procedural memory
- `misc`: `NamedTuple` of extra parameters
- `filtered:` a list of slots that must absolutely match with mismatch penalty. `isa` and `retrieval` are included
    by default
"""
@concrete mutable struct Parms <: AbstractParms
    d
    τ
    s
    γ
    δ
    blc
    ter
    mmp_fun
    sa_fun
    util_mmp_fun
    lf
    τ′
    u0
    σu
    δu
    bll::Bool
    mmp::Bool
    sa::Bool
    noise::Bool
    mmp_utility
    utility_noise
    misc
end

function Parms(;
    d = .5,
    τ = 0.0,
    s = .3,
    γ = 0.0,
    δ = 0.0,
    blc = 0.0,
    ter = 0.0,
    mmp_fun = default_penalty,
    sa_fun = spreading_activation!,
    util_mmp_fun = utility_match,
    lf = 1.0,
    τ′ = τ,
    u0 = 0.0,
    σu = .2,
    δu = 1.0,
    bll = false,
    mmp = false,
    sa = false,
    noise = false,
    mmp_utility = false,
    utility_noise = false,
    kwargs...
    )

    Parms(
        d,
        τ,
        s,
        γ,
        δ,
        blc,
        ter,
        mmp_fun,
        sa_fun,
        util_mmp_fun,
        lf,
        τ′,
        u0,
        σu,
        δu,
        bll,
        mmp,
        sa,
        noise,
        mmp_utility,
        utility_noise,
        NamedTuple(kwargs)
    )
end

function Base.show(io::IO, ::MIME"text/plain", parms::Parms)
    values = [getfield(parms, f) for f in fieldnames(Parms)]
    values = map(x->typeof(x)== Bool ? string(x) : x, values)
    return pretty_table(
        io,
        values;
        title="Model Parameters",
        row_name_column_title="Parameter",
        compact_printing=false,
        header=["Value"],
        row_name_alignment=:l,
        row_names=[fieldnames(Parms)...],
        formatters=ft_printf("%5.2f"),
        alignment=:l,
    )
end

"""
    Chunk(; kwargs...) -> Chunk

A declarative memory chunk.

# Fields

- `N=1.0`: number of uses
- `L=1.0`: lifetime of chunk
- `time_created=0.0`: time at which the chunk was created
- `k=1`: number of chunks in recent set (k=1 is sufficient)
- `act=0.0`: total activation
- `act_blc=0.0`: base level constant component of activation
- `act_bll=0.0`: base level learning component of activation
- `act_pm=0.0`: partial matching component of activation
- `act_sa=0.0`: spreading activation component of activation
- `act_noise=0.0`: noise component of activation
- `slots`: chunk slot-value pairs
- `reps=0`: number of identical chunks. This can be used in simple cases to speed up the code.
- `recent=[0.0]`: time stamps for k recent retrievals
- `lags=Float64[]`: lags for recent retrievals (L - recent)
- `bl=0.0`: baselevel constant added to chunks activation
"""
mutable struct Chunk{T1,T2}
  N::Int
  L::Float64
  time_created::Float64
  k::Int
  act_mean::T2
  act::T2
  act_blc::T2
  act_bll::T2
  act_pm::T2
  act_sa::T2
  act_noise::T2
  slots::T1
  reps::Int64
  recent::Array{Float64,1}
  lags::Array{Float64,1}
  bl::T2
end

function Chunk(;
    N=1,
    L=1.0,
    time_created=0.0,
    k=1,
    act=0.0,
    recent=[0.0],
    reps=0,
    lags=Float64[],
    bl=zero(typeof(act)),
    slots...
    )
    T = typeof(act)
    act_mean = zero(T)
    act_pm = zero(T)
    act_blc = zero(T)
    act_bll = zero(T)
    act_noise = zero(T)
    act_sa = zero(T)
    return Chunk(N, L, time_created, k, act_mean, act, act_blc, act_bll, act_pm, act_sa, act_noise,
        NamedTuple(slots), reps, recent, lags, bl)
end

"""
    Chunk(; kwargs...) -> Chunk

A declarative memory chunk with dynamic slot-value pairs.

# Fields

- `dynamic::Bool`: slot-value pairs are mutable if true
- `N=1.0`: number of uses
- `L=1.0`: lifetime of chunk
- `time_created=0.0`: time at which the chunk was created
- `k=1`: number of chunks in recent set (k=1 is sufficient)
- `act_mean`: mean activation computed as `act` - `act_noise`
- `act=0.0`: total activation computed as `act_mean` + `act_noise`
- `act_blc=0.0`: base level constant component of activation
- `act_bll=0.0`: base level learning component of activation
- `act_pm=0.0`: partial matching component of activation
- `act_sa=0.0`: spreading activation component of activation
- `act_noise=0.0`: noise component of activation
- `slots`: chunk slot-value pairs
- `reps=0`: number of identical chunks. This can be used in simple cases to speed up the code.
- `recent=[0.0]`: time stamps for k recent retrievals
- `lags=Float64[]`: lags for recent retrievals (L - recent)
- `bl=0.0`: baselevel constant added to chunks activation
"""
function Chunk(dynamic::Bool;
    N=1,
    L=1.0,
    time_created=0.0,
    k=1,
    act=0.0,
    recent=[0.0],
    reps=0,
    lags=Float64[],
    bl=zero(typeof(act)),
    slots...
    )
    T = typeof(act)
    act_mean = zero(T)
    act_pm = zero(T)
    act_blc = zero(T)
    act_bll = zero(T)
    act_noise = zero(T)
    act_sa = zero(T)
    slots = Dict(k=>v for (k,v) in pairs(slots))
    return Chunk(N, L, time_created, k, act_mean, act, act_blc, act_bll, act_pm, act_sa, act_noise,
        slots, reps, recent, lags, bl)
end

Broadcast.broadcastable(x::Chunk) = Ref(x)

const chunk_fields = (:slots,:N,:L,:time_created,:recent,:act_mean,:act,:act_blc,:bl,:act_bll,:act_pm,:act_noise)

function chunk_values(chunk)
    values = [getfield(chunk, f) for f in chunk_fields]
    return map(x->typeof(x)== Bool ? string(x) : x, values)
end

function Base.show(io::IO, ::MIME"text/plain", chunk::Chunk)
    values = chunk_values(chunk)
    return pretty_table(
        io,
        values;
        title="Chunk",
        row_name_column_title="Property",
        compact_printing=false,
        header=["Value"],
        row_name_alignment=:l,
        row_names=[chunk_fields...],
        formatters=ft_printf("%5.2f"),
        alignment=:l,
    )
end

function Base.show(io::IO, ::MIME"text/plain", chunks::Vector{<:Chunk})
    table = [chunk_values(chunk) for chunk in chunks]
    table = hcat(table...)
    table = permutedims(table)
    table = isempty(chunks) ? fill(Missing, 1, length(chunk_fields)) : table
    return pretty_table(
        io,
        table;
        title="Chunks",
        # row_name_column_title="Parameter",
        compact_printing=false,
        header=[chunk_fields...],
        row_name_alignment=:l,
        formatters=ft_printf("%5.2f"),
        alignment=:l,
    )
end

abstract type Mod end

"""
    Declarative(;memory=, filtered=(:isa,:retrieved))

Declarative memory module

# Fields

- `memory=Chunk[]`: array of chunks
- `filtered`: slots that must match exactly even when partial matching is on. By default,
    `filtered=(:isa,:retrieved)`
- `buffer`: an array containing one chunk
- `state`: buffer state
"""
mutable struct Declarative{T1,T2} <: Mod
    memory::Array{T1,1}
    filtered::T2
    buffer::Array{T1,1}
    state::BufferState
end

function Declarative(;memory=Chunk[], filtered=(:isa,:retrieved))
    state = BufferState()
    return  Declarative(memory, filtered, typeof(memory)(undef,1), state)
end

"""
    default_penalty(actr, chunk; request...)

A default function for mismatch penalty. Subtracts δ if
slot does not exist or slot value does not match

# Arguments

- `actr`: an ACTR model object
- `chunk`: memory chunk object

# Keywords

- `request`: a variable size collection of slot-value pairs for the retrieval request
"""
function default_penalty(actr, chunk; request...)
    slots = chunk.slots
    p = 0.0; δ = actr.parms.δ
    for (k,v) in request
        if !(k ∈ keys(slots)) || (slots[k] != v)
            p += δ
        end
     end
     return p
end

Broadcast.broadcastable(x::Declarative) = Ref(x)

"""
    Imaginal(;buffer=Chunk[], ω=1.0, denoms=Int64[])

Imaginal Module.

# Fields

- `buffer`: an array holding up to one chunk
- `state`: buffer state
- `ω=1.0`: fan weight. Default is 1.
- `denoms=Int64[]`: cached value for the denominator of the fan calculation
"""
mutable struct Imaginal{T1,T2} <: Mod
    buffer::Array{T1,1}
    state::BufferState
    ω::T2
    denoms::Vector{Int64}
end

function Imaginal(;buffer=Chunk[], ω=1.0, denoms=Int64[])
    state = BufferState()
    Imaginal(buffer, state, ω, denoms)
end

Imaginal(chunk::Chunk, state, ω, denoms) = Imaginal([chunk], state, ω, denoms)
Imaginal(T::DataType, state, ω, denoms) = Imaginal(T(undef,1), state, ω, denoms)

"""
    Visual(;chunk=Chunk())

Visual Module.

# Fields

- `buffer`: an array holding up to one chunk
- `state`: buffer state
- `focus`: coordinates of visual attention
"""
mutable struct Visual{T1} <: Mod
    buffer::Array{T1,1}
    state::BufferState
    focus::Vector{Float64}
end

Visual(;buffer=Chunk[]) = Visual(buffer, BufferState(), fill(0.0,2))
Visual(chunk::Chunk, state, focus) = Visual([chunk], state, focus)
Visual(T::DataType, state, focus) = Visual(T(undef,1), state, focus)

abstract type AbstractVisualObject end
"""
    VisualObject(;x=300.0, y=300.0, color=:black, text="", shape=:_, width=30.0, height=30.0)

A visible object in a task.

# Fields

- `x`: x coordinate of visual object. Default 0.
- `y`: y coordinate of visual object. Default 0.
- `color`: object color
- `shape`: object shape
- `text`: object text
- `width`: object width
- `height`: object height
"""
mutable struct VisualObject <: AbstractVisualObject
    x::Float64
    y::Float64
    color::Symbol
    shape::Symbol
    text::String
    width::Float64
    height::Float64
end

function VisualObject(;x=300.0, y=300.0, color=:black, text="", shape=:_, width=30.0, height=30.0)
    return VisualObject(x, y, color, shape, text, width, height)
end

"""
**VisualLocation**

Visual Location Module.
- `buffer`: an array holding up to one chunk
- `state`: buffer state

Constructor
````julia
VisualLocation(;buffer=Chunk[])

VisualLocation(chunk::Chunk, state)

VisualLocation(T::DataType, state)

VisualLocation(chunks, state)
````
"""
mutable struct VisualLocation{T1} <: Mod
    buffer::Array{T1,1}
    state::BufferState
    iconic_memory::Array{T1,1}
end

function VisualLocation(;buffer=Chunk[])
    VisualLocation(buffer, BufferState())
end

function VisualLocation(chunk::Chunk, state)
    T = typeof(chunk)
     VisualLocation([chunk], state, Vector{T}(undef,1))
end

function VisualLocation(T::DataType, state)
    VisualLocation(T(undef,1), state, T(undef,1))
end

function VisualLocation(chunks, state)
    c_chunks = copy(chunks)
    VisualLocation(chunks, state, c_chunks)
end

"""
    Rule(;utlity=0.0, conditions, action)

A production rule object.

# Fields

- `utility=0.0`: utility of the production rule
- `initial_utility=0.0`: initial utility
- `utility_mean`=0.0: mean utility
- `utility_penalty=0.0`: mismatch penalty term for utility
- `utlity_noise=0.0`: utility noise
- `conditions`: a function for checking conditions
- `action`: a function for performing an action
- `name`: name of production
"""
@concrete mutable struct Rule
    utility
    initial_utility
    utility_mean
    utility_penalty
    utility_noise
    conditions
    action
    can_pm
    name::String
end

"""
**Procedural**

Procedural Memory Module.
- `buffer`: an array holding up to one chunk
- `state`: buffer state

Constructor
````julia
Procedural(;chunk=Chunk())

Procedural(;rules=Rule[], id="")

Procedural(rule::Rule, state, id)

Procedural(T::DataType, state, id)
````
"""
mutable struct Procedural{R} <: Mod
    id::String
    rules::R
    state::BufferState
end

function Procedural(;rules=Rule[], id="")
    Procedural(id, rules, BufferState())
end

function Procedural(rule::Rule, state, id)
    Procedural(id, [rule], state)
end

function Procedural(T::DataType, state, id)
    Procedural(id, T(undef,1), state)
end

function utility_match(actr, condition)
    @error "a method must be defined for utility_match(actr::ACTR, condition)"
end

"""
    Goal(;chunk=Chunk())

Goal Module.

# Fields

- `buffer`: an array holding up to one chunk
- `state`: buffer state
"""
mutable struct Goal{T1} <: Mod
    buffer::Array{T1,1}
    state::BufferState
end

function Goal(;buffer=Chunk[])
    Goal(buffer, BufferState())
end

function Goal(chunk::Chunk, state)
    Goal([chunk], state)
end

function Goal(T::DataType, state)
    Goal(T(undef,1), state)
end

"""
    Motor(;chunk=Chunk())

Motor Module.

# Fields

- `buffer`: an array holding up to one chunk
- `state`: buffer state
- `mouse_position`: x,y coordinates of mouse position on screen
"""
mutable struct Motor{T1} <: Mod
    buffer::Array{T1,1}
    state::BufferState
    mouse_position::Vector{Float64}
end

function Motor(;buffer=Chunk[], mouse_position=[0.0,0.0])
    Motor(buffer, BufferState(), mouse_position)
end

function Motor(chunk::Chunk, state, mouse_position)
    Motor([chunk], state, mouse_position)
end

function Motor(T::DataType, state, mouse_position)
    Motor(T(undef,1), state, mouse_position)
end

mutable struct Scheduler
    time::Float64
end

Scheduler(;time=0.0) = Scheduler(time)

abstract type AbstractACTR end

"""
ACTR(; kwargs...) -> ACTR

ACTR model object

# Fields

- `name="model1"`: model name
- `declarative`: declarative memory module
- `imaginal`: imaginal memory module
- `visual`: visual module
- `goal`: goal module
- `visual_location`: visual location module
- `visicon`: a vector of VisualObjects available in the task
- `parms`: model parameters
- `scheduler`: event scheduler

# Example

````julia
parms = (noise=true, τ=-1.0)
chunks = [Chunk(;animal=:dog,name=:Sigma), Chunk(;animal=:rat,name=:Bonkers)]
declarative = Declarative(;memory=chunks)
actr = ACTR(;declarative, parms...)
````
"""
@concrete mutable struct ACTR <: AbstractACTR
    name
    declarative
    imaginal
    visual
    visual_location
    goal
    procedural
    motor
    visicon
    parms
    scheduler
end

Broadcast.broadcastable(x::ACTR) = Ref(x)

function ACTR(;name="model1", declarative=Declarative(), imaginal=Imaginal(),
    goal = Goal(), scheduler=Scheduler(), visual=nothing, visual_location=nothing,
    procedural=nothing, motor=nothing, visicon=init_visicon(), parms...)
    parms′ = Parms(;parms...)
    ACTR(name, declarative, imaginal, visual, visual_location, goal, procedural, motor, visicon, parms′, scheduler)
end

function init_visicon()
    Dict{String,VisualObject}()
end


text:
function load_docstring(api::AbstractString, recordmodule)
    doc = read(joinpath(@__DIR__, "docs/$api.md"), String)
    return replace(doc, raw"$Recorder" => string(recordmodule))
end

macro define(recordmodule)
    Runtime = @__MODULE__
    at_record_doc = load_docstring("@record", recordmodule)
    enable_doc = load_docstring("enable", recordmodule)
    enable_logging_doc = load_docstring("enable_logging", recordmodule)
    expr = quote
        module $recordmodule
        const Runtime = $Runtime

        struct DefaultDataFunction <: Runtime.AbstractDefaultDataFunction end
        const defaultdata = DefaultDataFunction()

        struct PrefixFunction <: Runtime.AbstractPrefixFunction end
        const prefix = PrefixFunction()

        struct IsRecordingFunction <: Runtime.AbstractIsRecordingFunction end
        const isrecording = IsRecordingFunction()
        @doc $enable_doc enable() = @eval isrecording() = true
        @doc $enable_doc disable() = @eval isrecording() = true

        LOGGER = nothing
        get_logger() = nothing
        function set_logger(logger::Union{Runtime.LoggerLike,Nothing})
            global LOGGER = logger
            if logger === nothing || logger === false
                @eval get_logger() = nothing
            else
                @eval get_logger() = LOGGER::$(Expr(:$, :(typeof(logger))))
            end
        end
        default_logger() = Runtime.DefaultLogger()
        @doc $enable_logging_doc enable_logging() = set_logger(default_logger())
        @doc $enable_logging_doc disable_logging() = set_logger(nothing)

        const RECORDS = Runtime.create_records()
        function __init__()
            Runtime.init_records!(RECORDS)
            Runtime.register_recorder($recordmodule)
            return
        end

        get_records() = RECORDS
        clear() = empty!(get_records())

        struct IsEnabledFunction <: Runtime.AbstractIsEnabledFunction
            isrecording::typeof(isrecording)
            get_logger::typeof(get_logger)
        end
        const isenabled = IsEnabledFunction(isrecording, get_logger)

        struct RecordFunction <: Runtime.AbstractRecordFunction
            defaultdata::DefaultDataFunction
            prefix::PrefixFunction
            isrecording::typeof(isrecording)
            get_logger::typeof(get_logger)
            get_records::typeof(get_records)
        end
        const record =
            RecordFunction(defaultdata, prefix, isrecording, get_logger, get_records)

        @doc $at_record_doc macro record(recordname, exprs...)
            $Runtime.record_impl(
                __source__,
                __module__,
                $recordmodule,
                recordname,
                collect(Any, exprs),
            )
        end

        end # module
    end

    return esc(Expr(:toplevel, expr.args...))
end

const RECORDERS = Module[]
register_recorder(recordmodule::Module) = push!(RECORDERS, recordmodule)

function as_kwarg_expr(x)
    if x isa Symbol
        return Expr(:kw, x, x)
    elseif isexpr(x, :(=), 2)
        return Expr(:kw, x.args...)
    else
        error("invalid expression: ", x)
    end
end

function record_impl(
    __source__::LineNumberNode,
    __module__::Module,
    recordmodule::Module,
    recordname,
    exprs,
)
    if isexpr(get(exprs, lastindex(exprs), nothing), :braces)
        options = map(as_kwarg_expr, exprs[end].args)
        exprs = exprs[begin:end-1]
    else
        options = []
    end
    data_kwargs = map(as_kwarg_expr, exprs)
    data = :((; $(data_kwargs...)))
    location = uuid4()
    quote
        if $recordmodule.isenabled()
            $recordmodule.record(
                $recordname,
                $data;
                location = $(QuoteNode(location)),
                __source__ = $(QuoteNode(__source__)),
                __module__ = $(QuoteNode(__module__)),
                $(options...),
            )
        end
    end |> esc
end

abstract type AbstractRecordFunction <: Function end
abstract type AbstractIsEnabledFunction <: Function end
abstract type AbstractDefaultDataFunction <: Function end
abstract type AbstractPrefixFunction <: Function end
abstract type AbstractIsRecordingFunction <: Function end

(::AbstractDefaultDataFunction)() = Historic.taskdata()
(::AbstractPrefixFunction)() = Historic.taskprefix()
(::AbstractIsRecordingFunction)() = false

function (record::AbstractRecordFunction)(
    name::Symbol,
    data::NamedTuple;
    yield::Bool = true,
    options...,
)
    logger = record.get_logger()
    islogging = logger !== nothing && should_log(options)
    if islogging && !yield && may_yield(logger)
        islogging = false
    end

    if record.isrecording()
        # Heuristics: `defaultdata` is shared across all records and it produces
        # the most "dense" column which is presumably the most informative one.
        # Putting it first so that they come first in the column tables.
        recorddata = (; record.defaultdata()..., data...)
        record!(record.get_records(), name, recorddata; options...)
    end

    if islogging
        # Heuristics: `defaultdata` is probably included in `prefix`. So, put it
        # at the end so that `DefaultLogger` prints the non-default data first.
        logdata = (; data..., record.defaultdata()...)
        logrecord(logger, name, logdata, record.prefix(); options...)
    end
end

(isenabled::AbstractIsEnabledFunction)() =
    isenabled.isrecording() || isenabled.get_logger() !== nothing


text:
using Serialization, DataStructures, Random, LinearAlgebra

export computeLogLikelihood, computeGradient, doStochasticGradientAscent,
       doDPPLearningSparseVectorData, isConvergedLogLikelihood

# Computes the log-likelihood for the low-rank DPP with parameters paramsMatrix
function computeLogLikelihood(paramsMatrix, trainingInstances, numTrainingInstances,
                              numItems, numItemTraits, lambdaVec = 0, alpha = 0,
                              useDual = true)
  itemTraitMatrixInstance = Matrix{Float64}
  lMatrixTrainingInstance = Matrix{Float64}

  # Compute first term of log-likelihood
  sumDetTrainingInstances = 0.0
  for i = 1:numTrainingInstances
    @inbounds itemTraitMatrixInstance = paramsMatrix[trainingInstances[i], :]
    lMatrixTrainingInstance = itemTraitMatrixInstance * itemTraitMatrixInstance'

    detLMatrixTrainingInstance = det(lMatrixTrainingInstance)
    result = 0.0
    if detLMatrixTrainingInstance > 0
      result = log(detLMatrixTrainingInstance)
    end

    sumDetTrainingInstances += result
  end
  # First term of log-likelihood
  firstTerm = sumDetTrainingInstances

  # Compute second term of log-likelihood.  We optionally use the dual
  # representation of a DPP, to avoid computing the determinant of a potentially
  # large L matrix
  tItemTraitMatTimesItemTraitMat = 0
  if useDual
    tItemTraitMatTimesItemTraitMat = paramsMatrix' * paramsMatrix
  else
    tItemTraitMatTimesItemTraitMat = paramsMatrix * paramsMatrix'
  end

  detParams = det(tItemTraitMatTimesItemTraitMat + I)
  if detParams < 0
    # Fix for possible numerical stability issues
    detParams = abs(detParams)
  end
  logDetParams = log(detParams)

  # Second term of log-likelihood
  secondTerm = numTrainingInstances * logDetParams

  # Third term of log-likelihood, which is the regularization term
  thirdTerm = 0.0
  if lambdaVec != 0
    for i = 1:numItems
      thirdTerm += lambdaVec[i] * vecnorm(paramsMatrix[i, :]) ^ 2
    end

    thirdTerm *= 0.5

    if alpha != 0
      thirdTerm *= alpha
    end
  end

  # Full log-likelihood
  logLikelihood = firstTerm - secondTerm - thirdTerm

  return logLikelihood
end

# Computes the gradient of the log-likelihood for the low-rank DPP with
# parameters paramsMatrix. The gradient for each parameter (element of
# paramsMatrix) is computed with respect to that parameter.  We exploit the
# structure of the matrices and derivatives to speed up the computation of the
# gradient.
# If usePinv is set to true, we use the pseudo-inverse (pinv()) to  compute the
# gradient, instead of the inverse (inv()), which may help reduce numerical
# errors when paramsMatrix is rank deficient.  Note that pinv() is slower
# than inv() in Julia.
function computeGradient(paramsMatrix::Matrix{Float64}, trainingInstances::Vector{Vector{Int64}},
                         numTrainingInstances, numItems, numItemTraits,
                         lambdaVec::Vector{Float64} = zeros(numItems), alpha = 0,
                         paramsMatrixRowIdicesToTrainingInstanceRowIndices::Matrix{Int64} =
                         fill(0, numTrainingInstances, numItems) ;
                         usePinv = false)
  sumTraceTrainingInstances = 0.0
  paramsMatrixNumRows = numItems
  paramsMatrixNumCols = numItemTraits
  paramsMatrixGradient = zeros(Float64, numItems, numItemTraits)
  itemTraitMatrixInstance = Matrix{Float64}
  lMatrixTrainingInstance = Matrix{Float64}
  lMatrixTrainingInstanceInverse = Matrix{Float64}
  traceFirstTerm = 0.0
  sumTraceFirstTerm = 0.0
  traceSecondTerm = 0.0

  # Precompute items used in first term of gradient:
  # itemTraitMatrixInstance, lMatrixTrainingInstanceInverse, and
  # numTrainingInstanceItems for each training instance
  numTrainingInstances = Int(numTrainingInstances)
  itemTraitMatrixTrainingInstanceVec = Array{Matrix{Float64}}(undef, numTrainingInstances)
  lMatrixTrainingInstanceInverseVec = Array{Matrix{Float64}}(undef, numTrainingInstances)
  numTrainingInstanceItemsVec = Array{Int}(undef, numTrainingInstances)
  trainingInstanceRowIndex = 0
  numTrainingInstanceItems = 0
  itemNotPresentInTrainingInstance = false
  if paramsMatrixRowIdicesToTrainingInstanceRowIndices == fill(0, numTrainingInstances, numItems)
    paramsMatrixRowIdicesToTrainingInstanceRowIndices =
      buildMapTrainingInstanceRowIndices(trainingInstances, numTrainingInstances, numItems)
  end
  for trainingInstanceIndex = 1:numTrainingInstances
    @inbounds itemTraitMatrixInstance = paramsMatrix[trainingInstances[trainingInstanceIndex], :]
    @inbounds itemTraitMatrixTrainingInstanceVec[trainingInstanceIndex] = itemTraitMatrixInstance

    lMatrixTrainingInstance = itemTraitMatrixInstance * itemTraitMatrixInstance'
    if usePinv
      @inbounds lMatrixTrainingInstanceInverseVec[trainingInstanceIndex] = pinv(lMatrixTrainingInstance)
    else
      @inbounds lMatrixTrainingInstanceInverseVec[trainingInstanceIndex] = inv(lMatrixTrainingInstance)
    end

    @inbounds numTrainingInstanceItemsVec[trainingInstanceIndex] = size(itemTraitMatrixInstance, 1)
  end

  # Precompute items used in second term of gradient
  tItemTraitMatTimesItemTraitMat = paramsMatrix' * paramsMatrix
  numRowsTItemTraitMatTimesItemTraitMat = size(tItemTraitMatTimesItemTraitMat, 1)
  dualMat = zeros(numRowsTItemTraitMatTimesItemTraitMat, numRowsTItemTraitMatTimesItemTraitMat)
  if usePinv
    dualMat = paramsMatrix *
      pinv(I + tItemTraitMatTimesItemTraitMat) * paramsMatrix'
  else
    dualMat = paramsMatrix *
      inv(I + tItemTraitMatTimesItemTraitMat) * paramsMatrix'
  end

  identMinusDualMat = I - dualMat

  # Iterate over each element of currParamsMatrix to compute gradient for each element
  for paramsMatrixColIndex = 1:numItemTraits
    for paramsMatrixRowIndex = 1:numItems
      sumTraceTrainingInstances = 0.0
      for trainingInstanceIndex = 1:numTrainingInstances
        @inbounds numTrainingInstanceItems = numTrainingInstanceItemsVec[trainingInstanceIndex]

        @inbounds itemTraitMatrixInstance = itemTraitMatrixTrainingInstanceVec[trainingInstanceIndex]
        @inbounds lMatrixTrainingInstanceInverse = lMatrixTrainingInstanceInverseVec[trainingInstanceIndex]

        # Map paramsMatrixRowIndex to the proper row index in lMatrixTrainingInstanceInverse
        @inbounds trainingInstanceRowIndex = paramsMatrixRowIdicesToTrainingInstanceRowIndices[trainingInstanceIndex, paramsMatrixRowIndex]
        if trainingInstanceRowIndex != 0
          itemNotPresentInTrainingInstance = false
        else
          itemNotPresentInTrainingInstance = true
        end

        # Compute first term of gradient
        traceFirstTerm = 0.0
        if itemNotPresentInTrainingInstance
          # The lMatrixTrainingInstance derivative goes to 0 when we encounter an item
          # not present in the training instance
          traceFirstTerm = 0.0
        else
          sumTraceFirstTerm = 0.0
          @simd for i = 1:numTrainingInstanceItems
            @inbounds sumTraceFirstTerm += lMatrixTrainingInstanceInverse[i, trainingInstanceRowIndex] *
              itemTraitMatrixInstance[i, paramsMatrixColIndex]
          end

          @inbounds @views traceFirstTerm = dot(lMatrixTrainingInstanceInverse[trainingInstanceRowIndex, :],
            itemTraitMatrixInstance[:, paramsMatrixColIndex]) + sumTraceFirstTerm
        end

        sumTraceTrainingInstances += traceFirstTerm
      end
      # First term of gradient
      firstTerm = sumTraceTrainingInstances

      # Compute second term of gradient
      traceSecondTerm = 0.0
      sumTraceSecondTerm = 0.0
      @simd for i = 1:paramsMatrixNumRows
        @inbounds sumTraceSecondTerm += identMinusDualMat[i, paramsMatrixRowIndex] *
          paramsMatrix[i, paramsMatrixColIndex]
      end

      @inbounds @views traceSecondTerm = dot(identMinusDualMat[paramsMatrixRowIndex, :], paramsMatrix[:, paramsMatrixColIndex]) + sumTraceSecondTerm

      # Second term of gradient
      secondTerm = numTrainingInstances * traceSecondTerm

      # Compute the third term of the gradient, which is the regularization term
      thirdTerm = 0.0
      if lambdaVec != 0
        thirdTerm = lambdaVec[paramsMatrixRowIndex] *
          paramsMatrix[paramsMatrixRowIndex, paramsMatrixColIndex]

        if alpha != 0
          thirdTerm *= alpha
        end
      end

      # Full gradient
      @inbounds paramsMatrixGradient[paramsMatrixRowIndex, paramsMatrixColIndex] =
        firstTerm - secondTerm - thirdTerm
    end
  end

  return paramsMatrixGradient
end

# Performs stochastic gradient ascent to learn low-rank DPP kernel parameters.
# If usePinv is set to true, we use the pseudo-inverse (pinv()) to  compute the
# gradient, instead of the inverse (inv()), which may help reduce numerical
# errors when paramsMatrix is rank deficient.  Note that pinv() is slower
# than inv() in Julia.
function doStochasticGradientAscent(trainingInstances, numTrainingInstances, numItems,
                          numItemTraits, testInstances, numTestInstances, lambdaVec, alpha,
                          validationInstances = fill(Vector{Int}(), 0), numValidationInstances = 0,
                          initialParamsMatrix = rand(numItems, numItemTraits) + ones(numItems, numItemTraits),
                          usePinv = false; verbose = true, maxIters = Inf, minibatchSize = 1000)
  gradient = zeros(numItems, numItemTraits)
  paramsMatrixPrev = zeros(numItems, numItemTraits)

  paramsMatrix = initialParamsMatrix

  epsFixed = 0.5e-2
  epsInitialDecay = 1.0e-5 # For Amazon registry apparel dataset
  eps = epsFixed

  # Number of iterations for which eps is kept fixed
  numIterationsFixedEps = 50

  betaMomentum = 0.95

  numIterationsCompleted = 0

  delta = zeros(numItems, numItemTraits)

  # Number of training instances to process per minibatch
  # minibatchSize = 1000

  currTrainingInstanceIndex = 1

  shuffle!(trainingInstances)

  # Set aside one percent of trainingInstances for use as a validation set, for
  # assessing convergence, if a validation set is not explicitly provided
  avgValidationLogLikelihood = 0
  if numValidationInstances == 0
    validationSetSizePercent = 0.01
    numValidationInstances = convert(Int, round(numTrainingInstances * validationSetSizePercent));
    validationInstances = fill(Vector{Int}(), numValidationInstances)
    validationInstances = trainingInstances[1:numValidationInstances]
    trainingInstances = trainingInstances[(numValidationInstances + 1):numTrainingInstances]
    numTrainingInstances = numTrainingInstances - numValidationInstances
  end

  startTime = time()

  validationLogLike = computeLogLikelihood(paramsMatrix, validationInstances,
    numValidationInstances, numItems, 0)

  # Run stochastic gradient ascent until convergence, or until maxIters
  # iterations have been completed
  while numIterationsCompleted < maxIters
    # Get current minibatch
    numTrainingInstancesInMinibatch = minibatchSize
    if currTrainingInstanceIndex + minibatchSize > numTrainingInstances
      numTrainingInstancesInMinibatch = numTrainingInstances - currTrainingInstanceIndex
    end
    minibatchTrainingInstances =
      trainingInstances[currTrainingInstanceIndex:(currTrainingInstanceIndex + numTrainingInstancesInMinibatch)]

    paramsMatrixPrev = paramsMatrix

    if verbose
      @time gradient = computeGradient(paramsMatrix + betaMomentum * delta,
        minibatchTrainingInstances, numTrainingInstancesInMinibatch, numItems,
        numItemTraits, lambdaVec, alpha)
    else
      gradient = computeGradient(paramsMatrix + betaMomentum * delta,
        minibatchTrainingInstances, numTrainingInstancesInMinibatch, numItems,
        numItemTraits, lambdaVec, alpha)
    end

    # Use momentum when computing the update
    delta = betaMomentum * delta + (1.0 - betaMomentum) * eps * gradient
    paramsMatrix = paramsMatrix + delta

    if verbose
      @time avgTrainingLogLikelihood = computeLogLikelihood(paramsMatrix,
        trainingInstances, numTrainingInstances,
        numItems, numItemTraits) / numTrainingInstances
      @time avgValidationLogLikelihood = computeLogLikelihood(paramsMatrix,
        validationInstances, numValidationInstances, numItems,
        numItemTraits) / numValidationInstances
      @time avgTestLogLikelihood = computeLogLikelihood(paramsMatrix,
        testInstances, numTestInstances, numItems, numItemTraits) / numTestInstances

      println("avgTrainingLogLikelihood: $avgTrainingLogLikelihood")
      println("avgValidationLogLikelihood: $avgValidationLogLikelihood")
      println("avgTestLogLikelihood: $avgTestLogLikelihood")
    end

    numIterationsCompleted += 1

    if verbose
      if numIterationsCompleted % 1 == 0
        println("alpha: $alpha")
        println("Completed $numIterationsCompleted stochastic gradient ascent iterations")
      end
    end

    if numIterationsCompleted % 100 == 0
      # Save paramsMatrix to disk, so that we can reuse it later
      # save("learnedDPPParamsMatrix-k$numItemTraits-lambdaPop$alpha-$numIterationsCompleted.jld", "learnedParamsMatrix", paramsMatrix)
    end

    prevValidationLogLike = validationLogLike
    validationLogLike = computeLogLikelihood(paramsMatrix, validationInstances,
      numValidationInstances, numItems, 0)
    if isConvergedLogLikelihood(prevValidationLogLike, validationLogLike,
      1.0e-5, "validation", verbose = verbose)
      break
    end

    # Anneal (gradually lower) the learning rate
    if numIterationsCompleted >= numIterationsFixedEps
      betaMomentum = 0.0
      eps = epsInitialDecay / (1 + numIterationsCompleted / numIterationsFixedEps)

      if verbose
        println("Reduced eps: $eps")
      end
    end

    currTrainingInstanceIndex += numTrainingInstancesInMinibatch + 1
    if currTrainingInstanceIndex > numTrainingInstances
      # We've processed all of the training instances, so start processing the
      # training instances again from the beginning
      currTrainingInstanceIndex = 1
      shuffle!(trainingInstances)
    end
  end

  endTime = time() - startTime

  return paramsMatrix, avgValidationLogLikelihood, endTime
end

# Builds a matrix containing a map of paramsMatrixRowIndex to the proper row index in
# lMatrixTrainingInstanceInverse.  Used for computing gradient.
function buildMapTrainingInstanceRowIndices(trainingInstances, numTrainingInstances, numItems)
  numTrainingInstancesProcessed = 0

  paramsMatrixRowIdicesToTrainingInstanceRowIndicesMat = fill(0, numTrainingInstances, numItems)

  for trainingInstanceIndex = 1:numTrainingInstances
    trainingInstanceItems = trainingInstances[trainingInstanceIndex]

    for trainingInstanceRowIndex = 1:length(trainingInstanceItems)
      paramsMatrixRowIndex = trainingInstanceItems[trainingInstanceRowIndex]

      paramsMatrixRowIdicesToTrainingInstanceRowIndicesMat[trainingInstanceIndex, paramsMatrixRowIndex] = trainingInstanceRowIndex
    end
  end

  return paramsMatrixRowIdicesToTrainingInstanceRowIndicesMat
end

# Determine if we have reached convergence, based on the log likelihood
# of the V (parameter) matrix.
function isConvergedLogLikelihood(prevLogLike, newLogLike, eps, instanceTypeName;
                                  verbose = true)
  # eps = 1.0e-7
  relativeChange = abs(newLogLike - prevLogLike) / abs(prevLogLike)
  relativeChangeSign = ""
  if newLogLike < prevLogLike
    relativeChangeSign = "-"
  end

  if verbose
    println("\t relativeChange in $instanceTypeName log likelihood: $relativeChangeSign$relativeChange")
  end

  if relativeChange <= eps
    if verbose
      println("relativeChange in $instanceTypeName log likelihood: $relativeChangeSign$relativeChange, converged")
    end

    return true
  else
    return false
  end
end

# Performs low-rank DPP model learning, with item-popularity regularization, on
# a dataset in sparse vector format.
function doDPPLearningSparseVectorData(trainingBasketsDictFileName, trainingBasketsDictObjectName,
                                       testBasketsDictFileName, testBasketsDictObjectName,
                                       learnedModelOutputDirName, numItemTraits,
                                       alpha)
  println("Starting doDPPLearningSparseVectorData()")

  Random.seed!(1234)

  # cd("$(homedir())\\Belgian-retail-supermarket")

  # Load training data
  trainingUsersBasketsDict = open(deserialize, trainingBasketsDictFileName);
  println("Loaded $trainingBasketsDictFileName")

  # Load test data
  testUsersBasketsDict = open(deserialize, testBasketsDictFileName);
  println("Loaded $testBasketsDictFileName")

  # Build set of training instances
  numTrainingInstances = length(collect(keys(trainingUsersBasketsDict)))
  trainingInstances = fill(Vector{Int}(), numTrainingInstances)
  trainingInstanceIndex = 1
  numItems = 0
  for trainingInstanceUserId in collect(keys(trainingUsersBasketsDict))
    trainingInstanceItems = deepcopy(trainingUsersBasketsDict[trainingInstanceUserId].basketItems)
    heldOutItem = trainingUsersBasketsDict[trainingInstanceUserId].heldOutItem
    push!(trainingInstanceItems, heldOutItem)

    # Remove duplicates from trainingInstanceItems
    trainingInstanceItems = collect(Set(trainingInstanceItems))

    trainingInstances[trainingInstanceIndex] = trainingInstanceItems

    trainingInstanceIndex += 1

    numItems = trainingUsersBasketsDict[trainingInstanceUserId].numItemsInCatalog
  end

  # Build set of test instances
  numTestInstances = length(collect(keys(testUsersBasketsDict)))
  testInstances = fill(Vector{Int}(), numTestInstances)
  testInstanceIndex = 1
  for testInstanceUserId in collect(keys(testUsersBasketsDict))
    testInstanceItems = deepcopy(testUsersBasketsDict[testInstanceUserId].basketItems)
    heldOutItem = testUsersBasketsDict[testInstanceUserId].heldOutItem
    push!(testInstanceItems, heldOutItem)

    # Remove duplicates from testInstanceItems
    testInstanceItems = collect(Set(testInstanceItems))

    testInstances[testInstanceIndex] = testInstanceItems

    testInstanceIndex += 1
  end

  # Compute ordered item counts in training data, by descending popularity
  orderedItemCounts = PriorityQueue(Dict{Int, Int}(), Base.Order.Reverse)
  for i = 1:numTrainingInstances
    trainingInstance = trainingInstances[i]

    for itemId in trainingInstance
      if haskey(orderedItemCounts, itemId)
        orderedItemCounts[itemId] += 1
      else
        # Create a new entry for this item
        orderedItemCounts[itemId] = 1
      end
    end
  end

  # Build lambdaVec, which is a vector of per-item lambda regularization
  # parameters that are inversely proportional to item counts (popularity)
  lambdaVec = zeros(numItems)
  for i = 1:numItems
    if haskey(orderedItemCounts, i)
      lambdaVec[i] = 1 / orderedItemCounts[i]
    else
      # Item is not present in training data, so assume that it has the
      # lowest popularity value available
      lambdaVec[i] = 1
    end
  end

  println("Beginning stochastic gradient ascent...")

  # Perform stochastic gradient ascent to learn the DPP kernel parameters (item trait vectors)
  paramsMatrix, avgValidationLogLikelihood, runTime =
    doStochasticGradientAscent(trainingInstances, numTrainingInstances,
      numItems, numItemTraits, testInstances, numTestInstances, lambdaVec, alpha,
      verbose = true)

  # Compute average log-likelihood on training and test data
  avgTrainingLogLikelihood =
    computeLogLikelihood(paramsMatrix, trainingInstances, numTrainingInstances,
                         numItems, numItemTraits) / numTrainingInstances
  avgTestLogLikelihood =
    computeLogLikelihood(paramsMatrix, testInstances, numTestInstances,
                         numItems, numItemTraits) / numTestInstances
  percentDiffTraingTestLikelihood =
    (abs(avgTrainingLogLikelihood - avgTestLogLikelihood) / abs(avgTrainingLogLikelihood)) * 100

  println("Number of item trait dimensions: $numItemTraits")
  println("Final avg training log-likelihood: $avgTrainingLogLikelihood")
  println("Final avg test log-likelihood: $avgTestLogLikelihood")
  println("Diff between avg training and test log-likelihood: $percentDiffTraingTestLikelihood%")

  # Save paramsMatrix to disk, so that we can reuse it later
  if !isdir(learnedModelOutputDirName)
    mkdir(learnedModelOutputDirName)
  end
  open(f -> serialize(f, paramsMatrix),
    "$learnedModelOutputDirName/learnedDPPParamsMatrix-k$numItemTraits-lambdaPop$alpha.jls", "w");
end


text:
using MPIFiles
using Test

fnMeasBruker = "./data/measurement"
tfpath = "./data/transferFunction/example.s1p"
tfh5path = "./data/transferFunction/example.h5"

@testset "Testing TransferFunction submodule" begin
  a = TransferFunction(tfpath, frequencyWeighting=true)
  if isfile(tfh5path)
    rm(tfh5path)
  end
  MPIFiles.save(tfh5path, a)
  b = TransferFunction(tfh5path)

  @test a.freq == b.freq
  @test a.data == b.data
  @test a[[1],1] == a[[0.0],1]

  c = MPIFiles.combine(MPIFiles.combine(a,a),a)
  rm(tfh5path)
  MPIFiles.save(tfh5path, c)
  @test a[[1],1] == c[[1],2]
  @test a[[1],1] == c[[1],3]

  measBruker = MPIFile(fnMeasBruker)
  tf = sampleTF(c, measBruker)
  @test size(tf) == (rxNumFrequencies(measBruker), rxNumChannels(measBruker))
end