text:
using BinaryBuilder

# Collection of sources required to build HDF5
name = "HDF5"
version = v"1.12.0"

sources = [
    FileSource("https://files.pythonhosted.org/packages/d5/f9/676c6a5c13806289da6177c538ce772e3e5b04ea10d76e6e72e9f0d042de/h5py-3.1.0-cp39-cp39-macosx_10_9_x86_64.whl", "cb74df83709d6d03d11e60b9480812f58da34f194beafa8c8314dbbeeedfe0a6"),

    FileSource("https://files.pythonhosted.org/packages/40/1b/dd36e8aa1b7c9d82a1f0aaabece4393797ded0ff613437dfb8c0780b33a9/h5py-3.1.0-cp39-cp39-manylinux1_x86_64.whl", "80c623be10479e81b64fa713b7ed4c0bbe9f02e8e7d2a2e5382336087b615ce4"),

    ArchiveSource("http://repo.msys2.org/mingw/i686/mingw-w64-i686-hdf5-1.12.0-2-any.pkg.tar.zst", "d9ade0d0fddfdeca3ea9de00b066e330e1573c547609a12b81c6a080b2c19f3e", unpack_target = "i686-w64-mingw32"),
    ArchiveSource("http://repo.msys2.org/mingw/i686/mingw-w64-i686-szip-2.1.1-2-any.pkg.tar.xz", "58b5efe1420a2bfd6e92cf94112d29b03ec588f54f4a995a1b26034076f0d369", unpack_target = "i686-w64-mingw32"),
    ArchiveSource("http://repo.msys2.org/mingw/i686/mingw-w64-i686-zlib-1.2.11-7-any.pkg.tar.xz", "addf6c52134027407640f1cbdf4efc5b64430f3a286cb4e4c4f5dbb44ce55a42", unpack_target = "i686-w64-mingw32"),
    # We need some special compiler support libraries from mingw for i686
    ArchiveSource("http://repo.msys2.org/mingw/i686/mingw-w64-i686-gcc-libs-10.2.0-5-any.pkg.tar.zst", "e03a63b24695951a1e80def754d6bd128744eaf1e562308ded31e989636e7651", unpack_target = "i686-w64-mingw32"),

    ArchiveSource("http://repo.msys2.org/mingw/x86_64/mingw-w64-x86_64-hdf5-1.12.0-2-any.pkg.tar.zst", "549462ad99a079ff725ac4bd1f662d3594515320ea324a7263a647578b258d86", unpack_target = "x86_64-w64-mingw32"),
    ArchiveSource("http://repo.msys2.org/mingw/x86_64/mingw-w64-x86_64-szip-2.1.1-2-any.pkg.tar.xz", "ec8fe26370b0673c4b91f5ccf3404907dc7c24cb9d75c7b8830aa93a7c13ace7", unpack_target = "x86_64-w64-mingw32"),
    ArchiveSource("http://repo.msys2.org/mingw/x86_64/mingw-w64-x86_64-zlib-1.2.11-7-any.pkg.tar.xz", "1decf05b8ae6ab10ddc9035929014837c18dd76da825329023da835aec53cec2", unpack_target = "x86_64-w64-mingw32"),

    # Can't use conda-forge on other platforms since it links too many libraries, but apparently on aarch64 is fine
    ArchiveSource("https://anaconda.org/conda-forge/hdf5/1.12.0/download/linux-aarch64/hdf5-1.12.0-nompi_h1022a3e_102.tar.bz2", "605aff906fd0fca9a52da6ad9b48607fab5cb26e2615d3827a1f318d6e103c4a", unpack_target = "aarch64-linux-gnu"),
]

# Bash recipe for building across all platforms
script = raw"""
cd ${WORKSPACE}/srcdir/
mkdir -p ${libdir} ${includedir}

if [[ ${target} == *mingw* ]]; then
    cd ${target}/mingw${nbits}

    rm -f bin/{*_cpp*,*fortran*,*f90*} # we do not need these
    mv bin/*.dll ${libdir}
    mv include/* ${includedir}

    install_license share/doc/hdf5/COPYING
elif [[ ${target} == aarch64-* ]]; then
    cd ${target}

    rm -f lib/{*_cpp*,*_fortran*} # we do not need these
    mv lib/* ${libdir}
    mv include/* ${includedir}

    install_license info/licenses/COPYING
else
    if [[ ${target} == x86_64-linux-gnu ]]; then
        WHL_FILE="*h5py-*manylinux1_x86_64*.whl"
        LIBSDIR=h5py.libs
    elif [[ ${target} == x86_64-apple-darwin* ]]; then
        WHL_FILE="*h5py-*macosx*.whl"
        LIBSDIR=h5py/.dylibs
    else
        echo "ERROR: Unsupported platform ${target}" >&2
        exit 1
    fi

    unzip ${WHL_FILE}
    mv ${LIBSDIR}/* ${libdir}

    # Use headers and license from aarch64-linux-gnu with the hope that they'll be fine
    mv aarch64-*/include/* ${includedir}
    install_license aarch64-*/info/licenses/COPYING
fi

# We need to be able to access `libhdf5` and `libhdf5_hl` directly, so symlink it from the hashed filename from manylinux pypi
if [[ ${target} == *86*linux* ]]; then
    libhdf5name=$(basename ${libdir}/libhdf5-*.${dlext}*)
    base="${libhdf5name%%.*}"
    ext="${libhdf5name#$base}"
    ln -s ${libhdf5name} ${libdir}/libhdf5${ext}

    libhdf5_hlname=$(basename ${libdir}/libhdf5_hl-*.${dlext}*)
    base="${libhdf5_hlname%%.*}"
    ext="${libhdf5_hlname#$base}"
    ln -s ${libhdf5_hlname} ${libdir}/libhdf5_hl${ext}
fi

"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Platform("x86_64", "linux"),
    Platform("aarch64", "linux"; libc="glibc"),
    Platform("x86_64", "macos"),
    Platform("x86_64", "windows"),
    Platform("i686", "windows"),
]

# The products that we will ensure are always built
products = [
    LibraryProduct("libhdf5", :libhdf5),
    LibraryProduct("libhdf5_hl", :libhdf5_hl),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("Zlib_jll"),
    Dependency("OpenSSL_jll"),
    Dependency("LibCURL_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: vchuravy/Yggdrasil

max_issues_repo_name: vchuravy/Yggdrasil

max_forks_repo_name: vchuravy/Yggdrasil


text:
module Bolt

export CosmoParams, AbstractCosmoParams
export Background, AbstractBackground
export IonizationHistory, AbstractIonizationHistory
export Peebles
export Hierarchy, boltsolve, BasicNewtonian
export source_grid, quadratic_k, cltt
export z2a, a2z, x2a, a2x, z2x, x2z

using Parameters
using Unitful, UnitfulAstro, NaturallyUnitful
using NLsolve
using OrdinaryDiffEq
using Interpolations
using OffsetArrays
using QuadGK
using ThreadPools
using ForwardDiff, DiffResults
using NumericalIntegration

using FFTW
import SpecialFunctions: lgamma, sphericalbesselj
import AbstractFFTs: fftfreq, Plan, plan_fft!, plan_ifft!
import LinearAlgebra: mul!, ldiv!

import PhysicalConstants.CODATA2018: ElectronMass, ProtonMass,
    FineStructureConstant, ThomsonCrossSection, NewtonianConstantOfGravitation

abstract type AbstractCosmoParams{T} end

@with_kw struct CosmoParams{T} <: AbstractCosmoParams{T} @deftype T
    h = 0.7  # hubble factor
    Ω_r = 5.042e-5  # radiation density
    Ω_b = 0.046  # baryon density
    Ω_m = 0.224  # matter density
    n = 1.0  # spectral index
    Y_p = 0.0  # primordial helium fraction, currently unused
    N_ν = 3.046 #effective number of relativisic species (PDG25 value)
end

include("util.jl")
include("fftlog.jl")
include("background.jl")
include("ionization.jl")
include("perturbations.jl")
include("spectra.jl")

end

ext: jl


text:
function SteeringBVP(f::LinearDynamics{Dx,Du}, j::TimePlusQuadraticControl{Du};
                     compile::Union{Val{false},Val{true}}=Val(false)) where {Dx,Du}
    compile === Val(true) ? SteeringBVP(f, j, EmptySteeringConstraints(), LinearQuadraticHelpers(f.A, f.B, f.c, j.R)) :
                            SteeringBVP(f, j, EmptySteeringConstraints(), EmptySteeringCache())
end

struct LinearQuadraticHelpers{FGinv<:Function,
                              FexpAt<:Function,
                              Fcdrift<:Function,
                              Fcost<:Function,
                              Fdcost<:Function,
                              Fddcost<:Function,
                              Fx<:Function,
                              Fu<:Function} <: SteeringCache
    Ginv::FGinv
    expAt::FexpAt
    cdrift::Fcdrift
    cost::Fcost
    dcost::Fdcost
    ddcost::Fddcost
    x::Fx
    u::Fu
    symbolic_exprs::Dict{String,Union{SymPy.Sym,Vector{SymPy.Sym},Matrix{SymPy.Sym}}}
end

function (bvp::LinearQuadraticSteering{Dx,Du,<:LinearQuadraticHelpers})(x0::StaticVector{Dx}, xf::StaticVector{Dx},
                                                                        c_max::T=eltype(x0)(1e6)) where {Dx,Du,T<:Number}    # TODO: handle c_max == Inf
    x0 == xf && return (cost=T(0), controls=BVPControl(T(0), x0, xf, bvp.cache.x, bvp.cache.u))
    t = optimal_time(bvp, x0, xf, c_max)
    (cost=bvp.cache.cost(x0, xf, t), controls=BVPControl(t, x0, xf, bvp.cache.x, bvp.cache.u))
end

function LinearQuadraticHelpers(A_::AbstractMatrix, B_::AbstractMatrix, c_::AbstractVector, R_::AbstractMatrix)
    A, B, c, R = Array(A_), Array(B_), Vector(c_), Array(R_)
    Dx, Du = size(B)
    t, s = SymPy.symbols("t s", real=true)
    x = collect(SymPy.symbols(join(("x$i" for i in 1:Dx), " "), real=true))
    y = collect(SymPy.symbols(join(("y$i" for i in 1:Dx), " "), real=true))

    expAt = (A*t).exp()
    expAs = (A*s).exp()
    expAt_s = (A*(t - s)).exp()
    G = SymPy.integrate.(expAt*B*inv(R)*B'*expAt', t)
    Ginv = inv(G)
    cdrift = SymPy.integrate.(expAt, t)*c
    xbar = expAt*x + cdrift
    cost = t + (y - xbar)'*Ginv*(y - xbar)
    dcost = diff(cost, t)
    ddcost = diff(cost, t, 2)
    x_s = expAs*x + SymPy.integrate.(expAs, s)*c + SymPy.integrate.(expAs*B*inv(R)*B'*expAs', s)*expAt_s'*Ginv*(y-xbar)
    u_s = inv(R)*B'*expAt_s.transpose()*Ginv*(y-xbar)  # transpose needed here (and technically above) to avoid `Any`s

    symbolic_exprs = Dict{String,Union{SymPy.Sym,Vector{SymPy.Sym},Matrix{SymPy.Sym}}}(
        "Ginv" => Ginv,
        "expAt" => expAt,
        "cdrift" => cdrift,
        "cost" => cost,
        "dcost" => dcost,
        "ddcost" => ddcost,
        "x_s" => x_s,
        "u_s" => u_s,
        "t" => t,
        "s" => s,
        "x" => x,
        "y" => y
    )
    for (k, v) in symbolic_exprs
        symbolic_exprs[k] = SymPy.collect.(SymPy.expand.(v), t)
    end

    symbol_dict = merge(Dict(Symbol("x$i") => :(x[$i]) for i in 1:Dx),
                        Dict(Symbol("y$i") => :(y[$i]) for i in 1:Dx))
    sarray    = A_ isa StaticArray
    t_args    = :((t::T) where {T})
    xyt_args  = :((x::AbstractVector, y::AbstractVector, t::T) where {T})
    xyts_args = :((x::AbstractVector, y::AbstractVector, t::T, s) where {T})
    LinearQuadraticHelpers(
        code2func(sympy2code.(symbolic_exprs["Ginv"], Ref(symbol_dict)), t_args, sarray),
        code2func(sympy2code.(symbolic_exprs["expAt"], Ref(symbol_dict)), t_args, sarray),
        code2func(sympy2code.(symbolic_exprs["cdrift"], Ref(symbol_dict)), t_args, sarray),
        code2func(sympy2code.(symbolic_exprs["cost"], Ref(symbol_dict)), xyt_args, sarray),
        code2func(sympy2code.(symbolic_exprs["dcost"], Ref(symbol_dict)), xyt_args, sarray),
        code2func(sympy2code.(symbolic_exprs["ddcost"], Ref(symbol_dict)), xyt_args, sarray),
        code2func(sympy2code.(symbolic_exprs["x_s"], Ref(symbol_dict)), xyts_args, sarray),
        code2func(sympy2code.(symbolic_exprs["u_s"], Ref(symbol_dict)), xyts_args, sarray),
        symbolic_exprs
    )
end

function sympy2code(x, symbol_dict = Dict())
    code = foldl(replace, (".+" => " .+",
                           ".-" => " .-",
                           ".*" => " .*",
                           "./" => " ./",
                           ".^" => " .^"); init=SymPy.sympy.julia_code(x))
    expr = Meta.parse(code)
    MacroTools.postwalk(x -> x isa AbstractFloat ? :(T($x)) : get(symbol_dict, x, x), expr)
end

code2func(code, args, static_array = true) = eval(:($args -> $code))
function code2func(code::AbstractVector, args, static_array = true)
    N = length(code)
    if static_array
        body = :(SVector{$N}($(code...)))
    else
        body = :([$(code...)])
    end
    eval(:($args -> $body))
end
function code2func(code::AbstractMatrix, args, static_array = true)
    M,N = size(code)
    if static_array
        body = :(SMatrix{$M,$N}($(code...)))
    else
        body = Expr(:vcat, (Expr(:row, code[i,:]...) for i in 1:M)...)
    end
    eval(:($args -> $body))
end

function optimal_time(bvp::SteeringBVP{D,C,EmptySteeringConstraints,<:LinearQuadraticHelpers},
                      x0::StaticVector{Dx},
                      xf::StaticVector{Dx},
                      t_max::T) where {Dx,Du,T<:Number,D<:LinearDynamics{Dx,Du},C<:TimePlusQuadraticControl{Du}}
    cost   = (s -> (Base.@_inline_meta; bvp.cache.cost(x0, xf, s)))    # closures and optimizers below both @inline-d
    dcost  = (s -> (Base.@_inline_meta; bvp.cache.dcost(x0, xf, s)))
    ddcost = (s -> (Base.@_inline_meta; bvp.cache.ddcost(x0, xf, s)))
    t = (T === Float64 ? newton(dcost, ddcost, t_max/100, t_max) :
                         bisection(dcost, t_max/100, t_max))
    t !== nothing ? t : golden_section(cost, t_max/100, t_max)
end


text:
push!(LOAD_PATH,"src")
using DalitzPlotAnalysis
using amplitudes_compass

event = let s = 1.5, s1 = 0.8
    [s,s1,0.3,0.1π,0.3,0.1π]
end
 # = [COMPASS_wave(i,event...) for e=1:size(mm,1), i=1:88];

cal = [COMPASS_wave(i,event...) for i in 1:88]
for i in 1:88
    println(cal[i])
end

@show Z(1,0,true,0,1,event[3:end]...)

@show ClebschGordon((2*[1,1,1,0,1,1])...)
@show WignerDϵ(true,1,0,0,event[3:4]...,0.0)
@show Wignerd(1,1,0,event[5])


text:
# packages -
using LinearAlgebra
using GLPK

# my codes -
include("Flux.jl")
include("Expa.jl")
include("Stoichiometric.jl")
include("Utility.jl")


text:
module Tetris

  using HTTP
  using JSON

  using WebIO
  using JSExpr

  using DataStructures
  using CodecZlib

  using Colors
  using Crayons

  using Nullables
  using UUIDs

  using Dates
  using REPL

  using Random

  include("abstracts.jl")
  include("consts.jl")
  include("api.jl")

  include("structs/index.jl")
  include("clocks/index.jl")

  include("renders/index.jl")
  include("tetrominoes/index.jl")

  include("functions/index.jl")
  include("actions/index.jl")

  include("raw_terminal.jl")
  include("setup.jl")

  macro tetris(cur_expr::Expr)
    cur_expr = quote
      setup()

      $(
        Expr(
          :macrocall,
          Symbol("@async"),
          esc(cur_expr)
        )
      )
    end

    return cur_expr
  end

  export @tetris

  function tetris_js(cur_scope::Scope, cur_string::JSString)
    is_repl && return
    evaljs(cur_scope, cur_string)
  end

  function __init__()
    global terminal

    global is_ijulia
    global is_juno
    global is_vscode
    global is_ide
    global is_repl

    terminal = REPL.Terminals.TTYTerminal(
      get(ENV, "TERM", Sys.iswindows() ? "" : "dumb"),
      stdin, stdout, stderr
    )

    is_ijulia = isdefined(Main, :IJulia) && Main.IJulia.inited
    is_juno = isdefined(Main, :Juno) && Main.Juno.isactive()
    is_vscode = isdefined(Main, :_vscodeserver)

    is_ide = is_ijulia || is_juno || is_vscode
    is_repl = !is_ide
  end

end


text:
# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.

using UnitCommitment, LinearAlgebra, JuMP
USE_GUROBI = (Base.find_package("Gurobi") != nothing)
USE_CBC = !USE_GUROBI
if USE_GUROBI
  using Gurobi
else
  using Cbc
end

NUM_THREADS = 4
LOG_LEVEL = 1

@testset "Model" begin
    @testset "Run" begin
        instance = UnitCommitment.read_benchmark("test/case14")
        #instance = UnitCommitment.read_benchmark("matpower/case3375wp/2017-02-01")
        #instance = UnitCommitment.read_benchmark("matpower/case1888rte/2017-02-01")
        for line in instance.lines, t in 1:4
            line.normal_flow_limit[t] = 10.0
        end
        #for formulation in [UnitCommitment.DefaultFormulation, UnitCommitment.TightFormulation]
        for formulation in [UnitCommitment.TightFormulation]
            @info string("Running test of ", formulation)
            if USE_CBC
              optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => LOG_LEVEL)
            end
            if USE_GUROBI
              optimizer = optimizer_with_attributes(Gurobi.Optimizer, "Threads" => NUM_THREADS)
            end
            model = build_model(instance=instance,
                                optimizer=optimizer,
                                variable_names=true,
                                formulation=formulation)

            JuMP.write_to_file(model.mip, "test.mps")

            # Optimize and retrieve solution
            UnitCommitment.optimize!(model)
            solution = get_solution(model)

            # Verify solution
            @test UnitCommitment.validate(instance, solution)

            # Reoptimize with fixed solution
            UnitCommitment.fix!(model, solution)
            UnitCommitment.optimize!(model)
            @test UnitCommitment.validate(instance, solution)

            #@show solution
        end # loop over components
    end # end testset Run
end # end test


text:
"""
`WAVChunk(id, data)` represents a RIFF chunk.
Symbol `id` is the four-character chunk ID.
"""
struct WAVChunk
    id::Symbol
    data::Vector{UInt8}
end

"""
A marker in a .wav file. `start_time` and `duration` are in samples.
"""
mutable struct WAVMarker
    label::String
    start_time::UInt32
    duration::UInt32
end

function read32(data0::UInt8, data1::UInt8, data2::UInt8, data3::UInt8)::UInt32
    UInt32(data0) | (UInt32(data1) << 8) | (UInt32(data2) << 16) | (UInt32(data3) << 24)
end

function write32(data::UInt32)::Vector{UInt8}
    [ UInt8((data>>0) & 0xff), UInt8((data>>8) & 0xff), UInt8((data>>16) & 0xff), UInt8((data>>24) & 0xff) ]
end

function write16(data::UInt16)::Vector{UInt8}
    [ UInt8((data>>0) & 0xff), UInt8((data>>8) & 0xff) ]
end

function read_adtl(markers::Dict{UInt32, WAVMarker}, adtl::Vector{UInt8})
    title  = String(adtl[1:4])
    size   = read32(adtl[5:8 ]...)
    cue_id = read32(adtl[9:12]...)

    # The adtl entry must have even length. Therefore, if the reported length
    # is odd, actually read an extra byte.
    if size % 2 == 1
        size += 1
    end

    marker = get(markers, cue_id, nothing)
    if marker == nothing
        marker = WAVMarker("", 0, 0)
        markers[cue_id] = marker
    end

    if title == "labl"
        i = 13;
        while (adtl[i] != 0)
            i += 1;
        end
        marker.label = String(adtl[13:(i-1)])
    elseif title == "ltxt"
        marker.duration = read32(adtl[13:16]...)
    end
    adtl[(size+5+4):end]
end

"""Read the contents of the LIST chunk to extract marker names and durations"""
function read_list(markers::Dict{UInt32, WAVMarker}, list::Vector{UInt8})
    if list[1:4] == b"adtl"
        adtl = list[5:end]
        while (length(adtl) >= 12)
            adtl = read_adtl(markers, adtl)
        end
    end
end

"""Read the contents of the cue chunk to extract marker start times"""
function read_cue(markers::Dict{UInt32, WAVMarker}, cue::Vector{UInt8})
    ncue = read32(cue[1:4]...)
    cue = cue[5:end]
    for i = 1:ncue
        cue_id = read32(cue[1:4]...)

        marker = get(markers, cue_id, nothing)
        if marker == nothing
            marker = WAVMarker("", 0, 0)
            markers[cue_id] = marker
        end
        marker.start_time = read32(cue[21:24]...)
        st = marker.start_time

        cue = cue[25:end]
    end
end

"""
    wav_cue_read(chunks::Vector{WAVChunk})

Takes a `Vector{WAVChunk}` (as returned by `wavread`) and returns
a `Vector{WAVMarker}`, where a `WAVMarker` is defined as:

```julia
mutable struct WAVMarker
    label::String
    start_time::UInt32
    duration::UInt32
end
```

Field values `start_time` and `duration` are in samples.

# Example
```julia
using WAV
x, fs, bits, in_chunks = wavread("in.wav")
markers = wav_cue_read(in_chunks)
```
"""
function wav_cue_read(chunks::Vector{WAVChunk})
    markers = Dict{UInt32, WAVMarker}()

    # See if list and cue chunks are present
    list_chunks = chunks[findall(c -> c.id == :LIST, chunks)]
    cue_chunks = chunks[findall(c -> c.id == Symbol("cue "), chunks)]

    for l in list_chunks
        read_list(markers, l.data)
    end

    for c in cue_chunks
        read_cue(markers, c.data)
    end

    markers
end

function write_cue(markers::Dict{UInt32, WAVMarker})
    cue = write32(UInt32(length(markers)))
    for (cue_id, marker) in markers
        cue =   [   cue;
                    write32(cue_id);
                    write32(marker.start_time);
                    b"data";
                    write32(UInt32(0));
                    write32(UInt32(0));
                    write32(marker.start_time)
                ]
    end
    cue
end

function write_marker_list(markers::Dict{UInt32, WAVMarker})
    list = b"adtl"

    # Create all the labl entries
    for (cue_id, marker) in markers
        labl = [write32(cue_id); codeunits(marker.label); 0x0]

        # The note and label entries must have an even number of bytes.
        # So, for the null terminated text in the label and note, we add a minimum of
        # one null terminator, but if that creates an odd number of bytes in the labl
        # or note entry, then add a second null terminator.
        if (length(labl) % 2) == 1
            labl = [labl; 0x0]
        end

        list = [list; b"labl"; write32(UInt32(length(labl))); labl]
    end

    # Create all the ltxt entries
    for (cue_id, marker) in markers
        ltxt =  [   write32(cue_id);
                    write32(marker.duration);
                    b"rgn ";
                    write16(UInt16(0)); # country
                    write16(UInt16(0)); # language
                    write16(UInt16(0)); # dialect
                    write16(UInt16(0))  # codepage
                ]
        list = [list; b"ltxt"; write32(UInt32(length(ltxt))); ltxt]
    end
    list
end

"""
    wav_info_write(tags::Dict{Symbol, String})::Vector{WAVChunk}

Converts a dictionary of INFO tags into a list of WAV chunks
appropriate for passing to `wavwrite`.

`tags` is a dictionary where the keys are symbols representing four-character
RIFF INFO tag IDs as specified in
https://www.sno.phy.queensu.ca/~phil/exiftool/TagNames/RIFF.html#Info
The values of the dictionary correspond to the tag data.
"""
function wav_info_write(tags::Dict{Symbol, String})
    info_data = b"INFO"

    # Create all the tag entries
    for t in keys(tags)
        tag = [codeunits(tags[t]); 0x0]

        # The tag entries must have an even number of bytes.
        # So, for the null terminated text in the tag, we add a minimum of
        # one null terminator, but if that creates an odd number of bytes in the tag
        # or note entry, then add a second null terminator.
        if (length(tag) % 2) == 1
            tag = [tag; 0x0]
        end

        info_data = [info_data; codeunits(String(t)); write32(UInt32(length(tag))); tag]
    end
    [WAVChunk(:LIST, info_data)]
end

function read_tag(tags::Dict{Symbol, String}, t::Vector{UInt8})
    tag_id = Symbol(String(t[1:4]))
    size   = read32(t[5:8 ]...)

    # The adtl entry must have even length. Therefore, if the reported length
    # is odd, actually read an extra byte.
    if size % 2 == 1
        size += 1
    end

    i = 9
    while t[i] != 0
        i += 1
    end
    tags[tag_id] = String(t[9:(i-1)])

    # Skip the null terminator
    t[(i+1):end]
end

"""
    wav_info_read(chunks::Vector{WAVChunk})::Dict{Symbol, String}

Given a list of chunks as returned by `wavread`, return a
Dict{Symbol, String} where the keys are symbols representing
four-character RIFF INFO tag IDs as specified in
https://www.sno.phy.queensu.ca/~phil/exiftool/TagNames/RIFF.html#Info
"""
function wav_info_read(chunks::Vector{WAVChunk})
    tags = Dict{Symbol, String}()

    list_chunks = chunks[findall(c -> c.id == :LIST, chunks)]
    for l in list_chunks
        list_data = l.data
        if list_data[1:4] == b"INFO"
            list_data = list_data[5:end]
            while (length(list_data) >= 12)
                list_data = read_tag(tags, list_data)
            end
        end
    end

    tags
end

"""
    wav_cue_write(markers::Dict{UInt32, WAVMarker})

Turns `WAVMarker`s into a `Vector{WAVChunk}` (as accepted by
`wavwrite`). The key for the dictionary is the ID of the marker to be
written to file.

Example:
```julia
out_chunks = wav_cue_write(markers)
wavwrite(x, "out.wav", Fs=fs, nbits=16, compression=WAVE_FORMAT_PCM, chunks=out_chunks)
```
"""
function wav_cue_write(markers::Dict{UInt32, WAVMarker})
    chunks = WAVChunk[]
    push!(chunks, WAVChunk(Symbol("cue "), write_cue(markers)))
    push!(chunks, WAVChunk(:LIST, write_marker_list(markers)))
    chunks
end


text:
show_name = "Stargate SG-1"
show_dir = joinpath(tempdir(), show_name)
episodes = joinpath.(show_dir, [joinpath("Season $i", "$show_name - s$(i)e$(e).mkv") for i ∈ lpad.(1:10, 2, '0'), e ∈ lpad.(1:20, 2, '0')])
mkpath.(dirname.(episodes))
touch.(episodes)

@testset "Discovery" begin
    @testset "Episodes" begin
        @test all(discover_episodes(dirname(episodes[1, 1])) .== Episode.(1:20, episodes[1, :]))
        @test all(discover_seasons(show_name, show_dir) .== [Season(i, Episode.(1:20, episodes[i, :])) for i ∈ 1:10])
    end
end


text:
using LAR

# test coarse surface approximation (r=1, R=3)
W,lar = toroidalsurface(1,3,10,20);
larview(W,lar.FV)
larview(W,lar.EV)
larview(W,lar.VV)
viewexploded(W,lar.FV)
viewexploded(W,lar.EV)
viewexploded(W,lar.VV)

# test non-manifold surface  (r=3, R=3)
V,lar = toroidalsurface(1,1,40,80);
larview(V,lar.FV)
viewexploded(V,lar.FV)

# test coarsest surface approximation (r=1, R=3)
W,lar = toroidalsurface(1,1,3,3);
larview(W,lar.FV)
larview(W,lar.EV)
viewexploded(W,lar.FV)
viewexploded(W,lar.EV)


text:
_NamedDimsArray_pullback(ȳ::AbstractArray) = (NoTangent(), ȳ, NoTangent())
_NamedDimsArray_pullback(ȳ::Tangent) = (NoTangent(), ȳ.data, NoTangent())
_NamedDimsArray_pullback(ȳ::AbstractThunk) = _NamedDimsArray_pullback(unthunk(ȳ))
function ChainRulesCore.rrule(::Type{<:NamedDimsArray}, values, names)
    return NamedDimsArray(values, names), _NamedDimsArray_pullback
end

function ChainRulesCore.rrule(T::Type{<:NamedDimsArray}, values)
    NamedDimsArray_values_pullback(ȳ) = _NamedDimsArray_pullback(ȳ)[1:2]
    return T(values), NamedDimsArray_values_pullback
end


text:
module aIFileIO

using CSV, DataFrames

include("aIFileFunctions.jl")
include("aIDFFunctions.jl")

export log_file, read_log, count_log_events, read_sourceFolder
export process_sourceFolder, count_CSV_files, count_TXT_files, set_sourcePath, show_sourcePath
export show_CSV_files, show_TXT_files

export get_file_symbols

end


text:
############################
# Specialized factorizations
############################

export id

"""
An interpolative decomposition.

For a matrix `A`, the interpolative decomposition `F` contains the matrices `B`
and `P` computed by `id()`. See the documentation of `id()` for details.

# References

\\cite{Cheng2005, Liberty2007}
"""
struct Interpolative{T} <: Factorization{T}
    B :: AbstractMatrix{T}
    P :: AbstractMatrix{T}
end

"""
    id(A, k, l)

Compute and return the interpolative decomposition of `A`: A ≈ B * P

Where:
* `B`'s columns are a subset of the columns of `A`
* some subset of `P`'s columns are the `k x k` identity, no entry of `P` exceeds magnitude 2, and
* ||B * P - A|| ≲ σ(A, k+1), the (`k+1`)st singular value of `A`.

# Arguments

`A`: Matrix to factorize

`k::Int`: Number of columns of A to return in B

`l::Int`: Length of random vectors to project onto

# Output

`(::Interpolative)`: interpolative decomposition.

# Implementation note

This is a hacky version of the algorithms described in \\cite{Liberty2007}
and \\cite{Cheng2005}. The former refers to the factorization (3.1) of the
latter.  However, it is not actually necessary to compute this
factorization in its entirely to compute an interpolative decomposition.

Instead, it suffices to find some permutation of the first k columns of Y =
R * A, extract the subset of A into B, then compute the P matrix as B\\A
which will automatically compute P using a suitable least-squares
algorithm.

The approximation we use here is to compute the column pivots of Y,
rather then use the true column pivots as would be computed by a column-
pivoted QR process.

# References

\\cite[Algorithm I]{Liberty2007}

```bibtex
@article{Cheng2005,
    author = {Cheng, H and Gimbutas, Z and Martinsson, P G and Rokhlin, V},
    doi = {10.1137/030602678},
    issn = {1064-8275},
    journal = {SIAM Journal on Scientific Computing},
    month = jan,
    number = {4},
    pages = {1389--1404},
    title = {On the Compression of Low Rank Matrices},
    volume = {26},
    year = {2005}
}
```
"""
function id(A, k::Int, l::Int)
    m, n = size(A)
    R = randn(l, m)
    Y = R * A #size l x n

    #Compute column pivots of first k columns of Y
    maxvals = map(j->maximum(abs.(view(Y, :, j))), 1:n)
    piv = sortperm(maxvals, rev=true)[1:k]

    B = A[:, piv]
    Interpolative(B, B\A)
end


text:
function setup()::Nothing
    if !isdefined(Main, :MicroCoverage_lock)
        @eval Main const MicroCoverage_lock = ReentrantLock()
        lock(Main.MicroCoverage_lock) do
            @eval Main const MicroCoverage_tracker = Vector{Bool}(undef, 0)
            @eval Main const MicroCoverage_tracker_linenumbernodes = Vector{LineNumberNode}(undef, 0)
        end
    end
    return nothing
end

function start(path::AbstractString)
    setup()
    if isfile(path)
        return _start_tracking_file(path)
    else
        throw(ArgumentError("Path must be an existing file"))
    end
end

function stop(; dump_coverage::Bool = false,
                dump_coverage_io::IO = stdout,
                min_padding::Integer = 1,
                max_padding::Integer = 24)
    setup()
    all_tracked_files = tracked_files()
    for x in all_tracked_files
        _restore_original_julia_file(x)
    end
    for x in all_tracked_files
        _write_coverage(x; dump_coverage = dump_coverage,
                           dump_coverage_io = dump_coverage_io,
                           min_padding = min_padding,
                           max_padding = max_padding)
    end
    return nothing
end

function clean(path::AbstractString = pwd())
    if isfile(path)
        return _clean_file(convert(String, realpath(path)))
    elseif isdir(path)
        return _clean_directory(convert(String, realpath(path)))
    else
        throw(ArgumentError("Path must be an existing file or directory"))
    end
end


text:
using Distributions
using RmathDist

distlist = [
            Cauchy(0,1),
            Cauchy(1e10,1e10),
            Cauchy(1e-10,1e-10)
            ]

for d in distlist
    for x in [-logspace(-300,300,1_000_000),logspace(-300,300,1_000_000)]
        test_rmath(pdf,d,x)
        if abs(x) < 1e140
            # rmath overflows for large x
            test_rmath(logpdf,d,x,1e5)
        end

        test_rmath(cdf,d,x)
        test_rmath(ccdf,d,x)
        test_rmath(logcdf,d,x)
        test_rmath(logccdf,d,x)
    end

    for p in logspace(-300,0,1_000_000)
        test_rmath(quantile,d,p,1e5)
        test_rmath(cquantile,d,p,1e5)
    end

    for lp in -logspace(-300,300,1_000_000)
        test_rmath(invlogcdf,d,lp,1e5)
        test_rmath(invlogccdf,d,lp,1e5)
    end
end


text:
# Run the model for the ensemble of solutions -
include("SolveBalances.jl")
include("Common.jl")
using PyPlot
using PyCall
@pyimport numpy as np

# Load the ensemble files from disk -
pc_array_full = readdlm("./data/pc_array_O1_O2.dat")
rank_array = readdlm("./data/rank_array_O1_O2.dat")

# Select the desired rank -
idx_rank = find(rank_array .<= 5.0)

# Setup time scale -
tStart = 0.0;
tStop = 25.0;
tStep = 0.1;
number_of_timesteps = 200
time_experimental = linspace(tStart,tStop,number_of_timesteps)

# initialize data array -
distance_array_set = Set{Array{Float64,2}}()

data_dictionary = Dict()
initial_condition_array = zeros(19)
initial_condition_array[1] = 0.0    # 1 zymosan
initial_condition_array[2] = 1.9      # 2 C4
initial_condition_array[3] = 0.322    # 3 C2
initial_condition_array[8] = 7.57     # 4 C3
initial_condition_array[14] = 0.195   # 5 C5
initial_condition_array[17] = 2.23    # 6 Factor H
initial_condition_array[18] = 0.417   # 7 C4BP

initial_condition_array[19] = 0.55     # 19 C3a (ic initial)
initial_condition_array[15] = 0.0006   # 15 C5a (ic initial)

data_dictionary["INITIAL_CONDITION_ARRAY"] = initial_condition_array;
number_parameters = 0

for (index,rank_index_value) in enumerate(idx_rank)

  if (mod(index,10)==0)

    # Grab the parameter set from the cache -
    parameter_array = pc_array_full[:,rank_index_value]
    data_dictionary["PARAMETER_ARRAY"] = parameter_array

    number_parameters = length(parameter_array)

    # Run the model -
    (t_base,x_base) = SolveBalances(tStart,tStop,tStep,data_dictionary)

    # Initialize the distance array -
    distance_array = zeros(length(parameter_array),length(parameter_array))

    # For each parameter combination, we need to measured the distance between the systems states -
    for (parameter_index_outer,parameter_value_outer) in enumerate(parameter_array)
      for (parameter_index_inner,parameter_value_inner) in enumerate(parameter_array)

        # Update the parameters -
        copy_of_parameter_array = copy(parameter_array)
        copy_of_parameter_array[parameter_index_outer] = parameter_value_outer*(1.10)
        copy_of_parameter_array[parameter_index_inner] = parameter_value_inner*(1.10)

        # Update the data dictionary -
        data_dictionary["PARAMETER_ARRAY"] = copy_of_parameter_array

        # Run the model -
        (t_perturbed,x_pertubed) = SolveBalances(tStart,tStop,tStep,data_dictionary)

        # measure the distance -
        local_distance = norm(norm(x_base - x_pertubed))

        # add this value to the distance_array -
        distance_array[parameter_index_outer,parameter_index_inner] = local_distance

        @show (parameter_index_outer,parameter_index_inner,index)
      end # end inner -
    end # end outer -

    # When I get here, I should have a full distance array -
    push!(distance_array_set,distance_array)

  end # end the mod -
end # end rank iteration -

# write -
write_array_set_to_disk(copy(distance_array_set),"./clustergram_data","DArray-control")

# calculate the array -
mean_array = calculate_mean_array_from_set(copy(distance_array_set),number_parameters,number_parameters)


text:
#=
Cairo drawing helper functions
=#

#=
Backgrounds
=#
function black_bg(ctx::CairoContext)
    set_source_rgb(ctx,0,0,0)
    paint(ctx)
    nothing
end

function white_bg(ctx::CairoContext)
    set_source_rgb(ctx,1,1,1)
    paint(ctx)
    nothing
end

#=
Draw Shapes
=#
function draw_circle(xi::Float64,xnew::Float64,yi::Float64,ynew::Float64,ctx::CairoContext)
    draw_circle(xi,xnew,yi,ynew,50.0,ctx)
end

function draw_circle(xi::Float64,xnew::Float64,yi::Float64,ynew::Float64,r::Float64,ctx::CairoContext)

    set_line_width(ctx, 9)
    set_source_rgb(ctx, 1.0, 1.0, 1.0)
    translate(ctx, (xnew-xi)/2, (ynew-yi)/2)
    arc(ctx, 0, 0, r, 0, 2 * pi)
    stroke_preserve(ctx)
    set_source_rgb(ctx, 1.0, 1.0, 1.0)
    fill(ctx)
    nothing
end


text:
struct TopicModelError <: Exception
    msg::String
end

Base.showerror(io::IO, e::TopicModelError) = print(io, "TopicModelError: ", e.msg)

showdocs(model::TopicModel, doc::Document) = showdocs(model.corp, doc)
showdocs(model::TopicModel, docs::Vector{Document}) = showdocs(model.corp, docs)
showdocs(model::TopicModel, d::Integer) = showdocs(model.corp, d)
showdocs(model::TopicModel, doc_indices::Vector{<:Integer}) = showdocs(model.corp, doc_indices)
showdocs(model::TopicModel, doc_range::UnitRange{<:Integer}) = showdocs(model.corp, collect(doc_range))
showdocs(model::TopicModel) = showdocs(model.corp)

showtitles(model::TopicModel, doc::Document) = showtitles(model.corp, doc)
showtitles(model::TopicModel, docs::Vector{Document}) = showtitles(model.corp, docs)
showtitles(model::TopicModel, d::Integer) = showtitles(model.corp, d)
showtitles(model::TopicModel, doc_indices::Vector{<:Integer}) = showtitles(model.corp, doc_indices)
showtitles(model::TopicModel, doc_range::UnitRange{<:Integer}) = showtitles(model.corp, collect(doc_range))
showtitles(model::TopicModel) = showtitles(model.corp)

getvocab(model::TopicModel) = getvocab(model.corp)
getusers(model::TopicModel) = getusers(model.corp)

### Display output for TopicModel objects.
Base.show(io::IO, model::LDA) = print(io, "Latent Dirichlet allocation model with $(model.K) topics.")
Base.show(io::IO, model::fLDA) = print(io, "Filtered latent Dirichlet allocation model with $(model.K) topics.")
Base.show(io::IO, model::CTM) = print(io, "Correlated topic model with $(model.K) topics.")
Base.show(io::IO, model::fCTM) = print(io, "Filtered correlated topic model with $(model.K) topics.")
Base.show(io::IO, model::CTPF) = print(io, "Collaborative topic Poisson factorization model with $(model.K) topics.")
Base.show(io::IO, model::gpuLDA) = print(io, "GPU accelerated latent Dirichlet allocation model with $(model.K) topics.")
Base.show(io::IO, model::gpuCTM) = print(io, "GPU accelerated correlated topic model with $(model.K) topics.")
Base.show(io::IO, model::gpuCTPF) = print(io, "GPU accelerated collaborative topic Poisson factorization model with $(model.K) topics.")

function check_model(model::LDA)
	"Check LDA parameters."

	check_corp(model.corp)
	isequal(collect(1:model.V), sort(collect(keys(model.corp.vocab))))				|| throw(TopicModelError("Corpus vocab keys must form unit range of length V."))
	isequal(model.M, length(model.corp))											|| throw(TopicModelError("M must equal the number of documents in the corpus."))
	isequal(model.N, [length(doc.terms) for doc in model.corp])						|| throw(TopicModelError("N must contain document lengths."))
	isequal(model.C, [sum(doc.counts) for doc in model.corp])						|| throw(TopicModelError("C must contain sums of document counts."))
	isequal(length(model.alpha), model.K)											|| throw(TopicModelError("alpha must be of length K."))
	all(isfinite.(model.alpha))														|| throw(TopicModelError("alpha must be finite."))
	all(model.alpha .> 0)															|| throw(TopicModelError("alpha must be positive."))
	isequal(size(model.beta), (model.K, model.V))									|| throw(TopicModelError("beta must be of size (K, V)."))
	(isstochastic(model.beta, dims=2) | isempty(model.beta))						|| throw(TopicModelError("beta must be a right stochastic matrix."))
	isequal(size(model.beta_old), (model.K, model.V))								|| throw(TopicModelError("beta_old must be of size (K, V)."))
	(isstochastic(model.beta_old, dims=2) | isempty(model.beta_old))				|| throw(TopicModelError("beta_old must be a right stochastic matrix."))
	isequal(model.beta_temp, zeros(model.K, model.V))								|| throw(TopicModelError("beta_temp must be a zero matrix of size (K, V)."))
	isequal(length(model.Elogtheta), model.M)										|| throw(TopicModelError("Elogtheta must be of length M."))
	all(Bool[isequal(length(model.Elogtheta[d]), model.K) for d in 1:model.M])		|| throw(TopicModelError("Elogtheta must contain vectors of length K."))
	all(Bool[all(isfinite.(model.Elogtheta[d])) for d in 1:model.M])				|| throw(TopicModelError("Elogtheta must be finite."))
	all(Bool[all(model.Elogtheta[d] .<= 0) for d in 1:model.M])						|| throw(TopicModelError("Elogtheta must be nonpositive."))
	isequal(length(model.Elogtheta_old), model.M)									|| throw(TopicModelError("Elogtheta_old must be of length M."))
	all(Bool[isequal(length(model.Elogtheta_old[d]), model.K) for d in 1:model.M])	|| throw(TopicModelError("Elogtheta_old must contain vectors of length K."))
	all(Bool[all(isfinite.(model.Elogtheta_old[d])) for d in 1:model.M])			|| throw(TopicModelError("Elogtheta_old must be finite."))
	all(Bool[all(model.Elogtheta_old[d] .<= 0) for d in 1:model.M])					|| throw(TopicModelError("Elogtheta_old must be nonpositive."))
	isequal(length(model.gamma), model.M)											|| throw(TopicModelError("gamma must be of length M."))
	all(Bool[isequal(length(model.gamma[d]), model.K) for d in 1:model.M])			|| throw(TopicModelError("gamma must contain vectors of length K."))
	all(Bool[all(isfinite.(model.gamma[d])) for d in 1:model.M])					|| throw(TopicModelError("gamma must be finite."))
	all(Bool[all(model.gamma[d] .> 0) for d in 1:model.M])							|| throw(TopicModelError("gamma must be positive."))
	isfinite(model.elbo)															|| throw(TopicModelError("elbo must be finite."))
	nothing
end

function check_model(model::fLDA)
	"Check fLDA parameters."

	check_corp(model.corp)
	isequal(collect(1:model.V), sort(collect(keys(model.corp.vocab))))				|| throw(TopicModelError("Corpus vocab keys must form unit range of length V."))
	isequal(model.M, length(model.corp))											|| throw(TopicModelError("M must equal the number of documents in the corpus."))
	isequal(model.N, [length(doc.terms) for doc in model.corp])						|| throw(TopicModelError("N must contain document lengths."))
	isequal(model.C, [sum(doc.counts) for doc in model.corp])						|| throw(TopicModelError("C must contain sums of document counts."))
	(0 <= model.eta <= 1)															|| throw(TopicModelError("eta must belong to the interval [0,1]."))
	isequal(length(model.alpha), model.K)											|| throw(TopicModelError("alpha must be of length K."))
	all(isfinite.(model.alpha))														|| throw(TopicModelError("alpha must be finite."))
	all(model.alpha .> 0)															|| throw(TopicModelError("alpha must be positive."))
	isequal(length(model.kappa), model.V)											|| throw(TopicModelError("kappa must be of length V"))
	(isprobvec(model.kappa) | isempty(model.kappa))									|| throw(TopicModelError("kappa must be a probability vector."))
	isequal(length(model.kappa_old), model.V)										|| throw(TopicModelError("kappa_old must be of length V."))
	(isprobvec(model.kappa_old) | isempty(model.kappa_old))							|| throw(TopicModelError("kappa_old must be a probability vector."))
	isequal(model.kappa_temp, zeros(model.V))										|| throw(TopicModelError("kappa_temp must be a zero vector of length V."))
	isequal(size(model.beta), (model.K, model.V))									|| throw(TopicModelError("beta must be of size (K, V)."))
	(isstochastic(model.beta, dims=2) | isempty(model.beta))						|| throw(TopicModelError("beta must be a right stochastic matrix."))
	isequal(size(model.beta_old), (model.K, model.V))								|| throw(TopicModelError("beta_old must be of size (K, V)."))
	(isstochastic(model.beta_old, dims=2) | isempty(model.beta_old))				|| throw(TopicModelError("beta_old must be a right stochastic matrix."))
	isequal(model.beta_temp, zeros(model.K, model.V))								|| throw(TopicModelError("beta_temp must be a zero matrix of size (K, V)."))
	isequal(length(model.Elogtheta), model.M)										|| throw(TopicModelError("Elogtheta must be of length M."))
	all(Bool[isequal(length(model.Elogtheta[d]), model.K) for d in 1:model.M])		|| throw(TopicModelError("Elogtheta must contain vectors of length K."))
	all(Bool[all(isfinite.(model.Elogtheta[d])) for d in 1:model.M])				|| throw(TopicModelError("Elogtheta must be finite."))
	all(Bool[all(model.Elogtheta[d] .<= 0) for d in 1:model.M])						|| throw(TopicModelError("Elogtheta must be nonpositive."))
	isequal(length(model.Elogtheta_old), model.M)									|| throw(TopicModelError("Elogtheta_old must be of length M."))
	all(Bool[isequal(length(model.Elogtheta_old[d]), model.K) for d in 1:model.M])	|| throw(TopicModelError("Elogtheta_old must contain vectors of length K."))
	all(Bool[all(isfinite.(model.Elogtheta_old[d])) for d in 1:model.M])			|| throw(TopicModelError("Elogtheta_old must be finite."))
	all(Bool[all(model.Elogtheta_old[d] .<= 0) for d in 1:model.M])					|| throw(TopicModelError("Elogtheta_old must be nonpositive."))
	isequal(length(model.gamma), model.M)											|| throw(TopicModelError("gamma must be of length M."))
	all(Bool[isequal(length(model.gamma[d]), model.K) for d in 1:model.M])			|| throw(TopicModelError("gamma must contain vectors of length K."))
	all(Bool[all(isfinite.(model.gamma[d])) for d in 1:model.M])					|| throw(TopicModelError("gamma must be finite."))
	all(Bool[all(model.gamma[d] .> 0) for d in 1:model.M])							|| throw(TopicModelError("gamma must be positive."))
	isequal(length(model.tau), model.M)												|| throw(TopicModelError("tau must be of length M."))
	all(Bool[isequal(length(model.tau[d]), model.N[d]) for d in 1:model.M])			|| throw(TopicModelError("tau must contain vectors of lengths N."))
	all(Bool[all(0 .<= model.tau[d] .<= 1) for d in 1:model.M])						|| throw(TopicModelError("tau must belong to the interval [0,1]."))
	isfinite(model.elbo)															|| throw(TopicModelError("elbo must be finite."))
	nothing
end

function check_model(model::CTM)
	"Check CTM parameters."

	check_corp(model.corp)
	isequal(collect(1:model.V), sort(collect(keys(model.corp.vocab))))			|| throw(TopicModelError("Corpus vocab keys must form unit range of length V."))
	isequal(model.M, length(model.corp))										|| throw(TopicModelError("M must equal the number of documents in the corpus."))
	isequal(model.N, [length(doc.terms) for doc in model.corp])					|| throw(TopicModelError("N must contain document lengths."))
	isequal(model.C, [sum(doc.counts) for doc in model.corp])					|| throw(TopicModelError("C must contain sums of document counts."))
	all(isfinite.(model.mu))													|| throw(TopicModelError("mu must be finite."))
	isequal(size(model.sigma), (model.K, model.K))								|| throw(TopicModelError("sigma must be of size (K, K)."))
	isposdef(model.sigma)														|| throw(TopicModelError("sigma must be positive-definite."))
	isequal(size(model.invsigma), (model.K, model.K))							|| throw(TopicModelError("invsigma must be of size (K, K)."))
	isposdef(model.invsigma)													|| throw(TopicModelError("invsigma must be positive-definite."))
	isequal(size(model.beta), (model.K, model.V))								|| throw(TopicModelError("beta must be of size (K, V)."))
	(isstochastic(model.beta, dims=2) | isempty(model.beta))					|| throw(TopicModelError("beta must be a right stochastic matrix."))
	isequal(size(model.beta_old), (model.K, model.V))							|| throw(TopicModelError("beta_old must be of size (K, V)."))
	(isstochastic(model.beta_old, dims=2) | isempty(model.beta_old))			|| throw(TopicModelError("beta_old must be a right stochastic matrix."))
	isequal(model.beta_temp, zeros(model.K, model.V))							|| throw(TopicModelError("beta_temp must be a zero matrix of size (K, V)."))
	isequal(length(model.lambda), model.M)										|| throw(TopicModelError("lambda must be of length M."))
	all(Bool[isequal(length(model.lambda[d]), model.K) for d in 1:model.M])		|| throw(TopicModelError("lambda must contain vectors of length K."))
	all(Bool[all(isfinite.(model.lambda[d])) for d in 1:model.M])				|| throw(TopicModelError("lambda must be finite."))
	isequal(length(model.lambda_old), model.M)									|| throw(TopicModelError("lambda_old must be of length M."))
	all(Bool[isequal(length(model.lambda_old[d]), model.K) for d in 1:model.M])	|| throw(TopicModelError("lambda_old must contain vectors of length K."))
	all(Bool[all(isfinite.(model.lambda_old[d])) for d in 1:model.M])			|| throw(TopicModelError("lambda_old must be finite."))
	isequal(length(model.vsq), model.M)											|| throw(TopicModelError("vsq must be of length M."))
	all(Bool[isequal(length(model.vsq[d]), model.K) for d in 1:model.M])		|| throw(TopicModelError("vsq must contain vectors of length K."))
	all(Bool[all(isfinite.(model.vsq[d])) for d in 1:model.M])					|| throw(TopicModelError("vsq must be finite."))
	all(Bool[all(model.vsq[d] .> 0) for d in 1:model.M])						|| throw(TopicModelError("vsq must be positive."))
	isequal(length(model.logzeta), model.M)										|| throw(TopicModelError("logzeta must be of length M."))
	all(isfinite.(model.logzeta))												|| throw(TopicModelError("logzeta must be finite."))
	isfinite(model.elbo)														|| throw(TopicModelError("elbo must be finite."))
	nothing
end

function check_model(model::fCTM)
	"Check fCTM parameters."

	check_corp(model.corp)
	isequal(collect(1:model.V), sort(collect(keys(model.corp.vocab))))			|| throw(TopicModelError("Corpus vocab keys must form unit range of length V."))
	isequal(model.M, length(model.corp))										|| throw(TopicModelError("M must equal the number of documents in the corpus."))
	isequal(model.N, [length(doc.terms) for doc in model.corp])					|| throw(TopicModelError("N must contain document lengths."))
	isequal(model.C, [sum(doc.counts) for doc in model.corp])					|| throw(TopicModelError("C must contain sums of document counts."))
	(0 <= model.eta <= 1)														|| throw(TopicModelError("eta must belong to the interval [0,1]."))
	all(isfinite.(model.mu))													|| throw(TopicModelError("mu must be finite."))
	isequal(size(model.sigma), (model.K, model.K))								|| throw(TopicModelError("sigma must be of size (K, K)."))
	isposdef(model.sigma)														|| throw(TopicModelError("sigma must be positive-definite."))
	isequal(size(model.invsigma), (model.K, model.K))							|| throw(TopicModelError("invsigma must be of size (K, K)."))
	isposdef(model.invsigma)													|| throw(TopicModelError("invsigma must be positive-definite."))
	isequal(length(model.kappa), model.V)										|| throw(TopicModelError("kappa must be of length V."))
	(isprobvec(model.kappa) | isempty(model.kappa))								|| throw(TopicModelError("kappa must be a probability vector."))
	isequal(length(model.kappa_old), model.V)									|| throw(TopicModelError("kappa_old must be of length V."))
	(isprobvec(model.kappa_old) | isempty(model.kappa_old))						|| throw(TopicModelError("kappa_old must be a probability vector."))
	isequal(model.kappa_temp, zeros(model.V))									|| throw(TopicModelError("kappa_temp must be zero vector of length V."))
	isequal(size(model.beta), (model.K, model.V))								|| throw(TopicModelError("beta must be of size (K, V)."))
	(isstochastic(model.beta, dims=2) | isempty(model.beta))					|| throw(TopicModelError("beta must be a right stochastic matrix."))
	isequal(size(model.beta_old), (model.K, model.V))							|| throw(TopicModelError("beta_old must be of size (K, V)."))
	(isstochastic(model.beta_old, dims=2) | isempty(model.beta_old))			|| throw(TopicModelError("beta_old must be a right stochastic matrix."))
	isequal(model.beta_temp, zeros(model.K, model.V))							|| throw(TopicModelError("beta_temp must be a zero matrix of size (K, V)."))
	isequal(length(model.lambda), model.M)										|| throw(TopicModelError("lambda must be of length M."))
	all(Bool[isequal(length(model.lambda[d]), model.K) for d in 1:model.M])		|| throw(TopicModelError("lambda must contain vectors of length K."))
	all(Bool[all(isfinite.(model.lambda[d])) for d in 1:model.M])				|| throw(TopicModelError("lambda must be finite."))
	isequal(length(model.lambda_old), model.M)									|| throw(TopicModelError("lambda_old must be of length M."))
	all(Bool[isequal(length(model.lambda_old[d]), model.K) for d in 1:model.M])	|| throw(TopicModelError("lambda_old must contain vectors of length K."))
	all(Bool[all(isfinite.(model.lambda_old[d])) for d in 1:model.M])			|| throw(TopicModelError("lambda_old must be finite."))
	isequal(length(model.vsq), model.M)											|| throw(TopicModelError("vsq must be of length M."))
	all(Bool[isequal(length(model.vsq[d]), model.K) for d in 1:model.M])		|| throw(TopicModelError("vsq must contain vectors of length K."))
	all(Bool[all(isfinite.(model.vsq[d])) for d in 1:model.M])					|| throw(TopicModelError("vsq must be finite."))
	all(Bool[all(model.vsq[d] .> 0) for d in 1:model.M])						|| throw(TopicModelError("vsq must be positive."))
	isequal(length(model.logzeta), model.M)										|| throw(TopicModelError("logzeta must be of length M."))
	all(isfinite.(model.logzeta))												|| throw(TopicModelError("logzeta must be finite."))
	isequal(length(model.tau), model.M)											|| throw(TopicModelError("tau must be of length M."))
	all(Bool[isequal(length(model.tau[d]), model.N[d]) for d in 1:model.M])		|| throw(TopicModelError("tau must contain vectors of lengths N."))
	all(Bool[all(0 .<= model.tau[d] .<= 1) for d in 1:model.M])					|| throw(TopicModelError("tau must belong to the interval [0,1]."))
	isfinite(model.elbo)														|| throw(TopicModelError("elbo must be finite."))
	nothing
end

function check_model(model::CTPF)
	"Check CTPF parameters."

	check_corp(model.corp)
	isequal(collect(1:model.V), sort(collect(keys(model.corp.vocab))))			|| throw(TopicModelError("Corpus vocab keys must form unit range of length V."))
	isequal(collect(1:model.U), sort(collect(keys(model.corp.users))))			|| throw(TopicModelError("Corpus users keys must form unit range of length U."))
	isequal(model.M, length(model.corp))										|| throw(TopicModelError("M must be equal to the number of documents in the corpus."))
	isequal(model.N, [length(model.corp[d].terms) for d in 1:model.M])			|| throw(TopicModelError("N must contain document lengths"))
	isequal(model.C, [sum(model.corp[d].counts) for d in 1:model.M])			|| throw(TopicModelError("C must contain sums of document counts."))
	isequal(model.R, [length(model.corp[d].readers) for d in 1:model.M])		|| throw(TopicModelError("R must contain numbers of readers in documents."))
	model.a > 0																	|| throw(TopicModelError("a must be positive."))
	model.b > 0																	|| throw(TopicModelError("b must be positive."))
	model.c > 0																	|| throw(TopicModelError("c must be positive."))
	model.d > 0																	|| throw(TopicModelError("d must be positive."))
	model.e > 0																	|| throw(TopicModelError("e must be positive."))
	model.f > 0																	|| throw(TopicModelError("f must be positive."))
	model.g > 0																	|| throw(TopicModelError("g must be positive."))
	model.h > 0																	|| throw(TopicModelError("h must be positive."))
	isequal(size(model.alef), (model.K, model.V))								|| throw(TopicModelError("alef must be of size (K, V)."))
	all(isfinite.(model.alef))													|| throw(TopicModelError("alef must be finite."))
	all(model.alef .> 0)														|| throw(TopicModelError("alef must be positive."))
	isequal(size(model.alef_old), (model.K, model.V))							|| throw(TopicModelError("alef_old must be of size (K, V)."))
	all(isfinite.(model.alef_old))												|| throw(TopicModelError("alef_old must be finite."))
	all(model.alef_old .> 0)													|| throw(TopicModelError("alef_old must be positive."))
	isequal(model.alef_temp, fill(model.a, model.K, model.V))					|| throw(TopicModelError("alef_temp must be a fill a matrix of size (K, V)."))
	isequal(length(model.bet), model.K)											|| throw(TopicModelError("bet must be of length K"))
	all(isfinite.(model.bet))													|| throw(TopicModelError("bet must be finite."))
	all(model.bet .> 0)															|| throw(TopicModelError("bet must be positive."))
	isequal(length(model.bet_old), model.K)										|| throw(TopicModelError("bet_old must be of length K"))
	all(isfinite.(model.bet_old))												|| throw(TopicModelError("bet_old must be finite."))
	all(model.bet_old .> 0)														|| throw(TopicModelError("bet_old must be positive."))
	isequal(length(model.gimel), model.M)										|| throw(TopicModelError("gimel must be of length M."))
	all(Bool[isequal(length(model.gimel[d]), model.K) for d in 1:model.M])		|| throw(TopicModelError("gimel must contain vectors of length K."))
	all(Bool[all(isfinite.(model.gimel[d])) for d in 1:model.M])				|| throw(TopicModelError("gimel must be finite."))
	all(Bool[all(model.gimel[d] .> 0) for d in 1:model.M])						|| throw(TopicModelError("gimel must be positive."))
	isequal(length(model.gimel_old), model.M)									|| throw(TopicModelError("gimel_old must be of length M."))
	all(Bool[isequal(length(model.gimel_old[d]), model.K) for d in 1:model.M])	|| throw(TopicModelError("gimel_old must contain vectors of length K."))
	all(Bool[all(isfinite.(model.gimel_old[d])) for d in 1:model.M])			|| throw(TopicModelError("gimel_old must be finite."))
	all(Bool[all(model.gimel_old[d] .> 0) for d in 1:model.M])					|| throw(TopicModelError("gimel_old must be positive."))
	isequal(length(model.dalet), model.K)										|| throw(TopicModelError("dalet must be of length K."))
	all(isfinite.(model.dalet))													|| throw(TopicModelError("dalet must be finite."))
	all(model.dalet .> 0)														|| throw(TopicModelError("dalet must be positive."))
	isequal(length(model.dalet_old), model.K)									|| throw(TopicModelError("dalet_old must be of length K."))
	all(isfinite.(model.dalet_old))												|| throw(TopicModelError("dalet_old must be finite."))
	all(model.dalet_old .> 0)													|| throw(TopicModelError("dalet_old must be positive."))
	isequal(size(model.he), (model.K, model.U))									|| throw(TopicModelError("he must be of size (K, U)"))
	all(isfinite.(model.he))													|| throw(TopicModelError("he must be finite."))
	all(model.he .> 0)															|| throw(TopicModelError("he must be positive."))
	isequal(size(model.he_old), (model.K, model.U))								|| throw(TopicModelError("he_old must be of size (K, U)"))
	all(isfinite.(model.he_old))												|| throw(TopicModelError("he_old must be finite."))
	all(model.he_old .> 0)														|| throw(TopicModelError("he_old must be positive."))
	isequal(model.he_temp, fill(model.e, model.K, model.U))						|| throw(TopicModelError("he_temp must be a fill e matrix of size (K, U)."))
	isequal(length(model.vav), model.K)											|| throw(TopicModelError("vav must be of length K."))
	all(isfinite.(model.vav))													|| throw(TopicModelError("vav must be finite."))
	all(model.vav .> 0)															|| throw(TopicModelError("vav must be positive."))
	isequal(length(model.vav_old), model.K)										|| throw(TopicModelError("vav_old must be of length K."))
	all(isfinite.(model.vav_old))												|| throw(TopicModelError("vav_old must be finite."))
	all(model.vav_old .> 0)														|| throw(TopicModelError("vav_old must be positive."))
	isequal(length(model.zayin), model.M)										|| throw(TopicModelError("zayin must be of length M."))
	all(Bool[isequal(length(model.zayin[d]), model.K) for d in 1:model.M])		|| throw(TopicModelError("zayin must contain vectors of length K."))
	all(Bool[all(isfinite.(model.zayin[d])) for d in 1:model.M])				|| throw(TopicModelError("zayin must be finite."))
	all(Bool[all(model.zayin[d] .> 0) for d in 1:model.M])						|| throw(TopicModelError("zayin must be positive."))
	isequal(length(model.zayin_old), model.M)									|| throw(TopicModelError("zayin_old must be of length M."))
	all(Bool[isequal(length(model.zayin_old[d]), model.K) for d in 1:model.M])	|| throw(TopicModelError("zayin_old must contain vectors of length K."))
	all(Bool[all(isfinite.(model.zayin_old[d])) for d in 1:model.M])			|| throw(TopicModelError("zayin_old must be finite."))
	all(Bool[all(model.zayin_old[d] .> 0) for d in 1:model.M])					|| throw(TopicModelError("zayin_old must be positive."))
	isequal(length(model.het), model.K)											|| throw(TopicModelError("het must be of length K."))
	all(isfinite.(model.het))													|| throw(TopicModelError("het must be finite."))
	all(model.het .> 0)															|| throw(TopicModelError("het must be positive."))
	isequal(length(model.het_old), model.K)										|| throw(TopicModelError("het_old must be of length K."))
	all(isfinite.(model.het_old))												|| throw(TopicModelError("het_old must be finite."))
	all(model.het_old .> 0)														|| throw(TopicModelError("het_old must be positive."))
	isfinite(model.elbo)														|| throw(TopicModelError("elbo must be finite."))
	nothing
end

function check_model(model::gpuLDA)
	"Check gpuLDA parameters."

	check_corp(model.corp)
	isequal(collect(1:model.V), sort(collect(keys(model.corp.vocab))))							|| throw(TopicModelError("Corpus vocab keys must form unit range of length V."))
	isequal(model.M, length(model.corp))														|| throw(TopicModelError("M must equal the number of documents in the corpus."))
	isequal(model.N, [length(doc.terms) for doc in model.corp])									|| throw(TopicModelError("N must contain document lengths."))
	isequal(model.C, [sum(doc.counts) for doc in model.corp])									|| throw(TopicModelError("C must contain sums of document counts."))
	isequal(length(model.alpha), model.K)														|| throw(TopicModelError("alpha must be of length K."))
	all(isfinite.(model.alpha))																	|| throw(TopicModelError("alpha must be finite."))
	all(model.alpha .> 0)																		|| throw(TopicModelError("alpha must be positive."))
	isequal(size(model.beta), (model.K, model.V))												|| throw(TopicModelError("beta must be of size (K, V)."))
	(isstochastic(model.beta, dims=2) | isempty(model.beta))									|| throw(TopicModelError("beta must be a right stochastic matrix."))
	isequal(length(model.Elogtheta), model.M)													|| throw(TopicModelError("Elogtheta must be of length M."))
	all(Bool[isequal(length(model.Elogtheta[d]), model.K) for d in 1:model.M])					|| throw(TopicModelError("Elogtheta must contain vectors of length K."))
	all(Bool[all(isfinite.(model.Elogtheta[d])) for d in 1:model.M])							|| throw(TopicModelError("Elogtheta must be finite."))
	all(Bool[all(model.Elogtheta[d] .<= 0) for d in 1:model.M])									|| throw(TopicModelError("Elogtheta must be nonpositive."))
	isequal(length(model.gamma), model.M)														|| throw(TopicModelError("gamma must be of length M."))
	all(Bool[isequal(length(model.gamma[d]), model.K) for d in 1:model.M])						|| throw(TopicModelError("gamma must contain vectors of length K."))
	all(Bool[all(isfinite.(model.gamma[d])) for d in 1:model.M])								|| throw(TopicModelError("gamma must be finite."))
	all(Bool[all(model.gamma[d] .> 0) for d in 1:model.M])										|| throw(TopicModelError("gamma must be positive."))
	isequal(length(model.phi), model.M)															|| throw(TopicModelError("phi must be of length M."))
	all(Bool[isequal(size(model.phi[d]), (model.K, model.N[d])) for d in 1:model.M])			|| throw(TopicModelError("phi must contain matrices of size (K, N)."))
	all(Bool[isstochastic(model.phi[d], dims=1) | isempty(model.phi[d]) for d in 1:model.M])	|| throw(TopicModelError("phi must contain left stochastic matrices."))
	isfinite(model.elbo)																		|| throw(TopicModelError("elbo must be finite."))
	nothing
end

function check_model(model::gpuCTM)
	"Check gpuCTM parameters."

	check_corp(model.corp)
	isequal(collect(1:model.V), sort(collect(keys(model.corp.vocab))))							|| throw(TopicModelError("Corpus vocab keys must form unit range of length V."))
	isequal(model.M, length(model.corp))														|| throw(TopicModelError("M must equal the number of documents in the corpus."))
	isequal(model.N, [length(doc.terms) for doc in model.corp])									|| throw(TopicModelError("N must contain document lengths."))
	isequal(model.C, [sum(doc.counts) for doc in model.corp])									|| throw(TopicModelError("C must contain sums of document counts."))
	all(isfinite.(model.mu))																	|| throw(TopicModelError("mu must be finite."))
	isequal(size(model.sigma), (model.K, model.K))												|| throw(TopicModelError("sigma must be of size (K, K)."))
	isposdef(model.sigma)																		|| throw(TopicModelError("sigma must be positive-definite."))
	isequal(size(model.invsigma), (model.K, model.K))											|| throw(TopicModelError("invsigma must be of size (K, K)."))
	isposdef(model.invsigma)																	|| throw(TopicModelError("invsigma must be positive-definite."))
	isequal(size(model.beta), (model.K, model.V))												|| throw(TopicModelError("beta must be of size (K, V)."))
	(isstochastic(model.beta, dims=2) | isempty(model.beta))									|| throw(TopicModelError("beta must be a right stochastic matrix."))
	isequal(length(model.lambda), model.M)														|| throw(TopicModelError("lambda must be of length M."))
	all(Bool[isequal(length(model.lambda[d]), model.K) for d in 1:model.M])						|| throw(TopicModelError("lambda must contain vectors of length K."))
	all(Bool[all(isfinite.(model.lambda[d])) for d in 1:model.M])								|| throw(TopicModelError("lambda must be finite."))
	isequal(length(model.vsq), model.M)															|| throw(TopicModelError("vsq must be of length M."))
	all(Bool[isequal(length(model.vsq[d]), model.K) for d in 1:model.M])						|| throw(TopicModelError("vsq must contain vectors of length K."))
	all(Bool[all(isfinite.(model.vsq[d])) for d in 1:model.M])									|| throw(TopicModelError("vsq must be finite."))
	all(Bool[all(model.vsq[d] .> 0) for d in 1:model.M])										|| throw(TopicModelError("vsq must be positive."))
	isequal(length(model.logzeta), model.M)														|| throw(TopicModelError("logzeta must be of length M."))
	all(isfinite.(model.logzeta))																|| throw(TopicModelError("logzeta must be finite."))
	isequal(length(model.phi), model.M)															|| throw(TopicModelError("phi must be of length M."))
	all(Bool[isequal(size(model.phi[d]), (model.K, model.N[d])) for d in 1:model.M])			|| throw(TopicModelError("phi must contain matrices of size (K, N)."))
	all(Bool[isstochastic(model.phi[d], dims=1) | isempty(model.phi[d]) for d in 1:model.M])	|| throw(TopicModelError("phi must contain left stochastic matrices."))
	isfinite(model.elbo)																		|| throw(TopicModelError("elbo must be finite."))
	nothing
end

function check_model(model::gpuCTPF)
	"Check gpuCTPF parameters."

	check_corp(model.corp)
	isequal(collect(1:model.V), sort(collect(keys(model.corp.vocab))))					|| throw(TopicModelError("Corpus vocab keys must form unit range of length V."))
	isequal(collect(1:model.U), sort(collect(keys(model.corp.users))))					|| throw(TopicModelError("Corpus users keys must form unit range of length U."))
	isequal(model.M, length(model.corp))												|| throw(TopicModelError("M must be equal to the number of documents in the corpus."))
	isequal(model.N, [length(model.corp[d].terms) for d in 1:model.M])					|| throw(TopicModelError("N must contain document lengths"))
	isequal(model.C, [sum(model.corp[d].counts) for d in 1:model.M])					|| throw(TopicModelError("C must contain sums of document counts."))
	isequal(model.R, [length(model.corp[d].readers) for d in 1:model.M])				|| throw(TopicModelError("R must contain numbers of readers in documents."))
	model.a > 0																			|| throw(TopicModelError("a must be positive."))
	model.b > 0																			|| throw(TopicModelError("b must be positive."))
	model.c > 0																			|| throw(TopicModelError("c must be positive."))
	model.d > 0																			|| throw(TopicModelError("d must be positive."))
	model.e > 0																			|| throw(TopicModelError("e must be positive."))
	model.f > 0																			|| throw(TopicModelError("f must be positive."))
	model.g > 0																			|| throw(TopicModelError("g must be positive."))
	model.h > 0																			|| throw(TopicModelError("h must be positive."))
	isequal(size(model.alef), (model.K, model.V))										|| throw(TopicModelError("alef must be of size (K, V)."))
	all(isfinite.(model.alef))															|| throw(TopicModelError("alef must be finite."))
	all(model.alef .> 0)																|| throw(TopicModelError("alef must be positive."))
	isequal(length(model.bet), model.K)													|| throw(TopicModelError("bet must be of length K"))
	all(isfinite.(model.bet))															|| throw(TopicModelError("bet must be finite."))
	all(model.bet .> 0)																	|| throw(TopicModelError("bet must be positive."))
	isequal(length(model.gimel), model.M)												|| throw(TopicModelError("gimel must be of length M."))
	all(Bool[isequal(length(model.gimel[d]), model.K) for d in 1:model.M])				|| throw(TopicModelError("gimel must contain vectors of length K."))
	all(Bool[all(isfinite.(model.gimel[d])) for d in 1:model.M])						|| throw(TopicModelError("gimel must be finite."))
	all(Bool[all(model.gimel[d] .> 0) for d in 1:model.M])								|| throw(TopicModelError("gimel must be positive."))
	isequal(length(model.dalet), model.K)												|| throw(TopicModelError("dalet must be of length K."))
	all(isfinite.(model.dalet))															|| throw(TopicModelError("dalet must be finite."))
	all(model.dalet .> 0)																|| throw(TopicModelError("dalet must be positive."))
	isequal(size(model.he), (model.K, model.U))											|| throw(TopicModelError("he must be of size (K, U)"))
	all(isfinite.(model.he))															|| throw(TopicModelError("he must be finite."))
	all(model.he .> 0)																	|| throw(TopicModelError("he must be positive."))
	isequal(length(model.vav), model.K)													|| throw(TopicModelError("vav must be of length K."))
	all(isfinite.(model.vav))															|| throw(TopicModelError("vav must be finite."))
	all(model.vav .> 0)																	|| throw(TopicModelError("vav must be positive."))
	isequal(length(model.zayin), model.M)												|| throw(TopicModelError("zayin must be of length M."))
	all(Bool[isequal(length(model.zayin[d]), model.K) for d in 1:model.M])				|| throw(TopicModelError("zayin must contain vectors of length K."))
	all(Bool[all(isfinite.(model.zayin[d])) for d in 1:model.M])						|| throw(TopicModelError("zayin must be finite."))
	all(Bool[all(model.zayin[d] .> 0) for d in 1:model.M])								|| throw(TopicModelError("zayin must be positive."))
	isequal(length(model.het), model.K)													|| throw(TopicModelError("het must be of length K."))
	all(isfinite.(model.het))															|| throw(TopicModelError("het must be finite."))
	all(model.het .> 0)																	|| throw(TopicModelError("het must be positive."))
	isequal(length(model.phi), model.M)													|| throw(TopicModelError("phi must be of length M."))
	all(Bool[isequal(size(model.phi[d]), (model.K, model.N[d])) for d in 1:model.M])	|| throw(TopicModelError("phi must contain matrices of size (K, N)."))
	all(Bool[isstochastic(model.phi[d], dims=1) for d in 1:model.M])					|| throw(TopicModelError("phi must contain left stochastic matrices."))
	isequal(length(model.xi), model.M)													|| throw(TopicModelError("phi must be of length M."))
	all(Bool[isequal(size(model.xi[d]), (2model.K, model.R[d])) for d in 1:model.M])	|| throw(TopicModelError("xi must contain matrices of size (2K, R)."))
	all(Bool[isstochastic(model.xi[d], dims=1) for d in 1:model.M])						|| throw(TopicModelError("xi must contain left stochastic matrices."))
	isfinite(model.elbo)																|| throw(TopicModelError("elbo must be finite"))
	nothing
end

function update_buffer!(model::gpuLDA)
	"Update gpuLDA model data in GPU RAM."

	terms = vcat([doc.terms for doc in model.corp]...) .- 1
	terms_sortperm = sortperm(terms) .- 1
	counts = vcat([doc.counts for doc in model.corp]...)

	J = zeros(Int, model.V)
	for j in terms
		J[j+1] += 1
	end

	N_cumsum = cumsum([0; model.N])
	J_cumsum = cumsum([0; J])

	model.terms_buffer = cl.Buffer(Int, model.context, (:r, :copy), hostbuf=terms)
	model.terms_sortperm_buffer = cl.Buffer(Int, model.context, (:r, :copy), hostbuf=terms_sortperm)
	model.counts_buffer = cl.Buffer(Int, model.context, (:r, :copy), hostbuf=counts)

	model.N_cumsum_buffer = cl.Buffer(Int, model.context, (:r, :copy), hostbuf=N_cumsum)
	model.J_cumsum_buffer = cl.Buffer(Int, model.context, (:r, :copy), hostbuf=J_cumsum)

	model.alpha_buffer = cl.Buffer(Float32, model.context, (:rw, :copy), hostbuf=model.alpha)
	model.beta_buffer = cl.Buffer(Float32, model.context, (:rw, :copy), hostbuf=model.beta)
	model.Elogtheta_buffer = cl.Buffer(Float32, model.context, (:rw, :copy), hostbuf=hcat(model.Elogtheta...))
	model.Elogtheta_sum_buffer = cl.Buffer(Float32, model.context, (:rw, :copy), hostbuf=zeros(Float32, model.K))
	model.Elogtheta_dist_buffer = cl.Buffer(Float32, model.context, (:rw, :copy), hostbuf=zeros(Float32, model.M))
	model.gamma_buffer = cl.Buffer(Float32, model.context, :rw, model.K * model.M)
	model.phi_buffer = cl.Buffer(Float32, model.context, :rw, model.K * sum(model.N))
end

function update_buffer!(model::gpuCTM)
	"Update gpuCTM model data in GPU RAM."

	terms = vcat([doc.terms for doc in model.corp]...) .- 1
	terms_sortperm = sortperm(terms) .- 1
	counts = vcat([doc.counts for doc in model.corp]...)

	J = zeros(Int, model.V)
	for j in terms
		J[j+1] += 1
	end

	N_cumsum = cumsum([0; model.N])
	J_cumsum = cumsum([0; J])

	model.C_buffer = cl.Buffer(Int, model.context, (:r, :copy), hostbuf=model.C)
	model.terms_buffer = cl.Buffer(Int, model.context, (:r, :copy), hostbuf=terms)
	model.terms_sortperm_buffer = cl.Buffer(Int, model.context, (:r, :copy), hostbuf=terms_sortperm)
	model.counts_buffer = cl.Buffer(Int, model.context, (:r, :copy), hostbuf=counts)

	model.N_cumsum_buffer = cl.Buffer(Int, model.context, (:r, :copy), hostbuf=N_cumsum)
	model.J_cumsum_buffer = cl.Buffer(Int, model.context, (:r, :copy), hostbuf=J_cumsum)

	model.sigma_buffer = cl.Buffer(Float32, model.context, (:rw, :copy), hostbuf=Matrix(model.sigma))
	model.invsigma_buffer = cl.Buffer(Float32, model.context, (:r, :copy), hostbuf=Matrix(model.invsigma))
	model.mu_buffer = cl.Buffer(Float32, model.context, (:rw, :copy), hostbuf=model.mu)
	model.beta_buffer = cl.Buffer(Float32, model.context, (:rw, :copy), hostbuf=model.beta)
	model.lambda_buffer = cl.Buffer(Float32, model.context, (:rw, :copy), hostbuf=hcat(model.lambda...))
	model.lambda_dist_buffer = cl.Buffer(Float32, model.context, (:rw, :copy), hostbuf=zeros(Float32, model.M))
	model.vsq_buffer = cl.Buffer(Float32, model.context, (:rw, :copy), hostbuf=hcat(model.vsq...))
	model.logzeta_buffer = cl.Buffer(Float32, model.context, :rw, model.M)
	model.phi_buffer = cl.Buffer(Float32, model.context, :rw, model.K * sum(model.N))
	model.phi_count_buffer = cl.Buffer(Float32, model.context, :rw, model.K * model.M)

	cl.set_arg!(model.lambda_kernel, 15, cl.LocalMem(Float32, model.K))
	cl.set_arg!(model.lambda_kernel, 16, cl.LocalMem(Float32, model.K))
	cl.set_arg!(model.lambda_kernel, 17, cl.LocalMem(Float32, model.K^2))
end

function update_buffer!(model::gpuCTPF)
	"Update gpuCTPF model data in GPU RAM."

	terms = vcat([doc.terms for doc in model.corp]...) .- 1
	terms_sortperm = sortperm(terms) .- 1
	counts = vcat([doc.counts for doc in model.corp]...)

	readers = vcat([doc.readers for doc in model.corp]...) .- 1
	readers_sortperm = sortperm(readers) .- 1
	ratings = vcat([doc.ratings for doc in model.corp]...)

	J = zeros(Int, model.V)
	for j in terms
		J[j+1] += 1
	end

	Y = zeros(Int, model.U)
	for r in readers
		Y[r+1] += 1
	end

	N_cumsum = cumsum([0; model.N])
	J_cumsum = cumsum([0; J])
	R_cumsum = cumsum([0; model.R])
	Y_cumsum = cumsum([0; Y])

	model.terms_buffer = cl.Buffer(Int, model.context, (:r, :copy), hostbuf=terms)
	model.terms_sortperm_buffer = cl.Buffer(Int, model.context, (:r, :copy), hostbuf=terms_sortperm)
	model.counts_buffer = cl.Buffer(Int, model.context, (:r, :copy), hostbuf=counts)
	model.readers_buffer = cl.Buffer(Int, model.context, (:r, :copy), hostbuf=[readers; 0])
	model.ratings_buffer = cl.Buffer(Int, model.context, (:r, :copy), hostbuf=[ratings; 0])
	model.readers_sortperm_buffer = cl.Buffer(Int, model.context, (:r, :copy), hostbuf=[readers_sortperm; 0])
	model.N_cumsum_buffer = cl.Buffer(Int, model.context, (:r, :copy), hostbuf=N_cumsum)
	model.J_cumsum_buffer = cl.Buffer(Int, model.context, (:r, :copy), hostbuf=J_cumsum)
	model.R_cumsum_buffer = cl.Buffer(Int, model.context, (:r, :copy), hostbuf=R_cumsum)
	model.Y_cumsum_buffer = cl.Buffer(Int, model.context, (:r, :copy), hostbuf=Y_cumsum)

	model.alef_buffer = cl.Buffer(Float32, model.context, (:rw, :copy), hostbuf=model.alef)
	model.bet_buffer = cl.Buffer(Float32, model.context, (:rw, :copy), hostbuf=model.bet)
	model.vav_buffer = cl.Buffer(Float32, model.context, (:rw, :copy), hostbuf=model.vav)
	model.gimel_buffer = cl.Buffer(Float32, model.context, (:rw, :copy), hostbuf=hcat(model.gimel...))
	model.zayin_buffer = cl.Buffer(Float32, model.context, (:rw, :copy), hostbuf=hcat(model.zayin...))
	model.dalet_buffer = cl.Buffer(Float32, model.context, (:rw, :copy), hostbuf=model.dalet)
	model.het_buffer = cl.Buffer(Float32, model.context, (:rw, :copy), hostbuf=model.het)
	model.phi_buffer = cl.Buffer(Float32, model.context, :rw, model.K * sum(model.N))

	if model.U > 0
		model.he_buffer = cl.Buffer(Float32, model.context, (:rw, :copy), hostbuf=model.he)
	else
		model.he_buffer = cl.Buffer(Float32, model.context, (:rw, :copy), hostbuf=Float32[0])
	end

	if sum(model.R) > 0
		model.xi_buffer = cl.Buffer(Float32, model.context, :rw, 2 * model.K * sum(model.R))
	else
		model.xi_buffer = cl.Buffer(Float32, model.context, :rw, 1)
	end
end

function update_host!(model::TopicModel)
	nothing
end

function update_host!(model::gpuLDA)
	"Update gpuLDA model data in CPU RAM."

	N_cumsum = zeros(Int, model.M + 1)
	for d in 1:model.M
		N_cumsum[d+1] = N_cumsum[d] + model.N[d]
	end

	model.beta = reshape(cl.read(model.queue, model.beta_buffer), model.K, model.V)
	Elogtheta_host = reshape(cl.read(model.queue, model.Elogtheta_buffer), model.K, model.M)
	model.Elogtheta = [Elogtheta_host[:,d] for d in 1:model.M]
	model.Elogtheta_sum = cl.read(model.queue, model.Elogtheta_sum_buffer)
	model.Elogtheta_dist = cl.read(model.queue, model.Elogtheta_dist_buffer)[1:model.M]
	gamma_host = reshape(cl.read(model.queue, model.gamma_buffer), model.K, model.M)
	model.gamma = [gamma_host[:,d] for d in 1:model.M]
	phi_host = reshape(cl.read(model.queue, model.phi_buffer), model.K, sum(model.N))
	model.phi = [phi_host[:,N_cumsum[d]+1:N_cumsum[d+1]] for d in 1:model.M]
end

function update_host!(model::gpuCTM)
	"Update gpuCTM model data in CPU RAM."

	N_cumsum = zeros(Int, model.M + 1)
	for d in 1:model.M
		N_cumsum[d+1] = N_cumsum[d] + model.N[d]
	end

	model.mu = cl.read(model.queue, model.mu_buffer)
	model.sigma = Symmetric(reshape(cl.read(model.queue, model.sigma_buffer), model.K, model.K))
	model.beta = reshape(cl.read(model.queue, model.beta_buffer), model.K, model.V)
	lambda_host = reshape(cl.read(model.queue, model.lambda_buffer), model.K, model.M)
	model.lambda = [lambda_host[:,d] for d in 1:model.M]
	model.lambda_dist = cl.read(model.queue,  model.lambda_dist_buffer)[1:model.M]
	vsq_host = reshape(cl.read(model.queue, model.vsq_buffer), model.K, model.M)
	model.vsq = [vsq_host[:,d] for d in 1:model.M]
	model.logzeta = cl.read(model.queue, model.logzeta_buffer)[1:model.M]
	phi_host = reshape(cl.read(model.queue, model.phi_buffer), model.K, sum(model.N))
	model.phi = [phi_host[:,N_cumsum[d]+1:N_cumsum[d+1]] for d in 1:model.M]
end

function update_host!(model::gpuCTPF)
	"Update gpuCTPF model data in CPU RAM."

	N_cumsum = zeros(Int, model.M + 1)
	for d in 1:model.M
		N_cumsum[d+1] = N_cumsum[d] + model.N[d]
	end

	R_cumsum = zeros(Int, model.M + 1)
	for d in 1:model.M
		R_cumsum[d+1] = R_cumsum[d] + model.R[d]
	end

	model.alef = reshape(cl.read(model.queue, model.alef_buffer), model.K, model.V)
	model.bet = cl.read(model.queue, model.bet_buffer)
	model.vav = cl.read(model.queue, model.vav_buffer)
	gimel_host = reshape(cl.read(model.queue, model.gimel_buffer), model.K, model.M)
	model.gimel = [gimel_host[:,d] for d in 1:model.M]
	zayin_host = reshape(cl.read(model.queue, model.zayin_buffer), model.K, model.M)
	model.zayin = [zayin_host[:,d] for d in 1:model.M]
	model.dalet = cl.read(model.queue, model.dalet_buffer)
	model.het = cl.read(model.queue, model.het_buffer)
	phi_host = reshape(cl.read(model.queue, model.phi_buffer), model.K, sum(model.N))
	model.phi = [phi_host[:,N_cumsum[d]+1:N_cumsum[d+1]] for d in 1:model.M]

	if model.U > 0
		model.he = reshape(cl.read(model.queue, model.he_buffer), model.K, model.U)
	end

	if sum(model.R) > 0
		xi_host = reshape(cl.read(model.queue, model.xi_buffer), 2 * model.K, sum(model.R))
		model.xi = [xi_host[:,R_cumsum[d]+1:R_cumsum[d+1]] for d in 1:model.M]
	end
end

function check_elbo!(model::TopicModel, checkelbo::Real, printelbo::Bool, k::Int, tol::Real)
	"Check and print value of delta_elbo."
	"If abs(delta_elbo) < tol, terminate algorithm."

	if k % checkelbo == 0
		update_host!(model)
		delta_elbo = -(model.elbo - update_elbo!(model))
		printelbo && println(k, " ∆elbo: ", round(delta_elbo, digits=3))

		if delta_elbo < tol
			return true
		end
	end
	false
end

function gendoc(model::Union{LDA, gpuLDA, fLDA}, laplace_smooth::Real=0.0)
	"Generate artificial document from LDA or gpuLDA generative model."
	"laplace_smooth governs the amount of Laplace smoothing applied to the topic-term distribution."

	(laplace_smooth >= 0) || throw(ArgumentError("laplace_smooth parameter must be nonnegative."))

	C = rand(Poisson(mean(model.C)))
	theta = rand(Dirichlet(model.alpha))

	topic_dist = Categorical(theta)
	vocab_dist = [Categorical((model.beta[i,:] .+ laplace_smooth) / (1 + laplace_smooth * model.V)) for i in 1:model.K]

	term_count = Dict{Int, Int}()
	for _ in 1:C
		z = rand(topic_dist)
		w = rand(vocab_dist[z])
		haskey(term_count, w) ? term_count[w] += 1 : term_count[w] = 1
	end

	doc = Document(terms=collect(keys(term_count)), counts=collect(values(term_count)))
	return doc
end

function gendoc(model::Union{CTM, gpuCTM, fCTM}, laplace_smooth::Real=0.0)
	"Generate artificial document from CTM or gpuCTM generative model."
	"laplace_smooth governs the amount of Laplace smoothing applied to the topic-term distribution."

	(laplace_smooth >= 0) || throw(ArgumentError("laplace_smooth parameter must be nonnegative."))

	C = rand(Poisson(mean(model.C)))
	theta = additive_logistic(rand(MvNormal(model.mu, model.sigma)))

	topic_dist = Categorical(theta)
	vocab_dist = [Categorical((model.beta[i,:] .+ laplace_smooth) / (1 + laplace_smooth * model.V)) for i in 1:model.K]

	term_count = Dict{Int, Int}()
	for _ in 1:C
		z = rand(topicdist)
		w = rand(vocab_dist[z])
		haskey(term_count, w) ? term_count[w] += 1 : term_count[w] = 1
	end

	doc = Document(terms=collect(keys(term_count)), counts=collect(values(term_count)))
	return doc
end

function gencorp(model::TopicModel, M::Integer; laplace_smooth::Real=0.0)
	"Generate artificial corpus using specified generative model."
	"laplace_smooth governs the amount of Laplace smoothing applied to the topic-term distribution."

	(M > 0)					|| throw(ArgumentError("corp_size parameter must be a positive integer."))
	(laplace_smooth >= 0)	|| throw(ArgumentError("laplace_smooth parameter must be nonnegative."))

	corp = Corpus(vocab=model.corp.vocab, users=model.corp.users)
	corp.docs = [gendoc(model, laplace_smooth) for d in 1:M]
	return corp
end

function showtopics(model::TopicModel, V::Integer=15; topics::Union{<:Integer, Vector{<:Integer}, UnitRange{<:Integer}}=1:model.K, cols::Integer=4)
	"Display the top T terms for each topic."
	"topics parameter controls which topics are displayed."
	"cols parameter controls the number of topic columns displayed per line."

	(V > 0)									|| throw(ArgumentError("Number of displayed terms must be a positive integer."))
	checkbounds(Bool, 1:model.K, topics)	|| throw(ArgumentError("Some topic indices are outside range."))
	(cols > 0)								|| throw(ArgumentError("cols must be a positive integer."))
	V = min(V, model.V)
	cols = min(cols, length(topics))

	maxjspacings = [maximum([length(model.corp.vocab[j]) for j in topic[1:V]]) for topic in model.topics]
	topic_blocks = Iterators.partition(topics, cols)

	for (n, topic_block) in enumerate(topic_blocks)
		for j in 0:V
			for (k, i) in enumerate(topic_block)
				if j == 0
					jspacing = max(4, maxjspacings[i] - length("$i") - 2)
					k == cols ? print(Crayon(foreground=:yellow, bold=true), "topic $i") : print(Crayon(foreground=:yellow, bold=true), "topic $i" * " "^jspacing)
				else
					jspacing = max(6 + length("$i"), maxjspacings[i]) - length(model.corp.vocab[model.topics[i][j]]) + 4
					k == cols ? print(Crayon(foreground=:white, bold=false), model.corp.vocab[model.topics[i][j]]) : print(Crayon(foreground=:white, bold=false), model.corp.vocab[model.topics[i][j]] * " "^jspacing)
				end
			end
			println()
		end

		if n < length(topic_blocks)
			println()
		end
	end
end

function showlibs(model::Union{CTPF, gpuCTPF}, users::Vector{<:Integer})
	"Display the documents in a user(s) library."

	checkbounds(Bool, 1:model.U, users) || throw(ArgumentError("Some user indices are outside range."))

	for (n, u) in enumerate(users)
		if isempty(model.libs[u])
			continue
		end

		@juliadots "User $u\n"
		try
			if model.corp.users[u][1:5] != "#user"
				@juliadots model.corp.users[u] * "\n"
			end

		catch
			@juliadots model.corp.users[u] * "\n"
		end

		for d in model.libs[u]
			print(Crayon(foreground=:yellow, bold=true), " • ")
			isempty(model.corp[d].title) ? println(Crayon(foreground=:white, bold=false), "Document $d") : println(Crayon(foreground=:white, bold=false), "$(model.corp[d].title)")
		end

		if n < length(users)
			println()
		end
	end
end

showlibs(model::Union{CTPF, gpuCTPF}, user::Integer) = showlibs(model, [user])
showlibs(model::Union{CTPF, gpuCTPF}, user_range::UnitRange{<:Integer}) = showlibs(model, collect(user_range))
showlibs(model::Union{CTPF, gpuCTPF}) = showlibs(model, 1:length(model.libs))

function showdrecs(model::Union{CTPF, gpuCTPF}, docs::Vector{<:Integer}, U::Integer=15; cols::Integer=1)
	"Display the top U user recommendations for a document(s)."
	"cols parameter controls the number of topic columns displayed per line."

	checkbounds(Bool, 1:model.M, docs) 	|| throw(ArgumentError("Some document indices are outside range."))
	(U > 0) 							|| throw(ArgumentError("Number of displayed users must be a positive integer."))
	(cols > 0)							|| throw(ArgumentError("cols must be a positive integer."))
	U = min(U, model.U)

	for (n, d) in enumerate(docs)
		if isempty(model.drecs[d])
			continue
		end

		@juliadots "Document $d\n"
		if !isempty(model.corp[d].title)
			@juliadots model.corp[d].title * "\n"
		end

		usercols = collect(Iterators.partition(model.drecs[d][1:U], Int(ceil(U / cols))))
		rankcols = collect(Iterators.partition(1:U, Int(ceil(U / cols))))

		for i in 1:length(usercols[1])
			for j in 1:length(usercols)
				try
					uspacing = maximum([length(model.corp.users[u]) for u in usercols[j]]) - length(model.corp.users[usercols[j][i]]) + 4
					rspacing = maximum([length("$r") for r in rankcols[j]]) - length(string(rankcols[j][i]))
					print(Crayon(foreground=:yellow, bold=true), string(rankcols[j][i]) * ". " * " "^rspacing)
					j == length(usercols) ? print(Crayon(foreground=:white, bold=false), model.corp.users[usercols[j][i]]) : print(Crayon(foreground=:white, bold=false), users[usercols[j][i]] * " "^uspacing)

				catch
					nothing
				end
			end
			println()
		end

		if n < length(docs)
			println()
		end
	end
end

showdrecs(model::Union{CTPF, gpuCTPF}, doc::Integer, U::Integer=15; cols::Integer=1) = showdrecs(model, [doc], U, cols=cols)
showdrecs(model::Union{CTPF, gpuCTPF}, docs::UnitRange{<:Integer}, U::Integer=15; cols::Integer=1) = showdrecs(model, collect(docs), U, cols=cols)

function showurecs(model::Union{CTPF, gpuCTPF}, users::Vector{<:Integer}, M::Integer=15; cols::Integer=1)
	"# Show the top 'M' document recommendations for a user(s)."
	"If a document has no title, the document's index in the corpus will be shown instead."

	checkbounds(Bool, 1:model.U, users) || throw(ArgumentError("Some user indices are outside range."))
	(M > 0) 							|| throw(ArgumentError("Number of displayed documents must be a positive integer."))
	(cols > 0)							|| throw(ArgumentError("cols must be a positive integer."))
	M = min(M, model.M)

	for (n, u) in enumerate(users)
		if isempty(model.urecs[u])
			continue
		end

		@juliadots "User $u\n"
		try
			if model.corp.users[u][1:5] != "#user"
				@juliadots model.corp.users[u] * "\n"
			end

		catch
			@juliadots model.corp.users[u] * "\n"
		end

		docucols = collect(Iterators.partition(model.urecs[u][1:M], Int(ceil(M / cols))))
		rankcols = collect(Iterators.partition(1:M, Int(ceil(M / cols))))

		for i in 1:length(docucols[1])
			for j in 1:length(docucols)
				try
					!isempty(model.corp[docucols[j][i]].title) ? title = model.corp[docucols[j][i]].title : title = "Document $(docucols[j][i])"
					dspacing = maximum([max(4 + length("$(docucols[j][i])"), length(model.corp[d].title)) for d in docucols[j]]) - length(title) + 4
					rspacing = maximum([length("$r") for r in rankcols[j]]) - length(string(rankcols[j][i]))
					print(Crayon(foreground=:yellow, bold=true), string(rankcols[j][i]) * ". " * " "^rspacing)
					j == length(docucols) ? print(Crayon(foreground=:white, bold=false), title) : print(Crayon(foreground=:white, bold=false), title * " "^dspacing)

				catch
					nothing
				end
			end
			println()
		end

		if n < length(users)
			println()
		end
	end
end

showurecs(model::Union{CTPF, gpuCTPF}, user::Integer, M::Integer=15; cols::Integer=1) = showurecs(model, [user], M, cols=cols)
showurecs(model::Union{CTPF, gpuCTPF}, users::UnitRange{<:Integer}, M::Integer=15; cols::Integer=1) = showurecs(model, collect(users), M, cols=cols)

function predict(corp::Corpus, train_model::Union{LDA, gpuLDA}; iter::Integer=10, tol::Real=1/train_model.K^2)
	"Predict topic distributions for corpus of documents based on trained LDA model."

	check_corp(corp)
	check_model(train_model)
	(corp.vocab == train_model.corp.vocab)	|| throw(CorpusError("Predict Corpus and train_model Corpus must have identical vocabularies."))
	(tol .>= 0)								|| throw(ArgumentError("Tolerance parameter must be nonnegative."))
	(iter .>= 0)							|| throw(ArgumentError("Iteration parameter must be nonnegative."))

	model = LDA(corp, train_model.K)
	model.alpha = train_model.alpha
	model.beta = train_model.beta
	model.topics = train_model.topics

	for d in 1:model.M
		for v in 1:iter
			update_phi!(model, d)
			update_gamma!(model, d)
			update_Elogtheta!(model, d)
			if norm(model.Elogtheta[d] - model.Elogtheta_old[d]) < tol
				break
			end
		end
	end

	return model
end

function predict(corp::Corpus, train_model::fLDA; iter::Integer=10, tol::Real=1/train_model.K^2)
	"Predict topic distributions for corpus of documents based on trained fLDA model."

	check_corp(corp)
	check_model(train_model)
	(corp.vocab == train_model.corp.vocab)	|| throw(CorpusError("Predict Corpus and train_model Corpus must have identical vocabularies."))
	(tol .>= 0)								|| throw(ArgumentError("Tolerance parameter must be nonnegative."))
	(iter .>= 0)							|| throw(ArgumentError("Iteration parameter must be nonnegative."))

	model = fLDA(corp, train_model.K)
	model.alpha = train_model.alpha
	model.beta = train_model.beta
	model.topics = train_model.topics

	for d in 1:model.M
		for v in 1:iter
			update_phi!(model, d)
			update_tau!(model, d)
			update_gamma!(model, d)
			update_Elogtheta!(model, d)
			if norm(model.Elogtheta[d] - model.Elogtheta_old[d]) < vtol
				break
			end
		end
	end

	return model
end

function predict(corp::Corpus, train_model::Union{CTM, gpuCTM}; iter::Integer=10, tol::Real=1/train_model.K^2, niter::Integer=1000, ntol::Real=1/train_model.K^2)
	"Predict topic distributions for corpus of documents based on trained CTM model."

	check_corp(corp)
	check_model(train_model)
	(corp.vocab == train_model.corp.vocab)	|| throw(CorpusError("Predict Corpus and train_model Corpus must have identical vocabularies."))
	all([tol, ntol] .>= 0)					|| throw(ArgumentError("Tolerance parameters must be nonnegative."))
	all([iter, niter] .>= 0)				|| throw(ArgumentError("Iteration parameters must be nonnegative."))

	model = CTM(corp, train_model.K)
	model.mu = train_model.mu
	model.sigma = train_model.sigma
	model.invsigma = train_model.invsigma
	model.beta = train_model.beta
	model.topics = train_model.topics

	for d in 1:model.M
		for v in 1:iter
			update_phi!(model, d)
			update_logzeta!(model, d)
			update_vsq!(model, d, niter, ntol)
			update_lambda!(model, d, niter, ntol)
			if norm(model.lambda[d] - model.lambda_old[d]) < tol
				break
			end
		end
	end

	return model
end

function predict(corp::Corpus, train_model::fCTM; iter::Integer=10, tol::Real=1/train_model.K^2, niter::Integer=1000, ntol::Real=1/train_model.K^2)
	"Predict topic distributions for corpus of documents based on trained fCTM model."

	check_corp(corp)
	check_model(train_model)
	(corp.vocab == train_model.corp.vocab)	|| throw(CorpusError("Predict Corpus and train_model Corpus must have identical vocabularies."))
	all([tol, ntol] .>= 0)					|| throw(ArgumentError("Tolerance parameters must be nonnegative."))
	all([iter, niter] .>= 0)				|| throw(ArgumentError("Iteration parameters must be nonnegative."))

	model = fCTM(corp, train_model.K)
	model.mu = train_model.mu
	model.sigma = train_model.sigma
	model.invsigma = train_model.invsigma
	model.beta = train_model.beta
	model.topics = train_model.topics

	for d in 1:model.M
		for v in 1:iter
			update_phi!(model, d)
			update_tau!(model, d)
			update_logzeta!(model, d)
			update_lambda!(model, d, niter, ntol)
			update_vsq!(model, d, niter, ntol)
			if norm(model.lambda[d] - model.lambda_old[d]) < vtol
				break
			end
		end
	end

	return model
end

function topicdist(model::Union{LDA, gpuLDA, fLDA}, d::Integer)
	"Get the LDA topic distribution for a document as a probability vector."

	(d <= length(model.corp)) || throw(CorpusError("Document index outside corpus range."))

	topic_distribution = model.gamma[d] / sum(model.gamma[d])
	return topic_distribution
end

function topicdist(model::Union{CTM, gpuCTM, fCTM}, d::Integer)
	"Get the CTM topic distribution for a document as a probability vector."

	(d <= length(model.corp)) || throw(CorpusError("Document index outside corpus range."))

	topic_distribution = additive_logistic(model.lambda[d] + 0.5 * model.vsq[d])
	return topic_distribution
end

function topicdist(model::Union{CTPF, gpuCTPF}, d::Integer)
	"Get the CTPF topic distribution for a document as a probability vector."

	(d <= length(model.corp)) || throw(CorpusError("Document index outside corpus range."))

	topic_distribution = model.gimel[d] / sum(model.gimel[d])
	return topic_distribution
end

function topicdist(model::TopicModel, doc_indices::Vector{<:Integer})
	"Get TopicModel topic distributions for document(s) as a probability vector."

	issubset(doc_indices, 1:length(model.corp)) || throw(CorpusError("Some document indices outside corpus range."))

	topic_distributions = Vector{typeof(model.elbo)}[]
	for d in doc_indices
		push!(topic_distributions, topicdist(model, d))
	end

	return topic_distributions
end

topicdist(model::TopicModel, doc_range::UnitRange{<:Integer}) = topicdist(model, collect(doc_range))


text:
# Hermitian lattices
@attributes mutable struct HermLat{S, T, U, V, W} <: AbsLat{S}
  space::HermSpace{S, T, U, W}
  pmat::V
  gram::U
  rational_span::HermSpace{S, T, U, W}
  base_algebra::S
  involution::W
  automorphism_group_generators::Vector{U}
  automorphism_group_order::fmpz
  generators
  minimal_generators
  norm
  scale

  function HermLat{S, T, U, V, W}() where {S, T, U, V, W}
    z = new{S, T, U, V, W}()
    return z
  end

  function HermLat(E::S, G::U, P::V) where {S, U, V}
    @assert degree(E) == 2
    A = automorphisms(E)
    a = gen(E)
    if A[1](a) == a
      involution = A[2]
    else
      involution = A[1]
    end

    K = base_field(E)

    space = HermSpace(E, G)

    z = new{S, typeof(K), U, V, typeof(involution)}(space, P)
    z.base_algebra = E
    z.involution = involution
    return z
  end

  function HermLat(E::S, G::U) where {S, U}
    n = nrows(G)
    M = pseudo_matrix(identity_matrix(E, n))
    return HermLat(E, G, M)
  end
end


text:
import StatsBase

# function simulate(rng::AbstractRNG, cm::ConditionalModel{A,B,M,Argvals,EmptyNTtype}, N::Int) where {A,B,M,Argvals}
#     m = Model(cm)
#     cm0 = setReturn(m, nothing)(argvals(cm))
#     info = StructArray(simulate(rng, cm0, N))
#     vals = [predict(cm, pars) for pars in info]
#     return StructArray{Noted}((vals, info))
# end

# function simulate(cm::ConditionalModel{A,B,M,Argvals,EmptyNTtype}, N::Int) where {A,B,M,Argvals}
#     return simulate(GLOBAL_RNG, cm, N)
# end

# function simulate(rng::AbstractRNG, cm::ConditionalModel{A,B,M,Argvals,EmptyNTtype}) where {A,B,M,Argvals}
#     m = Model(cm)
#     cm0 = setReturn(m, nothing)(argvals(cm))
#     info = simulate(rng, cm0)
#     val = predict(cm, info)
#     return Noted(val, info)
# end

# function simulate(cm::ConditionalModel{A,B,M,Argvals,EmptyNTtype}) where {A,B,M,Argvals}
#     return simulate(GLOBAL_RNG, cm)
# end

using GeneralizedGenerated
using Random: GLOBAL_RNG
using NestedTuples
using SampleChains
using TupleVectors

EmptyNTtype = NamedTuple{(),Tuple{}} where T<:Tuple
export simulate

function simulate(rng::AbstractRNG, d::ConditionalModel, N::Int; trace_assignments=false)
    x = simulate(rng, d)
    T = typeof(x)
    ta = TupleVector(undef, x, N)
    @inbounds ta[1] = x

    for j in 2:N
        @inbounds ta[j] = simulate(rng, d; trace_assignments)
    end

    return ta
end

simulate(d::ConditionalModel, N::Int; trace_assignments=false) = simulate(GLOBAL_RNG, d, N; trace_assignments)

@inline function simulate(rng::AbstractRNG, c::ConditionalModel; trace_assignments=false)
    m = Model(c)
    return _simulate(getmoduletypencoding(m), m, argvals(c), Val(trace_assignments))(rng)
end

@inline function simulate(m::ConditionalModel; trace_assignments=false)
    simulate(GLOBAL_RNG, m; trace_assignments)
end

@inline function simulate(rng::AbstractRNG, m::Model; trace_assignments=false)
    return _simulate(getmoduletypencoding(m), m, NamedTuple())(rng)
end

simulate(m::Model; trace_assignments=false) = simulate(GLOBAL_RNG, m; trace_assignments)

sourceSimulate(m::Model; trace_assignments=false) = sourceSimulate(trace_assignments)(m)
sourceSimulate(jd::ConditionalModel; trace_assignments=false) = sourceSimulate(jd.model; trace_assignments)

export sourceSimulate
function sourceSimulate(trace_assignments=false)
    ta = trace_assignments
    function(_m::Model)
        pars = sort(sampled(_m))

        tracekeys = sort(trace_assignments ? parameters(_m) : sampled(_m))
        _traceName = Dict((k => Symbol(:_trace_, k) for k in tracekeys))

        function proc(_m, st::Assign)
            q = @q begin $(st.x) = $(st.rhs) end
            if trace_assignments
                _traceName[st.x] = Symbol(:_trace_, st.x)
                push!(q.args, :($(_traceName[st.x]) = $(st.x)))
            end
            q
        end

        proc(_m, st::Sample)  = @q begin
            $(_traceName[st.x]) = simulate(_rng, $(st.rhs); trace_assignments=$ta)
            $(st.x) = value($(_traceName[st.x]))
        end

        proc(_m, st::Return)  = :(_value = $(st.rhs))
        proc(_m, st::LineNumber) = nothing

        _traces = map(x -> Expr(:(=), x,_traceName[x]), tracekeys)

        wrap(kernel) = @q begin
            _rng -> begin
                _value = nothing
                $kernel
                _trace = $(Expr(:tuple, _traces...))
                return (trace = _trace, value = _value)
            end
        end

        buildSource(_m, proc, wrap) |> MacroTools.flatten
    end
end

# simulate(rng::AbstractRNG, d::Distribution; trace_assignments=false) = rand(rng, d)

# simulate(rng::AbstractRNG, d::iid{Int}; trace_assignments=false) = [simulate(rng, d.dist; trace_assignments) for j in 1:d.size]

trace(x::NamedTuple) = x.trace
trace(x) = x

using MeasureTheory: AbstractMeasure

simulate(μ::AbstractMeasure; trace_assignments=false) = simulate(Random.GLOBAL_RNG, μ; trace_assignments)

simulate(rng::AbstractRNG, μ::AbstractMeasure; trace_assignments=false) = rand(rng, μ)

@gg function _simulate(M::Type{<:TypeLevel}, _m::Model, _args, trace_assignments::Val{V}) where {V}
    trace_assignments = V
    body = type2model(_m) |> sourceSimulate(trace_assignments) |> loadvals(_args, NamedTuple())
    @under_global from_type(_unwrap_type(M)) @q let M
        $body
    end
end

@gg function _simulate(M::Type{<:TypeLevel}, _m::Model, _args::NamedTuple{()}, trace_assignments::Val{V}) where {V}
    trace_assignments = V
    body = type2model(_m) |> sourceSimulate(trace_assignments)
    @under_global from_type(_unwrap_type(M)) @q let M
        $body
    end
end