text:
"""
A collection of various stable way to do a linear regression

#### References

Sections 3 and 4 of "Numerically stable and accurate stochastic simulation
approaches for solving dynamic economic models."

Kenneth L Judd, Serguei Maliar, and Lilia Maliar.

"""

using Parameters
using MathProgBase

abstract type AbstractLinReg end
abstract type AbstractLSMethod <: AbstractLinReg end
abstract type AbstractLADMethod <: AbstractLinReg end

# --------- #
# Utilities #
# --------- #

function standardize(x::AbstractMatrix, intercept::Bool=false)
    _x = intercept ? view(x, :, 2:size(x, 2)) : view(x, :, :)
    μ = mean(_x, dims=1)
    σ = std(_x, dims=1)
    xn = (_x .- μ) ./ σ
    xn, vec(μ), vec(σ)
end

function standardize(y::AbstractVector, intercet=false)
    μ = mean(y)
    σ = std(y)
    yn = (y .- μ) ./ σ
    yn, μ, σ
end

function add_intercept(x::AbstractArray{T}) where T
    [ones(T, size(x, 1)) x]
end

# --- #
# OLS #
# --- #

@with_kw struct OLS <: AbstractLSMethod
    standardize::Bool = true
    intercept::Bool = true
    @assert !(standardize && !(intercept)) "must have intercept if standardizing"
end

# API methods
should_standardize(m::OLS) = m.standardize
should_add_intercept(m::OLS) = m.intercept
slopes(::OLS, x, y) = pinv(x'x)*x'y

# ------ #
# LSLdiv #
# ------ #

@with_kw struct LSLdiv <: AbstractLSMethod
    standardize::Bool = true
    intercept::Bool = true
    @assert !(standardize && !(intercept)) "must have intercept if standardizing"
end

# API methods
should_standardize(m::LSLdiv) = m.standardize
should_add_intercept(m::LSLdiv) = m.intercept
slopes(::LSLdiv, x, y) = x\y

# ----- #
# LSSVD #
# ----- #

@with_kw struct LSSVD <: AbstractLSMethod
    standardize::Bool = true
    intercept::Bool = true
    @assert !(standardize && !(intercept)) "must have intercept if standardizing"
end

# API methods
should_standardize(m::LSSVD) = m.standardize
should_add_intercept(m::LSSVD) = m.intercept
function slopes(m::LSSVD, x, y)
    U, S, V = svd(x, full=false)
    S_inv = diagm(0 => 1.0 ./ S)
    V*S_inv*U'y
end

# ----------- #
# RLSTikhonov #
# ----------- #

@with_kw struct RLSTikhonov <: AbstractLSMethod
    η::Float64 = -5.0
    standardize::Bool = true
    intercept::Bool = true
    @assert !(standardize && !(intercept)) "must have intercept if standardizing"
    @assert η < 0.0 "penalty parameter for RLSTikhonov must be negative"
end

const RLST = RLSTikhonov

# API methods
should_standardize(m::RLSTikhonov) = m.standardize
should_add_intercept(m::RLSTikhonov) = m.intercept
function slopes(m::RLSTikhonov, x, y)
    nobs = size(x, 1); nx = size(x, 2)
    pinv(x'x + (nobs/nx*10^(m.η))*I)*x'y
end

# ----- #
# LADPP #
# ----- #

@with_kw struct LADPP <: AbstractLADMethod
    standardize::Bool = true
    intercept::Bool = true
    @assert (standardize && intercept) "standardize and intercept must both be true"
end

"""
Solves the linear program

``min_{v⁺; v⁻; β} 1v⁺ + 1v⁻``

such that

``\\begin{cases}v⁺ - v⁻ + X \\beta &= y \\\\ v⁺ \\ge 0, \\quad v⁻ \\ge 0``

"""
function slopes(m::LADPP, x, y::AbstractMatrix)
    nobs = size(x, 1)
    nx = size(x, 2)
    ny = size(y, 2)

    # The linear programming coefficient vector, x, is ordered [v⁺; v⁻; β]
    # we will now construct the matrices for the upper and lower bounds. Note
    # that we artificailly say that β must be between -300 and 300
    lb = [fill(0.0, 2*nobs); fill(-300, nx)]
    ub = [fill(Inf, 2*nobs); fill(300, nx)]

    # objective
    c = [ones(2*nobs); zeros(nx)]

    # equality constraints
    A = [I -I x]

    # we need to do linprog equation by equation (1 column of y at a time),
    # so we will pre-allocate a coefficient matrix and fill it in
    β = Array(Float64, nx, ny)

    for i in 1:ny
        # equality constraint
        b = y[:, i]
        out = linprog(c, A, '=', b, lb, ub).sol
        β[:, i] = out[end-nx+1:end]

    end
    β
end

# ----- #
# LADDP #
# ----- #

@with_kw struct LADDP <: AbstractLADMethod
    standardize::Bool = true
    intercept::Bool = true
    @assert (standardize && intercept) "standardize and intercept must both be true"
end

"""
Solves the linear program

``min_{q} -y'q``

such that

``X'q = 0 \\quad -1 \\le q \\le 1``

"""
function slopes(m::LADDP, x, y::AbstractMatrix)
    nobs = size(x, 1)
    nx = size(x, 2)
    ny = size(y, 2)

    # The linear programming coefficient vector, x, is ordered [q]
    # we will now construct the matrices for the upper and lower bounds. Note
    # that we artificailly say that β must be between -300 and 300
    lb = fill(-1.0, nobs)
    ub = fill(1.0, nobs)

    # equality constraint
    b = zeros(nx)

    # equality constraints
    A = x'

    # we need to do linprog equation by equation (1 column of y at a time),
    # so we will pre-allocate a coefficient matrix and fill it in
    β = Array(Float64, nx, ny)

    for i in 1:ny
        # equality constraint
        c = -y[:, i]
        out = linprog(c, A, '=', b, lb, ub)
        β[:, i] = -out.attrs[:lambda][1:nx]

    end
    β
end

# ------ #
# RLADPP #
# ------ #

@with_kw struct RLADPP <: AbstractLADMethod
    η::Float64 = -5.0
    standardize::Bool = true
    intercept::Bool = true
    @assert η < 0.0 "penalty parameter for RLSTikhonov must be negative"
    @assert (standardize && intercept) "standardize and intercept must both be true"
end

"""
Solves the linear program

``min_{c:=[v⁺; v⁻; ψ⁺; ψ⁻]} 1v⁺ + 1v⁻ + 1ψ⁺+ 1ψ⁻``

such that

``\\begin{cases}v⁺ - v⁻ + Xψ⁺ - Xψ⁻ &= y \\\\ c &\\ge 0``

"""
function slopes(m::RLADPP, x, y::AbstractMatrix)
    nobs = size(x, 1)
    nx = size(x, 2)
    ny = size(y, 2)

    # The linear programming coefficient vector, x, is ordered [v⁺; v⁻; ψ⁺; ψ⁻]
    lb = fill(0.0, 2*nobs + 2*nx)
    ub = fill(Inf, 2*nobs + 2*nx)

    # objective
    c = [ones(2*nobs); fill(10^(m.η)*nobs/nx, 2*nx)]

    # equality constraints
    A = [I -I x -x]

    # we need to do linprog equation by equation (1 column of y at a time),
    # so we will pre-allocate a coefficient matrix and fill it in
    β = Array(Float64, nx, ny)

    for i in 1:ny
        # equality constraint
        b = y[:, i]
        out = linprog(c, A, '=', b, lb, ub).sol
        ψ⁺ = out[end-(2*nx)+1:end-nx]
        ψ⁻ = out[end-nx+1:end]
        β[:, i] = ψ⁺ - ψ⁻
    end
    β
end

# ------ #
# RLADDP #
# ------ #

@with_kw struct RLADDP <: AbstractLADMethod
    η::Float64 = -5.0
    standardize::Bool = true
    intercept::Bool = true
    @assert η < 0.0 "penalty parameter for RLSTikhonov must be negative"
    @assert (standardize && intercept) "standardize and intercept must both be true"
end

"""
Solves the linear program

``min_{q} -y'q``

such that

``X'q \\le η1 \\quad -X'q \\le η1, \\quad -1 \\le q \\le 1``

"""
function slopes(m::RLADDP, x, y::AbstractMatrix)
    nobs = size(x, 1)
    nx = size(x, 2)
    ny = size(y, 2)

    # The linear programming coefficient vector, x, is ordered [q]
    # we will now construct the matrices for the upper and lower bounds. Note
    # that we artificailly say that β must be between -300 and 300
    lb = fill(-1.0, nobs)
    ub = fill(1.0, nobs)

    # constraints equations
    A = [x'; -x']
    b = fill(10^(m.η)*nobs/nx, 2*nx)

    # we need to do linprog equation by equation (1 column of y at a time),
    # so we will pre-allocate a coefficient matrix and fill it in
    β = Array(Float64, nx, ny)

    for i in 1:ny
        # equality constraint
        c = -y[:, i]
        out = linprog(c, A, '<', b, lb, ub)
        ψ⁺ = -out.attrs[:lambda][1:nx]
        ψ⁻ = -out.attrs[:lambda][nx+1:2*nx]
        β[:, i] = ψ⁺ - ψ⁻
    end
    β

end

# ------ #
# RLSSVD #
# ------ #

@with_kw struct RLSSVD <: AbstractLSMethod
    κ::Float64=100_000.0
    standardize::Bool = true
    intercept::Bool = true
    @assert κ > 0 "κ must be positive"
    @assert !(standardize && !(intercept)) "must have intercept if standardizing"
end

# API methods
should_standardize(m::RLSSVD) = m.standardize
should_add_intercept(m::RLSSVD) = m.intercept
function slopes(m::RLSSVD, x, y)
    U, S, V = svd(x, full=false)

    # find number of sufficiently independent principal components
    condition_numbers = S[1] ./ S[2:end]
    r = findfirst(_1 -> _1 > m.κ, condition_numbers)
    r = r == 0 ? length(S) : r

    # extract these columns of V, U, S
    Vr = V[:, 1:r]
    Ur = U[:, 1:r]
    Sr_inv = diagm(0 => 1.0 ./ S[1:r])

    # compute estimate
    Vr*Sr_inv*Ur'y
end

# -------- #
# Generics #
# -------- #

function slopes(m::AbstractLADMethod, x, y::AbstractVector)
    βmat = slopes(m, x, y[:, :])
    βmat[:, 1]
end

should_standardize(::AbstractLinReg) = true
should_add_intercept(::AbstractLinReg) = true

_reg_nonorm_noint(m::AbstractLinReg, x, y) = slopes(m, x, y)
_reg_nonorm_int(m::AbstractLinReg, x, y) = slopes(m, add_intercept(x), y)

function _reg_norm_noint(m::AbstractLinReg, x, y)
    xn, μx, σx = standardize(x)
    yn, μy, σy = standardize(y)

    # Compute standardized coefficients, then undo normalization
    β = slopes(m, xn, yn)
    # update slopes
    for j in 1:size(y, 2)
        for i in 1:size(β, 1)
            β[i, j] *= σy[j]/σx[i]
        end
    end
    β
end

function _reg_norm_int(m::AbstractLinReg, x, y::AbstractVector)
    xn, μx, σx = standardize(x)
    yn, μy, σy = standardize(y)

    # Compute standardized coefficients, then undo normalization
    β1 = slopes(m, xn, yn)
    for i in eachindex(β1)
        β1[i] *= σy/σx[i]
    end
    β0 = μy - dot(μx, β1)

    β = similar(β1, length(β1) + 1)
    β[1] = β0
    β[2:end] = β1
    β
end

function _reg_norm_int(m::AbstractLinReg, x, y::AbstractMatrix)
    xn, μx, σx = standardize(x)
    yn, μy, σy = standardize(y)

    # Compute standardized coefficients, then undo normalization
    β1 = slopes(m, xn, yn)
    β0 = similar(β1, 1, size(y, 2))

    # update slopes
    for j in 1:size(y, 2)
        for i in 1:size(β1, 1)
            β1[i, j] *= σy[j]/σx[i]
        end
    end

    # compute intercept
    for j in 1:size(y, 2)
        β0[j] = μy[j] - dot(β1[:, j], μx)
    end

    Float64[β0; β1]
end

regress(m::AbstractLinReg, x, y) =
    (should_standardize(m)  && should_add_intercept(m))  ? _reg_norm_int(m, x, y) :
    (!should_standardize(m) && should_add_intercept(m))  ? _reg_nonorm_int(m, x, y) :
    (should_standardize(m)  && !should_add_intercept(m)) ? _reg_norm_noint(m, x, y) :
                                                         _reg_nonorm_noint(m, x, y)


text:
# ship image processing

using RoME, Distributions
using Images, ImageView, ImageDraw
using FileIO
using FixedPointNumbers
using Gadfly
using RoMEPlotting, KernelDensityEstimatePlotting
using StatsBase
using Colors
using ColorTypes

datadir = "/home/dehann/software/multibeam/data/dataset-1/"
datafile = datadir*"1468970745714629_raw.png"
datafile = datadir*"1468970485675274_raw.png"
datafile = datadir*"1468969086135548_raw.png"

img = load(datafile)
imshow(img)

rangebins, azibins = size(img)

# calibrate to SI
ranges = 2.25+9.0/512.0*(0:rangebins-1)
azis = (14.4 - 0.3*(0.5:1:azibins-0.5)) * pi/180.0

#   Factor graphing

N = 75
fg = initfg()

addVariable!(fg, :x0, Pose2, N=N)
addFactor!(fg, [:x0], PriorPose2(MvNormal(zeros(3), 0.001^2*eye(3))))

doautoinit!(fg, :x0)

prevsym = :null
prevsymm = :null

for idx in 1:96

lsym = Symbol("l$(idx)")
addVariable!(fg, lsym, Point2, N=N)

llsym = Symbol("l$(idx+1000)")
# addVariable!(fg, llsym, Point2, N=N)

bss = IIF.AliasingScalarSampler(collect(ranges),  Float64.(img[:,idx]), SNRfloor=0.70)
tr = Pose2Point2BearingRange(Normal(azis[idx], 0.005), bss)
addFactor!(fg, [:x0;lsym], tr)
# addFactor!(fg, [:x0;lsym;llsym], tr, multihypo=[1.0;0.5;0.5])

if prevsym != :null
  pp2 = Point2Point2(MvNormal(zeros(2), 0.1^2*eye(2)))
  addFactor!(fg, [prevsym,lsym], pp2)

  # pp22 = Point2Point2(MvNormal(zeros(2), 0.5^2*eye(2)))
  # addFactor!(fg, [prevsymm,llsym], pp22)
end
prevsym = lsym
prevsymm = llsym

end

# the measurement point
# addVariable!(fg, :l48, Point2)
# bss = IIF.AliasingScalarSampler(collect(ranges),  Float64.(img[:,48]), SNRfloor=0.7)
# tr = Pose2Point2BearingRange(Normal(0.0, 0.005), bss)
# addFactor!(fg, [:x0;:l48], tr)

# @async writeGraphPdf(fg)
#
# pts = IIF.approxConv(fg, :x0l48f1, :l48)

# p1 = plotKDE(kde!(pts), axis=[-1.0 10.0;
#                               -1.0 1.0], levels=8) # , c=["red"]
#

getSolverParams(fg).N = N
@time solveTree!(fg)

drawPosesLandms(fg, spscale=0.1, drawhist=false, meanmax=:max)

drawPosesLandms(fg, spscale=0.1, drawhist=false, meanmax=:mean, xmin=5.3,xmax=5.7, ymin=-2.0,ymax=2.0)

plotKDE(fg, [Symbol("l$i") for i in 0:2:60], levels=1)

# savejld(fg)
# loadjld(fg)

imgc = convert(Array{RGB{N0f8}},img)
# imshow(imgc)
# imgdd = RGB.(imgd,imgd,imgd)

# xy = KDE.getKDEMax(getBelief(fg, :l1))

# imgc = ImageDraw.draw(imgc, CirclePointRadius(48, 256, 2), RGB{N0f8}(1.0,0.0,0.0))

for ll in ls(fg)[2]
  xy = KDE.getKDEMax(getBelief(fg, ll))
  the = atan2(-xy[2],xy[1]) * 180.0/pi
  ran = norm(xy)

  u = round(Int, (ran - 2.25)/(9.0/512.0))
  v = round(Int, (the + 14.4)/(0.3))
  @show v
  v = v .< 2 ? 2 : v
  v = 95 .< v ? 95 : v
  @show v, u

  imgc = ImageDraw.draw(imgc, CirclePointRadius(v, u, 2), RGB{N0f8}(1.0,0.0,0.0))

end

imgc

## DEV WORK

plot(x=ranges, y=img[:,48], Geom.line)

imgf = Float64.(img[:,48])
imgf .-= quantile(imgf, 0.5)
imgf[imgf .< 0.0] = 0.0
imgf ./= norm(imgf)

l1 = layer(x=ranges, y=100*imgf, Geom.line)
plot(l1)

wim = StatsBase.ProbabilityWeights(imgf)

smpls = zeros(100)
StatsBase.alias_sample!(ranges,wim,smpls)

l2 = layer(x=smpls, Geom.histogram(bincount=70), Theme(default_color=colorant"orange"))

plot(l1,l2)

pts = rand(bss, 100)

l3 = layer(x=pts, Geom.histogram(bincount=70), Theme(default_color=colorant"green"))

plot(l1,l2,l3)

#

ext: jl


text:
using SelfPropelledModel
using NaturalES
using LinearAlgebra

transpv(v)=collect(v')

function mysim(P)
    N=160
    sv=sample_particles(N,()->500.0.*(rand(2).-0.5),()->(rand(2).-0.5))
    frames=[]
    for i in 1:170
        X=foldl(hcat,[[p.x...] for p in sv])
        push!(frames,X)
        spm_step!(sv,P[1],P[2],P[3],P[4])
    end
    transpv.(frames)
end

function metrics(f,nf)
    Δ=nf.-f
    Np,Nc=size(Δ)
    m0=0.0
    m1=0.0
    m2=0.0
    m3=0.0
    m4=0.0
    for i in 3:Np-2, j in i-2:i+2
        s=normalize(Δ[i,:])⋅normalize(Δ[j,:])
        m0+=1
        m1+=s
        m2+=s*s
        m3+=s*s*s
        m4+=s*s*s*s
    end
    M1=m1/m0
    M2=m2/m0
    M3=m3/m0
    M4=m4/m0

    (M1,
    M2-M1*M1,
    M3-3M1*M2+2M1*M1*M1
    )
end

function metrics(frames)
    m=metrics(frames[1],frames[2])
    N=1
    for i in 2:length(frames)-1
        m=m.+metrics(frames[i],frames[i+1])
        N+=1
    end
    m./N
end
mysol=zeros(4)
Nmysol=0
function cost(x)
    global mysol.+=x
    global Nmysol+=1
    target=(0.20677718821137364, 0.5559917303075799, -0.14723289814743612)
    s=sqrt(sum(abs2,metrics(mysim(x)).-target))
    println(s," -> ",mysol./Nmysol)
    s
end

A=mysim([0.4,0.5,1.0,1.5])
metrics(A)

separable_nes(cost,[0.5,0.5,0.5,0.5],0.1)

x=mysol./mysol[3]
x[4]=1.5

frames=mysim(x)
using PyPlot
pts=PermutedDimsArray(cat(frames...,dims=3),(1,3,2))|>collect
pygui(true)
figure(figsize=(5,5))
for i in 1:160
    plot(pts[i,:,1],pts[i,:,2],color="black",lw=0.5)
end
xlim((-250,250))
ylim((-250,250))
tight_layout()


text:
function getUserData(id::AbstractString, db=getDB())
    data = DBInterface.execute(db, """
        SELECT * FROM userdata WHERE id = '$id'
    """) |> DataFrame
    if length(data[!, 1]) == 0
        return (false, missing)
    else
        return (true, getUserData(data))
    end
end
function getUserData(data::DataFrame)
    pos = Position(
      data.x[1],
      data.y[1]
    )
    inv = JSON3.read(data.inv[1], Vector{InventoryItem})
    Player(
        pos,
        data.char[1],
        data.emp[1],
        inv,
        data.rank[1],
        data.id[1]
    )
end
function getUserData(id::String, cache::PlayerCache, force=false)
    if id in keys(cache) && !force
        cache[id]
    else
        ply = getUserData(id)
        cache[ply.id] = ply
        ply
    end
end

function saveUserData(ply::Player, db=getDB())
    DBInterface.execute(db, """
    UPDATE userdata
    SET x = $(ply.pos.x), y = $(ply.pos.y), char = $(ply.char), emp = $(ply.emp), inv = '$(JSON3.write(ply.inv))', rank = $(ply.rank)
    WHERE id = '$(ply.id)'
    """)
end
saveUserData(id::AbstractString, cache::PlayerCache, db=getDB()) = saveUserData(getUserData(id, cache), db)

function createUser(id::String, char::Int)
    pos = Position(
        rand(200:500),
        rand(300:400)
    )
    ply = Player(
        pos,
        char,
        0,
        [],
        0,
        id
    )
    createUser(ply)
    ply
end
function createUser(ply::Player, db=getDB())
    DBInterface.execute(db, """
    INSERT INTO userdata (x, y, char, emp, inv, rank, id)
    VALUES ($(ply.pos.x), $(ply.pos.y), $(ply.char), $(ply.emp), '$(JSON3.write(ply.inv))', $(ply.rank), '$(ply.id)')
    """)
end


text:
using RemoteGtkREPL, Pkg #need Pkg for REPLMode
RemoteGtkREPL.estalbish_connection(ARGS)

#= RemoteGtkREPL.remotecall_fetch(include_string, GtkREPLWorker.gtkrepl, Main,
    "$(GtkREPLWorker.remote_mod).add_remote_console_cb($(GtkREPLWorker.id), $(GtkREPLWorker.port))"
) =#

# put important things in a module for safety
#= module GtkREPLWorker

    using RemoteGtkREPL, Sockets

    gtkrepl_port = parse(Int,ARGS[1])
    global const id = parse(Int, ARGS[2]) # console/worker id
    global const remote_mod = ARGS[3]    # the module calling us as a String

    port, server = RemoteGtkREPL.start_server()

    global const gtkrepl = connect(gtkrepl_port)

    RemoteGtkREPL.init(gtkrepl, id)
end =#

#ploting stuff
#= function gadfly()
    @eval begin
        RemoteGtkREPL.gadfly()

        export figure
        import Base: show, display

        show(io::IO, p::Gadfly.Plot) = write(io,"Gadfly.Plot(...)")
        function display(p::Gadfly.Plot)
            remotecall_fetch(display, GtkREPLWorker.gtkrepl, p)
            nothing
        end

        figure() = remotecall_fetch(RemoteGtkREPL.eval_command_remotely, GtkREPLWorker.gtkrepl,"figure()", Main)
        figure(i::Integer) = remotecall_fetch(RemoteGtkREPL.eval_command_remotely, GtkREPLWorker.gtkrepl,"figure($i)", Main)
    end
end =#

#start Gadfly by default
#using Gadfly
#gadfly()


text:
# Syntactic structure of a SQL query.

# Base type.

"""
A part of a SQL query.
"""
abstract type AbstractSQLClause
end

function dissect(scr::Symbol, ClauseType::Type{<:AbstractSQLClause}, pats::Vector{Any})
    scr_core = gensym(:scr_core)
    ex = Expr(:&&, :($scr_core isa $ClauseType), Any[dissect(scr_core, pat) for pat in pats]...)
    :($scr isa SQLClause && (local $scr_core = $scr[]; $ex))
end

# Opaque wrapper that serves as a specialization barrier.

"""
An opaque wrapper over an arbitrary SQL clause.
"""
struct SQLClause <: AbstractSQLClause
    core::AbstractSQLClause

    SQLClause(@nospecialize core::AbstractSQLClause) =
        new(core)
end

Base.getindex(c::SQLClause) =
    c.core

Base.convert(::Type{SQLClause}, c::SQLClause) =
    c

Base.convert(::Type{SQLClause}, @nospecialize c::AbstractSQLClause) =
    SQLClause(c)

Base.convert(::Type{SQLClause}, obj) =
    convert(SQLClause, convert(AbstractSQLClause, obj)::AbstractSQLClause)

(c::AbstractSQLClause)(c′) =
    c(convert(SQLClause, c′))

(c::AbstractSQLClause)(c′::SQLClause) =
    rebase(c, c′)

rebase(c::SQLClause, c′) =
    convert(SQLClause, rebase(c[], c′))

rebase(::Nothing, c′) =
    c′

Base.:(==)(@nospecialize(::AbstractSQLClause), @nospecialize(::AbstractSQLClause)) =
    false

Base.:(==)(c1::SQLClause, c2::SQLClause) =
    c1[] == c2[]

Base.:(==)(c1::SQLClause, @nospecialize(c2::AbstractSQLClause)) =
    c1[] == c2

Base.:(==)(@nospecialize(c1::AbstractSQLClause), c2::SQLClause) =
    c1 == c2[]

@generated function Base.:(==)(c1::C, c2::C) where {C <: AbstractSQLClause}
    exs = Expr[]
    for f in fieldnames(c1)
        push!(exs, :(isequal(c1.$(f), c2.$(f))))
    end
    Expr(:||, :(c1 === c2), Expr(:&&, exs...))
end

Base.hash(c::SQLClause, h::UInt) =
    hash(c[], h)

@generated function Base.hash(c::AbstractSQLClause, h::UInt)
    ex = :(h + $(hash(c)))
    for f in fieldnames(c)
        ex = :(hash(c.$(f), $ex))
    end
    ex
end

# Pretty-printing.

Base.show(io::IO, c::AbstractSQLClause) =
    print(io, quoteof(c, limit = true))

Base.show(io::IO, ::MIME"text/plain", c::AbstractSQLClause) =
    pprint(io, c)

function PrettyPrinting.quoteof(c::SQLClause; limit::Bool = false, unwrap::Bool = false)
    ctx = QuoteContext(limit = limit)
    ex = quoteof(c[], ctx)
    if unwrap
        ex = Expr(:ref, ex)
    end
    ex
end

PrettyPrinting.quoteof(c::AbstractSQLClause; limit::Bool = false) =
    quoteof(convert(SQLClause, c), limit = limit, unwrap = true)

PrettyPrinting.quoteof(c::SQLClause, ctx::QuoteContext) =
    if !ctx.limit
        quoteof(c[], ctx)
    else
        :…
    end

PrettyPrinting.quoteof(cs::Vector{SQLClause}, ctx::QuoteContext) =
    if isempty(cs)
        Any[]
    elseif !ctx.limit
        Any[quoteof(c, ctx) for c in cs]
    else
        Any[:…]
    end

PrettyPrinting.quoteof(names::Vector{Symbol}, ctx::QuoteContext) =
    if isempty(names)
        Any[]
    elseif !ctx.limit
        Any[QuoteNode(name) for name in names]
    else
        Any[:…]
    end

# Concrete clause types.

include("clauses/aggregate.jl")
include("clauses/as.jl")
include("clauses/case.jl")
include("clauses/cte.jl")
include("clauses/from.jl")
include("clauses/function.jl")
include("clauses/group.jl")
include("clauses/having.jl")
include("clauses/identifier.jl")
include("clauses/join.jl")
include("clauses/keyword.jl")
include("clauses/limit.jl")
include("clauses/literal.jl")
include("clauses/operator.jl")
include("clauses/order.jl")
include("clauses/partition.jl")
include("clauses/select.jl")
include("clauses/sort.jl")
include("clauses/union.jl")
include("clauses/variable.jl")
include("clauses/where.jl")
include("clauses/window.jl")
include("clauses/with.jl")


text:
__precompile__(true)
module UnitfulHartree
using Unitful
import Unitful

@unit mₑ "mₑ" ElectronMass 9.1093835611e-31*u"kg" false
@unit e₀ "eₒ" ElementaryCharge 1.602176620898e-19*u"C" false
@unit kₑ "kₑ" CoulombForceConstant 1/(4π)u"ϵ0" false
@unit ħ "ħ" ReducedPlanckConstant Unitful.ħ false
@unit a₀ "a₀" BohrRadius 1ħ^2/(1kₑ*mₑ*e₀^2) false
@unit Eₕ "Eₕ" HartreeEnergy mₑ*e₀^4*kₑ^2/(1ħ^2) true
@unit Ry "Ry" RydbergEnergy 0.5Eₕ true
@unit rₑ "rₑ" ClassicalElectronRadius 1e₀^2*kₑ/(1mₑ*Unitful.c^2) false
const α = 1e₀^2*1kₑ/(1Unitful.c*ħ)
const mₚ = 1836.15mₑ
const μ_b = e₀*ħ/(2mₑ)
const ϵ₀ = 1/(4π*kₑ)

# Some gymnastics required here because if we precompile, we cannot add to
# Unitful.basefactors at compile time and expect the changes to persist to runtime.
const localunits = Unitful.basefactors
function __init__()
    merge!(Unitful.basefactors, localunits)
    Unitful.register(UnitfulHartree)
end

end # module

max_stars_repo_name: mdavezac/UnitfulHartree

max_issues_repo_name: mdavezac/UnitfulHartree

max_forks_repo_name: mdavezac/UnitfulHartree


text:
"""
`prettyprint(name::ASCIIString)`

        prints simple a hexagonal layout view of the rectangular grid layout of the beeswax
        program "name" for more comfortable viewing.

        Equivalent to 'prettyprint(name,1)'.

"""
function prettyprint(name::AbstractString)
    prettyprint(name,1)
end

"""
`prettyprint(name::ASCIIString,style::Int)`

        prints either a simple hexagonal layout of the beeswax program "name" or a full
        hexagonal layout with grid overlay to make following the code easier.

    Available styles

 *   1     simple hexagonal layout.
 *   2     advanced hexagonal layout with overlay.

"""
function prettyprint(name::AbstractString,style::Int)
    prog=open(name)
    code=readlines(prog)
    close(prog)
    rows=length(code)
    cols=0
    for i=1:rows
        code[i]=chomp(code[i])
        cols=maximum([cols,strwidth(code[i])])
    end

    if style == 1
        for row=1:rows
            print(" "^(rows-row))
            code[row]=replace(code[row]," ","•")
            code[row]=rpad(code[row],cols,"•")
            for i in code[row]
                print("$i ")
            end
            println()
        end
    elseif style == 2
        for row=1:rows
            print(" "^(2*(rows-row)))
            code[row]=replace(code[row]," ","•")
            code[row]=rpad(code[row],cols,"•")
            for i=1:endof(code[row])
                try
                    i < endof(code[row]) ? print(code[row][i]," — ") : print(code[row][i])
                catch
                end
            end
            println()
            row<rows ? print(" "^(2*(rows-row)-2)) : nothing
            row<rows ? print(" / \\"^(cols-1)*" /\n") : nothing
        end
    end
end


text:
"""
Resolve citations in bibliography, and format bibnodes and citations.
"""
abstract type ResolveCitations <: Builder.DocumentPipeline end
Selectors.order(::Type{ResolveCitations}) = 5.5
Selectors.matcher(::Type{ResolveCitations}) = true
Selectors.strict(::Type{ResolveCitations}) = false
function Selectors.runner(::Type{ResolveCitations}, doc::Documents.Document)
    Utilities.log(doc, "resolving citations.")
    resolve_citations(doc)
end

function resolve_citations(doc::Documents.Document)
	for (src, page) in doc.internal.pages
		empty!(page.globals.meta)
		for element in page.elements
			resolve_citations(page.mapping[element], page, doc)
		end
	end
end

function resolve_citations(element, page, doc)
	walk_and_replace(page.globals.meta, element) do el
		resolve_citation(el, page.globals.meta, page, doc)
	end
end

resolve_citation(other, meta, page, doc) = nothing
function resolve_citation(cite::CitationNode, meta, page, doc)
	return cite.bibnode.database[cite.key][2]
end


text:
using Test
using Manifolds
using Manifolds:
    _derivative,
    _derivative!,
    differential,
    differential!,
    gradient,
    gradient!,
    _gradient,
    _gradient!,
    _jacobian,
    _jacobian!

import Manifolds: gradient

struct TestRiemannianBackend <: AbstractRiemannianDiffBackend end
function Manifolds.gradient(::AbstractManifold, f, p, ::TestRiemannianBackend)
    return collect(1.0:length(p))
end

using FiniteDifferences, FiniteDiff
using LinearAlgebra: Diagonal, dot

@testset "Differentiation backend" begin
    fd51 = Manifolds.FiniteDifferencesBackend()
    @testset "default_differential_backend" begin
        #ForwardDiff is loaded first in utils.
        @test default_differential_backend() === Manifolds.ForwardDiffBackend()

        @test length(fd51.method.grid) == 5
        # check method order
        @test typeof(fd51.method).parameters[2] == 1
        fd71 = Manifolds.FiniteDifferencesBackend(central_fdm(7, 1))
        @test set_default_differential_backend!(fd71) == fd71
        @test default_differential_backend() == fd71
    end

    using ForwardDiff

    fwd_diff = Manifolds.ForwardDiffBackend()
    @testset "ForwardDiff" begin
        @test default_differential_backend() isa Manifolds.FiniteDifferencesBackend

        @test set_default_differential_backend!(fwd_diff) == fwd_diff
        @test default_differential_backend() == fwd_diff
        @test set_default_differential_backend!(fd51) isa Manifolds.FiniteDifferencesBackend
        @test default_differential_backend() isa Manifolds.FiniteDifferencesBackend

        set_default_differential_backend!(fwd_diff)
        @test default_differential_backend() == fwd_diff
        set_default_differential_backend!(fd51)
    end

    using FiniteDiff

    finite_diff = Manifolds.FiniteDiffBackend()
    @testset "FiniteDiff" begin
        @test default_differential_backend() isa Manifolds.FiniteDifferencesBackend

        @test set_default_differential_backend!(finite_diff) == finite_diff
        @test default_differential_backend() == finite_diff
        @test set_default_differential_backend!(fd51) isa Manifolds.FiniteDifferencesBackend
        @test default_differential_backend() isa Manifolds.FiniteDifferencesBackend

        set_default_differential_backend!(finite_diff)
        @test default_differential_backend() == finite_diff
        set_default_differential_backend!(fd51)
    end

    using ReverseDiff

    reverse_diff = Manifolds.ReverseDiffBackend()
    @testset "ReverseDiff" begin
        @test default_differential_backend() isa Manifolds.FiniteDifferencesBackend

        @test set_default_differential_backend!(reverse_diff) == reverse_diff
        @test default_differential_backend() == reverse_diff
        @test set_default_differential_backend!(fd51) isa Manifolds.FiniteDifferencesBackend
        @test default_differential_backend() isa Manifolds.FiniteDifferencesBackend

        set_default_differential_backend!(reverse_diff)
        @test default_differential_backend() == reverse_diff
        set_default_differential_backend!(fd51)
    end

    using Zygote
    zygote_diff = Manifolds.ZygoteDiffBackend()

    @testset "gradient" begin
        set_default_differential_backend!(fd51)
        r2 = Euclidean(2)

        c1(t) = [sin(t), cos(t)]
        f1(x) = x[1] + x[2]^2
        function f1!(y, x)
            y .= x[1] + x[2]^2
            return y
        end
        f2(x) = 3 * x[1] * x[2] + x[2]^3

        @testset "Inference" begin
            X = [-1.0, -1.0]
            @test (@inferred _derivative(c1, 0.0, Manifolds.ForwardDiffBackend())) ≈
                  [1.0, 0.0]
            @test (@inferred _derivative!(c1, X, 0.0, Manifolds.ForwardDiffBackend())) === X
            @test X ≈ [1.0, 0.0]

            @test (@inferred _derivative(c1, 0.0, finite_diff)) ≈ [1.0, 0.0]
            @test (@inferred _gradient(f1, [1.0, -1.0], finite_diff)) ≈ [1.0, -2.0]
        end

        @testset for backend in [fd51, fwd_diff, finite_diff]
            set_default_differential_backend!(backend)
            @test _derivative(c1, 0.0) ≈ [1.0, 0.0]
            X = [-1.0, -1.0]
            @test _derivative!(c1, X, 0.0) === X
            @test isapprox(X, [1.0, 0.0])
        end
        @testset for backend in [fd51, fwd_diff, finite_diff, reverse_diff, zygote_diff]
            set_default_differential_backend!(backend)
            X = [-1.0, -1.0]
            @test _gradient(f1, [1.0, -1.0]) ≈ [1.0, -2.0]
            @test _gradient!(f1, X, [1.0, -1.0]) === X
            @test X ≈ [1.0, -2.0]
        end
        @testset for backend in [finite_diff]
            set_default_differential_backend!(backend)
            X = [-0.0 -0.0]
            @test _jacobian(f1, [1.0, -1.0]) ≈ [1.0 -2.0]
            # The following seems not to worf for :central, but it does for forward
            fdf = Manifolds.FiniteDiffBackend(Val(:forward))
            @test_broken _jacobian!(f1!, X, [1.0, -1.0], fdf) === X
            @test_broken X ≈ [1.0 -2.0]
        end
        set_default_differential_backend!(Manifolds.NoneDiffBackend())
        @testset for backend in [fd51, Manifolds.ForwardDiffBackend()]
            @test _derivative(c1, 0.0, backend) ≈ [1.0, 0.0]
            @test _gradient(f1, [1.0, -1.0], backend) ≈ [1.0, -2.0]
        end

        set_default_differential_backend!(fd51)
    end
end

rb_onb_default = TangentDiffBackend(
    default_differential_backend(),
    Manifolds.ExponentialRetraction(),
    Manifolds.LogarithmicInverseRetraction(),
    DefaultOrthonormalBasis(),
)

rb_onb_fd51 = TangentDiffBackend(Manifolds.FiniteDifferencesBackend())

rb_onb_fwd_diff = TangentDiffBackend(Manifolds.ForwardDiffBackend())

rb_onb_finite_diff = TangentDiffBackend(Manifolds.FiniteDiffBackend())

rb_onb_default2 = TangentDiffBackend(
    default_differential_backend();
    basis=CachedBasis(
        DefaultOrthonormalBasis(),
        [[0.0, -1.0, 0.0], [sqrt(2) / 2, 0.0, -sqrt(2) / 2]],
    ),
)

rb_proj = Manifolds.RiemannianProjectionBackend(default_differential_backend())

@testset "Riemannian differentials" begin
    s2 = Sphere(2)
    p = [0.0, 0.0, 1.0]
    q = [1.0, 0.0, 0.0]
    c1(t) = geodesic(s2, q, p, t)

    Xval = [-sqrt(2) / 2, 0.0, sqrt(2) / 2]
    @test isapprox(s2, c1(π / 4), differential(s2, c1, π / 4, rb_onb_default), Xval)
    X = similar(p)
    differential!(s2, c1, X, π / 4, rb_onb_default)
    @test isapprox(s2, c1(π / 4), X, Xval)

    @testset for backend in [rb_onb_fd51, rb_onb_fwd_diff, rb_onb_finite_diff]
        @test isapprox(s2, c1(π / 4), differential(s2, c1, π / 4, backend), Xval)
        X = similar(p)
        differential!(s2, c1, X, π / 4, backend)
        @test isapprox(s2, c1(π / 4), X, Xval)
    end
end

@testset "Riemannian gradients" begin
    s2 = Sphere(2)
    f1(p) = p[1]

    q = [sqrt(2) / 2, 0, sqrt(2) / 2]
    X = similar(q)
    for backend in [rb_onb_default, rb_proj]
        @test isapprox(s2, q, gradient(s2, f1, q, backend), [0.5, 0.0, -0.5])
        @test gradient!(s2, f1, X, q, backend) === X
        @test isapprox(s2, q, X, [0.5, 0.0, -0.5])
    end
    X = similar(q)
    for backend in [rb_onb_default, rb_proj]
        gradient!(s2, f1, X, q, backend)
        @test isapprox(s2, q, X, [0.5, 0.0, -0.5])
    end

    # Test the gradient fallback
    @test gradient(s2, f1, q, TestRiemannianBackend()) == [1.0, 2.0, 3.0]
    X = similar(q)
    @test gradient!(s2, f1, X, q, TestRiemannianBackend()) === X
    @test X == [1.0, 2.0, 3.0]
end

@testset "Default Errors for the ODEExponentialRetraction" begin
    @test_throws DomainError ODEExponentialRetraction(
        ProjectionRetraction(),
        CachedBasis(DefaultOrthogonalBasis(), []), # not yet supported
    )
    @test_throws DomainError ODEExponentialRetraction(
        ExponentialRetraction(), # not possible
        DefaultOrthogonalBasis(),
    )
    @test_throws DomainError ODEExponentialRetraction(
        ExponentialRetraction(), # not possible
        CachedBasis(DefaultOrthogonalBasis(), []), # combination not possible
    )
end


text:
mutable struct InputSettings
    T::Vector{Int64}
    n_items::Int64
    n_groups::Int64
    groups::Vector{String}
    irt_model::String
    irt_parametrization::String
    irt_D::Float64
    enemy_sets_var::Vector{String}
    friend_sets_var::Vector{String}
    item_use_min::Vector{Int64}
    item_use_max::Vector{Int64}
    length_min::Vector{Int64}
    length_max::Vector{Int64}
    length_weight::Vector{Float64}
    expected_score_var::Vector{String}
    expected_score_pts::Vector{Vector{Float64}}
    expected_score_min::Vector{Vector{Float64}}
    expected_score_max::Vector{Vector{Float64}}
    mean_vars::Vector{Vector{String}}
    mean_vars_min::Vector{Vector{Float64}}
    mean_vars_max::Vector{Vector{Float64}}
    obj_type::String
    obj_points::Vector{Vector{Float64}}
    obj_targets::Vector{Vector{Float64}}
    obj_aux_int::Int64
    obj_aux_float::Float64
    categories::Vector{String}
    InputSettings(
        T,
        n_items,
        n_groups,
        groups,
        irt_model,
        irt_parametrization,
        irt_D,
        enemy_sets_var,
        friend_sets_var,
        item_use_min,
        item_use_max,
        length_min,
        length_max,
        length_weight,
        expected_score_var,
        expected_score_pts,
        expected_score_min,
        expected_score_max,
        mean_vars,
        mean_vars_min,
        mean_vars_max,
        obj_type,
        obj_points,
        obj_targets,
        obj_aux_int,
        obj_aux_float,
        categories,
    ) = new(
        T,
        n_items,
        n_groups,
        groups,
        irt_model,
        irt_parametrization,
        irt_D,
        enemy_sets_var,
        friend_sets_var,
        item_use_min,
        item_use_max,
        length_min,
        length_max,
        length_weight,
        expected_score_var,
        expected_score_pts,
        expected_score_min,
        expected_score_max,
        mean_vars,
        mean_vars_min,
        mean_vars_max,
        obj_type,
        obj_points,
        obj_targets,
        obj_aux_int,
        obj_aux_float,
        categories,
    )
    function InputSettings(;
        T = Int64[],
        n_items = zero(Int64),
        n_groups = zero(Int64),
        groups = String[],
        irt_model = "",
        irt_parametrization = "at-b",
        irt_D = 1.0,
        enemy_sets_var = String[],
        friend_sets_var = String[],
        item_use_min = Int64[],
        item_use_max = Int64[],
        length_min = Int64[],
        length_max = Int64[],
        length_weight = Float64[],
        expected_score_var = String[],
        expected_score_pts = Float64[],
        expected_score_min = Float64[],
        expected_score_max = Float64[],
        mean_vars = Vector{Vector{String}}(undef, 0),
        mean_vars_min = Vector{Vector{Float64}}(undef, 0),
        mean_vars_max = Vector{Vector{Float64}}(undef, 0),
        obj_type = "",
        obj_points = Vector{Vector{Float64}}(undef, 0),
        obj_targets = Vector{Vector{Float64}}(undef, 0),
        obj_aux_int = zero(Int64),
        obj_aux_float = zero(Float64),
        categories = String[],
    )
        InputSettings(
            T,
            n_items,
            n_groups,
            groups,
            irt_model,
            irt_parametrization,
            irt_D,
            enemy_sets_var,
            friend_sets_var,
            item_use_min,
            item_use_max,
            length_min,
            length_max,
            length_weight,
            expected_score_var,
            expected_score_pts,
            expected_score_min,
            expected_score_max,
            mean_vars,
            mean_vars_min,
            mean_vars_max,
            obj_type,
            obj_points,
            obj_targets,
            obj_aux_int,
            obj_aux_float,
            categories,
        )
    end
end


text:
function show{P}(io::IO, x::ArbFloat{P})
    if isexact(x)
      s = string(midpoint(x))
    else
      s = string(x)
    end
    print(io, s)
end

show{P}(x::ArbFloat{P}) = show(Base.STDOUT, x)

function showcompact{P}(io::IO, x::ArbFloat{P})
    s = stringCompact(x)
    print(io, s)
end

function showall{P}(io::IO, x::ArbFloat{P})
    s = stringAll(x)
    print(io, s)
end

function showallcompact{P}(io::IO, x::ArbFloat{P})
    s = stringAllCompact(x)
    print(io, s)
end

function showsmart{P}(io::IO, x::ArbFloat{P})
    s = smartstring(x)
    print(io, s)
end
# showsmart is not a Base show function, it needs explict version without io parameter
showsmart{P}(x::ArbFloat{P}) = showsmart(Base.STDOUT, x)

function showmany{P,N}(io::IO, x::NTuple{N,ArbFloat{P}}, stringformer::Function)
    if N==0
       print(io,"()")
       return nothing
    elseif N==1
       print(io,string("( ",stringformer(x[1]),", )"))
       return nothing
    end

    ss = Vector{String}(N)
    for i in 1:N
      ss[i] = stringformer(x[i])
    end

    println(io,string("( ", ss[1], ","));
    for s in ss[2:end-1]
      println(io, string("  ", s, ","))
    end
    println(io,string("  ", ss[end], " )"))
end

showmany{P,N}(x::NTuple{N,ArbFloat{P}}, stringformer::Function) =
   showmany(Base.STDOUT,x,stringformer)

function showmany{P}(io::IO, x::Vector{ArbFloat{P}}, stringformer::Function)
    n = length(x)

    if n==0
       print(io,"[]")
       return nothing
    elseif n==1
       print(io,string("[ ",stringformer(x[1])," ]"))
       return nothing
    end

    ss = Vector{String}(n)
    for i in 1:n
      ss[i] = stringformer(x[i])
    end

    println(io,string("[ ", ss[1], ","));
    for s in ss[2:end-1]
      println(io, string("  ", s, ","))
    end
    println(io,string("  ", ss[end], " ]"))
end

showmany{P}(x::Vector{ArbFloat{P}}, stringformer::Function) =
    showmany(Base.STDOUT,x,stringformer)

show{P,N}(io::IO, x::NTuple{N,ArbFloat{P}}) = showmany(io, x, string)
showall{P,N}(io::IO, x::NTuple{N,ArbFloat{P}}) = showmany(io, x, stringAll)
showcompact{P,N}(io::IO, x::NTuple{N,ArbFloat{P}}) = showmany(io, x, stringCompact)
# showsmart is not a Base show function, it needs explict version without io parameter
showsmart{P,N}(io::IO, x::NTuple{N,ArbFloat{P}}) = showmany(io, x, smartstring)
showsmart{P,N}(x::NTuple{N,ArbFloat{P}}) = showmany(Base.STDOUT, x, smartstring)

show{P}(io::IO, x::Vector{ArbFloat{P}}) = showmany(io, x, string)
showall{P}(io::IO, x::Vector{ArbFloat{P}}) = showmany(io, x, stringAll)
showcompact{P}(io::IO, x::Vector{ArbFloat{P}}) = showmany(io, x, stringCompact)
# showsmart is not a Base show function, it needs explict version without io parameter
showsmart{P}(io::IO, x::Vector{ArbFloat{P}}) = showmany(io, x, smartstring)
showsmart{P}(x::Vector{ArbFloat{P}}) = showmany(Base.STDOUT, x, smartstring)

show{P,N}(io::IO, x::Vararg{ArbFloat{P},N}) = showmany(io, x, string)
showall{P,N}(io::IO, x::Vararg{ArbFloat{P},N}) = showmany(io, x, stringAll)
showcompact{P,N}(io::IO, x::Vararg{ArbFloat{P},N}) = showmany(io, x, stringCompact)
# showsmart is not a Base show function, it needs explict version without io parameter
showsmart{P,N}(io::IO, x::Vararg{ArbFloat{P},N}) = showmany(io, x, smartstring)
showsmart{P,N}(x::Vararg{ArbFloat{P},N}) = showmany(Base.STDOUT, x, smartstring)


text:
using LinearAlgebra
"""
This code is loosely based on a routine originally copyright Federal Reserve Bank of Atlanta
and written by Iskander Karibzhanov.
"""
mutable struct KalmanFilter{S<:Number}
    T::AbstractMatrix{S}
    R::AbstractMatrix{S}
    C::AbstractVector{S}
    Q::AbstractMatrix{S}
    Z::AbstractMatrix{S}
    D::AbstractVector{S}
    E::AbstractMatrix{S}
    s_t::AbstractVector{S} # s_{t|t-1} or s_{t|t}
    P_t::AbstractMatrix{S} # P_{t|t-1} or P_{t|t}
    loglh_t::U where {U<:Number}     # P(y_t | y_{1:t})
    converged::Bool
    PZV::AbstractMatrix{S}
end

"""
```
KalmanFilter(T, R, C, Q, Z, D, E, [s_0, P_0])
```

Outer constructor for the `KalmanFilter` type.
"""
function KalmanFilter(T::AbstractMatrix{S}, R::AbstractMatrix{S}, C::AbstractVector{S}, Q::AbstractMatrix{S},
                      Z::AbstractMatrix{S}, D::AbstractVector{S}, E::AbstractMatrix{S},
                      s_0::AbstractVector{S} = Vector{S}(undef, 0),
                      P_0::AbstractMatrix{S} = Matrix{S}(undef, 0, 0),
                      converged::Bool = false,
                      PZV::AbstractMatrix{S} = Matrix{S}(undef, 0, 0)) where {S<:Real}
    if isempty(s_0) || isempty(P_0)
        if issparse(T) # can't call eigvals on a sparse matrix
            s_0, P_0 = init_stationary_states(Matrix(T), R, C, Q)
        else
            s_0, P_0 = init_stationary_states(T, R, C, Q)
        end
    end

    return KalmanFilter(T, R, C, Q, Z, D, E, s_0, P_0, NaN, converged, PZV)
end

"""
```
init_stationary_states(T, R, C, Q)
```

Compute the initial state `s_0` and state covariance matrix `P_0` under the
stationarity condition:

```
s_0  =  (I - T) \\ C
P_0 = reshape(I - kron(T, T)) \\ vec(R*Q*R'), Ns, Ns)
```

where:

- `kron(T, T)` is a matrix of dimension `Ns^2` x `Ns^2`, the Kronecker
  product of `T`
- `vec(R*Q*R')` is the `Ns^2` x 1 column vector constructed by stacking the
  `Ns` columns of `R*Q*R'`

All eigenvalues of `T` are inside the unit circle when the state space model
is stationary. When the preceding formula cannot be applied, the initial state
vector estimate is set to `C` and its covariance matrix is given by `1e6 * I`.
"""
function init_stationary_states(T::AbstractMatrix{S}, R::AbstractMatrix{S}, C::AbstractVector{S},
                                Q::AbstractMatrix{S};
                                check_is_stationary::Bool = true) where {S<:Real}
    if check_is_stationary
        e = eigvals(T)
        if all(abs.(e) .< 1)
            s_0 = (UniformScaling(1) - T)\C
            P_0 = solve_discrete_lyapunov(T, R*Q*R')
        else
            Ns = size(T, 1)
            s_0 = C
            P_0 = 1e6 * Matrix(1.0I, Ns, Ns)
        end
    else
        s_0 = (Matrix{eltype(T)}(I, size(T)...) - T)\C
        P_0 = solve_discrete_lyapunov(T, R*Q*R')
    end
    return s_0, P_0
end

function init_stationary_states(T::TrackedMatrix{S}, R::TrackedMatrix{S}, C::TrackedVector{S},
                                Q::TrackedMatrix{S};
                                check_is_stationary::Bool = true) where {S<:Real}
    mat_type = return_tracker_parameter_type(T)
    if check_is_stationary
        e = eigvals(T)
        if all(abs.(e) .< 1)
            s_0 = (UniformScaling(1) - T)\C
            P_0 = solve_discrete_lyapunov(T, R*Q*R')
        else
            Ns = size(T, 1)
            s_0 = C
            P_0 = 1e6 * Matrix(1.0I, Ns, Ns)
        end
    else
        s_0 = (Matrix{mat_type}(I, size(T)...) - T)\C
        P_0 = solve_discrete_lyapunov(T, R*Q*R')
    end
    return s_0, P_0
end

"""
```
kalman_filter(y, T, R, C, Q, Z, D, E, s_0 = AbstractVector(), P_0 = Matrix();
    outputs = [:loglh, :pred, :filt], Nt0 = 0)

kalman_filter(regime_indices, y, Ts, Rs, Cs, Qs, Zs, Ds, Es,
    s_0 = AbstractVector(), P_0 = Matrix(); outputs = [:loglh, :pred, :filt],
    Nt0 = 0)
```

This function implements the Kalman filter for the following state-space model:

```
s_{t+1} = C + T*s_t + R*ϵ_t    (transition equation)
y_t     = D + Z*s_t + u_t      (measurement equation)

ϵ_t ∼ N(0, Q)
u_t ∼ N(0, E)
Cov(ϵ_t, u_t) = 0
```

### Inputs

- `y`: `Ny` x `Nt` matrix containing data `y_1, ... , y_T`
- `s_0`: optional `Ns` x 1 initial state vector
- `P_0`: optional `Ns` x `Ns` initial state covariance matrix

**Method 1 only:**

- `T`: `Ns` x `Ns` state transition matrix
- `R`: `Ns` x `Ne` matrix in the transition equation mapping shocks to states
- `C`: `Ns` x 1 constant vector in the transition equation
- `Q`: `Ne` x `Ne` matrix of shock covariances
- `Z`: `Ny` x `Ns` matrix in the measurement equation mapping states to
  observables
- `D`: `Ny` x 1 constant vector in the measurement equation
- `E`: `Ny` x `Ny` matrix of measurement error covariances

**Method 2 only:**

- `regime_indices`: `AbstractVector{UnitRange{Int}}` of length `n_regimes`, where
  `regime_indices[i]` indicates the time periods `t` in regime `i`
- `Ts`: `AbstractVector{Matrix{S}}` of `T` matrices for each regime
- `Rs`
- `Cs`
- `Qs`
- `Zs`
- `Ds`
- `Es`

where:

- `Nt`: number of time periods for which we have data
- `Ns`: number of states
- `Ne`: number of shocks
- `Ny`: number of observables

### Keyword Arguments

- `outputs`: some subset of `[:loglh, :pred, :filt]` specifying which outputs to
  compute and return. There will always be the same number of return values,
  but, for example, `s_pred` and `P_pred` will be returned as empty arrays if
  `:pred` is not in `outputs`
- `Nt0`: number of presample periods to omit from all return values

### Outputs

- `loglh`: length `Nt` vector of conditional log-likelihoods P(y_t | y_{1:t-1})
- `s_pred`: `Ns` x `Nt` matrix of one-step predicted state vectors s_{t|t-1}
- `P_pred`: `Ns` x `Ns` x `Nt` array of mean squared errors P_{t|t-1} of
  predicted state vectors
- `s_filt`: `Ns` x `Nt` matrix of filtered state vectors s_{t|t}
- `P_filt`: `Ns` x `Ns` x `Nt` matrix containing mean squared errors P_{t|t} of
  filtered state vectors
- `s_0`: `Ns` x 1 initial state vector. This may have been reassigned to the
  last presample state vector if `Nt0 > 0`
- `P_0`: `Ns` x `Ns` initial state covariance matrix. This may have been
  reassigned to the last presample state covariance if `Nt0 > 0`
- `s_T`: `Ns` x 1 final filtered state `s_{T|T}`
- `P_T`: `Ns` x `Ns` final filtered state covariance matrix `P_{T|T}`

### Notes

When `s_0` and `P_0` are omitted, they are computed using
`init_stationary_states`.
"""
function kalman_filter(regime_indices::Vector{UnitRange{Int}}, y::AbstractArray,
                       Ts::Vector{<:AbstractMatrix{S}}, Rs::Vector{<:AbstractMatrix{S}},
                       Cs::Vector{<:AbstractVector{S}}, Qs::Vector{<:AbstractMatrix{S}},
                       Zs::Vector{<:AbstractMatrix{S}}, Ds::Vector{<:AbstractVector{S}},
                       Es::Vector{<:AbstractMatrix{S}}, s_0::AbstractVector{S} = Vector{S}(undef, 0),
                       P_0::AbstractMatrix{S} = Matrix{S}(undef, 0, 0);
                       outputs::AbstractVector{Symbol} = [:loglh, :pred, :filt],
                       Nt0::Int = 0, tol::AbstractFloat = 0.0) where {S<:Real}

    # Determine outputs
    return_loglh = :loglh in outputs
    return_pred  = :pred  in outputs
    return_filt  = :filt  in outputs

    # Dimensions
    Ns = size(Ts[1], 1) # number of states
    Nt = size(y, 2)     # number of periods of data

    # Initialize inputs and outputs
    k = KalmanFilter(Ts[1], Rs[1], Cs[1], Qs[1], Zs[1], Ds[1], Es[1], s_0, P_0)

    mynan = convert(S, NaN)
    s_pred = return_pred  ? fill(mynan, Ns, Nt)     : Matrix{S}(undef, 0, 0)
    P_pred = return_pred  ? fill(mynan, Ns, Ns, Nt) : Array{S, 3}(undef, 0, 0, 0)
    s_filt = return_filt  ? fill(mynan, Ns, Nt)     : Matrix{S}(undef, 0, 0)
    P_filt = return_filt  ? fill(mynan, Ns, Ns, Nt) : Array{S, 3}(undef, 0, 0, 0)
    loglh  = return_loglh ? fill(mynan, Nt)         : Vector{S}(undef, 0)

    # Populate s_0 and P_0
    s_0 = k.s_t
    P_0 = k.P_t

    # Iterate through regimes
    s_t = s_0
    P_t = P_0

    for i = 1:length(regime_indices)
        ts = regime_indices[i]

        loglh_i, s_pred_i, P_pred_i, s_filt_i, P_filt_i, s_0, P_0, s_t, P_t =
            kalman_filter(y[:, ts], Ts[i], Rs[i], Cs[i], Qs[i], Zs[i], Ds[i], Es[i],
                          s_t, P_t; outputs = outputs, Nt0 = 0, tol = tol)

        if return_loglh
            loglh[ts] = loglh_i
        end
        if return_pred
            s_pred[:,    ts] = s_pred_i
            P_pred[:, :, ts] = P_pred_i
        end
        if return_filt
            s_filt[:,    ts] = s_filt_i
            P_filt[:, :, ts] = P_filt_i
        end
    end

    # Populate s_T and P_T
    s_T, P_T = s_t, P_t

    # Remove presample periods
    loglh, s_pred, P_pred, s_filt, P_filt =
        remove_presample!(Nt0, loglh, s_pred, P_pred, s_filt, P_filt; outputs = outputs)

    return loglh, s_pred, P_pred, s_filt, P_filt, s_0, P_0, s_T, P_T
end

function kalman_filter(y::AbstractArray, T::AbstractMatrix{S}, R::AbstractMatrix{S},
                       C::AbstractVector{S}, Q::AbstractMatrix{S}, Z::AbstractMatrix{S},
                       D::AbstractVector{S}, E::AbstractMatrix{S},
                       s_0::AbstractVector{S} = Vector{S}(undef, 0),
                       P_0::AbstractMatrix{S} = Matrix{S}(undef, 0, 0);
                       outputs::AbstractVector{Symbol} = [:loglh, :pred, :filt],
                       Nt0::Int = 0, tol::AbstractFloat = 0.0) where {S<:Real}

    # Determine outputs
    return_loglh = :loglh in outputs
    return_pred  = :pred  in outputs
    return_filt  = :filt  in outputs

    # Dimensions
    Ns = size(T, 1) # number of states
    Nt = size(y, 2) # number of periods of data

    # Initialize inputs and outputs
    k = KalmanFilter(T, R, C, Q, Z, D, E, s_0, P_0)

    mynan = convert(S, NaN)
    loglh  = return_loglh ? fill(mynan, Nt)         : Vector{S}(undef, 0)
    s_pred = return_pred  ? fill(mynan, Ns, Nt)     : Matrix{S}(undef, 0, 0)
    P_pred = return_pred  ? fill(mynan, Ns, Ns, Nt) : Array{S, 3}(undef, 0, 0, 0)
    s_filt = return_filt  ? fill(mynan, Ns, Nt)     : Matrix{S}(undef, 0, 0)
    P_filt = return_filt  ? fill(mynan, Ns, Ns, Nt) : Array{S, 3}(undef, 0, 0, 0)

    # Populate initial states
    s_0 = k.s_t
    P_0 = k.P_t

    # Loop through periods t
    RQR′ = R * Q * R'
    for t = 1:Nt
        # Forecast
        forecast!(k, RQR′) # calculate RQR′ once for efficiency
        if return_pred
            s_pred[:,    t] = k.s_t
            P_pred[:, :, t] = k.P_t
        end

        # Update and compute log-likelihood
        update!(k, y[:, t]; return_loglh = return_loglh, tol = tol)
        if return_filt
            s_filt[:,    t] = k.s_t
            P_filt[:, :, t] = k.P_t
        end
        if return_loglh
            loglh[t]        = k.loglh_t
        end

        # Update s_0 and P_0 if Nt0 > 0
        if t == Nt0
            s_0 = k.s_t
            P_0 = k.P_t
        end
    end

    # Populate final states
    s_T = k.s_t
    P_T = k.P_t

    # Remove presample periods
    loglh, s_pred, P_pred, s_filt, P_filt =
        remove_presample!(Nt0, loglh, s_pred, P_pred, s_filt, P_filt; outputs = outputs)

    return loglh, s_pred, P_pred, s_filt, P_filt, s_0, P_0, s_T, P_T
end

function kalman_likelihood(regime_indices::Vector{UnitRange{Int}}, y::AbstractArray,
                           Ts::Vector{<:AbstractMatrix{S}}, Rs::Vector{<:AbstractMatrix{S}},
                           Cs::Vector{<:AbstractVector{S}}, Qs::Vector{<:AbstractMatrix{S}},
                           Zs::Vector{<:AbstractMatrix{S}}, Ds::Vector{<:AbstractVector{S}},
                           Es::Vector{<:AbstractMatrix{S}}, s_0::AbstractVector{S} = Vector{S}(undef, 0),
                           P_0::AbstractMatrix{S} = Matrix{S}(undef, 0, 0);
                           Nt0::Int = 0, tol::AbstractFloat = 0.0) where {S<:Real}
    # Dimensions
    Nt = size(y, 2) # number of periods of data

    @assert first(regime_indices[1]) == 1
    @assert last(regime_indices[end]) == Nt

    # Initialize inputs and outputs
    k = KalmanFilter(Ts[1], Rs[1], Cs[1], Qs[1], Zs[1], Ds[1], Es[1], s_0, P_0)

    mynan = convert(S, NaN)
    loglh = fill(mynan, Nt)

    # Iterate through regimes
    s_t = k.s_t
    P_t = k.P_t

    if length(regime_indices) == 1
        loglh = kalman_likelihood(y, Ts[1], Rs[1], Cs[1], Qs[1],
                                  Zs[1], Ds[1], Es[1], s_t, P_t;
                                  Nt0 = 0, tol = tol, switching = false)
    else
        for i = 1:length(regime_indices)
            ts = regime_indices[i]
            loglh[ts], s_t, P_t = kalman_likelihood(y[:, ts], Ts[i], Rs[i], Cs[i], Qs[i],
                                                    Zs[i], Ds[i], Es[i], s_t, P_t;
                                                    Nt0 = 0, tol = tol, switching = true)
        end
    end

    # Remove presample periods
    loglh = remove_presample!(Nt0, loglh)

    return loglh
end

function kalman_likelihood(y::AbstractArray, T::AbstractMatrix{S}, R::AbstractMatrix{S}, C::AbstractVector{S},
                           Q::AbstractMatrix{S}, Z::AbstractMatrix{S}, D::AbstractVector{S}, E::AbstractMatrix{S},
                           s_0::AbstractVector{S} = Vector{S}(undef, 0),
                           P_0::AbstractMatrix{S} = Matrix{S}(undef, 0, 0);
                           Nt0::Int = 0, tol::AbstractFloat = 0.0,
                           switching::Bool = true) where {S<:Real}
    # Dimensions
    Nt = size(y, 2) # number of periods of data

    # Initialize inputs and outputs
    k = KalmanFilter(T, R, C, Q, Z, D, E, s_0, P_0)

    mynan = convert(S, NaN)
    loglh = fill(mynan, Nt)

    # Loop through periods t
    RQR′ = R * Q * R'
    for t = 1:Nt
        # Forecast
        forecast!(k, RQR′) # calculate RQR′ once for efficiency

        # Update and compute log-likelihood
        update!(k, y[:, t]; return_loglh = true, tol = tol)
        loglh[t] = k.loglh_t
    end

    # Remove presample periods
    loglh = remove_presample!(Nt0, loglh)

    if switching
        return loglh, k.s_t, k.P_t
    else
        return loglh
    end
end

function kalman_likelihood(y::AbstractArray, T::TrackedMatrix{S}, R::TrackedMatrix{S},
                           C::TrackedVector{S}, Q::TrackedMatrix{S}, Z::TrackedMatrix{S},
                           D::TrackedVector{S}, E::TrackedMatrix{S},
                           s_0::TrackedVector{S} = param(Vector{S}(undef, 0)),
                           P_0::TrackedMatrix{S} = param(Matrix{S}(undef, 0, 0));
                           Nt0::Int = 0, tol::AbstractFloat = 0.0) where {S<:Real}
    # Dimensions
    Nt = size(y, 2) # number of periods of data

    # Initialize inputs and outputs
    k = KalmanFilter(T, R, C, Q, Z, D, E, s_0, P_0)

    mynan = convert(S, NaN)
    loglh = fill(Tracker.TrackedReal{S}(mynan), Nt)

    # Loop through periods t
    RQR′ = R * Q * R'
    for t = 1:Nt
        # Forecast
        forecast!(k, RQR′) # calculate RQR′ once for efficiency

        # Update and compute log-likelihood
        update!(k, y[:, t]; return_loglh = true, tol = tol)
        loglh[t] = k.loglh_t
    end

    # Remove presample periods
    loglh = remove_presample!(Nt0, loglh)

    return loglh
end

"""
```
forecast!(k::KalmanFilter)
```

Compute the one-step-ahead states s_{t|t-1} and state covariances P_{t|t-1} and
assign to `k`.
"""
function forecast!(k::KalmanFilter{S}, RQR′::AbstractMatrix = k.R * k.Q * k.R') where {S<:Real}
    T, C           = k.T, k.C
    s_filt, P_filt = k.s_t, k.P_t

    k.s_t = T*s_filt + C         # s_{t|t-1} = T*s_{t-1|t-1} + C
    k.P_t = T*P_filt*T' + RQR′   # P_{t|t-1} = Var s_{t|t-1} = T*P_{t-1|t-1}*T' + R*Q*R'
    return nothing
end

"""
```
update!(k::KalmanFilter{S}, y_obs)
```

Compute the filtered states s_{t|t} and state covariances P_{t|t}, and the
log-likelihood P(y_t | y_{1:t-1}) and assign to `k`.
"""
function update!(k::KalmanFilter{S}, y_obs::AbstractArray;
                 return_loglh::Bool = true, tol::AbstractFloat = 0.0) where {S<:Real}

    # Keep rows of measurement equation corresponding to non-missing observables
    nonmissing = .![ismissing(x) ? true : isnan(x) for x in y_obs]

    y_obs = y_obs[nonmissing]
    Z     = k.Z[nonmissing, :]
    D     = k.D[nonmissing]
    E     = k.E[nonmissing, nonmissing]

    Ny    = length(y_obs)

    s_pred = k.s_t
    P_pred = k.P_t

    y_pred = Z*s_pred + D             # y_{t|t-1} = Z*s_{t|t-1} + D

    V_pred     = Z*P_pred*Z' + E      # V_{t|t-1} = Var y_{t|t-1} = Z*P_{t|t-1}*Z' + E
    V_pred     = (V_pred + V_pred')/2 # V_pred should be symmetric; this guarantees symmetry and divides by 2 so entries aren't double
    V_pred_inv = inv(V_pred)
    dy         = y_obs - y_pred       # dy = y_t - y_{t|t-1} (prediction error)

    if !k.converged
        PZV = P_pred'*Z'*V_pred_inv
    end

    if size(k.PZV) == size(PZV)
        if (tol > 0.0) && (maximum(abs.(PZV - k.PZV)) < tol)
            k.converged = true
        end
    end
    k.PZV = PZV

    k.s_t = s_pred + PZV*dy       # s_{t|t} = s_{t|t-1} + P_{t|t-1}'*Z'/V_{t|t-1}*dy
    k.P_t = P_pred - PZV*Z*P_pred # P_{t|t} = P_{t|t-1} - P_{t|t-1}'*Z'/V_{t|t-1}*Z*P_{t|t-1}

    if return_loglh
        # p(y_t | y_{1:t-1})
        k.loglh_t = -(Ny*log(2π) + log(det(V_pred)) + dy'*V_pred_inv*dy)/2
    end
    return nothing
end


text:
using Test,SphericalHarmonicModes
using SphericalHarmonicModes: ModeMissingError, SHModeRange,
check_if_lm_range_is_valid, flip

using Aqua
@testset "project quality" begin
    Aqua.test_all(SphericalHarmonicModes)
end

promote_typeof(x, y) = mapreduce(typeof, promote_type, (x, y))

@testset "constructors" begin

	@testset "Special Ranges" begin
	    @testset "ZeroTo" begin
	        r = ZeroTo(3)
	        @test first(r) == 0
	        @test last(r) == 3
	        @test !isempty(r)

	        r = ZeroTo(UInt8(3))
			@test first(r) == 0
	        @test last(r) == 3
	        @test !isempty(r)

            for rin in Any[0:2, ZeroTo(2), ZeroTo{true}(2), ZeroTo{false}(2)]
                r = ZeroTo(rin)
                @test first(r) == 0
                @test last(r) == 2
                r = ZeroTo{false}(rin)
                @test first(r) == 0
                @test last(r) == 2
            end
            for rin in Any[FullRange(0), ToZero(0), SingleValuedRange(0)]
                r = ZeroTo(rin)
                @test r isa ZeroTo
                @test first(r) == last(r) == 0
                r = ZeroTo{false}(rin)
                @test r isa ZeroTo{false}
                @test first(r) == last(r) == 0
            end
            @test_throws ArgumentError ZeroTo(1:2)
            @test_throws Exception ZeroTo{true}(ZeroTo(1))
            for r in Any[ToZero(1), FullRange(1), SingleValuedRange(1)]
                @test_throws Exception ZeroTo(r)
                @test_throws Exception ZeroTo{false}(r)
                @test_throws Exception ZeroTo{true}(r)
            end
	    end
	    @testset "FullRange" begin
	        r = FullRange(3)
	        @test first(r) == -3
	        @test last(r) == 3
	        @test !isempty(r)

            r = FullRange(-1:1)
            @test first(r) == -1
            @test last(r) == 1
            @test_throws ArgumentError FullRange(0:1)

            for rin in Any[-2:2, FullRange(2), FullRange{false}(2), FullRange{true}(2)]
                r = FullRange(rin)
                @test r isa FullRange
                @test first(r) == -2
                @test last(r) == 2
                r = FullRange{false}(rin)
                @test r isa FullRange{false}
                @test first(r) == -2
                @test last(r) == 2
            end
            for rin in Any[ZeroTo(0), ToZero(0), SingleValuedRange(0)]
                r = FullRange(rin)
                @test r isa FullRange
                @test first(r) == last(r) == 0
                r = FullRange{false}(rin)
                @test r isa FullRange{false}
                @test first(r) == last(r) == 0
            end
            @test_throws ArgumentError FullRange(1:2)
            @test_throws Exception FullRange{true}(FullRange(1))
            for r in Any[ToZero(1), ZeroTo(1), SingleValuedRange(1)]
                @test_throws Exception FullRange(r)
                @test_throws Exception FullRange{false}(r)
                @test_throws Exception FullRange{true}(r)
            end
	    end
	    @testset "ToZero" begin
	        r = ToZero(3)
	        @test first(r) == -3
	        @test last(r) == 0
	        @test !isempty(r)

            r = ToZero(-2:0)
            @test first(r) == -2
            @test last(r) == 0
            @test_throws ArgumentError ToZero(-2:-1)

            for rin in Any[-2:0, ToZero(2), ToZero{true}(2), ToZero{false}(2)]
                r = ToZero(rin)
                @test r isa ToZero
                @test first(r) == -2
                @test last(r) == 0
                r = ToZero{false}(rin)
                @test r isa ToZero{false}
                @test first(r) == -2
                @test last(r) == 0
            end
            @test_throws ArgumentError ToZero(1:2)
            @test_throws Exception ToZero{true}(ToZero(1))
            for rin in Any[ZeroTo(0), FullRange(0), SingleValuedRange(0)]
                r = ToZero(rin)
                @test r isa ToZero
                @test first(r) == last(r) == 0
                r = ToZero{false}(rin)
                @test r isa ToZero{false}
                @test first(r) == last(r) == 0
            end
            for r in Any[ZeroTo(1), FullRange(1), SingleValuedRange(1)]
                @test_throws Exception ToZero(r)
                @test_throws Exception ToZero{false}(r)
                @test_throws Exception ToZero{true}(r)
            end
	    end
	    @testset "SingleValuedRange" begin
	    	n = 3
	        r = SingleValuedRange(n)
	        @test r == n:n
	        @test first(r) == n
	        @test last(r) == n
	        @test length(r) == 1
	        @test !isempty(r)

            r = SingleValuedRange(n:n)
            @test first(r) == n
            @test last(r) == n
            @test_throws ArgumentError SingleValuedRange(1:2)

            for rin in Any[1:1, SingleValuedRange(1)]
                r = SingleValuedRange(rin)
                @test first(r) == 1
                @test last(r) == 1
            end
            @test_throws ArgumentError SingleValuedRange(1:2)
            for rin in Any[ZeroTo(0), FullRange(0), ToZero(0)]
                r = SingleValuedRange(rin)
                @test first(r) == last(r) == 0
            end
            for r in Any[ZeroTo(1), FullRange(1), ToZero(1)]
                @test_throws Exception SingleValuedRange(r)
            end
	    end
        @testset "promotion" begin
            @test promote_typeof(ZeroTo(1), ZeroTo(1)) == typeof(ZeroTo(1))
            @test promote_typeof(ZeroTo(1), SingleValuedRange(1)) == UnitRange{Int}
            @test promote_typeof(FullRange(1), SingleValuedRange(1)) == UnitRange{Int}
            @test promote_typeof(FullRange(1), ZeroTo(1)) == UnitRange{Int}
            @test promote_typeof(SingleValuedRange(1), SingleValuedRange(2)) == SingleValuedRange
        end
	end

	@testset "LM" begin
		@test eltype(LM(1:2)) == Tuple{Int,Int}
		@test eltype(LM(big(1):big(2))) == Tuple{BigInt,Int}

		@test LM{UnitRange{Int},UnitRange{Int}}(2:3,2:2) == LM(2:3,2:2)
		@test LM(1:3,2:2) == LM(2:3,2:2)
		@test LM(1:3) == LM(1:3, -3:3)
		@test LM(1:3, FullRange) == LM(1:3, -3:3)
		@test LM(1:3, ZeroTo) == LM(1:3, 0:3)
		@test LM(1:3, ToZero) == LM(1:3, -3:0)
		@test LM(ZeroTo(3)) == LM(0:3, -3:3)
		@test LM(ZeroTo(3), 1:2) == LM(0:3, 1:2)
		@test LM(ZeroTo(3), FullRange) == LM(0:3, -3:3)
		@test LM(ZeroTo(3), ZeroTo) == LM(0:3, 0:3)
		@test LM(ZeroTo(3), ToZero) == LM(0:3, -3:0)

		@test LM(SingleValuedRange(3), 1:2) == LM(3:3, 1:2)
        @test LM(3, 1:2) == LM(3:3, 1:2)
        @test LM(3) == LM(3:3, -3:3)
        @test LM(3, 2) == LM(3:3, 2:2)

		@test_throws ArgumentError LM(1:0)
		@test_throws ArgumentError LM(1:1, 1:0)
		@test_throws ArgumentError LM(-1:1)
		@test_throws ArgumentError LM(0:1, 1:3)
		@test_throws ArgumentError LM(0:1, 3:3)
		@test_throws ArgumentError LM(0:1, -3:0)
		@test_throws ArgumentError LM(SingleValuedRange(1), 3:4)
		@test_throws ArgumentError LM(SingleValuedRange(1), ZeroTo(3))
		@test_throws ArgumentError LM(SingleValuedRange(1), ToZero(3))
		@test_throws ArgumentError LM(SingleValuedRange(1), FullRange(3))

		@testset "type-stability" begin
		    for LT in (:ZeroTo, :SingleValuedRange),
		    	MT in (:ZeroTo, :ToZero, :FullRange)
		    	@eval @test LM($LT(3), $MT) isa LM{<:$LT, <:$MT}
		    end
		    @test LM(ZeroTo(3), SingleValuedRange(2)) isa LM{UnitRange{Int}, SingleValuedRange}
		    @test LM(SingleValuedRange(3), SingleValuedRange(2)) isa LM{SingleValuedRange, SingleValuedRange}
		end

        @testset "promotion" begin
            @test promote_typeof(LM(1,1), LM(2,2)) == LM{SingleValuedRange, SingleValuedRange}
            @test promote_typeof(LM(1), LM(1)) == LM{SingleValuedRange, FullRange{true}}
            @test promote_typeof(LM(1:1), LM(1)) == LM{UnitRange{Int64}, FullRange{true}}
            @test promote_typeof(LM(1:1, ZeroTo), LM(1, ZeroTo)) == LM{UnitRange{Int}, ZeroTo{true}}
            @test promote_typeof(LM(1:1, ZeroTo), LM(1, ZeroTo(1))) == LM{UnitRange{Int}, UnitRange{Int}}
            @test promote_typeof(LM(1:1, 1:1), LM(1,1)) == LM{UnitRange{Int}, UnitRange{Int}}

            iterslist = Any[LM(1), LM(1:2), LM(1:2, 0:0), LM(1:1, ZeroTo), LM(1:2, ZeroTo(1)),
                        LM(1, ZeroTo), LM(1:1, FullRange), LM(1, FullRange), LM(ZeroTo(3), ZeroTo),
                        LM(ZeroTo(2)), LM(ZeroTo(2), FullRange(1))]

            for x in iterslist, y in iterslist
                @test isconcretetype(promote_typeof(x, y))
            end
        end
        @testset "conversion" begin
            for lr in Any[0:1, ZeroTo(1), SingleValuedRange(1)], mr in Any[0:1, ZeroTo(1), ToZero(1), SingleValuedRange(1)]
                convert(LM{UnitRange{Int}, UnitRange{Int}}, LM(lr, mr)) isa LM{UnitRange{Int}, UnitRange{Int}}
            end
            for LT in [UnitRange{Int}, ZeroTo{false}]
                for MT in [ZeroTo{false}, ZeroTo{true}]
                    @test convert(LM{LT, MT}, LM(0:1, 0:1)) isa LM{LT, MT}
                end
                for MT in [UnitRange{Int}, ToZero{false}, ToZero{true}]
                    @test convert(LM{LT, MT}, LM(0:1, -1:0)) isa LM{LT, MT}
                end
                for MT in [FullRange{false}, FullRange{true}]
                    @test convert(LM{LT, MT}, LM(0:1, -1:1)) isa LM{LT, MT}
                end
                @test convert(LM{LT, SingleValuedRange}, LM(0:1, 0:0)) isa LM{LT, SingleValuedRange}
            end
            for LT in [UnitRange{Int}, SingleValuedRange]
                for MT in [ZeroTo{false}, ZeroTo{true}]
                    @test convert(LM{LT, MT}, LM(1:1, 0:1)) isa LM{LT, MT}
                end
                for MT in [ToZero{false}, ToZero{true}]
                    @test convert(LM{LT, MT}, LM(1:1, -1:0)) isa LM{LT, MT}
                end
                for MT in [FullRange{false}, FullRange{true}]
                    @test convert(LM{LT, MT}, LM(1:1, -1:1)) isa LM{LT, MT}
                end
                @test convert(LM{LT, SingleValuedRange}, LM(1:1, 1:1)) isa LM{LT, SingleValuedRange}
            end
        end
	end

	@testset "ML" begin
		@test eltype(ML(1:2)) == Tuple{Int,Int}
		@test eltype(ML(big(1):big(2))) == Tuple{BigInt,Int}
		@test ML{UnitRange{Int},UnitRange{Int}}(2:3,2:2) == ML(2:3,2:2)
		@test ML(1:3,2:2) == ML(2:3,2:2)
		@test ML(1:3) == ML(1:3, -3:3)
		@test ML(1:3, FullRange) == ML(1:3, -3:3)
		@test ML(1:3, ZeroTo) == ML(1:3, 0:3)
		@test ML(1:3, ToZero) == ML(1:3, -3:0)
		@test ML(ZeroTo(3)) == ML(0:3, -3:3)
		@test ML(ZeroTo(3), 1:2) == ML(0:3, 1:2)
		@test ML(ZeroTo(3), FullRange) == ML(0:3, -3:3)
		@test ML(ZeroTo(3), ZeroTo) == ML(0:3, 0:3)
		@test ML(ZeroTo(3), ToZero) == ML(0:3, -3:0)

		@test ML(SingleValuedRange(3), 1:2) == ML(3:3, 1:2)

		@test_throws ArgumentError ML(1:0)
		@test_throws ArgumentError ML(1:1, 1:0)
		@test_throws ArgumentError ML(-1:1)
		@test_throws ArgumentError ML(0:1, 1:3)
		@test_throws ArgumentError ML(0:1, 3:3)
		@test_throws ArgumentError ML(0:1, -3:0)
		@test_throws ArgumentError ML(SingleValuedRange(1), 3:4)
		@test_throws ArgumentError ML(SingleValuedRange(1), ZeroTo(3))
		@test_throws ArgumentError ML(SingleValuedRange(1), ToZero(3))
		@test_throws ArgumentError ML(SingleValuedRange(1), FullRange(3))

		@testset "type-stability" begin
		    for LT in (:ZeroTo, :SingleValuedRange),
		    	MT in (:ZeroTo, :ToZero, :FullRange)
		    	@eval @test ML($LT(3), $MT) isa ML{<:$LT, <:$MT}
		    end
		    @test ML(ZeroTo(3), SingleValuedRange(2)) isa ML{UnitRange{Int}, SingleValuedRange}
		    @test ML(SingleValuedRange(3), SingleValuedRange(2)) isa ML{SingleValuedRange, SingleValuedRange}
		end

        @testset "promotion" begin
            @test promote_typeof(ML(1,1), ML(2,2)) == ML{SingleValuedRange, SingleValuedRange}
            @test promote_typeof(ML(1), ML(1)) == ML{SingleValuedRange, FullRange{true}}
            @test promote_typeof(ML(1:1), ML(1)) == ML{UnitRange{Int64}, FullRange{true}}
            @test promote_typeof(ML(1:1, ZeroTo), ML(1, ZeroTo)) == ML{UnitRange{Int}, ZeroTo{true}}
            @test promote_typeof(ML(1:1, ZeroTo), ML(1, ZeroTo(1))) == ML{UnitRange{Int}, UnitRange{Int}}
            @test promote_typeof(ML(1:1, 1:1), ML(1,1)) == ML{UnitRange{Int}, UnitRange{Int}}

            iterslist = Any[ML(1), ML(1:2), ML(1:2, 0:0), ML(1:1, ZeroTo), ML(1:2, ZeroTo(1)),
                        ML(1, ZeroTo), ML(1:1, FullRange), ML(1, FullRange), ML(ZeroTo(3), ZeroTo),
                        ML(ZeroTo(2)), ML(ZeroTo(2), FullRange(1))]

            for x in iterslist, y in iterslist
                @test isconcretetype(promote_typeof(x, y))
            end
        end
        @testset "conversion" begin
            for lr in Any[0:1, ZeroTo(1), SingleValuedRange(1)], mr in Any[0:1, ZeroTo(1), ToZero(1), SingleValuedRange(1)]
                convert(ML{UnitRange{Int}, UnitRange{Int}}, ML(lr, mr)) isa ML{UnitRange{Int}, UnitRange{Int}}
            end
            for LT in [UnitRange{Int}, ZeroTo{false}]
                for MT in [ZeroTo{false}, ZeroTo{true}]
                    @test convert(ML{LT, MT}, ML(0:1, 0:1)) isa ML{LT, MT}
                end
                for MT in [UnitRange{Int}, ToZero{false}, ToZero{true}]
                    @test convert(ML{LT, MT}, ML(0:1, -1:0)) isa ML{LT, MT}
                end
                for MT in [FullRange{false}, FullRange{true}]
                    @test convert(ML{LT, MT}, ML(0:1, -1:1)) isa ML{LT, MT}
                end
                @test convert(ML{LT, SingleValuedRange}, ML(0:1, 0:0)) isa ML{LT, SingleValuedRange}
            end
            for LT in [UnitRange{Int}, SingleValuedRange]
                for MT in [ZeroTo{false}, ZeroTo{true}]
                    @test convert(ML{LT, MT}, ML(1:1, 0:1)) isa ML{LT, MT}
                end
                for MT in [ToZero{false}, ToZero{true}]
                    @test convert(ML{LT, MT}, ML(1:1, -1:0)) isa ML{LT, MT}
                end
                for MT in [FullRange{false}, FullRange{true}]
                    @test convert(ML{LT, MT}, ML(1:1, -1:1)) isa ML{LT, MT}
                end
                @test convert(ML{LT, SingleValuedRange}, ML(1:1, 1:1)) isa ML{LT, SingleValuedRange}
            end
        end
	end

	@testset "L2L1Triangle " begin
		Δl_max = 3
		l_range = 1:10
		m = L2L1Triangle(l_range, Δl_max)
		@test m == L2L1Triangle(l_range,LM(0:Δl_max,0:0))
		@test m == L2L1Triangle(l_range,ML(0:Δl_max,0:0))
		@test L2L1Triangle(m) === m
		@test L2L1Triangle(L1L2Triangle(m)) === m

		l2_range = 5:8
		mr = L2L1Triangle(l_range, Δl_max, l2_range)
		@test mr == L2L1Triangle(2:10,3,5:8)
		mr = L2L1Triangle(extrema(l_range)..., Δl_max, l2_range)
		@test mr == L2L1Triangle(2:10,3,5:8)

		l_range = 10:10
		l2_range = 1:8
		mr = L2L1Triangle(l_range, Δl_max, l2_range)
		@test mr == L2L1Triangle(10:10,3,7:8)

		l2_range = 1:20
		mr = L2L1Triangle(l_range, Δl_max, l2_range)
		@test mr == L2L1Triangle(10:10,3,7:13)

		l_range = 1:10
		l2_range = 1:3
		mr = L2L1Triangle(l_range, Δl_max, l2_range)
		@test mr == L2L1Triangle(1:6,3,1:3)

		@test eltype(m) == Tuple{Int,Int}

		@testset "conversion" begin
			@test convert(L1L2Triangle, m) === L1L2Triangle(m)
			@test convert(L2L1Triangle, m) === m
		end
	end

	@testset "L1L2Triangle " begin
		Δl_max = 3
		l_range = 1:10
		m = L1L2Triangle(l_range, Δl_max)
		@test m == L1L2Triangle(l_range,LM(0:Δl_max,0:0))
		@test m == L1L2Triangle(l_range,ML(0:Δl_max,0:0))

		@test L1L2Triangle(m) === m
		@test L1L2Triangle(L2L1Triangle(m)) === m

		l2_range = 5:8
		mr = L1L2Triangle(l_range, Δl_max, l2_range)
		@test mr == L1L2Triangle(2:10,3,5:8)
		mr = L1L2Triangle(extrema(l_range)..., Δl_max, l2_range)
		@test mr == L1L2Triangle(2:10,3,5:8)

		l_range = 10:10
		l2_range = 1:8
		mr = L1L2Triangle(l_range, Δl_max, l2_range)
		@test mr == L1L2Triangle(10:10,3,7:8)

		l2_range = 1:20
		mr = L1L2Triangle(l_range, Δl_max, l2_range)
		@test mr == L1L2Triangle(10:10,3,7:13)

		l_range = 1:10
		l2_range = 1:3
		mr = L1L2Triangle(l_range, Δl_max, l2_range)
		@test mr == L1L2Triangle(1:6,3,1:3)

		@test eltype(m) == Tuple{Int,Int}

		@testset "conversion" begin
			@test convert(L2L1Triangle, m) === L2L1Triangle(m)
			@test convert(L1L2Triangle, m) === m
		end
	end

	@testset "keys" begin
		m = LM(1:2)
		@test keys(m) == Base.OneTo(8)
		@test eachindex(m) == Base.OneTo(8)

		m = ML(1:2)
		@test keys(m) == Base.OneTo(8)
		@test eachindex(m) == Base.OneTo(8)

		m = L2L1Triangle(1:2,1)
		@test keys(m) == Base.OneTo(6)
		@test eachindex(m) == Base.OneTo(6)

		m = L1L2Triangle(1:2,1)
		@test keys(m) == Base.OneTo(6)
		@test eachindex(m) == Base.OneTo(6)
	end
end

@testset "length" begin
	l_cutoff=5

	function iterated_length_template(mr,f)
		try
			return f(mr)
		catch e
			if e isa ArgumentError
				show(mr)
				println()
				return -1
			else
				throw(e)
			end
		end
	end

	function testlength(iterated_length, mr)
		@test begin
			res = length(mr) == iterated_length(mr) == length(collect(mr))
			if !res
				@show mr, typeof(mr)
			end
			res
		end
	end

	iterated_length(mr::ML) = iterated_length_template(mr,
								x->sum(length(m_range(x,l)) for l in l_range(x)))

	iterated_length(mr::LM) = iterated_length_template(mr,
								x->sum(length(l_range(x,m)) for m in m_range(x)))

	@testset "ML" begin
		for l_min=0:l_cutoff, l_max=l_min:l_cutoff

			for T in [FullRange, ZeroTo, ToZero]
				mr = ML(l_min:l_max, T)
				testlength(iterated_length, mr)

                mr = ML(l_min:l_max, T(0))
                testlength(iterated_length, mr)
                mr2 = ML(UnitRange(l_range(mr)), UnitRange(m_range(mr)))
                @test begin
                    res = length(mr) == length(mr2)
                    if !res
                        @show typeof(mr)
                    end
                    res
                end
                @test collect(mr) == collect(mr2)

                mr = ML(l_min:l_max, T(l_min))
                testlength(iterated_length, mr)
                mr2 = ML(UnitRange(l_range(mr)), UnitRange(m_range(mr)))
                @test length(mr) == length(mr2)
                @test collect(mr) == collect(mr2)

				mr = ML(ZeroTo(l_max), T)
				testlength(iterated_length, mr)

				mr = ML(SingleValuedRange(l_min), T)
				testlength(iterated_length, mr)
			end

			mr = ML(l_min:l_min, SingleValuedRange(l_min))
			testlength(iterated_length, mr)

			mr = ML(ZeroTo(l_min), SingleValuedRange(l_min))
			testlength(iterated_length, mr)

			mr = ML(SingleValuedRange(l_min), SingleValuedRange(l_min))
			testlength(iterated_length, mr)

			for m_min=-l_max:l_max, m_max=m_min:l_max

				l_min_trimmed = max(l_min,max(m_min,-m_max))
				mr = ML(l_min_trimmed:l_max, m_min:m_max)

				testlength(iterated_length, mr)
			end
		end
	end

	@testset "LM" begin
		for l_min=0:l_cutoff, l_max=l_min:l_cutoff

			for T in [FullRange, ZeroTo, ToZero]
				mr = LM(l_min:l_max, T)
				testlength(iterated_length, mr)

                mr = LM(l_min:l_max, T(0))
                testlength(iterated_length, mr)
                mr2 = LM(UnitRange(l_range(mr)), UnitRange(m_range(mr)))
                @test begin
                    res = length(mr) == length(mr2)
                    if !res
                        @show typeof(mr)
                    end
                    res
                end
                @test collect(mr) == collect(mr2)

                mr = LM(l_min:l_max, T(l_min))
                testlength(iterated_length, mr)
                mr2 = LM(UnitRange(l_range(mr)), UnitRange(m_range(mr)))
                @test length(mr) == length(mr2)
                @test collect(mr) == collect(mr2)

				mr = LM(ZeroTo(l_max), T)
				testlength(iterated_length, mr)

				mr = LM(SingleValuedRange(l_min), T)
				testlength(iterated_length, mr)
			end

			mr = LM(l_min:l_min, SingleValuedRange(l_min))
			testlength(iterated_length, mr)

			mr = LM(ZeroTo(l_min), SingleValuedRange(l_min))
			testlength(iterated_length, mr)

			mr = LM(SingleValuedRange(l_min), SingleValuedRange(l_min))
			testlength(iterated_length, mr)

			for m_min=-l_max:l_max, m_max=m_min:l_max

				l_min_trimmed = max(l_min,max(m_min,-m_max))
				mr = LM(l_min_trimmed:l_max, m_min:m_max)

				testlength(iterated_length, mr)
			end
		end
	end

	@testset "LM == ML" begin
	    for l_min=0:l_cutoff, l_max=l_min:l_cutoff

	    	for T in [FullRange, ZeroTo, ToZero]
		    	m1 = LM(l_min:l_max, T)
				m2 = ML(l_min:l_max, T)
				@test length(m1) == length(m2)

		    	m1 = LM(ZeroTo(l_max), T)
				m2 = ML(ZeroTo(l_max), T)
				@test length(m1) == length(m2)

				m1 = LM(SingleValuedRange(l_min), T)
				m2 = ML(SingleValuedRange(l_min), T)
				@test length(m1) == length(m2)
			end

			m1 = LM(l_min:l_min, SingleValuedRange(l_min))
			m2 = ML(l_min:l_min, SingleValuedRange(l_min))
			@test length(m1) == length(m2)

			m1 = LM(ZeroTo(l_min), SingleValuedRange(l_min))
			m2 = ML(ZeroTo(l_min), SingleValuedRange(l_min))
			@test length(m1) == length(m2)

			m1 = LM(SingleValuedRange(l_min), SingleValuedRange(l_min))
			m2 = ML(SingleValuedRange(l_min), SingleValuedRange(l_min))
			@test length(m1) == length(m2)

	    	for m_min=-l_max:l_max,m_max=m_min:l_max

		    	l_min_trimmed = max(l_min,max(m_min,-m_max))

				m1 = LM(l_min_trimmed:l_max, m_min:m_max)
				m2 = ML(l_min_trimmed:l_max, m_min:m_max)

				@test length(m1) == length(m2)
			end
		end
	end

	iterated_length(mr::SphericalHarmonicModes.AbstractTriangleIterator) =
		iterated_length_template(mr, x->sum(length(l2_range(x,l1)) for l1 in l1_range(x)))

	@testset "L2L1Triangle" begin
		@testset "default l2minmax" begin
			for l_min=0:l_cutoff,l_max=l_min:l_cutoff,Δl_max=0:l_cutoff
				mr = L2L1Triangle(l_min,l_max,Δl_max)
				@test length(mr) == iterated_length(mr)
			end
		end

		@testset "all" begin
			for l_min=0:l_cutoff,l_max=l_min:l_cutoff,Δl_max=0:l_cutoff,
				l2_min=max(l_min-Δl_max,0):l_max+Δl_max,l2_max=l2_min:l_max+Δl_max

				l_min_trimmed = max(l_min,l2_min - Δl_max)
				l_max_trimmed = min(l_max,l2_max + Δl_max)
				mr = L2L1Triangle(l_min_trimmed,l_max_trimmed,Δl_max,l2_min,l2_max)
				@test begin
					res = length(mr) == iterated_length(mr)
					if !res
						println(mr)
					end
					res
				end
			end
		end
	end

	@testset "L1L2Triangle" begin
		@testset "default l2minmax" begin
			for l_min=0:l_cutoff,l_max=l_min:l_cutoff,Δl_max=0:l_cutoff
				mr = L1L2Triangle(l_min, l_max, Δl_max)
				@test length(mr) == iterated_length(mr)
			end
		end

		@testset "all" begin
			for l_min=0:l_cutoff,l_max=l_min:l_cutoff,Δl_max=0:l_cutoff,
				l2_min=max(l_min-Δl_max,0):l_max+Δl_max,l2_max=l2_min:l_max+Δl_max

				l_min_trimmed = max(l_min,l2_min - Δl_max)
				l_max_trimmed = min(l_max,l2_max + Δl_max)
				mr = L1L2Triangle(l_min_trimmed,l_max_trimmed,Δl_max,l2_min,l2_max)
				@test begin
					res = length(mr) == iterated_length(mr)
					if !res
						println(mr)
					end
					res
				end
			end
		end
	end
end

@testset "iterate" begin
    @testset "LM" begin
        m = LM(1:2)
        mode,state = iterate(m)
        @test mode == (2,-2)
        mode,state = iterate(m, state)
        @test mode == (1,-1)

        @test iterate(m, (last(m),length(m) + 1)) === nothing
    end
    @testset "ML" begin
        m = ML(1:2)
        mode,state = iterate(m)
        @test mode == (1,-1)
        mode,state = iterate(m, state)
        mode,state = iterate(m, state)
        mode,state = iterate(m, state)
        @test mode == (2,-2)

        @test iterate(m, (last(m),length(m) + 1)) === nothing
    end
    @testset "L2L1Triangle" begin
        m = L2L1Triangle(1:2, 1)
        mode,state = iterate(m)
        @test mode == (0,1)
        mode,state = iterate(m, state)
        mode,state = iterate(m, state)
        mode,state = iterate(m, state)
        @test mode == (1,2)

        @test iterate(m, (last(m),length(m) + 1)) === nothing
    end
	@testset "L1L2Triangle" begin
		ml2l1 = L2L1Triangle(1:20, 20)
		ml1l2 = L1L2Triangle(1:20, 20)
		for (v1, v2) in zip(ml1l2, ml2l1)
			@test v1 == reverse(v2)
			@test v1 in ml1l2
			@test v2 in ml2l1
		end
	end
end

@testset "LM ML modes" begin
	l_cutoff = 5
	for l_min in 0:l_cutoff, l_max=l_min:l_cutoff

		for T in [FullRange, ZeroTo, ToZero]
			m1 = LM(l_min:l_max, T)
			m2 = ML(l_min:l_max, T)

			@test sort(collect(m1)) == sort(collect(m2))

			m1 = LM(ZeroTo(l_max), T)
			m2 = ML(ZeroTo(l_max), T)

			@test sort(collect(m1)) == sort(collect(m2))
		end

		for m_min=-l_max:l_max, m_max=m_min:l_max

			l_min_trimmed = max(l_min,max(m_min,-m_max))
			l_range = l_min_trimmed:l_max
			m_range = m_min:m_max

			m1 = LM(l_range,m_range)
			m2 = ML(l_range,m_range)

			@test sort(collect(m1)) == sort(collect(m2))
		end
	end
end

@testset "modeindex" begin

	modeindex2(m, s::Integer, t::Integer) = modeindex2(m, (s,t))
    modeindex2(m, (s,t)::Tuple) = findfirst(isequal((s,t)), m)

	function testmodeindex(mr)
		for (ind,(s,t)) in enumerate(mr)
			@test begin
                res = modeindex(mr,(s,t)) == modeindex(mr,s,t) == modeindex2(mr,s,t) == ind
                if !res
                    @show mr, typeof(mr), s, t
                end
                res
            end
		end
	end

	l_cutoff = 5
	@testset "LM" begin
		for l_min=0:l_cutoff, l_max=l_min:l_cutoff

			for MT in [ZeroTo, ToZero, FullRange]
				mr = LM(l_min:l_max, MT)
				testmodeindex(mr)

				mr = LM(ZeroTo(l_max), MT)
				testmodeindex(mr)

                mr = LM(ZeroTo(l_max), MT(l_max))
                testmodeindex(mr)

				mr = LM(SingleValuedRange(l_min), MT)
				testmodeindex(mr)
			end

            mr = LM(1:4, ZeroTo(4))
            @test modeindex(mr, (3,3)) == modeindex2(mr, (3,3))

			mr = LM(l_min:l_max, SingleValuedRange(l_min))
			testmodeindex(mr)

			mr = LM(ZeroTo(l_max), SingleValuedRange(l_min))
			testmodeindex(mr)

			mr = LM(SingleValuedRange(l_min), SingleValuedRange(l_min))
			testmodeindex(mr)

			for m_min=-l_max:l_max, m_max=m_min:l_max

				l_min_trimmed = max(l_min,max(m_min,-m_max))

				mr = LM(l_min_trimmed:l_max, m_min:m_max)
				testmodeindex(mr)
			end
		end

		@test_throws ModeMissingError modeindex(LM(1:1), (2,2))
	end

	@testset "ML" begin
		for l_min=0:l_cutoff, l_max=l_min:l_cutoff

			for MT in [ZeroTo, ToZero, FullRange]
				mr = ML(l_min:l_max, MT)
				testmodeindex(mr)

				mr = ML(ZeroTo(l_max), MT)
				testmodeindex(mr)

                mr = ML(ZeroTo(l_max), MT(l_max))
                testmodeindex(mr)

				mr = ML(SingleValuedRange(l_min), MT)
				testmodeindex(mr)
			end

			mr = ML(l_min:l_max, SingleValuedRange(l_min))
			testmodeindex(mr)

			mr = ML(ZeroTo(l_max), SingleValuedRange(l_min))
			testmodeindex(mr)

			mr = ML(SingleValuedRange(l_min), SingleValuedRange(l_min))
			testmodeindex(mr)

			for m_min=-l_max:l_max, m_max=m_min:l_max

				l_min_trimmed = max(l_min,max(m_min,-m_max))

				mr = ML(l_min_trimmed:l_max, m_min:m_max)
				testmodeindex(mr)
			end
		end

		@test_throws ModeMissingError modeindex(ML(1:1), (2,2))
	end

    @testset "vector indexing" begin
        mr = ML(0:4, ZeroTo)
        (l,m) = first(mr)
        mr2 = ML(l,m)
        @test modeindex(mr, mr2) == modeindex(mr, l, m) == 1
        @test modeindex(flip(mr), mr2) == modeindex(flip(mr), l, m)

        for MT in [ZeroTo, FullRange, ToZero]
            mr = ML(0:4, MT)
            ms = ML(0:3, MT)
            @test modeindex(mr, ms) == [modeindex(mr, i) for i in ms]
            ms = ML(0:3, MT(3))
            @test modeindex(mr, ms) == [modeindex(mr, i) for i in ms]
        end

        mr = LM(0:4, ZeroTo)
        (l, m) = first(mr)
        mr2 = LM(l,m)
        @test modeindex(mr, mr2) == modeindex(mr, l, m) == 1
        @test modeindex(mr, LM(l:l,m)) == 1:1
    end

	@testset "L2L1Triangle" begin
		for l_min=0:l_cutoff,l_max=l_min:l_cutoff,
			Δl_max=0:l_cutoff,
			s′min=max(0,l_min-Δl_max):l_max+Δl_max,s′max=s′min:l_max+Δl_max

			l_max_trimmed = min(l_max,s′max + Δl_max)
			l_min_trimmed = max(l_min,s′min - Δl_max)

			m3 = L2L1Triangle(l_min_trimmed,l_max_trimmed,Δl_max,s′min,s′max)
			for (ind,(s′,s)) in enumerate(m3)
				@test begin
					res = modeindex(m3,(s′,s)) == modeindex(m3,s′,s) == modeindex2(m3,s′,s) == ind
					if !res
						println(m3)
					end
					res
				end
			end
		end
	end

	@testset "L1L2Triangle" begin
		for l_min=0:l_cutoff,l_max=l_min:l_cutoff,
			Δl_max=0:l_cutoff,
			s′min=max(0,l_min-Δl_max):l_max+Δl_max,s′max=s′min:l_max+Δl_max

			l_max_trimmed = min(l_max,s′max + Δl_max)
			l_min_trimmed = max(l_min,s′min - Δl_max)

			m3 = L1L2Triangle(l_min_trimmed,l_max_trimmed,Δl_max,s′min,s′max)
			for (ind,(s′,s)) in enumerate(m3)
				@test begin
					res = modeindex(m3,(s′,s)) == modeindex(m3,s′,s) == modeindex2(m3,s′,s) == ind
					if !res
						println(m3)
					end
					res
				end
			end
		end
	end

    @testset "other iterators" begin
        r = collect((i,i) for i in 1:1000)
        for (ind, t) in enumerate(r)
            @test modeindex(r, t) == ind
            @test modeindex(r, t...) == ind
        end
        @test modeindex(r, (0,0)) === nothing
    end
end

@testset "firstlast" begin
	@testset "first" begin
		@testset "ML" begin
			m1 = ML(rand(1:5):rand(6:10))
			@test first(collect(m1)) == first(m1)
		end

		@testset "LM" begin
			m2 = LM(rand(1:5):rand(6:10))
			@test first(collect(m2)) == first(m2)
		end

		@testset "L2L1Triangle" begin
			m3 = L2L1Triangle(rand(1:3):rand(4:10),rand(1:5))
			@test first(collect(m3)) == first(m3)
		end

		@testset "L1L2Triangle" begin
			m3 = L1L2Triangle(rand(1:3):rand(4:10),rand(1:5))
			@test first(collect(m3)) == first(m3)
		end
	end
	@testset "last" begin
		@testset "ML" begin
			m1 = ML(rand(1:5):rand(6:10))
			@test last(collect(m1)) == last(m1)
		end

		@testset "LM" begin
			m2 = LM(rand(1:5):rand(6:10))
			@test last(collect(m2)) == last(m2)
		end

		@testset "L2L1Triangle" begin
			m3 = L2L1Triangle(rand(1:3):rand(4:10),rand(1:5))
			@test last(collect(m3)) == last(m3)
		end

		@testset "L1L2Triangle" begin
			m3 = L1L2Triangle(rand(1:3):rand(4:10),rand(1:5))
			@test last(collect(m3)) == last(m3)
		end
	end
end

@testset "firstindex & lastindex" begin
    @testset "firstindex" begin
        mr = LM(1:2)
        @test firstindex(mr) == 1
        mr = ML(1:2)
        @test firstindex(mr) == 1
        mr = L2L1Triangle(1:3,1)
        @test firstindex(mr) == 1
    end
    @testset "lastindex" begin
        mr = LM(1:2)
        @test lastindex(mr) == 8
        mr = ML(1:2)
        @test lastindex(mr) == 8
        mr = L2L1Triangle(1:3,1)
        @test lastindex(mr) == 9
    end
end

@testset "size & axes" begin
    @testset "size" begin
        mr = LM(1:2)
        @test size(mr) == (8,)
        @test size(mr,1) == 8
        @test size(mr,2) == 1
        mr = ML(1:2)
        @test size(mr) == (8,)
        @test size(mr,1) == 8
        @test size(mr,2) == 1
        mr = L2L1Triangle(1:3,1)
        @test size(mr) == (9,)
        @test size(mr,1) == 9
        @test size(mr,2) == 1
    end
    @testset "axes" begin
        mr = LM(1:2)
        @test axes(mr) == (Base.OneTo(8),)
        @test axes(mr,1) == Base.OneTo(8)
        @test axes(mr,2) == Base.OneTo(1)
        mr = ML(1:2)
        @test axes(mr) == (Base.OneTo(8),)
        @test axes(mr,1) == Base.OneTo(8)
        @test axes(mr,2) == Base.OneTo(1)
        mr = L2L1Triangle(1:3,1)
        @test axes(mr) == (Base.OneTo(9),)
        @test axes(mr,1) == Base.OneTo(9)
        @test axes(mr,2) == Base.OneTo(1)
    end
end

@testset "flip" begin
    @testset "LM" begin
		m = LM(0:1)
		@test ML(m) == ML(0:1)
        @test convert(ML, m) == ML(0:1)
		@test flip(m) == ML(0:1)
		@test LM(m) === m
        @test convert(LM, m) === m
        @test flip(SphericalHarmonicModes.basetype(LM(1:2))) == ML
    end
    @testset "ML" begin
		m = ML(0:1)
		@test LM(m) == LM(0:1)
        @test convert(LM, m) == LM(0:1)
		@test flip(m) == LM(0:1)
		@test ML(m) === m
        @test convert(ML, m) === m
        @test flip(SphericalHarmonicModes.basetype(ML(1:2))) == LM
    end
end

@testset "ordering" begin
    for l_range in Any[0:1, ZeroTo(1)], m_range in Any[0:1, FullRange, ZeroTo]
        lm = LM(0:1)
        ml = SphericalHarmonicModes.flip(lm)
        @test SphericalHarmonicModes.ofordering(lm, ml) == lm
        @test SphericalHarmonicModes.ofordering(lm, lm) == lm
        @test SphericalHarmonicModes.ofordering(ml, lm) == ml
        @test SphericalHarmonicModes.ofordering(ml, lm) == ml
    end

	m21 = L2L1Triangle(1:3, 3)
	m12 = L1L2Triangle(1:3, 3)
	@test SphericalHarmonicModes.ofordering(m21, m21) === m21
	@test SphericalHarmonicModes.ofordering(m21, m12) === m21
	@test SphericalHarmonicModes.ofordering(m12, m21) === m12
	@test SphericalHarmonicModes.ofordering(m12, m12) === m12
end

@testset "intersect" begin
	@testset "same l, different m" begin
		for T in [LM, ML]
			@eval begin
			    l1 = $T(3:3,0:2)
			    l2 = $T(3:3,1:3)
			    @test intersect(l1,l2) == $T(3:3,1:2)

			    l1 = $T(1:3)
			    l2 = $T(1:3,1:3)
			    @test intersect(l1,l2) == $T(1:3,1:3)

			    for MT in [ZeroTo, ToZero]
				    local l1 = $T(1:3, MT)
				    local l2 = $T(1:3, MT)
				    @test intersect(l1,l2) == intersect(l2,l1) == $T(1:3, MT)

				    local l1 = $T(1:3, MT)
				    local l2 = $T(1:3, FullRange)
				    @test intersect(l1,l2) == intersect(l2,l1) == $T(1:3, MT)
				end

				l1 = $T(1:3, ZeroTo)
			    l2 = $T(1:3, ToZero)
			    @test intersect(l1,l2) == $T(1:3, ZeroTo(0))
			    @test intersect(l2,l1) == $T(1:3, ZeroTo(0))
			end
		end
	end
	@testset "different l, same m" begin
		for T in [LM, ML]
			@eval begin
			    l1 = $T(1:3,1:1)
			    l2 = $T(1:1,1:1)
			    @test intersect(l1,l2) == $T(1:1,1:1)

			    l1 = $T(2:2)
			    l2 = $T(1:3,-2:2)
			    @test intersect(l1,l2) == $T(2:2,-2:2)

			    l1 = $T(ZeroTo(2))
			    l2 = $T(1:3,-2:2)
			    @test intersect(l1,l2) == $T(1:2,-2:2)

			    l1 = $T(ZeroTo(2))
			    l2 = $T(ZeroTo(3),-2:2)
			    @test intersect(l1,l2) == $T(ZeroTo(2),-2:2)
			end
		end
	end
	@testset "different l, different m" begin
		for T in [LM,ML]
			@eval begin
			    l1 = $T(2:3,2:2)
			    l2 = $T(0:2)
			    @test intersect(l1,l2) == $T(2:2,2:2)

			    l1 = $T(2:3)
			    l2 = $T(1:10,-1:10)
			    @test intersect(l1,l2) == $T(2:3,-1:3)

			    l1 = $T(2:3)
			    l2 = $T(9:10, 5:5)
			    @test intersect(l1,l2) === nothing

			    for MT in [FullRange, ZeroTo, ToZero]
				    local l1 = $T(2:3,MT)
				    local l2 = $T(1:4,MT)
			    	@test intersect(l1,l2) == $T(2:3,MT)

			    	local l1 = $T(ZeroTo(3),MT)
				    local l2 = $T(ZeroTo(4),MT)
			    	@test intersect(l1,l2) == $T(ZeroTo(3),MT)
			    end
			end
		end
	end

    @test intersect(FullRange(2), FullRange(1)) === FullRange(1)

    @test intersect(SingleValuedRange(1), -1:1) === SingleValuedRange(1)
    @test intersect(SingleValuedRange(1), FullRange(1)) === SingleValuedRange(1)
    @test intersect(SingleValuedRange(1), ZeroTo(1)) === SingleValuedRange(1)
    @test intersect(SingleValuedRange(1), 2:3) === nothing
    @test intersect(-1:1, SingleValuedRange(1)) === SingleValuedRange(1)
    @test intersect(FullRange(1), SingleValuedRange(1)) === SingleValuedRange(1)
    @test intersect(ZeroTo(1), SingleValuedRange(1)) === SingleValuedRange(1)
    @test intersect(2:3, SingleValuedRange(1)) === nothing
    @test intersect(SingleValuedRange(1), SingleValuedRange(1)) === SingleValuedRange(1)
    @test intersect(SingleValuedRange(1), SingleValuedRange(2)) === nothing
end

@testset "show" begin
    io = IOBuffer()

    function testshow(io, x)
    	show(io, x)
    	show(io, MIME"text/plain"(), x)
    end

    testshow(io, SingleValuedRange(2))
    testshow(io, ZeroTo(2))
    testshow(io, ToZero(2))
    testshow(io, FullRange(2))

    testshow(io, LM(1:2,1:1))
    testshow(io, ML(1:2,1:1))
    testshow(io, L2L1Triangle(1:2,1))
    testshow(io, L1L2Triangle(1:2,1))

    showerror(io, ModeMissingError(LM(1:1), (2,2)) )
    showerror(io, ModeMissingError(L2L1Triangle(1:1,1), (3,3)) )
end

@testset "to_indices" begin
    modes = ML(0:3)
    @test to_indices(zeros(), (modes,), (:,)) == (1:length(modes),)
end