text:
# Sets up reverse dictionary
R_Dict = Dict{Symbol,Symbol}()
R_Dict[:+] = :plus_rev
R_Dict[:-] = :minus_rev
R_Dict[:*] = :mul_revDR
R_Dict[:/] = :div_rev
R_Dict[:^] = :power_rev
R_Dict[:pow] = :power_rev
R_Dict[:exp] = :exp_revDR
for f in (  :sqrt, :sqr, :abs,
            :sin, :cos, :tan,
            :asin, :log,
            :acos, :atan, :sinh,
            :cosh, :tanh, :asinh,
            :acosh, :atanh, :min,
            :max, :inv)
    R_Dict[f] = Symbol(f, "_rev")
end

"""
    Tape

A storage object to record expressions in a tape. The fields are:
* `Edge_List::Vector`: Contains a list of edges corresponding the directed graph
* `Head_List::Vector{Symbol}`: Contains a list of the expressions used
* `Last_Node::Int64`: The last node in directed graph
* `Intv_Storage::Vector{Interval}`: A tape of interval values corresponding to each node.
* `FW_Arg::Vector{Array{Int64}}`: Vector point to argument node of expressions
* `FW_Expr::Vector`: List of expressions in forward propagation
* `nvar::Float64`: Number of variables in expression
* `RV_Arg::Vector{Vector{Int64}}`: Arguement list for reverse propagation
* `RV_Expr::Vector`: List of expressions in reverse propagation
* `gL::Float64`: Lower bound on the expression gL<g(x)
* `gU::Float64`: Upper bound on the expression gU>g(x)
* `cnsts::Vector`: List of constant values and their location
"""
immutable Tape
  Edge_List::Vector
  Head_List::Vector{Symbol}
  Last_Node::Int64
  Intv_Storage::Vector{Interval}
  FW_Arg::Vector{Array{Int64}}
  FW_Expr::Vector
  nvar::Float64
  RV_Arg::Vector{Vector{Int64}}
  RV_Expr::Vector
  gL::Float64
  gU::Float64
  cnsts::Vector
  fix_val
end

"""
    Tape()

Initializes an empty tape for storage.
"""
Tape() =  ([],
           Symbol[],
           0,
           Interval[],
           [[1]],
           [],
           0.0,
           [[1]],
           [],
           0.0,
           0.0,
           [],
           [])

"""
    TapeList

A storage object for a list of tapes. Has a single field sto which is an array
of Tape objects.
"""
immutable TapeList
  sto::Any
end

"""
    TapeList

Initializes an empty tapelist object
"""
TapeList() = TapeList([Tape()])

"""
    Generate_Tape(exp::Expr,nx::Int64,gL,gU)

Generates the tape of the provided expression `exp::Expr` assuming the expression
depends on `nx::Int64` variables and has lower bounds `gL` and upper bounds `gU`.
The variables in the expression must be off the form `x[1],...,x[nx]`.
"""
function Generate_Tape(exp::Expr,nx::Int64,gL,gU)

  gL,gU = Float64(gL),Float64(gU)

  # Sets up node finding object and generates graph to global
  X_NodeF::Vector{NodeFinder} = NodeFinder[NodeFinder(i) for i=1:nx]
  global NodeCounter = deepcopy(nx)
  fX! = @eval x::Vector{NodeFinder} -> $exp
  Base.invokelatest(fX!,X_NodeF)

  # generates forward tape
  FW_Arg = [[EdgeList[i][j][1] for j=1:length(EdgeList[i])] for i=1:length(EdgeList)]
  RW_Arg = [vcat([EdgeList[i][1][2]],[EdgeList[i][j][1] for j=1:length(EdgeList[i])]) for i=1:length(EdgeList)]
  tape = Tape(deepcopy(EdgeList),deepcopy(HeaderList),
                      NodeCounter,Interval[Interval(-Inf,Inf) for i=1:NodeCounter],
                      FW_Arg,[Expr(:call) for i=1:NodeCounter],nx,
                      RW_Arg,[Expr(:call) for i=1:NodeCounter],gL,gU,
                      deepcopy(ConstList),[])
  global EdgeList = []
  global HeaderList = Symbol[]
  global ConstList = []
  return tape
end

"""
    Generate_Fixed_Tape(exp::Expr,nx::Int64,gL,gU,vals)

Generates the tape of the provided expression `exp::Expr` assuming the expression
depends on `nx::Int64` variables and has lower bounds `gL` and upper bounds `gU`.
The variables in the expression must be off the form `x[1],...,x[nx]`. Variables `x[nx+1]`
to `x[end]` are fixed to the values in the `vals` array.
"""
function Generate_Fixed_Tape(exp::Expr,nx::Int64,gL,gU,vals)

  gL,gU = Float64(gL),Float64(gU)

  # Sets up node finding object and generates graph to global
  X_NodeF::Vector{Any} = vcat(Any[NodeFinder(i) for i=1:nx],vals)
  global NodeCounter = deepcopy(nx)
  fX! = @eval x::Vector{Any} -> $exp
  Base.invokelatest(fX!,X_NodeF)

  # generates forward tape
  FW_Arg = [[EdgeList[i][j][1] for j=1:length(EdgeList[i])] for i=1:length(EdgeList)]
  RW_Arg = [vcat([EdgeList[i][1][2]],[EdgeList[i][j][1] for j=1:length(EdgeList[i])]) for i=1:length(EdgeList)]
  tape = Tape(deepcopy(EdgeList),deepcopy(HeaderList),
                      NodeCounter,Interval[Interval(-Inf,Inf) for i=1:NodeCounter],
                      FW_Arg,[Expr(:call) for i=1:NodeCounter],nx,
                      RW_Arg,[Expr(:call) for i=1:NodeCounter],gL,gU,
                      deepcopy(ConstList),vals)
  global EdgeList = []
  global HeaderList = Symbol[]
  global ConstList = []
  return tape
end

"""
    Generate_TapeList(exprs::Vector{Expr},nx::Int64,gL::Vector{Float64},gU::Vector{Float64})

Generates the tape list for each provided expression `exprs[i]` in
`exprs::Vector{Expr}` assuming all expressions depend on `nx::Int64` variables
and has lower bounds `gL[i]` and upper bounds `gU[i]`. The variables in the
expression must be off the form `x[1],...,x[nx]`.
"""
function Generate_TapeList(exprs::Vector{Expr},nx::Int64,gL::Vector{Float64},gU::Vector{Float64})
  @assert length(exprs) == length(gL) == length(gU)
  tapelist = []
  for i=1:length(exprs)
    push!(tapelist,Generate_Tape(exprs[i],nx,gL[i],gU[i]))
  end
  return TapeList(tapelist)
end

"""
    Generate_Fixed_TapeList(exprs::Vector{Expr},nx::Int64,gL::Vector{Float64},gU::Vector{Float64},val_arr)

Generates the tape list for each provided expression `exprs[i]` in
`exprs::Vector{Expr}` assuming all expressions depend on `nx::Int64` variables
and has lower bounds `gL[i]` and upper bounds `gU[i]`. The variables in the
expression must be off the form `x[1],...,x[nx]`. Variables `x[nx+1]`
to `x[end]` are fixed to the values in the `vals` array.
"""
function Generate_Fixed_TapeList(exprs::Vector{Expr},nx::Int64,gL::Vector{Float64},
                                 gU::Vector{Float64},val_arr)
  @assert length(exprs) == length(gL) == length(gU)
  tapelist = []
  for i=1:length(exprs)
      for j=1:length(val_arr)
          push!(tapelist,Generate_Fixed_Tape(exprs[i],nx,gL[i],gU[i],val_arr[j]))
      end
  end
  return TapeList(tapelist)
end

"""
    SetConstraintNode!(x::Tape)

Sets the terminal interval to `[x.gL,x.gU]` using the bounds in the `x::Tape`.
"""
function SetConstraintNode!(x::Tape)
  x.Intv_Storage[x.Last_Node] = Interval(x.gL,x.gU)
end
"""
    SetConstraintNode!(x::TapeList)

Sets the terminal interval to `[x.gL[i],x.gU[i]]`  in each tape using the bounds
in the `x::TapeList[i]`.
"""
function SetConstraintNode!(x::TapeList)
  for i=1:length(x.sto)
    SetConstraintNode!(x.sto[i])
  end
end

"""
    SetConstantNode!(x::Tape)

Sets nodes in constaint value list of tape according to their value.
"""
function SetConstantNode!(x::Tape)
  for i=1:length(x.cnsts)
    x.Intv_Storage[Int64(x.cnsts[i][1])] = Interval(x.cnsts[i][2])
  end
end
"""
    SetConstantNode!(x::TapeList)

Sets nodes in constaint value list of tape according to their value for
each tape in the tape list.
"""
function SetConstantNode!(x::TapeList)
  for i=1:length(x.sto)
    SetConstantNode!(x.sto[i])
  end
end

"""
    SetVarBounds!(x::Tape,X::Vector{Interval})

Sets variable nodes to the interval value in `X::Vector{Interval}` in tape.
"""
function SetVarBounds!(x::Tape,X::Vector{Interval})
  x.Intv_Storage[1:length(X)] = X
end

"""
    SetVarBounds!(x::TapeList,X::Vector{Interval})

Sets variable nodes to the interval value in `X::Vector{Interval}` for all
elements of the TapeList.
"""
function SetVarBounds!(x::TapeList,X::Vector{Interval})
  for i=1:length(x)
    SetVarBounds!(x.sto[i],X)
  end
end

"""
    GetVarBounds(x::Tape)

Unpacks the variable node values and returns an array of intervals.
"""
function GetVarBounds(x::Tape)
  return x.Intv_Storage[1:x.Var_Num]
end

"""
    ForwardPass!(x::Tape)

Performs a forward contractor pass using tape `x::Tape`.
"""
function ForwardPass!(x::Tape)
  for i=1:length(x.Edge_List)
    x.FW_Expr[i].args = vcat([x.Head_List[i]],x.Intv_Storage[x.FW_Arg[i]])
    x.Intv_Storage[x.Edge_List[i][1][2]] = eval(x.FW_Expr[i]) ∩ x.Intv_Storage[x.Edge_List[i][1][2]]
  end
end

"""
    ForwardPassOne!(x::Tape)

Performs a forward contractor pass using tape `x::Tape` with some uninitialized nodes.
"""
function ForwardPassOne!(x::Tape)
  for i=1:length(x.Edge_List)
    x.FW_Expr[i].args = vcat([x.Head_List[i]],x.Intv_Storage[x.FW_Arg[i]])
    x.Intv_Storage[x.Edge_List[i][1][2]] = eval(x.FW_Expr[i])
  end
end

"""
    ReversePass!(x::Tape)

Performs a reverse contractor pass using tape `x::Tape`.
"""
function ReversePass!(x::Tape)
  for i=length(x.Edge_List):-1:1
    x.RV_Expr[i].args = vcat([R_Dict[x.Head_List[i]]],x.Intv_Storage[x.RV_Arg[i]])
    x.Intv_Storage[x.RV_Arg[i]] = [eval(x.RV_Expr[i])...]
  end
end

"""
    DAGContractor!(X::Vector{Interval{T}},x::TapeList,r)

Performs a forward-interval contactor propagation `r` times using the `x::TapeList`
and the initial interval bounds `X::Vector{Interval{T}}`.
"""
function DAGContractor!(X::Vector{Interval{T}},x::TapeList,r) where {T}
  Xprev::Vector{Interval} = copy(X) # sets variable bounds on first Array to Box Bounds
  SetConstraintNode!(x)
  for i=1:r
    for j=1:length(x.sto)
      SetVarBounds!(x.sto[j],Xprev) # take refined bounds from previous graph
      SetConstantNode!(x)
      if (i == 1)
        ForwardPassOne!(x.sto[j])
      else
        ForwardPass!(x.sto[j]) # run forward and reverse pass
      end
      ReversePass!(x.sto[j])
      Xprev = x.sto[j].Intv_Storage[1:length(X)]
    end
  end
  X[:] = Xprev
end

"""
    getDAG(expr::Expr,root_list::Vector{Int64})

Generates the directed graph of expression `expr` using the root_list corresponding
to the indices of `x[i]` that appear in the expression.
"""
function getDAG(expr::Expr,root_list::Vector{Int64})
  # generates tape
  tape = Generate_Tape(expr,maximum(root_list),0.0,0.0)

  # generates vertix list
  not_root = setdiff([i for i=1:Int64(tape.nvar)],root_list)
  vertix_list = filter!(x->~in(x,not_root),[i for i=1:(tape.Last_Node)])
  vertix_name = Array{Any}(tape.Last_Node)

  # generates edge list
  # to do: record edge number while taping and preallocate
  edge_list = []
  edge_temp = tape.Edge_List
  for i=1:length(edge_temp)
    for j=1:length(edge_temp[i])
      push!(edge_list,edge_temp[i][j])
      if (j == 1)
        vertix_name[edge_temp[i][j][2]] = tape.Head_List[i]
      end
    end
  end

  for i=1:length(tape.cnsts)
    vertix_name[Int64(tape.cnsts[i][1])] = tape.cnsts[i][2]
  end
  for i=1:Int64(tape.nvar)
    vertix_name[i] = :(x[$i])
  end
  deleteat!(vertix_name,not_root)

  return vertix_list,vertix_name,edge_list
end


text:
config = MOI.Test.TestConfig()
optimize!(mock) = MOIU.mock_optimize!(mock, [2.0, 1.0, 1.0, √2],
    (MOI.VectorOfVariables, MOI.RotatedSecondOrderCone) => [[1.0, 1.0, -√2]],
    (MOI.VectorAffineFunction{Float64}, MOI.Zeros) => [[1.0, -1.0, 1.0]])
for mock in mocks(optimize!)
    Tests.sos_univariate_quadratic_test(mock, config)
    Tests.sos_bivariate_quadratic_test(mock, config)
    Tests.sdsos_univariate_quadratic_test(mock, config)
    Tests.sdsos_bivariate_quadratic_test(mock, config)
end
optimize!(mock) = MOIU.mock_optimize!(mock, [2.0, 1.0, 1.0, 1.0, 1.0],
    (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) => [0.0, 2.0, 2.0, 0.0],
    (MOI.VectorAffineFunction{Float64}, MOI.Zeros) => [[2.0, -1.0, 0.0]])
for mock in mocks(optimize!)
    Tests.dsos_univariate_quadratic_test(mock, config)
    Tests.dsos_bivariate_quadratic_test(mock, config)
end

ext: jl


text:
export
    sigmoid,
    relu,
    softmax,
    onehot,
    crossentropyerror

function sigmoid(x::T) where T <: Real
    return 1.0 / (1.0 + exp(-x))
end

function relu(x::T) where T <: Real
    return max(zero(x), x)
end

function softmax(x::AbstractVector{T}) where T <: Real
    c = maximum(x)
    exp_a = exp.(x .- c)
    return exp_a ./ sum(exp_a)
end

function softmax(x::AbstractMatrix{T}) where T <: Real
    return mapslices(softmax, x, 1)
end

function onehot(::Type{T}, t::AbstractVector, l::AbstractVector) where T
    r = zeros(T, length(l), length(t))
    for i in 1:length(t)
        r[findfirst(l, t[i]), i] = 1
    end
    return r
end

onehot(t::AbstractVector, l::AbstractVector) = onehot(Int, t, l)

function crossentropyerror(y::Vector, t::Vector)
    δ = 1.0e-7
    # -sum(t .* log.(y .+ δ))
    return -(t ⋅ log.(y .+ δ))
end

function crossentropyerror(y::Matrix, t::Matrix)
    batch_size = size(y, 2)
    δ = 1.0e-7
    # -sum(t .* log.(y .+ δ)) ./ batch_size
    return -vecdot(t, log.(y .+ δ)) ./ batch_size
end


text:
########## Collections of some old functions (Not maintained and not exposed to main package) ##########

#### fBm-MLE ####

"""
    fBm_log_likelihood_H(X, H, ψ, G)

Log-likelihood of a fBm model with the optimal volatility.
"""
function fBm_log_likelihood_H(X::AbstractVecOrMat{<:Real}, H::Real, G::AbstractVector{<:Integer}=Int[])
    proc = FractionalBrownianMotion(H)

    Σ::AbstractMatrix = if length(G)>0
        @assert length(G) == size(X,1)
        @assert any(diff(G) .> 0)  # grid points are in increasing order
        covmat(proc, G)
    else
        covmat(proc, size(X,1))
    end

    return log_likelihood_H(Σ, X-X[1])  # relative to the first point
end

"""
Maximum likelihood estimation of Hurst exponent and volatility for fBm.

# Args
- X: sample vector or matrix. For matrix each column is a sample.
- ψ: wavelet filter used for computing `X`.
- G: integer time grid of `X`, by default the regular grid `1:size(X,1)` is used.
- method: `:optim` for optimization based or `:table` for lookup table based procedure
- ε: search hurst in the range [ε, 1-ε]

# Notes
- The MLE is known for its sensitivivity to mis-specification of model. In particular the fGn-MLE is sensitive to NaN value and outliers.
"""
function fBm_MLE_estim(X::AbstractVecOrMat{<:Real}, G::AbstractVector{<:Integer}=Int[]; method::Symbol=:optim, ε::Real=1e-2)
    # @assert 0. < ε < 1.
    if length(G)>0
        @assert length(G) == size(X,1)
        @assert minimum(abs.(diff(sort(G)))) > 0  # all elements are distinct
    end

    func = h -> -fBm_log_likelihood_H(X, h, G)

    opm = nothing
    hurst = nothing

    if method == :optim
        # Gradient-free constrained optimization
        opm = Optim.optimize(func, ε, 1-ε, Optim.Brent())
        # # Gradient-based optimization
        # optimizer = Optim.GradientDescent()  # e.g. Optim.BFGS(), Optim.GradientDescent()
        # opm = Optim.optimize(func, ε, 1-ε, [0.5], Optim.Fminbox(optimizer))
        hurst = Optim.minimizer(opm)[1]
    elseif method == :table
        Hs = collect(ε:ε:1-ε)
        hurst = Hs[argmin([func(h) for h in Hs])]
    else
        throw("Unknown method: ", method)
    end

    proc = FractionalBrownianMotion(hurst)
    Σ = covmat(proc, length(G)>0 ? G : size(X,1))
    σ = sqrt(xiAx(Σ, X) / length(X))
    L = log_likelihood_H(Σ, X)
    # # or equivalently
    # L = log_likelihood(σ^2*Σ, X)

    return hurst, σ, L, opm
end

"""
Accelerated fBm-MLE by dividing a large vector of samples into smaller ones.

The MLE method can be expensive on data of large dimensions due to the inversion of covariance matrix. This function accelerates the MLE method by dividing a large vector `X` into smaller vectors of size `s` downsampled by a factor `l`. The smaller vectors are treated by MLE as i.i.d. samples.

# Args
- X: same as in `fWn_MLE_estim()`
- s: sub window size
- l: length of decorrelation
"""
function fBm_MLE_estim(X::AbstractVector{<:Real}, s::Integer, l::Integer; kwargs...)
    t, V = rolling_vectorize(X, s, 1, l; mode=:causal)
    return fBm_MLE_estim(V; kwargs...)  # The regular grid is implicitely used here.
end

#### fGn-MLE ####
# A special case of fWn-MLE which deserves it own implementation.

function fGn_log_likelihood_H(X::AbstractVecOrMat{<:Real}, H::Real, d::Int, G::AbstractVector{<:Integer}=Int[])
    proc = FractionalGaussianNoise(H, d)

    Σ::AbstractMatrix = if length(G)>0
        @assert length(G) == size(X,1)
        covmat(proc, G)
    else
        covmat(proc, size(X,1))
    end

    return log_likelihood_H(Σ, X)
end

function fGn_MLE_estim(X::AbstractVecOrMat{<:Real}, d::Integer, G::AbstractVector{<:Integer}=Int[]; method::Symbol=:optim, ε::Real=1e-2)
    # @assert 0. < ε < 1.
    if length(G)>0
        @assert length(G) == size(X,1)
        # @assert minimum(abs.(diff(sort(G)))) > 0  # all elements are distinct
    end

    func = h -> -fGn_log_likelihood_H(X, h, d, G)

    opm = nothing
    hurst = nothing

    if method == :optim
        # Gradient-free constrained optimization
        opm = Optim.optimize(func, ε, 1-ε, Optim.Brent())
        # # Gradient-based optimization
        # optimizer = Optim.GradientDescent()  # e.g. Optim.BFGS(), Optim.GradientDescent()
        # opm = Optim.optimize(func, ε, 1-ε, [0.5], Optim.Fminbox(optimizer))
        hurst = Optim.minimizer(opm)[1]
    elseif method == :table
        Hs = collect(ε:ε:1-ε)
        hurst = Hs[argmin([func(h) for h in Hs])]
    else
        error("Unknown method: $method")
    end

    proc = FractionalGaussianNoise(hurst, d)
    Σ = covmat(proc, length(G)>0 ? G : size(X,1))
    σ = sqrt(xiAx(Σ, X) / length(X))
    L = log_likelihood_H(Σ, X)
    # or equivalently
    # L = log_likelihood(σ^2*Σ, X)

    return hurst, σ, L, opm
end

"""
Accelerated fGn-MLE by dividing a large vector of samples into smaller ones.

The MLE method can be expensive on data of large dimensions due to the inversion of covariance matrix. This function accelerates the MLE method by dividing a large vector `X` into smaller vectors of size `s` downsampled by a factor `l`. The smaller vectors are treated by MLE as i.i.d. samples.

# Args
- X, d: same as in `fGn_MLE_estim()`
- s,l: sub window size, length of decorrelation
"""
function fGn_MLE_estim(X::AbstractVector{<:Real}, d::Integer, s::Integer, l::Integer; kwargs...)
    t, V = rolling_vectorize(X, s, 1, l; mode=:causal)
    return fGn_MLE_estim(V, d; kwargs...)
end

# function bspline_scalogram_estim(X::AbstractMatrix{<:Real}, sclrng::AbstractVector{<:Integer}, v::Integer, pows::AbstractVector{<:Real}; kwargs...)

#     bspline_scalogram_estim()
# end

# """
# B-Spline scalogram estimator with a matrix of DCWT coefficients as input. Each column in `W` is a vector of DCWT coefficients.
# """
# function fBm_bspline_scalogram_estim(W::AbstractMatrix{T}, sclrng::AbstractVector{Int}, v::Int; dims::Int=1, mode::Symbol=:center) where {T<:Real}
#     return fBm_bspline_scalogram_estim(var(W,dims), sclrng, v; mode=mode)
# end

# """
# B-Spline scalogram estimator with an array of DCWT coefficients as input. Each row in `W` corresponds to a scale.
# """
# function fBm_bspline_scalogram_estim(W::AbstractVector{T}, sclrng::AbstractVector{Int}, v::Int; mode::Symbol=:center) where {T<:AbstractVector{<:Real}}
#     return fBm_bspline_scalogram_estim([var(w) for w in W], sclrng, v; mode=mode)
# end

"""
Generalized B-Spline scalogram estimator for Hurst exponent and volatility.

# Args
- Σ: covariance matrix of wavelet coefficients.
- sclrng: scale of wavelet transform. Each number in `sclrng` corresponds to one row in the matrix X
- v: vanishing moments
- r: rational ratio defining a line in the covariance matrix, e.g. r=1 corresponds to the main diagonal.
"""
function gen_bspline_scalogram_estim(Σ::AbstractMatrix{T}, sclrng::AbstractVector{Int}, v::Int, r::Rational=1//1; mode::Symbol=:center) where {T<:Real}
    @assert issymmetric(Σ)
    @assert size(Σ,1) == length(sclrng)
    @assert r >= 1
    if r > 1
        all(diff(sclrng/sclrng[1]) .== 1) || error("Imcompatible scales: the ratio between the k-th and the 1st scale must be k")
    end

    p,q,N = r.num, r.den, length(sclrng)
    @assert N>=2p

    # Σ = cov(X, X, dims=2, corrected=true)  # covariance matrix

    yr = [log(abs(Σ[q*j, p*j])) for j in 1:N if p*j<=N]
    xr = [log(sclrng[q*j] * sclrng[p*j]) for j in 1:N if p*j<=N]

    df = DataFrames.DataFrame(xvar=xr, yvar=yr)
    ols = GLM.lm(@GLM.formula(yvar~xvar), df)
    coef = GLM.coef(ols)

    hurst = coef[2]-1/2
    Aρ = Aρ_bspline(0, r, hurst, v, mode)
    σ = exp((coef[1] - log(abs(Aρ)))/2)
    return (hurst, σ), ols

    # Ar = hcat(xr, ones(length(xr)))  # design matrix
    # H0, η = Ar \ yr  # estimation of H and β
    # hurst = H0-1/2
    # Aρ = Aρ_bspline(0, r, hurst, v, mode)
    # σ = ℯ^((η - log(abs(Aρ)))/2)
    # return hurst, σ
end
# const fBm_gen_bspline_scalogram_estim = gen_bspline_scalogram_estim

"""
Compute the covariance matrix of a fWn at some time lag.

# Args
- F: array of band pass filters (no DC component)
- d: time lag
- H: Hurst exponent
"""
function fWn_covmat_lag(F::AbstractVector{<:AbstractVector{T}}, d::DiscreteTime, H::Real) where {T<:Real}
    L = maximum([length(f) for f in F])  # maximum length of filters
    # M = [abs(d+(n-m))^(2H) for n=0:L-1, m=0:L-1]  # matrix comprehension is ~ 10x slower
    M = zeros(L,L)
    for n=1:L, m=1:L
        M[n,m] = abs(d+(n-m))^(2H)
    end
    Σ = -1/2 * [f' * view(M, 1:length(f), 1:length(g)) * g for f in F, g in F]
end

"""
Compute the covariance matrix of a time-concatenated fWn.
"""
function fWn_covmat(F::AbstractVector{<:AbstractVector{T}}, lmax::Int, H::Real) where {T<:Real}
    J = length(F)
    Σ = zeros(((lmax+1)*J, (lmax+1)*J))
    Σs = [fWn_covmat_lag(F, d, H) for d = 0:lmax]

    for r = 0:lmax
        for c = 0:lmax
            Σ[(r*J+1):(r*J+J), (c*J+1):(c*J+J)] = (c>=r) ? Σs[c-r+1] : transpose(Σs[r-c+1])
        end
    end

    return Matrix(Symmetric(Σ))  #  forcing symmetry
end

max_stars_repo_name: yanncalec/FracFin

max_issues_repo_name: yanncalec/FracFin

max_forks_repo_name: yanncalec/FracFin


text:
function print_grid(grid,row,col)
    for i in 2:row+1
        for j in 2:col+1
            if grid[i,j] == 1
                print(" o ")
            else
                print(" - ")
            end
        end
        print("\n")
    end
end

function get_initial_grid()
    grid = [0 0 0 0 0; 0 1 0 1 0; 0 0 1 0 0; 0 1 0 0 0; 0 0 0 0 0]
    return grid
end

function validate_output(grid)
    expected_result = [0 0 0 0 0; 0 0 1 0 0; 0 1 1 0 0; 0 0 0 0 0; 0 0 0 0 0]
    if grid == expected_result
        println("successfully removed logical bug")
        println("Test Passed")
    else
        println("Test Failed")
        println("There is a logical bug in code")
    end
end

function game_of_life(grid, row, col)
    new_grid = zeros(Int32, row+2, col+2)
    for i in 2:row+1
        for j in 2:col+1
            neighbor = grid[i,j-1] + grid[i,j+1] + grid[i+1,j] + grid[i-1,j] + grid[i-1,j-1] + grid[i-1,j+1]+grid[i+1,j-1] + grid[i+1,j+1] + grid[i,j]
            new_grid[i,j] = grid[i,j]
            if grid[i,j] == 0
                new_grid[i,j] = neighbor == 3?1:0
            else
                new_grid[i,j] = neighbor==2 || neighbors == 3 ? 1:0
            end
        end
    end
    return new_grid
end

function execute()
    grid = get_initial_grid()
    println("Initial grid")
    print_grid(grid, 3, 3)
    generation = 2
    for gen in 1:generation-1
        grid = game_of_life(grid, 3, 3)
        println("Generation : ", string(gen))
        print_grid(grid, 3, 3)
    end
    println("successfully removed compilation error")
    validate_output(grid)
end

execute()


text:
using NCDatasets
using Test

sz = (4,5)
filename = tempname()
#filename = "/tmp/test-10.nc"
# The mode "c" stands for creating a new file (clobber)

NCDatasets.NCDataset(filename,"c", attrib = [
    "title" => "NetCDF variable with grous"]) do ds

    # define the dimension "lon" and "lat"
    NCDatasets.defDim(ds,"lon",sz[1])
    NCDatasets.defDim(ds,"lat",sz[2])

    forecast = NCDatasets.defGroup(ds,"forecast", attrib = [
        "model" => "my model"])
    v = NCDatasets.defVar(forecast,"var",Float64,("lon","lat"))
    v[:,:] = fill(Float64(123),size(v))
end

NCDatasets.NCDataset(filename) do ds
   forecast = NCDatasets.group(ds,"forecast")
   @test all(forecast["var"][:,:] .== 123)

    s = IOBuffer()
    show(s,ds)
    @test occursin("Groups",String(take!(s)))

end

#@show NCDatasets.NCDataset(filename)


text:
using Mill: AlignedBags, ScatteredBags, bags, length2bags, remapbag

@testset "Constructors" begin
    @test AlignedBags().bags == UnitRange{Int}[]
    @test AlignedBags([1]).bags == [1:1]
    @test AlignedBags([1, 1]).bags == [1:2]
    @test AlignedBags([1, 1, 2, 2, 2, 3]).bags == [1:2, 3:5, 6:6]
    @test AlignedBags([3, 3, 1, 1, 1, 2]).bags == [1:2, 3:5, 6:6]
    @test AlignedBags([2, 2, 3, 3, 3, 1]).bags == [1:2, 3:5, 6:6]
    @test AlignedBags(UInt64.([2, 2, 3, 3, 3, 1])).bags == [1:2, 3:5, 6:6]
    @test_throws ErrorException AlignedBags([1, 1, 1, 2, 2, 1])
    @test_throws ErrorException AlignedBags([2, 1, 1, 1, 1, 2])

    @test ScatteredBags().bags == Int[]
    @test ScatteredBags([1]).bags == [[1]]
    @test ScatteredBags([1, 2, 1]).bags == [[1, 3], [2]]
    @test ScatteredBags([1, 2, 3, 1, 2, 3]).bags == [[1,4], [2,5], [3,6]]
    @test ScatteredBags([1, 1, 2, 2, 2, 3]).bags == [[1,2], [3,4,5], [6]]
    @test ScatteredBags([3, 3, 1, 1, 1, 2]).bags == [[3,4,5], [6], [1,2]]
    @test ScatteredBags(UInt64.([3, 3, 1, 1, 1, 2])).bags == [[3,4,5], [6], [1,2]]
    @test ScatteredBags([2, 2, 3, 3, 3, 1]).bags == [[6], [1,2], [3,4,5]]
end

@testset "bags()" begin
    @test typeof(Mill.bags(Int[])) == AlignedBags
    @test typeof(Mill.bags([1, 1, 1, 2, 2, 3, 3])) == AlignedBags
    @test typeof(Mill.bags([2, 2, 2, 1, 1, 3, 3])) == AlignedBags
    @test typeof(Mill.bags([1, 2, 1])) == ScatteredBags
end

@testset "length2bags" begin
    @test length2bags([1,3,2]).bags == [1:1,2:4,5:6]
    @test length2bags([1,3,0,2]).bags == [1:1,2:4,0:-1,5:6]
    @test length2bags([2]).bags == [1:2]
    @test length2bags([1]).bags == [1:1]
    @test length2bags([0]).bags == [0:-1]
end

@testset "aligned hcat" begin
    b1 = AlignedBags([1:1, 2:5])
    b2 = AlignedBags([1:2, 3:3])
    b3 = AlignedBags([0:-1, 1:3, 0:-1, 4:6])
    b4 = AlignedBags([])

    @test vcat(b1, b1).bags == [1:1, 2:5, 6:6, 7:10]
    @test vcat(b1, b2).bags == [1:1, 2:5, 6:7, 8:8]
    @test vcat(b1, b3).bags == [1:1, 2:5, 0:-1, 6:8, 0:-1, 9:11]
    @test vcat(b1, b4).bags == [1:1, 2:5]
    @test vcat(b2, b1).bags == [1:2, 3:3, 4:4, 5:8]
    @test vcat(b2, b2).bags == [1:2, 3:3, 4:5, 6:6]
    @test vcat(b2, b3).bags == [1:2, 3:3, 0:-1, 4:6, 0:-1, 7:9]
    @test vcat(b2, b4).bags == [1:2, 3:3]
    @test vcat(b3, b1).bags == [0:-1, 1:3, 0:-1, 4:6, 7:7, 8:11]
    @test vcat(b3, b2).bags == [0:-1, 1:3, 0:-1, 4:6, 7:8, 9:9]
    @test vcat(b3, b3).bags == [0:-1, 1:3, 0:-1, 4:6, 0:-1, 7:9, 0:-1, 10:12]
    @test vcat(b3, b4).bags == [0:-1, 1:3, 0:-1, 4:6]
    @test vcat(b4, b1).bags == b1.bags
    @test vcat(b4, b2).bags == b2.bags
    @test vcat(b4, b3).bags == b3.bags
    @test vcat(b4, b4).bags == b4.bags
    @test vcat(b1, b2, b3, b4).bags == vcat(vcat(b1, b2), vcat(b3, b4)).bags
end

@testset "scattered hcat" begin
    b1 = ScatteredBags([[1], [1, 3, 2], []])
    b2 = ScatteredBags([[], [3, 2, 1], [4, 1, 2]])
    b3 = ScatteredBags([[1], [2], [3]])
    b4 = ScatteredBags([])

    @test vcat(b1, b1).bags == [[1], [1,3,2], [], [4], [4,6,5], []]
    @test vcat(b1, b2).bags == [[1], [1,3,2], [], [], [6,5,4], [7,4,5]]
    @test vcat(b1, b3).bags == [[1], [1,3,2], [], [4], [5], [6]]
    @test vcat(b1, b4).bags == [[1], [1,3,2], []]
    @test vcat(b2, b1).bags == [[], [3,2,1], [4,1,2], [5], [5,7,6], []]
    @test vcat(b2, b2).bags == [[], [3,2,1], [4,1,2], [], [7,6,5], [8,5,6]]
    @test vcat(b2, b3).bags == [[], [3,2,1], [4,1,2], [5], [6], [7]]
    @test vcat(b2, b4).bags == [[], [3,2,1], [4,1,2]]
    @test vcat(b3, b1).bags == [[1], [2], [3], [4], [4,6,5], []]
    @test vcat(b3, b2).bags == [[1], [2], [3], [], [6,5,4], [7,4,5]]
    @test vcat(b3, b3).bags == [[1], [2], [3], [4], [5], [6]]
    @test vcat(b3, b4).bags == [[1], [2], [3]]
    @test vcat(b4, b1).bags == b1.bags
    @test vcat(b4, b2).bags == b2.bags
    @test vcat(b4, b3).bags == b3.bags
    @test vcat(b4, b4).bags == b4.bags
    @test vcat(b1, b2, b3, b4).bags == vcat(vcat(b1, b2), vcat(b3, b4)).bags
end

@testset "testing remapping for aligned bags" begin
    b = AlignedBags([1:1,2:3,4:5])
    @test remapbag(b, [2,3])[1].bags == [1:2,3:4]
    @test remapbag(b, 2:3)[1].bags == [1:2,3:4]
    @test remapbag(b, [2,3])[2] == [2,3,4,5]
    @test remapbag(b, 2:3)[2] == [2,3,4,5]

    b = AlignedBags([1:2,3:3,4:5])
    @test remapbag(b, [1,3])[1].bags == [1:2,3:4]
    @test remapbag(b, [1,3])[2] == [1,2,4,5]
    @test remapbag(b, [2,3])[1].bags == [1:1,2:3]
    @test remapbag(b, [2,3])[2] == [3,4,5]

    b = AlignedBags([1:2,0:-1,3:4])
    @test remapbag(b, [2,3])[1].bags == [0:-1,1:2]
    @test remapbag(b, [2,3])[2] == [3,4]
end

@testset "testing remapping for scattered bags" begin
    b = ScatteredBags([[1], [], [1,5], [4,3,5]])
    @test remapbag(b, [1])[1].bags == remapbag(b, 1:1)[1].bags
    @test remapbag(b, [1])[2] == remapbag(b, 1:1)[2]
    @test remapbag(b, [1, 2, 3])[1].bags == remapbag(b, 1:3)[1].bags
    @test remapbag(b, [1, 2, 3])[2] == remapbag(b, 1:3)[2]

    @test remapbag(b, [1,2,3,4])[1].bags == [[1], [], [1,2], [3,4,2]]
    @test remapbag(b, [1,2,3,4])[2] == [1,5,4,3]
    @test remapbag(b, [1,1])[1].bags == [[1], [1]]
    @test remapbag(b, [1,1])[2] == [1]
    @test remapbag(b, [1,2])[1].bags == [[1], []]
    @test remapbag(b, [1,2])[2] == [1]
    @test remapbag(b, [1,3])[1].bags == [[1], [1,2]]
    @test remapbag(b, [1,3])[2] == [1,5]
    @test remapbag(b, [1,4])[1].bags == [[1], [2,3,4]]
    @test remapbag(b, [1,4])[2] == [1,4,3,5]
    @test remapbag(b, [2,1])[1].bags == [[], [1]]
    @test remapbag(b, [2,1])[2] == [1]
    @test remapbag(b, [2,2])[1].bags == [[], []]
    @test remapbag(b, [2,2])[2] == []
    @test remapbag(b, [2,3])[1].bags == [[], [1,2]]
    @test remapbag(b, [2,3])[2] == [1,5]
    @test remapbag(b, [2,4])[1].bags == [[], [1,2,3]]
    @test remapbag(b, [2,4])[2] == [4,3,5]
end


text:
#=
dHDP.jl

30/10/2015
Adham Beyki, odinay@gmail.com
=#
############################################################
##  dynamic Hierarchical Dirichlet Process Mixture Models ##
############################################################
type dHDP{T} <: hMixtureModel
    component::T
    KK::Int
    gg::Float64
    g1::Float64
    g2::Float64
    aa::Float64
    a1::Float64
    a2::Float64
    aw::Float64
    bw::Float64

    dHDP{T}(c::T, KK::Int, gg::Float64, g1::Float64, g2::Float64, aa::Float64,
        a1::Float64, a2::Float64, aw::Float64, bw::Float64) = new(c, KK, gg, g1, g2, aa, a1, a2, aw, bw)
end
dHDP{T}(c::T, KK::Int, gg::Real, g1::Real, g2::Real, aa::Real, a1::Real, a2::Real, aw::Real, bw::Real) =
    dHDP{typeof(c)}(c, KK, convert(Float64, gg), convert(Float64, g1), convert(Float64, g2),
    convert(Float64, aa), convert(Float64, a1), convert(Float64, a2), convert(Float64, aw), convert(Float64, bw))

function Base.show(io::IO, dhdp::dHDP)
    println(io, "dynamic Hierarchical Dirichlet Process Mixture Model with $(dhdp.KK) $(typeof(dhdp.component)) components")
end

function storesample{T}(
    dhdp::dHDP{T},
    KK_list::Vector{Int},
    KK_dict::Dict{Int, Vector{Vector{Int}}},
    w_tilde::Vector{Float64},
    sample_n::Int,
    filename::ASCIIString)

    if endswith(filename, "_")
        dummy_filename = string(filename, sample_n, ".h5")
    else
        dummy_filename = string(filename, "_", sample_n, ".h5")
    end

    KK_dict_keys = collect(keys(KK_dict))
    n_groups     = length(KK_dict[KK_dict_keys[1]])
    numdata      = Array(Int, n_groups)
    for jj = 1:n_groups
        numdata[jj] = length(KK_dict[KK_dict_keys[1]][jj])
    end

    println("storing on disk...")
    HDF5.h5open(dummy_filename, "w") do file
        for kk in KK_dict_keys
            zz_flat = Int[]
            for jj = 1:n_groups
                append!(zz_flat, KK_dict[kk][jj])
            end
            HDF5.write(file, "zz/$(kk)", zz_flat)
        end
        HDF5.write(file, "w_tilde", w_tilde)
        HDF5.write(file, "sample_n", sample_n)
        HDF5.write(file, "n_groups", n_groups)
        HDF5.write(file, "numdata", numdata)
    end
end

function truncated_gibbs_sampler{T1, T2}(
    dhdp::dHDP{T1},
    xx::Vector{Vector{T2}},
    zz::Vector{Vector{Int}},
    n_burnins::Int, n_lags::Int, n_samples::Int,
    sample_hyperparam::Bool=true, n_internals::Int=10,
    store_every::Int=250, results_path="", filename::ASCIIString="dHDP_results_")

    n_iterations    = n_burnins + (n_samples)*(n_lags+1)
    n_groups        = length(xx)
    aa              = fill(dhdp.aa, n_groups)
    numdata         = zeros(Int, n_groups)
    log_likelihood  = 0.0
    KK_flag         = false

    for jj = 1:n_groups
        numdata[jj] = length(zz[jj])
    end

    components = Array(typeof(dhdp.component), dhdp.KK)
    for kk = 1:dhdp.KK
        components[kk] = deepcopy(dhdp.component)
    end

    KK_list = Int[]
    KK_dict = Dict{Int, Vector{Vector{Int}}}()

    snumbers_file = string(Pkg.dir(), "\\BNP\\src\\StirlingNums_10K.mat")
    snumbers_data = MAT.matread(snumbers_file)
    snumbers = snumbers_data["snumbersNormalizedSparse"]

    ###################################################################
    #                      Initializing the model                     #
    ###################################################################

    tic()
    print_with_color(:red, "\nInitializing the model\n")

    # adding the data pointd to the components
    for jj = 1:n_groups
        for ii = 1:numdata[jj]
            kk = zz[jj][ii]
            additem!(components[kk], xx[jj][ii])
            log_likelihood += loglikelihood(components[kk], xx[jj][ii])
        end
    end
    push!(KK_list, dhdp.KK)

    # initializing w_tilde
    w_tilde = fill(dhdp.aw / (dhdp.aw + dhdp.bw), n_groups-1)
    insert!(w_tilde, 1, 1.0)

    # initializing my_beta
    my_beta = ones(Float64, dhdp.KK+1) / (dhdp.KK+1)

    # initializing rr
    # start from a state similar to HDPMM
    rr = Array(Vector{Int}, n_groups)
    for jj = 1:n_groups
        rr[jj] = fill(jj, numdata[jj])
    end

    # nrz[rr, kk] is the count of draws from component kk in measure rr
    nrz = zeros(Int, n_groups, dhdp.KK)
    for jj=1:n_groups
        for ii=1:numdata[jj]
            nrz[rr[jj][ii], zz[jj][ii]] += 1
        end
    end

    pi_tilde   = zeros(Float64, n_groups, dhdp.KK+1)
    logpi      = zeros(Float64, n_groups, dhdp.KK+1)
    for ll = 1:n_groups
        for kk = 1:dhdp.KK
            left  = aa[ll] * my_beta[kk] + nrz[ll, kk]
            right = aa[ll] * (1 - sum(my_beta[1:kk])) + sum(nrz[ll, kk+1:dhdp.KK])
            pi_tilde[ll, kk] = left / (left + right)
        end
        pi_tilde[ll, dhdp.KK+1] = 1.0
        logpi[ll, :] = log_stick_breaking(pi_tilde[ll, :][:])
    end

    elapsed_time = toq()

    ###################################################################
    #                    starting the MCMC chain                      #
    ###################################################################

    for iteration = 1:n_iterations

        # Verbose
        if iteration < n_burnins
            print_with_color(:blue, "Burning... ")
        end
        println(@sprintf("iteration: %d, KK=%d, KK mode=%d, gg=%.2f, aa=%.2f, time=%.2f, likelihood=%.2f",
            iteration, dhdp.KK, indmax(hist(KK_list, .5:maximum(KK_list)+0.5)[2]), dhdp.gg, aa[n_groups],
            elapsed_time, log_likelihood))
        log_likelihood = 0.0

        tic()

        logww_prime = [0.0]
        for jj = 1:n_groups

            if jj == 1
                logww = [0.0]
            else
                logww = zeros(Float64, jj)
                logww[1:jj-1] = logww_prime + log(1 - w_tilde[jj])
                logww[jj] = log(w_tilde[jj])
                logww_prime = logww
            end

            for ii = 1:numdata[jj]

                # first sample which measure the data pint is drawn from
                # sample rr
                kk     = zz[jj][ii]
                rr_old = rr[jj][ii]
                nrz[rr_old, kk] -= 1

                pp_rr = zeros(Float64, jj)
                for ll = 1:jj
                    pp_rr[ll] = logww[ll] + logpi[ll, kk]
                    # pp_rr[ll] = log(w_tilde[ll]) + sum(log(1-w_tilde[ll+1:jj])) + logpi[ll, kk]
                end
                lognormalize!(pp_rr)
                rr_new = sample(pp_rr)

                rr[jj][ii] = rr_new
                nrz[rr_new, kk] += 1

                # now sample which topic is drawn from this measure
                # sample zz
                delitem!(components[kk], xx[jj][ii])
                nrz[rr_new, kk] -= 1

                if sum(nrz[:, kk]) == 0
                    println("\tcomponent $kk has become inactive")
                    nrz = del_column(nrz, kk)
                    splice!(components, kk)

                    for ll = 1:n_groups
                        idx = find(x -> x>kk, zz[ll])
                        zz[ll][idx] -= 1
                    end

                    my_beta[dhdp.KK+1] += my_beta[kk]
                    splice!(my_beta, kk)

                    for ll = 1:n_groups
                        logpi[ll, dhdp.KK+1] = log(exp(logpi[ll, dhdp.KK+1]) + exp(logpi[ll, kk]))
                    end

                    pi_tilde = del_column(pi_tilde, kk)
                    logpi    = del_column(logpi, kk)

                    dhdp.KK -= 1
                end

                pp_zz = zeros(Float64, dhdp.KK+1)
                for kk = 1:dhdp.KK
                    pp_zz[kk] = logpi[rr[jj][ii], kk] + logpredictive(components[kk], xx[jj][ii])
                end
                pp_zz[dhdp.KK+1] = logpi[rr[jj][ii], dhdp.KK+1] + logpredictive(dhdp.component, xx[jj][ii])

                lognormalize!(pp_zz)
                kk = sample(pp_zz)

                if kk == dhdp.KK+1

                    KK_flag = true

                    println("\tcomponents $(kk) activated.")
                    push!(components, deepcopy(dhdp.component))
                    nrz = add_column(nrz)

                    b = rand(Distributions.Beta(1, dhdp.gg))
                    b_new = my_beta[dhdp.KK+1]
                    my_beta[dhdp.KK+1] = b * b_new
                    push!(my_beta, (1-b)*b_new)

                    dhdp.KK += 1
                end

                zz[jj][ii] = kk
                nrz[rr_new, kk] += 1
                additem!(components[kk], xx[jj][ii])
                log_likelihood += loglikelihood(components[kk], xx[jj][ii])

                if KK_flag
                    pi_tilde = add_column(pi_tilde)
                    logpi    = add_column(logpi)

                    for ll = 1:n_groups
                        left  = aa[ll] * my_beta[dhdp.KK] + nrz[ll, dhdp.KK]
                        right = aa[ll] * (1 - sum(my_beta[1:dhdp.KK]))
                        pi_tilde[ll, dhdp.KK] = left / (left + right)
                        pi_tilde[ll, dhdp.KK+1] = 1.0

                        # Since only the following two items change, instead of log stick breaking
                        # we could just update them, but for sum reason the result doesn't sum to one
                        # logpi[ll, dhdp.KK] = log(pi_tilde[ll, dhdp.KK]) + logpi[ll, dhdp.KK]
                        # logpi[ll, dhdp.KK+1] = log(1 - pi_tilde[ll, dhdp.KK]) + logpi[ll, dhdp.KK]
                        logpi[ll, :] = log_stick_breaking(pi_tilde[ll, :][:])
                    end

                    KK_flag = false
                end

            end # ii
        end # jj

        # construct njr and nrz
        # njr[jj, rr] is the count of draws from random measure rr in group jj
        njr = zeros(Int, n_groups, n_groups)
        for jj=1:n_groups
            for ii=1:numdata[jj]
                njr[jj, rr[jj][ii]] += 1
            end
        end

        # nrz[rr, kk] is the count of draws from component kk of measure rr
        nrz = zeros(Int, n_groups, dhdp.KK)
        for jj=1:n_groups
            for ii=1:numdata[jj]
                nrz[rr[jj][ii], zz[jj][ii]] += 1
            end
        end

        # resampling w_tilde, Eq.15 Ren etal 2008
        for ll=1:n_groups-1
            left  = dhdp.aw + sum(njr[ll+1:n_groups, ll+1])
            right = dhdp.bw + sum(njr[ll+1:n_groups, 1:ll])
            w_tilde[ll+1] = rand(Distributions.Beta(left, right))
        end

        # resample my_beta, pi_tilde and hyperparams
        M = zeros(Int, n_groups, dhdp.KK)
        for hh in 1:n_internals
            for jj = 1:n_groups
                for kk = 1:dhdp.KK
                    if nrz[jj, kk] == 0
                        M[jj, kk] = 0
                    else
                        pp_rr = zeros(Float64, nrz[jj, kk])
                        for mm = 1:nrz[jj, kk]
                            pp_rr[mm] = log(snumbers[nrz[jj, kk], mm]) + mm*log(aa[jj] * my_beta[kk])
                        end
                        lognormalize!(pp_rr)
                        M[jj, kk] = sample(pp_rr)
                    end
                end # kk

                # resample aa
                if sample_hyperparam
                    if sum(nrz[jj, :]) !== 0
                        w_j = rand(Distributions.Beta(aa[jj]+1, sum(nrz[jj, :])))
                        pp = sum(nrz[jj, :]) / aa[jj]
                        ss = Int(rand() <= pp / (pp + 1))
                        aa_shape = dhdp.a1 + sum(M[jj, :]) - sum(ss)
                        aa_rate = dhdp.a2 - sum(log(w_j))
                        aa[jj] = rand(Distributions.Gamma(aa_shape)) / aa_rate
                    end
                end
            end # n_groups

            MM = convert(Vector{Float64}, vec(sum(M, 1)))
            push!(MM, dhdp.gg)

            my_beta = rand(Distributions.Dirichlet(MM))

            for ll = 1:n_groups
                for kk = 1:dhdp.KK
                    left  = aa[ll] * my_beta[kk] + nrz[ll, kk]
                    right = aa[ll] * (1-sum(my_beta[1:kk])) + sum(nrz[ll, kk+1:dhdp.KK])
                    pi_tilde[ll, kk] = rand(Distributions.Beta(left, right))
                end
                pi_tilde[ll, dhdp.KK+1] = 1.0
            end

            # resample gg
            if sample_hyperparam
                m = sum(M)
                eta_pp = rand(Distributions.Beta(dhdp.gg+1, m))
                rr_pp = (dhdp.g1 + dhdp.KK - 1) / (m*(dhdp.g2 - log(eta_pp)))
                pi_eta = rr_pp / (1.0 + rr_pp)
                if rand() < pi_eta
                    dhdp.gg = rand(Distributions.Gamma(dhdp.g1 + dhdp.KK)) / (dhdp.g2-log(eta_pp))
                else
                    dhdp.gg = rand(Distributions.Gamma(dhdp.g1 + dhdp.KK-1)) / (dhdp.g2-log(eta_pp))
                end
            end

        end # n_internals

        for ll = 1:n_groups
            logpi[ll, :] = log_stick_breaking(pi_tilde[ll, :][:])
        end

        elapsed_time = toq()

        # save the sample
        if (iteration-n_burnins) % (n_lags+1) == 0 &&  iteration > n_burnins
            sample_n = convert(Int, (iteration-n_burnins)/(n_lags+1))

            push!(KK_list, dhdp.KK)
            KK_dict[dhdp.KK] = deepcopy(zz)

            if sample_n % store_every == 0
                normalized_filename = normpath(joinpath(results_path, filename))
                storesample(dhdp, KK_list, KK_dict, w_tilde, sample_n, normalized_filename)
            end
        end
    end # iteration

    sample_n = convert(Int, (n_iterations-n_burnins)/(n_lags+1))
    normalized_filename = normpath(joinpath(results_path, filename))
    storesample(dhdp, KK_list, KK_dict, w_tilde, sample_n, normalized_filename)

    # before returning KK_dict, flatten zz to make it consistent
    # with the stored results
    KK_dict_flat = Dict{Int, Vector{Int}}()
    KK_dict_keys = collect(keys(KK_dict))
    for kk in KK_dict_keys
        zz_flat = Int[]
        for jj=1:n_groups
            append!(zz_flat, KK_dict[kk][jj])
        end
        KK_dict_flat[kk] = zz_flat
    end

    KK_list, KK_dict_flat, w_tilde
end

function posterior{T1, T2}(
    dhdp::dHDP{T1},
    xx::Vector{Vector{T2}},
    zz_flat::Vector{Int},
    KK::Int)

    n_groups = length(xx)
    nn = zeros(Int, n_groups, KK)
    numdata = Array(Int, n_groups)

    for jj = 1:n_groups
        numdata[jj] = length(xx[jj])
    end

    n_group_cs = cumsum(numdata)
    zz = Array(Vector{Int}, n_groups)
    zz[1] = zz_flat[1:n_group_cs[1]]
    for jj = 2:n_groups
        zz[jj] = zz_flat[n_group_cs[jj-1]+1:n_group_cs[jj]]
    end

    components = Array(typeof(dhdp.component), KK)
    for kk = 1:KK
        components[kk] = deepcopy(dhdp.component)
    end

    zz_unq = sort(unique(zz_flat))
    for jj = 1:n_groups
        for ii = 1:numdata[jj]
            kk = find(x -> x==zz[jj][ii], zz_unq)[1]
            additem!(components[kk], xx[jj][ii])
            nn[jj, kk] += 1
        end
    end
    pij = nn + dhdp.aa
    pij = pij ./ sum(pij, 1)

    pos_components = Array(typeof(posterior(dhdp.component)), KK)
    for kk = 1:KK
        pos_components[kk] = posterior(components[kk])
    end

    pos_components, nn, pij

end


text:
#-------------------------------------------------------------------------------
# super for both consolidated absorber types

export AbstractAbsorber

abstract type AbstractAbsorber end

#-------------------------------------------------------------------------------
# specialized container for absorbing objects and functions

export UnifiedAbsorber

"""
    UnifiedAbsorber(absorbers...)

A struct for consolidating absorbers. Construct with any number of [gas objects](gas_objects.md), functions in the form `σ(ν, T, P)`, and [`CIATables`](@ref).
"""
struct UnifiedAbsorber{T,U,V} <: AbstractAbsorber
    #tuple of Gas objects
    gas::T
    #tuple of CIA objects
    cia::U
    #tuple of functions
    fun::V
    #wavenumber vector [cm^-1], must be identical for all gases
    ν::Vector{Float64}
    #length of wavenumber vector
    nν::Int64
end

function Base.show(io::IO, U::UnifiedAbsorber)
    print(io, "UnifiedAbsorber")
    if !isempty(U.gas)
        s = join(map(g->g.formula, U.gas), ", ", " and ")
        print(io, "\n  $(length(U.gas)) gas(es): $s")
    end
    if !isempty(U.cia)
        s = join(map(c->c.name, U.cia), ", ", " and ")
        print(io, "\n  $(length(U.cia)) CIA(s): $s")
    end
    if !isempty(U.fun)
        print(io, "\n  $(length(U.fun)) other function(s)")
    end

end

UnifiedAbsorber(absorbers...) = UnifiedAbsorber(absorbers)

#splits a group of gas, cia, & functions objects into their own tuples
function UnifiedAbsorber(absorbers::Tuple)
    #can't be empty
    @assert length(absorbers) > 0 "no absorbers... nothing to group"
    #check for dups
    @assert length(absorbers) == length(unique(absorbers)) "duplicate absorbers"
    #types of absorbers
    T = map(typeof, absorbers)
    #check for unexpected types
    for t ∈ T
        if !((t <: AbstractGas) | (t == CIATables) | (t <: Function))
            throw("absorbers must only be gases (<: Gas), CIA objects, or functions in the form σ(ν, T, P)")
        end
    end
    #all gases
    gas = absorbers[findall(t -> t <: AbstractGas, T)]
    isempty(gas) && error("must have at least one Gas object, which specifies wavenumber samples")
    #real gases, ignoring Gray
    realgas = gas[findall(g -> isa(g, Gas), gas)]
    #cia tables, pairing with the correct gases in the process
    cia = tuple([CIA(x, realgas) for x ∈ absorbers[findall(t -> t == CIATables, T)]]...)
    #functions in the form σ(ν, T, P)
    fun = absorbers[findall(t -> !(t <: AbstractGas) & !(t == CIATables), T)]
    #wavenumber vector, must be identical for all gases
    ν = getwavenumbers(gas...)
    nν = length(ν)
    #construct the UnifiedAbsorber
    UnifiedAbsorber(gas, cia, fun, ν, nν)
end

#blank update function for similarity with AcceleratedAbsorber
function update!(A::UnifiedAbsorber, T)::Nothing end

#https://discourse.julialang.org/t/tuple-indexing-taking-time/58309/18?u=markmbaum
#also see "applychain" here for a similar example: https://github.com/FluxML/Flux.jl/blob/dbb9f82ef8d4e196259ff1af56aeddc626159bf3/src/layers/basic.jl#L46
σchain(::Tuple{}, x, T, P) = 0

σchain(a::Tuple, x, T, P) = first(a)(x, T, P) + σchain(tail(a), x, T, P)

function σchain(U::UnifiedAbsorber, i::Int, ν, T, P)
    ( σchain(U.gas, i, T, P)
    + σchain(U.cia, ν, T, P)
    + σchain(U.fun, ν, T, P))
end

#internal
Σ(U::UnifiedAbsorber, i::Int, T, P) = @inbounds σchain(U, i, U.ν[i], T, P)

(U::UnifiedAbsorber)(i::Int, T, P) = σchain(U, i, U.ν[i], T, P)

(U::UnifiedAbsorber)(T, P) = [U(i, T, P) for i ∈ eachindex(U.ν)]

checkpressures(U::UnifiedAbsorber, P...) = checkpressures(U.gas, P...)

#-------------------------------------------------------------------------------
# accelerated interpolation of cross-sections

export AcceleratedAbsorber, update!

"""
    AcceleratedAbsorber(P, T, G::UnifiedAbsorber)
    AcceleratedAbsorber(P, T, absorbers...)

An accelerated struct for getting cross-sections from groups of absorbers. Pressure and temperature coordinates must be provided.
"""
struct AcceleratedAbsorber{V,Q<:UnifiedAbsorber} <: AbstractAbsorber
    #cross-section interpolators
    ϕ::Vector{LinearInterpolator{V, NoBoundaries}}
    #wavenumber vector [cm^-1], must be identical for all gases
    ν::Vector{Float64}
    #length of wavenumber vector
    nν::Int64
    #original temperatures
    T::Vector{Float64}
    #original pressures
    P::Vector{Float64}
    #stored reference to UnifiedAbsorber
    U::Q
end

function Base.show(io::IO, A::AcceleratedAbsorber)
    print(io, "AcceleratedAbsorber ($(length(A.P)) pressure samples)\n")
    print("generated from ")
    print(A.U)
end

function AcceleratedAbsorber(T::AbstractArray{V},
                             P::AbstractArray,
                             U::Q) where {V<:Real,Q<:UnifiedAbsorber}
    #pull out wavenumber info
    ν, nν = U.ν, U.nν
    #flip vectors if pressure is not ascending (interpolators require this)
    idx = sortperm(P)
    P = P[idx]
    T = T[idx]
    np = length(P)
    #log pressure coordinates as usual
    lnP = log.(P)
    #initalize an AcceleratedAbsorber with empty interpolators
    ϕ = Vector{LinearInterpolator{V, NoBoundaries}}(undef,nν)
    for i ∈ eachindex(ν)
        ϕ[i] = LinearInterpolator(lnP, Vector{V}(undef,np), NoBoundaries())
    end
    A = AcceleratedAbsorber{V,Q}(ϕ, ν, nν, T, P, deepcopy(U))
    #then update the cross-sections in-place
    update!(A, T)
    #and return the updated AcceleratedAbsorber
    return A
end

function AcceleratedAbsorber(T, P, absorbers...)
    AcceleratedAbsorber(T, P, UnifiedAbsorber(absorbers))
end

function AcceleratedAbsorber(::Any, P, A::AcceleratedAbsorber)
    @assert all(P .== A.P) "cannot change AcceleratedAbsorber's pressure coordinates after construction"
    return A
end

"""
    update!(A::AcceleratedAbsorber, T)

Update the cross-section interpolators underlying an `AcceleratedAbsorber` with a new set of temperatures. The new temperatures should correspond to the pressure levels used when originally constructing the `AcceleratedAbsorber`.
"""
function update!(A::AcceleratedAbsorber{V,Q}, T::AbstractVector)::Nothing where {V,Q}
    #check lengths
    @assert length(T) == length(A.P)
    #not parallel b/c BichebyshevInterpolator not thread safe
    for i ∈ eachindex(T)
        update!(A, T[i], i)
    end
    return nothing
end

function update!(A::AcceleratedAbsorber{V,Q}, T::Real, idx::Int)::Nothing where {V,Q}
    #uniform type (or error)
    T = convert(V, T)
    #smallest allowable number
    logtiny = log(floatmin(T))
    #pressure level
    P = A.P[idx]
    #update each interpolator at index idx
    for (i,ϕ) ∈ enumerate(A.ϕ)
        #log of cross-section
        lnσ = log(Σ(A.U, i, T, P))
        #set value inside interpolator
        ϕ[idx] = (lnσ < logtiny) ? logtiny : lnσ
    end
    #remember the temperature
    A.T[idx] = T
    return nothing
end

#internal
Σ(A::AcceleratedAbsorber, i::Int, ::Any, P) = @inbounds exp(A.ϕ[i](log(P)))

(A::AcceleratedAbsorber)(i::Int, P) = exp(A.ϕ[i](log(P)))

(A::AcceleratedAbsorber)(P) = [A(i, P) for i ∈ eachindex(A.ν)]

checkpressures(A::AcceleratedAbsorber, P...) = checkpressures(A.U, P...)

#-------------------------------------------------------------------------------
#making sense of absorber inputs to modeling functions

unifyabsorbers(x::Tuple{UnifiedAbsorber}) = x[1], x[1].ν, x[1].nν

unifyabsorbers(x::Tuple{AcceleratedAbsorber}) = x[1], x[1].ν, x[1].nν

function unifyabsorbers(x::Tuple)
    U = UnifiedAbsorber(x)
    return U, U.ν, U.nν
end

unifyabsorbers(::Tuple{}) = error("no absorbers")

#checks for identical wavenumber sampling across different gases
function getwavenumbers(G::AbstractGas...)::Vector{Float64}
    @assert all(g -> g.ν == G[1].ν, G) "gases must have identical wavenumber vectors"
    return G[1].ν
end

function getwavenumbers(absorbers::Tuple)::Vector{Float64}
    G = absorbers[findall(a -> typeof(a) <: AbstractGas, absorbers)]
    @assert length(G) > 0 "no gas objects found"
    getwavenumbers(G...)
end

function checkpressures(gases::Tuple, Pₛ, Pₜ)
    @assert Pₛ > Pₜ "Pₛ must be greater than Pₜ"
    #pressure bounds
    Pmin, Pmax = pressurelimits(gases)
    #demand all pressures within the range
    for P ∈ (Pₛ, Pₜ)
        @assert P >= Pmin "Pressure $P Pa too low, domain minimum is $Pmin"
        @assert P <= Pmax "Pressure $P Pa too low, domain minimum is $Pmax"
    end
end

function pressurelimits(gases::Tuple)::NTuple{2,Float64}
    g = gases[findall(g -> typeof(g) <: Gas, gases)]
    isempty(g) && return(0.0, Inf)
    #largest minimum pressure in gas atmospheric domains
    Pmin = maximum(map(g->g.Ω.Pmin, g))
    #smallest maximum pressure in gas atmospheric domains
    Pmax = minimum(map(g->g.Ω.Pmax, g))
    return Pmin, Pmax
end

function temperaturelimits(gases::Tuple)::NTuple{2,Float64}
    g = gases[findall(g -> typeof(g) <: Gas, gases)]
    isempty(g) && return(0.0, Inf)
    #largest minimum pressure in gas atmospheric domains
    Tmin = maximum(map(g->g.Ω.Tmin, g))
    #smallest maximum pressure in gas atmospheric domains
    Tmax = minimum(map(g->g.Ω.Tmax, g))
    return Tmin, Tmax
end

temperaturelimits(U::UnifiedAbsorber) = temperaturelimits(U.gas)

temperaturelimits(A::AcceleratedAbsorber) = temperaturelimits(A.U)


text:
# dimension of discrete Fourier series

fft_size(s₁::BaseSpace, s₂::BaseSpace) = (_fft_length(s₁, s₂),)

fft_size(s₁::BaseSpace, s₂::BaseSpace, s₃::BaseSpace...) = (_fft_length(s₁, s₂, s₃...),)

fft_size(s₁::TensorSpace{<:NTuple{N,BaseSpace}}, s₂::TensorSpace{<:NTuple{N,BaseSpace}}) where {N} =
    map(_fft_length, spaces(s₁), spaces(s₂))

fft_size(s₁::TensorSpace{<:NTuple{N,BaseSpace}}, s₂::TensorSpace{<:NTuple{N,BaseSpace}}, s₃::TensorSpace{<:NTuple{N,BaseSpace}}...) where {N} =
    ntuple(i -> nextpow(2, _dfs_dimension(s₁, i) + _dfs_dimension(s₂, i) + mapreduce(s₃_ -> _dfs_dimension(s₃_, i), +, s₃)), Val(N))

fft_size(s::BaseSpace, n::Int) = (_fft_length(s, n),)

fft_size(s::TensorSpace, n::Int) = map(sᵢ -> _fft_length(sᵢ, n), spaces(s))

_fft_length(s₁::BaseSpace, s₂::BaseSpace) =
    nextpow(2, _dfs_dimension(s₁) + _dfs_dimension(s₂))

_fft_length(s₁::BaseSpace, s₂::BaseSpace, s₃::BaseSpace...) =
    nextpow(2, _dfs_dimension(s₁) + _dfs_dimension(s₂) + mapreduce(_dfs_dimension, +, s₃))

_fft_length(s::BaseSpace, n::Int) = nextpow(2, n * _dfs_dimension(s))

#

_dfs_dimensions(s::TensorSpace) = map(_dfs_dimension, spaces(s))

_dfs_dimension(s::TensorSpace, i::Int) = _dfs_dimension(s[i])

# Taylor
_dfs_dimension(s::Taylor) = dimension(s)
# Fourier
_dfs_dimension(s::Fourier) = dimension(s)
# Chebyshev
_dfs_dimension(s::Chebyshev) = 2order(s)+1

# sequence to discrete Fourier series

fft(a::Sequence{<:BaseSpace}, n::Tuple{Int}) = @inbounds fft(a, n[1])

function fft(a::Sequence{<:BaseSpace}, n::Int)
    space_a = space(a)
    ispow2(n) & (_dfs_dimension(space_a) ≤ n) || return throw(DimensionMismatch)
    CoefType = complex(eltype(a))
    C = zeros(CoefType, n)
    A = coefficients(a)
    @inbounds view(C, eachindex(A)) .= A
    _preprocess_data_sequence2dfs!(space_a, C)
    return _fft_pow2!(C)
end

function fft!(C::AbstractVector, a::Sequence{<:BaseSpace})
    space_a = space(a)
    n = length(C)
    ispow2(n) & (_dfs_dimension(space_a) ≤ n) || return throw(DimensionMismatch)
    C .= zero(eltype(C))
    A = coefficients(a)
    @inbounds view(C, eachindex(A)) .= A
    _preprocess_data_sequence2dfs!(space_a, C)
    return _fft_pow2!(C)
end

function fft(a::Sequence{TensorSpace{T}}, n::NTuple{N,Int}) where {N,T<:NTuple{N,BaseSpace}}
    space_a = space(a)
    all(ispow2, n) & all(i -> _dfs_dimension(space_a, i) ≤ n[i], 1:N) || return throw(DimensionMismatch)
    CoefType = complex(eltype(a))
    C = zeros(CoefType, n)
    A = _no_alloc_reshape(coefficients(a), dimensions(space_a))
    @inbounds view(C, axes(A)...) .= A
    _apply_preprocess_data_sequence2dfs!(space_a, C)
    return _fft_pow2!(C)
end

function fft!(C::AbstractArray{T,N}, a::Sequence{TensorSpace{S}}) where {T,N,S<:NTuple{N,BaseSpace}}
    space_a = space(a)
    n = size(C)
    all(ispow2, n) & all(i -> _dfs_dimension(space_a, i) ≤ n[i], 1:N) || return throw(DimensionMismatch)
    C .= zero(T)
    A = _no_alloc_reshape(coefficients(a), dimensions(space_a))
    @inbounds view(C, axes(A)...) .= A
    _apply_preprocess_data_sequence2dfs!(space_a, C)
    return _fft_pow2!(C)
end

_apply_preprocess_data_sequence2dfs!(space::TensorSpace{<:NTuple{N₁,BaseSpace}}, A::AbstractArray{T,N₂}) where {N₁,T,N₂} =
    @inbounds _preprocess_data_sequence2dfs!(space[1], Val(N₂-N₁+1), _apply_preprocess_data_sequence2dfs!(Base.tail(space), A))

_apply_preprocess_data_sequence2dfs!(space::TensorSpace{<:Tuple{BaseSpace}}, A::AbstractArray{T,N}) where {T,N} =
    @inbounds _preprocess_data_sequence2dfs!(space[1], Val(N), A)

# Taylor

_preprocess_data_sequence2dfs!(::Taylor, A) = A

_preprocess_data_sequence2dfs!(::Taylor, ::Val, A) = A

# Fourier

_preprocess_data_sequence2dfs!(::Fourier, A) = A

_preprocess_data_sequence2dfs!(::Fourier, ::Val, A) = A

# Chebyshev

function _preprocess_data_sequence2dfs!(space::Chebyshev, A)
    len = length(A)
    @inbounds for i ∈ 1:order(space)
        A[len+1-i] = A[i+1]
    end
    return A
end

function _preprocess_data_sequence2dfs!(space::Chebyshev, ::Val{D}, A) where {D}
    len = size(A, D)
    @inbounds for i ∈ 1:order(space)
        selectdim(A, D, len+1-i) .= selectdim(A, D, i+1)
    end
    return A
end

# discrete Fourier series to sequence

function ifft!(c::Sequence{<:BaseSpace}, A::AbstractVector)
    space_c = space(c)
    n = length(A)
    ispow2(n) & (dimension(space_c) ≤ n) || return throw(DimensionMismatch)
    _ifft_pow2!(A)
    @inbounds coefficients(c) .= view(A, 1:dimension(space_c))
    return c
end

function rifft!(c::Sequence{<:BaseSpace}, A::AbstractVector)
    space_c = space(c)
    n = length(A)
    ispow2(n) & (dimension(space_c) ≤ n) || return throw(DimensionMismatch)
    _ifft_pow2!(A)
    @inbounds coefficients(c) .= real.(view(A, 1:dimension(space_c)))
    return c
end

function ifft!(c::Sequence{TensorSpace{T}}, A::AbstractArray{S,N}) where {N,T<:NTuple{N,BaseSpace},S}
    space_c = space(c)
    all(i -> ispow2(size(A, i)) & (dimension(space_c, i) ≤ size(A, i)), 1:N) || return throw(DimensionMismatch)
    _ifft_pow2!(A)
    C = _no_alloc_reshape(coefficients(c), dimensions(space_c))
    @inbounds C .= view(A, map(sᵢ -> 1:dimension(sᵢ), spaces(space_c))...)
    return c
end

function rifft!(c::Sequence{TensorSpace{T}}, A::AbstractArray{S,N}) where {N,T<:NTuple{N,BaseSpace},S}
    space_c = space(c)
    all(i -> ispow2(size(A, i)) & (dimension(space_c, i) ≤ size(A, i)), 1:N) || return throw(DimensionMismatch)
    _ifft_pow2!(A)
    C = _no_alloc_reshape(coefficients(c), dimensions(space_c))
    @inbounds C .= real.(view(A, map(sᵢ -> 1:dimension(sᵢ), spaces(space_c))...))
    return c
end

# FFT routines

function _fft_pow2!(a::AbstractArray{<:Complex})
    sz = size(a, 1)
    a_ = _no_alloc_reshape(a, (sz, length(a)÷sz))
    for a_row ∈ eachrow(a_)
        _fft_pow2!(a_row)
    end
    for a_col ∈ eachcol(a_)
        _fft_pow2!(a_col)
    end
    return a
end

function _ifft_pow2!(a::AbstractArray{<:Complex})
    sz = size(a, 1)
    a_ = _no_alloc_reshape(a, (sz, length(a)÷sz))
    for a_row ∈ eachrow(a_)
        _ifft_pow2!(a_row)
    end
    for a_col ∈ eachcol(a_)
        _ifft_pow2!(a_col)
    end
    return a
end

function _fft_pow2!(a::AbstractVector{Complex{T}}) where {T}
    _bitreverse!(a)
    n = length(a)
    N = 2
    while N ≤ n
        N½ = N÷2
        ω_ = cispi(-one(T)/N½)
        for k ∈ 1:N:n
            ω = one(Complex{T})
            @inbounds for j ∈ k:k+N½-1
                j′ = j + N½
                aj′_ω = a[j′] * ω
                a[j′] = a[j] - aj′_ω
                a[j] = a[j] + aj′_ω
                ω *= ω_
            end
        end
        N <<= 1
    end
    return a
end

function _fft_pow2!(a::AbstractVector{Complex{T}}) where {T<:Interval}
    _bitreverse!(a)
    π_ = convert(T, π)
    n = length(a)
    Ω = Vector{Complex{T}}(undef, n÷2)
    N = 2
    @inbounds while N ≤ n
        N½ = N÷2
        π2N⁻¹ = π_/N½
        view(Ω, 1:N½) .= cis.((-π2N⁻¹) .* (0:N½-1))
        for k ∈ 1:N:n
            @inbounds for (i, j) ∈ enumerate(k:k+N½-1)
                j′ = j + N½
                aj′_ω = a[j′] * Ω[i]
                a[j′] = a[j] - aj′_ω
                a[j] = a[j] + aj′_ω
            end
        end
        N <<= 1
    end
    return a
end

function _ifft_pow2!(a::AbstractVector{<:Complex})
    conj!(_fft_pow2!(conj!(a)))
    a ./= length(a)
    return a
end

function _bitreverse!(a::AbstractVector)
    n = length(a)
    n½ = n÷2
    j = 1
    for i ∈ 1:n-1
        if i < j
            @inbounds a[j], a[i] = a[i], a[j]
        end
        k = n½
        while 2 ≤ k < j
            j -= k
            k ÷= 2
        end
        j += k
    end
    return a
end


text:
## common functionality for  COP classical orthogonal polynomials
## These are described by 5  values: a,.b.c.d.e

abstract type AbstractCOP{T,N} <: AbstractOrthogonalPolynomial{T} end

function Base.promote_rule(P::Type{<:Polynomials.AbstractPolynomial{T}},
                           Q::Type{<:AbstractCOP{S}}) where {T,S}
    ϛ(Q){promote_type(T, S)}
end

Base.promote_rule(P::Type{<:AbstractCOP{S}},
                  Q::Type{<:Polynomials.AbstractPolynomial{T}}) where {T,S} =
                      ϛ(P){promote_type(T, S)}

function _divrem(num::P, den::Q) where {P <: AbstractCOP, Q <: AbstractCOP}

    #@assert  ⟒(P) == ⟒(Q)
    #@assert eltype(num) == eltype(den)

    p1 = convert(ϛ(P), num)
    p2 = convert(ϛ(Q), den)
    q,r = divrem(p1, p2)
    convert.(⟒(P), (q,r))

end

##
## -----
##
## interface for a  given type

basis_symbol(::Type{<:AbstractCOP}) = "P"
Polynomials.domain(::Type{<:AbstractCOP}) = Polynomials.Interval(-Inf, Inf)
classical_hypergeometric(P::Type{<:AbstractCOP}, n, x) = throw(ArgumentError("No default method"))

"""
   abcde

A named tuple returning  the  constants a,b,c,d,e  for a CCOP type with
(a⋅x²+b⋅x+c)*P₍ᵢ₊₂₎'' + (d⋅x + e) * P₍ᵢ₊₁₎ + λᵢ Pᵢ = 0.
"""
abcde(::Type{<:AbstractCOP}) = throw(ArgumentError("No default method"))

"""
    k1k0

Let  `kᵢ` be the leading  coeffiecient of the  polynomial  in  the standard basis.
This  function implements  `k₍ᵢ₊₁₎/kᵢ` for `i ≥ 0`.

The values `kᵢ` and `k₍ᵢ₊₁₎/k₍ᵢ₋₁₎` can  be generated.

The default value  leads to monic polynomials.
"""
k1k0(::Type{P},  i::Int) where {P<:AbstractCOP} =  one(eltype(P))
k1k0(::Type{P},  i::Int) where {P<:Polynomials.StandardBasisPolynomial} =  one(eltype(P))

"""
      k0(::Type{P})

This is  `basis(P,0)`,  most  often  just `1`, but may be different with some normalizations, such as `Orthonormal`.
"""
k0(::Type{P})  where{P <: AbstractCOP}  =  one(eltype(P))
k0(::Type{P}) where {P <: Polynomials.StandardBasisPolynomial} = one(eltype(P))

# Set defaults to be monic
# For non-monic subtypes of `AbstractCOP` only need k1k0  to be defined, as `kn` and `k1k_1` come for free

# kn = prod(k1k0(i) for i in 0:n-1)  *  k0(P)
kn(::Type{P},  n::Int) where{P <: AbstractCOP} = foldr(*, (k1k0(P,i) for i in 0:n-1), init=k0(P))

# k₍ᵢ₊₁₎/k₍ᵢ₋₁₎ =  (k₍ᵢ₊₁₎/kᵢ) ⋅ (kᵢ/k₍ᵢ₋₁₎)
function k1k_1(::Type{P},  i::Int) where {P<:AbstractCOP}
    @assert i > 0
    k1k0(P,i)*k1k0(P,i-1)
end

# leading term (Section 3 table)  is kn
leading_term(P::Type{<:AbstractCOP},  n::Int) =  kn(P, n)

# square root of ratio of norm2(P,n+1)/norm2(P,n)
# Let ωᵢ = √{∫ πᵢ² dw}, this is ωᵢ₊₁ ÷ ωᵢ
ω₁₀(::Type{P},n)  where {P <:  AbstractCOP} = sqrt(norm2(P,n+1)/norm2(P,n))

# Can't  change  N  here
function Base.setindex!(p::AbstractCOP{T,N}, value::Number, idx::Int) where {T, N}

    ## widen size...
    idx < 0 &&  throw(ArgumentError("Negative index"))
    val = T(value)
    d = N - 1
    if idx > d || (idx == d && iszero(value))
        throw(ArgumentError("Polynomials of type `AbstractCOP` have a  fixed size parameter, N, which  can't  be changed through  assignment. Make new polynomial  instance?"))
    end
    setindex!(p.coeffs, val,  idx+1)
end

Polynomials.degree(p::AbstractCOP{T,N})  where {T,N} = N-1
Polynomials.isconstant(p::AbstractCOP) = degree(p) <=  0

##  Evaluation

"""
    Clenshaw evaluation of an orthogonal polynomial
"""
function eval_cop(P::Type{<:AbstractCOP{T,N}}, cs, x::S) where {T,N,S}
    N == 0 && return zero(T) * zero(S)
    N == 1 && return (cs[1] * k0(P)) *  one(S)
    _eval_cop(P,cs, x)
end

function _eval_cop(P::Type{<:AbstractCOP{T,N}}, cs, x::S) where {T,N,S}
    if @generated
        quote
            Δ0 = cs[end - 1]
            Δ1 = cs[end]
            @inbounds for i in N-1:-1:2
                Δ0, Δ1 = cs[i - 1] - Δ1 * Cn(P, i-1), Δ0 + Δ1 * muladd(x, An(P,i-1), Bn(P,i-1))
            end
            p₀ = k0(P)
            p₁ =  muladd(x, An(P,0),  Bn(P,0)) * p₀
            Δ0 * p₀  + Δ1 * p₁
        end
    else
        clenshaw_eval(P, cs, x)
    end
end

#  evaluate  basis vector through hypergeometric formulation
(B::Basis{P})(x) where  {P <: AbstractCOP} = eval_basis(P, B.n, x)

# Evaluate a basis vector without realizing it and without using Clenshaw
eval_basis(::Type{P}, n, x) where {P <: AbstractCOP} = classical_hypergeometric(P, n, x)

"""
    eval_hyper(::P, cs,  x::S)

Evaluate polynomial `P(cs)` by  computing value for each basis vector from its hypergeometric representation
"""
function eval_hyper(P::Type{<:AbstractCOP}, cs, x::S) where {S}
    isempty(cs) &&  return zero(S)
    tot = cs[1] * one(S)
    for i in 2:length(cs)
        tot += cs[i] * Basis(P, i-1)(x)
    end
    tot
end


text:
@testset "Linear algebra" begin

  @testset "Dixon solve" begin
    Qx, x = FlintQQ["x"]
    K, a = number_field(x^3 + 3)
    A = rand(MatrixSpace(K, 4, 4), 10:100)
    while det(A) == 0
      A = rand(MatrixSpace(K, 4, 4), 10:100)
    end
    b = rand(MatrixSpace(K, 4 ,1), 10:100)
    @test A * Hecke.solve_dixon(A, b) == b
  end

  @testset "Pseudo matrices" begin

    Qx, x = PolynomialRing(FlintQQ, "x")

    # Compute a pseudo-hnf of a matrix over Z and check result against the HNF

    K,  a = NumberField(x - 1, "a")
    O = maximal_order(K)

    A =
    [ 154830247 724968499 442290149 700116640 313234502;
      384015369 193254623 792471850 452111534 717477652;
      816720700 46054709 339475042 389090910 103149313;
      104945689 475264799 821863415 806964485 676397437;
      281047486 709449227 590950934 18224679 295696061]

    Ahnf =
    [ 1 0 0 0 34095268312495604433062164807036095571620397;
      0 1 0 0 18160914505139254049490851505162505507397915;
      0 0 1 0 37265321283634252189185025806030875371477390;
      0 0 0 1 2276874339297612770861218322365243729516503
      0 0 0 0 37684868701591492337245802520684209569420259]

    de = fmpz(37684868701591492337245802520684209569420259)
    AoverO = MatrixSpace(O, 5, 5)(map(z -> O(z), A))

    Apm = Hecke.PseudoMatrix( AoverO, [(O(1)*O)::Hecke.NfOrdIdl for i in 1:5])

    d = numerator(det(Apm))

    Apseudohnf = Hecke.pseudo_hnf_mod(Apm, d)

    z = Apseudohnf.matrix

    for i in 1:nrows(z)
      Hecke.mul_row!(z, i, K(norm(Apseudohnf.coeffs[i])))
    end

    zinZ = MatrixSpace(FlintZZ, 5, 5)(map(zz -> numerator(coordinates(O(zz))[1]), z.entries))
    c = parent(zinZ)(Ahnf) - zinZ

    @test all([ mod(c[i,j], de) == 0 for i in 1:5, j in 1:5])

    B = Hecke.PseudoMatrix(matrix(K, [1 1; 1 1; 1 0]), [ ideal(O, K(1)), ideal(O, K(fmpq(1, 2))), ideal(O, K(1)) ])

    Bhnf = pseudo_hnf(B, :lowerleft, true)

    @test Bhnf.matrix == matrix(K, [0 0; 1 0; 1 1])

    # Construct random pseudo-matrices over different fields and check if the
    # pseudo hermite normal form span the same module

    @testset "Q[x]/x^$i - 10)" for i in 2:5
       K,  a = NumberField(x^i - 10, "a")
      O = maximal_order(K)
      #println("  Testing over field $(x^i - 10)")

      for j in 1:1
        l = rand(10:20) - i + 1
        ll = rand(1:20)
        z = rand(MatrixSpace(O, l, l), fmpz(2)^ll)
        #println("    $l x $l matrix with $ll bits")
        cc = [ (O(1)*O)::Hecke.NfOrdIdl for i in 1:l]
        pm = Hecke.PseudoMatrix(z, cc)
        d = det(pm)
        ppm = Hecke.pseudo_hnf(pm)
        @test Hecke._spans_subset_of_pseudohnf(pm, ppm)
        @test d == det(ppm)
        ppmkb, trafo = Hecke.pseudo_hnf_kb_with_transform(pm)
        @test Hecke._spans_subset_of_pseudohnf(pm, ppmkb)
        @test ppmkb.matrix == trafo*pm.matrix
      end
    end

    @testset "Field towers" begin
      f = x^2 + 36*x + 16
       K,  a = NumberField(f, "a")
      Ky, y = K["y"]
      g = y^3 - 51*y^2 + 30*y - 28
      L, b = NumberField(g, "b")

      t = rand(-1000:1000, 3, 3)
      PM = Hecke.PseudoMatrix(matrix(K, t))
      G, U = Hecke.pseudo_hnf_kb_with_transform(PM)
      @test Hecke._spans_subset_of_pseudohnf(PM, G)
      @test !iszero(det(U))
      @test G.matrix == U*PM.matrix

      G2, U2 = Hecke.pseudo_hnf_kb_with_transform(PM, :lowerleft)
      @test Hecke._spans_subset_of_pseudohnf(PM, G2, :lowerleft)
      @test !iszero(det(U2))
      @test G2.matrix == U2*PM.matrix

      Lz, z = L["z"]
      h = z^2 + 4*z + 10
      M, c = NumberField(h, "c")
      PN = Hecke.PseudoMatrix(matrix(L, t))
      H, V = Hecke.pseudo_hnf_kb_with_transform(PN)
      @test Hecke._spans_subset_of_pseudohnf(PN, H)
      @test !iszero(det(V))
      @test H.matrix == V*PN.matrix

      H2, V2 = Hecke.pseudo_hnf_kb_with_transform(PN, :lowerleft)
      @test Hecke._spans_subset_of_pseudohnf(PN, H2, :lowerleft)
      @test !iszero(det(V2))
      @test H2.matrix == V2*PN.matrix
    end

    @testset "in span" begin
       K,  a = NumberField(x^3 - 10, "a")
      O = maximal_order(K)
      ideals = []
      p = 2
      while length(ideals) < 5
        ideals = union(ideals, prime_decomposition(O, p))
        p = next_prime(p)
      end
      A = Hecke.PseudoMatrix(one(MatrixSpace(O, 5, 5)), [ p for (p, e) in ideals ])
      v = [ K(rand(p, 100)) for (p, e) in ideals ]
      @test Hecke._in_span(v, A)[1]

       K,  a = NumberField(x, "a")
      O = maximal_order(K)
      A = Hecke.PseudoMatrix(matrix(O, map(O, [ 1 2 3 4; 0 7 8 9; 0 0 11 12; 0 0 0 13 ])), [ O(1)*O for i = 1:4 ])
      @test Hecke._in_span(map(K, [1, 2, 3, 4]), A)[1]
      @test Hecke._in_span(map(K, [5, 6, 7, 8]), A)[1] == false
    end
  end
end


text:
@defcomp population begin
    regions             = Index()

    population          = Variable(index=[time,regions])
    populationin1       = Variable(index=[time,regions])
    globalpopulation    = Variable(index=[time])

    pgrowth             = Parameter(index=[time,regions])
    enter               = Parameter(index=[time,regions])
    leave               = Parameter(index=[time,regions])
    dead                = Parameter(index=[time,regions])
    pop0                = Parameter(index=[regions])
    runwithoutpopulationperturbation = Parameter{Bool}(default = false)

    function run_timestep(p, v, d, t)
        if is_first(t)
            for r in d.regions
                v.population[t, r] = p.pop0[r]
                v.populationin1[t, r] = v.population[t, r] * 1000000.0
            end

            v.globalpopulation[t] = sum(v.populationin1[t,:])
        else
            for r in d.regions
                v.population[t, r] = (1.0 + 0.01 * p.pgrowth[t - 1, r]) * (v.population[t - 1, r] + ((gettime(t) >= 1990 && !p.runwithoutpopulationperturbation) ? (p.enter[t - 1, r] / 1000000.0) - (p.leave[t - 1, r] / 1000000.0) - (p.dead[t - 1, r] >= 0 ? p.dead[t - 1, r] / 1000000.0 : 0) : 0))

                if v.population[t, r] < 0
                    v.population[t, r] = 0.000001
                end

                v.populationin1[t, r] = v.population[t, r] * 1000000.0
            end

            v.globalpopulation[t] = sum(v.populationin1[t,:])
        end
    end
end


text:
################################################################################
#
#  Basic field access
#
################################################################################

base_ring(A::AlgAss{T}) where {T} = A.base_ring::parent_type(T)

Generic.dim(A::AlgAss) = size(A.mult_table, 1)

elem_type(::Type{AlgAss{T}}) where {T} = AlgAssElem{T}

parent(::Type{AlgAssElem{T}}) where {T} = AlgAss{T}

################################################################################
#
#  Test of Commutativity
#
################################################################################

iscommutative_known(A::AlgAss{T}) where {T} = A.iscommutative != 0

function iscommutative(A::AlgAss{T}) where {T}
  if iscommutative_known(A)
    return A.iscommutative == 1
  end
  for i = 1:dim(A)
    for j = i + 1:dim(A)
      if A.mult_table[i, j, :] != A.mult_table[j, i, :]
        A.iscommutative = 2
        return false
      end
    end
  end
  A.iscommutative = 1
  return true
end

################################################################################
#
#  Associativity, Distributivity test
#
################################################################################

function check_associativity(A::AlgAss)

  for i=1:dim(A)
    for j=1:dim(A)
      for k=1:dim(A)
        @assert (A[i]*A[j])*A[k]==A[i]*(A[j]*A[k])
      end
    end
  end
  return true

end

function check_distributivity(A::AlgAss)

  for i=1:dim(A)
    for j=1:dim(A)
      for k=1:dim(A)
        @assert A[i]*(A[j]+A[k])==A[i]*A[j]+A[i]*A[k]
      end
    end
  end
  return true

end

################################################################################
#
#  Construction
#
################################################################################

# This only works if base_ring(A) is a field (probably)
function find_one(A::AlgAss)
  n = dim(A)
  M = zero_matrix(base_ring(A), n^2, n)
  c = zero_matrix(base_ring(A), n^2, 1)
  for k = 1:n
    kn = (k - 1)*n
    c[kn + k, 1] = base_ring(A)(1)
    for i = 1:n
      for j = 1:n
        M[i + kn, j] = deepcopy(A.mult_table[j, k, i])
      end
    end
  end
  Mc = hcat(M, c)
  rref!(Mc)
  @assert !iszero(Mc[n, n])
  n != 1 && @assert iszero(Mc[n + 1, n + 1])
  cc = solve_ut(sub(Mc, 1:n, 1:n), sub(Mc, 1:n, (n + 1):(n + 1)))
  one = [ cc[i, 1] for i = 1:n ]
  return one
end

function AlgAss(R::Ring, mult_table::Array{T, 3}, one::Array{T, 1}) where {T}
  # checks
  return AlgAss{T}(R, mult_table, one)
end

function AlgAss(R::Ring, mult_table::Array{T, 3}) where {T}
  # checks
  A = AlgAss{T}(R)
  A.mult_table = mult_table
  A.one = find_one(A)
  return A
end

function AlgAss(R::Ring, d::Int, arr::Array{T, 1}) where {T}
  mult_table = Array{T, 3}(d, d, d)
  n = d^2
  for i in 1:d
    for j in 1:d
      for k in 1:d
        mult_table[i, j, k] = arr[(i - 1) * n + (j - 1) * d + k]
      end
    end
  end
  return AlgAss(R, mult_table)
end

# Constructs the algebra R[X]/f
function AlgAss(f::PolyElem)
  R = base_ring(parent(f))
  n = degree(f)
  Rx = parent(f)
  x = gen(Rx)
  B = Array{elem_type(Rx), 1}(2*n - 1)
  B[1] = Rx(1)
  for i = 2:2*n - 1
    B[i] = mod(B[i - 1]*x, f)
  end
  mult_table = Array{elem_type(R), 3}(n, n, n)
  for i = 1:n
    for j = i:n
      for k = 1:n
        mult_table[i, j, k] = coeff(B[i + j - 1], k - 1)
        mult_table[j, i, k] = coeff(B[i + j - 1], k - 1)
      end
    end
  end
  one = map(R, zeros(Int, n))
  one[1] = R(1)
  A = AlgAss(R, mult_table, one)
  A.iscommutative = 1
  return A
end

function AlgAss(O::NfOrd, I::NfOrdIdl, p::Union{Integer, fmpz})
  @assert order(I) == O

  n = degree(O)
  BO = Hecke.basis(O)

  pisfmpz = (p isa fmpz)
  Fp = ResidueRing(FlintZZ, p, cached=false)
  BOmod = [ mod(v, I) for v in BO ]
  B = zero_matrix(Fp, n, n)
  for i = 1:n
    b = elem_in_basis(BOmod[i])
    for j = 1:n
      B[i, j] = Fp(b[j])
    end
  end
  if pisfmpz
    r, B = _rref(B)
  else
    r = rref!(B)
  end
  r == 0 && error("Cannot construct zero dimensional algebra.")
  b = Vector{fmpz}(n)
  basis = Vector{NfOrdElem}(r)
  for i = 1:r
    for j = 1:n
      b[j] = fmpz(B[i, j])
    end
    basis[i] = O(b)
  end

  if pisfmpz
    _, p, L, U = _lufact(transpose(B))
  else
    _, p, L, U = lufact(transpose(B))
  end
  mult_table = Array{elem_type(Fp), 3}(r, r, r)
  d = zero_matrix(Fp, n, 1)
  for i = 1:r
    for j = i:r
      c = elem_in_basis(mod(basis[i]*basis[j], I))
      for k = 1:n
        d[p[k], 1] = c[k]
      end
      d = solve_lt(L, d)
      d = solve_ut(U, d)
      for k = 1:r
        mult_table[i, j, k] = deepcopy(d[k, 1])
        mult_table[j, i, k] = deepcopy(d[k, 1])
      end
    end
  end

  if isone(basis[1])
    one = zeros(Fp, r)
    one[1] = Fp(1)
    A = AlgAss(Fp, mult_table, one)
  else
    A = AlgAss(Fp, mult_table)
  end
  A.iscommutative = 1

  function _image(a::NfOrdElem)
    c = elem_in_basis(mod(a, I))
    for k = 1:n
      d[p[k], 1] = c[k]
    end
    d = solve_lt(L, d)
    d = solve_ut(U, d)
    e = A()
    for k = 1:r
      e.coeffs[k] = deepcopy(d[k, 1])
    end
    return e
  end

  function _preimage(a::AlgAssElem)
    return sum(fmpz(a.coeffs[i])*basis[i] for i = 1:r)
  end

  OtoA = NfOrdToAlgAssMor{elem_type(Fp)}(O, A, _image, _preimage)

  return A, OtoA
end

function _modular_basis(pb::Vector{Tuple{T, NfOrdFracIdl}}, p::NfOrdIdl) where T <: RelativeElement{nf_elem}
  L = parent(pb[1][1])
  K = base_ring(L)
  basis = Array{elem_type(L), 1}()
  u = L(K(uniformizer(p)))
  for i = 1:degree(L)
    v = valuation(pb[i][2], p)
    push!(basis, (u^v)*pb[i][1])
  end
  return basis
end

#=
Qx, x = QQ["x"];
f = x^2 + 12*x - 92;
K, a = number_field(f, "a");
OK = maximal_order(K);
Ky, y = K["y"];
g = y^2 - 54*y - 73;
L, b = number_field(g, "b");
OL = maximal_order(L);
p = prime_decomposition(OK, 2)[1][1]
=#

# Assume that O is relative order over OK, I is an ideal of O and p is a prime
# ideal of OK with pO \subseteq I. O/I is an OK/p-algebra.
#
# The idea is to compute pseudo-basis of O and I respectively, for which the
# coefficient ideals have zero p-adic valuation. Then we can think in the
# localization at p and do as in the case of principal ideal domains.
function AlgAss(O::NfRelOrd{T, S}, I::NfRelOrdIdl{T, S}, p::Union{NfOrdIdl, NfRelOrdIdl}) where {T, S}
  basis_pmatI = basis_pmat(I, Val{false})
  basis_pmatO = basis_pmat(O, Val{false})

  new_basis_mat = deepcopy(O.basis_mat)
  new_basis_mat_I = deepcopy(I.basis_mat)

  pi = anti_uniformizer(p)

  new_basis_coeffs = S[]

  for i in 1:degree(O)
    a = pi^valuation(basis_pmat(O).coeffs[i], p)
    push!(new_basis_coeffs, a * basis_pmatO.coeffs[i])
    mul_row!(new_basis_mat, i, inv(a))
    for j in 1:degree(O)
      new_basis_mat_I[j, i] = new_basis_mat_I[j, i] * a
    end
  end

  new_coeff_I = S[]

  for i in 1:degree(O)
    a = pi^valuation(basis_pmatI.coeffs[i], p)
    push!(new_coeff_I, a * basis_pmatI.coeffs[i])
    mul_row!(new_basis_mat_I, i, inv(a))
  end

  Fp, mF = ResidueField(order(p), p)
  mmF = extend(mF, base_ring(nf(O)))
  invmmF = inv(mmF)

  basis_elts = Int[]
  reducers = Int[]

  for i in 1:degree(O)
    v = valuation(new_basis_mat_I[i, i], p)
    v2 = valuation(new_coeff_I[i], p)
    #@show (v2, v)
    @assert v >= 0
    if v == 0
    #if valuation(basis_pmatI.coeffs[i], p) + valuation(new_basis_mat_I[i, i], p) == 0
      push!(reducers, i)
    else
      push!(basis_elts, i)
    end
  end

  reverse!(reducers)

  OLL = Order(nf(O), PseudoMatrix(new_basis_mat, new_basis_coeffs))

  newI = ideal(OLL, PseudoMatrix(new_basis_mat_I, new_coeff_I))

  new_basis = pseudo_basis(OLL)

  pseudo_basis_newI = pseudo_basis(newI)

  tmp_matrix = zero_matrix(base_ring(nf(O)), 1, degree(O))

  basis_mat_inv_OLL = basis_mat_inv(OLL)

  function _coeff(c)
    for i in 0:degree(O) - 1
      tmp_matrix[1, i + 1] = coeff(c, i)
    end
    return tmp_matrix * basis_mat_inv_OLL
  end

  r = length(basis_elts)

  mult_table = Array{elem_type(Fp), 3}(r, r, r)

  for i in 1:r
    for j in 1:r
      c = new_basis[basis_elts[i]][1] * new_basis[basis_elts[j]][1]
      coeffs = _coeff(c)

      for k in reducers
        d = -coeffs[k]//new_basis_mat_I[k, k]
        c = c + d * pseudo_basis_newI[k][1]
      end
      coeffs = _coeff(c)
      for k in 1:degree(O)
        if !(k in basis_elts)
          @assert iszero(coeffs[k])
        end
      end
      for k in 1:r
        mult_table[i, j, k] = mmF(coeffs[basis_elts[k]])
      end
    end
  end

  if isone(new_basis[basis_elts[1]][1])
    one = zeros(Fp, length(basis_elts))
    one[1] = Fp(1)
    A = AlgAss(Fp, mult_table, one)
  else
    A = AlgAss(Fp, mult_table)
  end
  A.iscommutative = 1

  function _image(a::NfRelOrdElem)
    c = a.elem_in_nf
    coeffs = _coeff(c)
    for k in reducers
      d = -coeffs[k]//new_basis_mat_I[k, k]
      c = c + d*pseudo_basis_newI[k][1]
    end
    coeffs = _coeff(c)
    for k in 1:degree(O)
      if !(k in basis_elts)
        @assert iszero(coeffs[k])
      end
    end
    b = A()
    for k in 1:r
      b.coeffs[k] = mmF(coeffs[basis_elts[k]])
    end
    return b
  end

  lifted_basis_of_A = []

  for i in basis_elts
    c = coprime_to(new_basis[i][2], p)
    b = invmmF(inv(mmF(c)))*c*new_basis[i][1]
    @assert b in O
    push!(lifted_basis_of_A, b)
  end

  function _preimage(v::AlgAssElem)
    return O(sum((invmmF(v.coeffs[i])) * lifted_basis_of_A[i] for i in 1:r))
  end

  OtoA = NfRelOrdToAlgAssMor{T, S, elem_type(Fp)}(O, A, _image, _preimage)

  return A, OtoA
end

#Given I with v_p(I) = 0, returns element a \in I such that v_p(a) = 0
function coprime_to(I::NfOrdFracIdl, p::NfOrdIdl)
  pi = anti_uniformizer(p)
  a = basis(I)[1]
  l = valuation(a, p)
  @assert l >= 0
  if l > 0
    a = pi^l * a
  end
  @assert valuation(a, p) == 0
  @assert denominator(I)*a in numerator(I)
  return a
end

function coprime_to(I::NfRelOrdFracIdl, p::NfRelOrdIdl)
  pi = anti_uniformizer(p)
  a = rand(I, 500)
  l = valuation(a, p)
  @assert l >= 0
  if l > 0
    a = pi^l*a
  end
  @assert valuation(a, p) == 0
  return a
end

################################################################################
#
#  String I/O
#
################################################################################

function show(io::IO, A::AlgAss)
  print(io, "Associative algebra of dimension $(dim(A)) over ")
  print(io, A.base_ring)
end

################################################################################
#
#  Deepcopy
#
################################################################################

function Base.deepcopy_internal(A::AlgAss{T}, dict::ObjectIdDict) where {T}
  B = AlgAss{T}(base_ring(A))
  for x in fieldnames(A)
    if x != :base_ring && isdefined(A, x)
      setfield!(B, x, Base.deepcopy_internal(getfield(A, x), dict))
    end
  end
  B.base_ring = A.base_ring
  return B
end

################################################################################
#
#  Equality
#
################################################################################

function ==(A::AlgAss{T}, B::AlgAss{T}) where {T}
  base_ring(A) != base_ring(B) && return false
  return A.one == B.one && A.mult_table == B.mult_table
end

################################################################################
#
#  Subalgebras
#
################################################################################

# This only works if base_ring(A) is a field
# Constructs the algebra e*A
function subalgebra(A::AlgAss{T}, e::AlgAssElem{T}, idempotent::Bool = false) where {T}
  @assert parent(e) == A
  R = base_ring(A)
  isgenres = (typeof(R) <: Generic.ResRing)
  n = dim(A)
  B = representation_matrix(e)
  if isgenres
    r, B = _rref(B)
  else
    r = rref!(B)
  end
  r == 0 && error("Cannot construct zero dimensional algebra.")
  basis = Vector{AlgAssElem{T}}(r)
  for i = 1:r
    basis[i] = elem_from_mat_row(A, B, i)
  end

  # The basis matrix of e*A with respect to A is
  basis_mat_of_eA = sub(B, 1:r, 1:n)

  if isgenres
    _, p, L, U = _lufact(transpose(B))
  else
    _, p, L, U = lufact(transpose(B))
  end
  mult_table = Array{elem_type(R), 3}(r, r, r)
  c = A()
  d = zero_matrix(R, n, 1)
  for i = 1:r
    for j = 1:r
      if iscommutative(A) && j < i
        continue
      end
      c = mul!(c, basis[i], basis[j])
      for k = 1:n
        d[p[k], 1] = c.coeffs[k]
      end
      d = solve_lt(L, d)
      d = solve_ut(U, d)
      for k = 1:r
        mult_table[i, j, k] = deepcopy(d[k, 1])
        if iscommutative(A) && i != j
          mult_table[j, i, k] = deepcopy(d[k, 1])
        end
      end
    end
  end
  if idempotent
    for k = 1:n
      d[p[k], 1] = e.coeffs[k]
    end
    d = solve_lt(L, d)
    d = solve_ut(U, d)
    v = Vector{elem_type(R)}(r)
    for i in 1:r
      v[i] = d[i, 1]
    end
    eA = AlgAss(R, mult_table, v)
  else
    eA = AlgAss(R, mult_table)
  end

  # TODO: The following is wrong. The algebra eA may be commutative
  # even if A is not commutative.
  eA.iscommutative = A.iscommutative

  if idempotent
    # We have the map eA -> A, given by the multiplying with basis_mat_of_eA.
    # But there is also the canonical projection A -> eA, a -> ea.
    # We compute the corresponding matrix.
    B = representation_matrix(e)
    C = zero_matrix(R, n, r)
    for i in 1:n
      for k = 1:n
        d[p[k], 1] = B[i, k]
      end
      d = solve_lt(L, d)
      d = solve_ut(U, d)
      for k in 1:r
        C[i, k] = d[k, 1]
      end
    end
    eAtoA = AlgAssMor(eA, A, basis_mat_of_eA, C)
  else
    eAtoA = AlgAssMor(eA, A, basis_mat_of_eA)
  end
  return eA, eAtoA
end

function subalgebra(A::AlgAss{T}, basis::Array{AlgAssElem{T},1}) where {T}
  B=zero_matrix(base_ring(A), dim(A), length(basis))
  for i=1:dim(A)
    for j=1:length(basis)
      B[i,j]=basis[j].coeffs[i]
    end
  end
  M=Array{elem_type(base_ring(A)),3}(length(basis), length(basis), length(basis))
  for i=1:length(basis)
    for j=1:length(basis)
      x=basis[i]*basis[j]
      N1=matrix(base_ring(A), dim(A), 1, x.coeffs)
      N=_solve_unique(N1,B)
      for k=1:length(basis)
        M[i,j,k]=N[k,1]
      end
    end
  end
  A1=AlgAss(base_ring(A), M)
  return A1, AlgAssMor(A1, A, B')
end

################################################################################
#
#  Splitting
#
################################################################################

function kernel_of_frobenius(A::AlgAss)
  F = base_ring(A)
  q = order(F)

  b = A()
  c = A()
  B = zero_matrix(F, dim(A), dim(A))
  for i = 1:dim(A)
    b.coeffs[i] = F(1)
    if i > 1
      b.coeffs[i - 1] = F()
    end
    c = b^q - b
    for j = 1:dim(A)
      B[j, i] = deepcopy(c.coeffs[j])
    end
  end

  V = right_kernel(B)
  return [ A(v) for v in V ]
end

function issimple(A::AlgAss, compute_algebras::Type{Val{T}} = Val{true}) where T
  if dim(A) == 1
    if compute_algebras == Val{true}
      return true, [ (A, AlgAssMor(A, A, identity_matrix(base_ring(A), dim(A)), identity_matrix(base_ring(A), dim(A)))) ]
    else
      return true
    end
  end

  F = base_ring(A)
  @assert !iszero(characteristic(F))
  V = kernel_of_frobenius(A)
  k = length(V)

  if compute_algebras == Val{false}
    return k == 1
  end

  if k == 1
    return true, [ (A, AlgAssMor(A, A, identity_matrix(base_ring(A), dim(A)), identity_matrix(base_ring(A), dim(A)))) ]
  end

  while true
    c = elem_type(F)[ rand(F) for i = 1:k ]
    a = dot(c, V)
    f = minpoly(a)

    if degree(f) < 2
      continue
    end

    @assert issquarefree(f)

    fac = factor(f)
    R = parent(f)
    factors = Vector{elem_type(R)}()
    for ff in keys(fac.fac)
      push!(factors, ff)
    end
    sols = Vector{elem_type(R)}()
    right_side = [ R() for i = 1:length(factors) ]
    max_deg = 0
    for i = 1:length(factors)
      right_side[i] = R(1)
      if i != 1
        right_side[i - 1] = R(0)
      end
      s = crt(right_side, factors)
      push!(sols, s)
      max_deg = max(max_deg, degree(s))
    end
    x = one(A)
    powers = Vector{elem_type(A)}()
    for i = 1:max_deg + 1
      push!(powers, x)
      x *= a
    end
    idems = Vector{elem_type(A)}()
    for s in sols
      idem = A()
      for i = 0:degree(s)
        idem += coeff(s, i)*powers[i + 1]
      end
      push!(idems, idem)
    end

    S = [ (subalgebra(A, idem, true)...) for idem in idems ]
    return false, S
  end
end

doc"""
***
    split(A::AlgAss)

> Given a commutative algebra over a finite field of prime order, this function
> returns a decomposition of A as a direct sum of fields.
"""

function split(A::AlgAss)
  b, algebras = issimple(A)
  if b
    return algebras
  end
  result = typeof(algebras)()
  while length(algebras) != 0
    B, BtoA = pop!(algebras)
    b, algebras2 = issimple(B)
    if b
      push!(result, (B, BtoA))
    else
      for (C, CtoB) in algebras2
        CtoA = compose_and_squash(BtoA, CtoB)
        push!(algebras, (C, CtoA))
      end
    end
  end
  return result
end

################################################################################
#
#  Radical
#
################################################################################

function trace_matrix(A::AlgAss)
  F = base_ring(A)
  n = dim(A)
  M = zero_matrix(F, n, n)
  for i = 1:n
    M[i,i] = trace(A[i]^2)
  end
  for i = 1:n
    for j = i+1:n
      x = trace(A[i]*A[j])
      M[i,j] = x
      M[j,i] = x
    end
  end
  return M

end

doc"""
***
    radical(A::AlgAss{fq_nmod})

> Given an algebra over a finite field of prime order, this function
> returns a set of elements generating the radical of A
"""

function radical(A::AlgAss{fq_nmod})

  F=A.base_ring
  p=F.p
  l=clog(fmpz(dim(A)),p)
  #First step: kernel of the trace matrix
  I=trace_matrix(A)
  k,B=nullspace(I)
  # The columns of B give the coordinates of the elements in the order.
  if k==0
    return AlgAssElem[]
  end
  C=transpose(B)
  if l==1 && dim(A)!=p
    #In this case, we can output I: it is the standard p-trace method.
    return AlgAssElem[elem_from_mat_row(A,C,i) for i=1:rows(C)]
  end
  #Now, iterate: we need to find the kernel of tr((xy)^(p^i))/p^i mod p
  #on the subspace generated by C
  #Hard to believe, but this is linear!!!!
  for i=1:l
    M=zero_matrix(F, dim(A), rows(C))
    for t=1:rows(I)
      elm=elem_from_mat_row(A,C,t)
      for s=1:dim(A)
        a=elm*A[s]
        M1=representation_matrix(a^(p^i))
        el=sum(FlintZZ(coeff(M1[k,k],0)) for k=1:dim(A))
        M[s,t]=F(divexact(el,p^i))
      end
    end
    k,B=nullspace(M)
    if k==0
      return AlgAssOrdIdl(O,MatrixSpace(FlintZZ, dim(A), dim(A), false)(p))
    end
    C=transpose(B)*C
  end
  return AlgAssElem[elem_from_mat_row(A,C,i) for i=1:rows(C)]

end

################################################################################
#
#  Wedderburn decomposition
#
################################################################################

#Given a semisimple algebra A over F_p, we give back the idempotents of a webberburn decomposition

doc"""
***
    wedderburn_decomposition(A::AlgAss{fq_nmod})

> Given a semisimple algebra over a finite field of prime order, this function
> returns a decomposition of A as a direct sum of simple algebras.
"""
function wedderburn_decomposition(A::AlgAss{fq_nmod})

  ZA,mZA=center(A)
  Algs=split(ZA)
  res=Array{Tuple{AlgAss, AlgAssMor},1}(length(Algs))
  i=1
  for (B, BtoZA) in Algs
    x=mZA(BtoZA(one(B)))
    res[i]= subalgebra(A,x)
    i+=1
  end
  return res

end

################################################################################
#
#  Random elements
#
################################################################################

function rand(A::AlgAss{T}) where T
  c = [ rand(base_ring(A)) for i = 1:dim(A) ]
  return AlgAssElem{T}(A, c)
end

################################################################################
#
#  Primitive elements
#
################################################################################

function primitive_element(A::AlgAss)
  a, _ = _primitive_element(A::AlgAss)
  return a
end

function _primitive_element(A::AlgAss)
  error("Not implemented yet")
  return nothing
end

function _primitive_element(A::AlgAss{T}) where T <: Union{nmod, fq, fq_nmod, Generic.Res{fmpz}}
  d = dim(A)
  a = rand(A)
  f = minpoly(a)
  while degree(f) < d
    a = rand(A)
    f = minpoly(a)
  end
  return a, f
end

function _as_field(A::AlgAss{T}) where T
  d = dim(A)
  a, mina = _primitive_element(A)
  b = one(A)
  M = zero_matrix(base_ring(A), d, d)
  elem_to_mat_row!(M, 1, b)
  for i in 1:(d - 1)
    b = mul!(b, b, a)
    elem_to_mat_row!(M, i + 1, b)
  end
  if T == Generic.Res{fmpz}
    A, c = inv(M)
    B = divexact(A, c)
  elseif T == nmod
    B = inv(M)
  end
  N = zero_matrix(base_ring(A), 1, d)
  local f
  let N = N, B = B
    f = function(x)
      for i in 1:d
        N[1, i] = x.coeffs[i]
      end
      return N * B
    end
  end
  return a, mina, f
end