text:
struct Field
    name::Symbol
    interpolation::Interpolation
    dim::Int
end

abstract type AbstractDofHandler end

"""
    DofHandler(grid::Grid)

Construct a `DofHandler` based on the grid `grid`.
"""
struct DofHandler{dim,C,T} <: AbstractDofHandler
    field_names::Vector{Symbol}
    field_dims::Vector{Int}
    # TODO: field_interpolations can probably be better typed: We should at least require
    #       all the interpolations to have the same dimension and reference shape
    field_interpolations::Vector{Interpolation}
    bc_values::Vector{BCValues{T}} # TODO: BcValues is created/handeld by the constrainthandler, so this can be removed
    cell_dofs::Vector{Int}
    cell_dofs_offset::Vector{Int}
    closed::ScalarWrapper{Bool}
    grid::Grid{dim,C,T}
    ndofs::ScalarWrapper{Int}
end

function DofHandler(grid::Grid)
    isconcretetype(getcelltype(grid)) || error("Grid includes different celltypes. Use MixedDofHandler instead of DofHandler")
    DofHandler(Symbol[], Int[], Interpolation[], BCValues{Float64}[], Int[], Int[], ScalarWrapper(false), grid, Ferrite.ScalarWrapper(-1))
end

function Base.show(io::IO, ::MIME"text/plain", dh::DofHandler)
    println(io, "DofHandler")
    println(io, "  Fields:")
    for i in 1:nfields(dh)
        println(io, "    ", repr(dh.field_names[i]), ", interpolation: ", dh.field_interpolations[i],", dim: ", dh.field_dims[i])
    end
    if !isclosed(dh)
        print(io, "  Not closed!")
    else
        println(io, "  Dofs per cell: ", ndofs_per_cell(dh))
        print(io, "  Total dofs: ", ndofs(dh))
    end
end

ndofs(dh::AbstractDofHandler) = dh.ndofs[]
ndofs_per_cell(dh::AbstractDofHandler, cell::Int=1) = dh.cell_dofs_offset[cell+1] - dh.cell_dofs_offset[cell]
isclosed(dh::AbstractDofHandler) = dh.closed[]
nfields(dh::AbstractDofHandler) = length(dh.field_names)
getfieldnames(dh::AbstractDofHandler) = dh.field_names
ndim(dh::AbstractDofHandler, field_name::Symbol) = dh.field_dims[find_field(dh, field_name)]
function find_field(dh::DofHandler, field_name::Symbol)
    j = findfirst(i->i == field_name, dh.field_names)
    j == 0 && error("did not find field $field_name")
    return j
end

# Calculate the offset to the first local dof of a field
function field_offset(dh::DofHandler, field_name::Symbol)
    offset = 0
    for i in 1:find_field(dh, field_name)-1
        offset += getnbasefunctions(dh.field_interpolations[i])::Int * dh.field_dims[i]
    end
    return offset
end

getfieldinterpolation(dh::DofHandler, field_idx::Int) = dh.field_interpolations[field_idx]
getfielddim(dh::DofHandler, field_idx::Int) = dh.field_dims[field_idx]
getbcvalue(dh::DofHandler, field_idx::Int) = dh.bc_values[field_idx]

function getfielddim(dh::DofHandler, name::Symbol)
    field_pos = findfirst(i->i == name, getfieldnames(dh))
    field_pos === nothing && error("did not find field $name")
    return dh.field_dims[field_pos]
end

"""
    dof_range(dh:DofHandler, field_name)

Return the local dof range for `field_name`. Example:

```jldoctest
julia> grid = generate_grid(Triangle, (3, 3));

julia> dh = DofHandler(grid); push!(dh, :u, 3); push!(dh, :p, 1); close!(dh);

julia> dof_range(dh, :u)
1:9

julia> dof_range(dh, :p)
10:12
```
"""
function dof_range(dh::DofHandler, field_name::Symbol)
    f = find_field(dh, field_name)
    offset = field_offset(dh, field_name)
    n_field_dofs = getnbasefunctions(dh.field_interpolations[f])::Int * dh.field_dims[f]
    return (offset+1):(offset+n_field_dofs)
end

function Base.push!(dh::DofHandler, name::Symbol, dim::Int, ip::Interpolation=default_interpolation(getcelltype(dh.grid)))
    @assert !isclosed(dh)
    @assert !in(name, dh.field_names)
    push!(dh.field_names, name)
    push!(dh.field_dims, dim)
    push!(dh.field_interpolations, ip)
    push!(dh.bc_values, BCValues(ip, default_interpolation(getcelltype(dh.grid))))
    return dh
end

# sort and return true (was already sorted) or false (if we had to sort)
function sortedge(edge::Tuple{Int,Int})
    a, b = edge
    a < b ? (return (edge, true)) : (return ((b, a), false))
end

sortface(face::Tuple{Int,Int}) = minmax(face[1], face[2])
function sortface(face::Tuple{Int,Int,Int})
    a, b, c = face
    b, c = minmax(b, c)
    a, c = minmax(a, c)
    a, b = minmax(a, b)
    return (a, b, c)
end

function sortface(face::Tuple{Int,Int,Int,Int})
    a, b, c, d = face
    c, d = minmax(c, d)
    b, d = minmax(b, d)
    a, d = minmax(a, d)
    b, c = minmax(b, c)
    a, c = minmax(a, c)
    a, b = minmax(a, b)
    return (a, b, c)
end

function close!(dh::DofHandler)
    dh, _, _, _ = __close!(dh)
    return dh
end

# close the DofHandler and distribute all the dofs
function __close!(dh::DofHandler{dim}) where {dim}
    @assert !isclosed(dh)

    # `vertexdict` keeps track of the visited vertices. We store the global vertex
    # number and the first dof we added to that vertex.
    vertexdicts = [Dict{Int,Int}() for _ in 1:nfields(dh)]

    # `edgedict` keeps track of the visited edges, this will only be used for a 3D problem
    # An edge is determined from two vertices, but we also need to store the direction
    # of the first edge we encounter and add dofs too. When we encounter the same edge
    # the next time we check if the direction is the same, otherwise we reuse the dofs
    # in the reverse order
    edgedicts = [Dict{Tuple{Int,Int},Tuple{Int,Bool}}() for _ in 1:nfields(dh)]

    # `facedict` keeps track of the visited faces. We only need to store the first dof we
    # added to the face; if we encounter the same face again we *always* reverse the order
    # In 2D a face (i.e. a line) is uniquely determined by 2 vertices, and in 3D a
    # face (i.e. a surface) is uniquely determined by 3 vertices.
    facedicts = [Dict{NTuple{dim,Int},Int}() for _ in 1:nfields(dh)]

    # celldofs are never shared between different cells so there is no need
    # for a `celldict` to keep track of which cells we have added dofs too.

    # We create the `InterpolationInfo` structs with precomputed information for each
    # interpolation since that allows having the cell loop as the outermost loop,
    # and the interpolation loop inside without using a function barrier
    interpolation_infos = InterpolationInfo[]
    for interpolation in dh.field_interpolations
        # push!(dh.interpolation_info, InterpolationInfo(interpolation))
        push!(interpolation_infos, InterpolationInfo(interpolation))
    end

    # not implemented yet: more than one facedof per face in 3D
    dim == 3 && @assert(!any(x->x.nfacedofs > 1, interpolation_infos))

    nextdof = 1 # next free dof to distribute
    push!(dh.cell_dofs_offset, 1) # dofs for the first cell start at 1

    # loop over all the cells, and distribute dofs for all the fields
    for (ci, cell) in enumerate(getcells(dh.grid))
        @debug println("cell #$ci")
        for fi in 1:nfields(dh)
            interpolation_info = interpolation_infos[fi]
            @debug println("  field: $(dh.field_names[fi])")
            if interpolation_info.nvertexdofs > 0
                for vertex in vertices(cell)
                    @debug println("    vertex#$vertex")
                    token = Base.ht_keyindex2!(vertexdicts[fi], vertex)
                    if token > 0 # haskey(vertexdicts[fi], vertex) # reuse dofs
                        reuse_dof = vertexdicts[fi].vals[token] # vertexdicts[fi][vertex]
                        for d in 1:dh.field_dims[fi]
                            @debug println("      reusing dof #$(reuse_dof + (d-1))")
                            push!(dh.cell_dofs, reuse_dof + (d-1))
                        end
                    else # token <= 0, distribute new dofs
                        for vertexdof in 1:interpolation_info.nvertexdofs
                            Base._setindex!(vertexdicts[fi], nextdof, vertex, -token) # vertexdicts[fi][vertex] = nextdof
                            for d in 1:dh.field_dims[fi]
                                @debug println("      adding dof#$nextdof")
                                push!(dh.cell_dofs, nextdof)
                                nextdof += 1
                            end
                        end
                    end
                end # vertex loop
            end
            if dim == 3 # edges only in 3D
                if interpolation_info.nedgedofs > 0
                    for edge in edges(cell)
                        sedge, dir = sortedge(edge)
                        @debug println("    edge#$sedge dir: $(dir)")
                        token = Base.ht_keyindex2!(edgedicts[fi], sedge)
                        if token > 0 # haskey(edgedicts[fi], sedge), reuse dofs
                            startdof, olddir = edgedicts[fi].vals[token] # edgedicts[fi][sedge] # first dof for this edge (if dir == true)
                            for edgedof in (dir == olddir ? (1:interpolation_info.nedgedofs) : (interpolation_info.nedgedofs:-1:1))
                                for d in 1:dh.field_dims[fi]
                                    reuse_dof = startdof + (d-1) + (edgedof-1)*dh.field_dims[fi]
                                    @debug println("      reusing dof#$(reuse_dof)")
                                    push!(dh.cell_dofs, reuse_dof)
                                end
                            end
                        else # token <= 0, distribute new dofs
                            Base._setindex!(edgedicts[fi], (nextdof, dir), sedge, -token) # edgedicts[fi][sedge] = (nextdof, dir),  store only the first dof for the edge
                            for edgedof in 1:interpolation_info.nedgedofs
                                for d in 1:dh.field_dims[fi]
                                    @debug println("      adding dof#$nextdof")
                                    push!(dh.cell_dofs, nextdof)
                                    nextdof += 1
                                end
                            end
                        end
                    end # edge loop
                end
            end
            if interpolation_info.nfacedofs > 0 && (interpolation_info.dim == dim)
                for face in faces(cell)
                    sface = sortface(face) # TODO: faces(cell) may as well just return the sorted list
                    @debug println("    face#$sface")
                    token = Base.ht_keyindex2!(facedicts[fi], sface)
                    if token > 0 # haskey(facedicts[fi], sface), reuse dofs
                        startdof = facedicts[fi].vals[token] # facedicts[fi][sface]
                        for facedof in interpolation_info.nfacedofs:-1:1 # always reverse (YOLO)
                            for d in 1:dh.field_dims[fi]
                                reuse_dof = startdof + (d-1) + (facedof-1)*dh.field_dims[fi]
                                @debug println("      reusing dof#$(reuse_dof)")
                                push!(dh.cell_dofs, reuse_dof)
                            end
                        end
                    else # distribute new dofs
                        Base._setindex!(facedicts[fi], nextdof, sface, -token)# facedicts[fi][sface] = nextdof,  store the first dof for this face
                        for facedof in 1:interpolation_info.nfacedofs
                            for d in 1:dh.field_dims[fi]
                                @debug println("      adding dof#$nextdof")
                                push!(dh.cell_dofs, nextdof)
                                nextdof += 1
                            end
                        end
                    end
                end # face loop
            end
            if interpolation_info.ncelldofs > 0 # always distribute new dofs for cell
                @debug println("    cell#$ci")
                for celldof in 1:interpolation_info.ncelldofs
                    for d in 1:dh.field_dims[fi]
                        @debug println("      adding dof#$nextdof")
                        push!(dh.cell_dofs, nextdof)
                        nextdof += 1
                    end
                end # cell loop
            end
        end # field loop
        # push! the first index of the next cell to the offset vector
        push!(dh.cell_dofs_offset, length(dh.cell_dofs)+1)
    end # cell loop
    dh.ndofs[] = maximum(dh.cell_dofs)
    dh.closed[] = true

    return dh, vertexdicts, edgedicts, facedicts

end

function celldofs!(global_dofs::Vector{Int}, dh::DofHandler, i::Int)
    @assert isclosed(dh)
    @assert length(global_dofs) == ndofs_per_cell(dh, i)
    unsafe_copyto!(global_dofs, 1, dh.cell_dofs, dh.cell_dofs_offset[i], length(global_dofs))
    return global_dofs
end

function cellnodes!(global_nodes::Vector{Int}, grid::Grid{dim,C}, i::Int) where {dim,C}
    @assert length(global_nodes) == nnodes(C)
    for j in 1:nnodes(C) # Currently assuming that DofHandler only has one celltype
        global_nodes[j] = grid.cells[i].nodes[j]
    end
    return global_nodes
end

function cellcoords!(global_coords::Vector{Vec{dim,T}}, grid::Grid{dim,C}, i::Int) where {dim,C,T}
    @assert length(global_coords) == nnodes(C)
    for j in 1:nnodes(C) # Currently assuming that DofHandler only has one celltype
        nodeid = grid.cells[i].nodes[j]
        global_coords[j] = grid.nodes[nodeid].x
    end
    return global_coords
end
cellcoords!(global_coords::Vector{<:Vec}, dh::DofHandler, i::Int) = cellcoords!(global_coords, dh.grid, i)

function celldofs(dh::DofHandler, i::Int)
    @assert isclosed(dh)
    n = ndofs_per_cell(dh, i)
    global_dofs = zeros(Int, n)
    unsafe_copyto!(global_dofs, 1, dh.cell_dofs, dh.cell_dofs_offset[i], n)
    return global_dofs
end

# Creates a sparsity pattern from the dofs in a DofHandler.
# Returns a sparse matrix with the correct storage pattern.
"""
    create_sparsity_pattern(dh::DofHandler)

Create the sparsity pattern corresponding to the degree of freedom
numbering in the [`DofHandler`](@ref). Return a `SparseMatrixCSC`
with stored values in the correct places.

See the [Sparsity Pattern](@ref) section of the manual.
"""
@inline create_sparsity_pattern(dh::DofHandler) = _create_sparsity_pattern(dh, false)

"""
    create_symmetric_sparsity_pattern(dh::DofHandler)

Create the symmetric sparsity pattern corresponding to the degree of freedom
numbering in the [`DofHandler`](@ref) by only considering the upper
triangle of the matrix. Return a `Symmetric{SparseMatrixCSC}`.

See the [Sparsity Pattern](@ref) section of the manual.
"""
@inline create_symmetric_sparsity_pattern(dh::DofHandler) = Symmetric(_create_sparsity_pattern(dh, true), :U)

function _create_sparsity_pattern(dh::DofHandler, sym::Bool)
    ncells = getncells(dh.grid)
    n = ndofs_per_cell(dh)
    N = sym ? div(n*(n+1), 2) * ncells : n^2 * ncells
    N += ndofs(dh) # always add the diagonal elements
    I = Int[]; resize!(I, N)
    J = Int[]; resize!(J, N)
    global_dofs = zeros(Int, n)
    cnt = 0
    for element_id in 1:ncells
        celldofs!(global_dofs, dh, element_id)
        @inbounds for j in 1:n, i in 1:n
            dofi = global_dofs[i]
            dofj = global_dofs[j]
            sym && (dofi > dofj && continue)
            cnt += 1
            if cnt > length(J)
                resize!(I, trunc(Int, length(I) * 1.5))
                resize!(J, trunc(Int, length(J) * 1.5))
            end
            I[cnt] = dofi
            J[cnt] = dofj

        end
    end
    @inbounds for d in 1:ndofs(dh)
        cnt += 1
        if cnt > length(J)
            resize!(I, trunc(Int, length(I) + ndofs(dh)))
            resize!(J, trunc(Int, length(J) + ndofs(dh)))
        end
        I[cnt] = d
        J[cnt] = d
    end
    resize!(I, cnt)
    resize!(J, cnt)
    V = zeros(length(I))
    K = sparse(I, J, V)
    return K
end

# dof renumbering
"""
    renumber!(dh::DofHandler, perm)

Renumber the degrees of freedom in the DofHandler according to the
permuation `perm`.

!!! warning
    Remember to do renumbering *before* adding boundary conditions,
    otherwise the mapping for the dofs will be wrong.
"""
function renumber!(dh::AbstractDofHandler, perm::AbstractVector{<:Integer})
    @assert isperm(perm) && length(perm) == ndofs(dh)
    cell_dofs = dh.cell_dofs
    for i in eachindex(cell_dofs)
        cell_dofs[i] = perm[cell_dofs[i]]
    end
    return dh
end

function WriteVTK.vtk_grid(filename::AbstractString, dh::AbstractDofHandler; compress::Bool=true)
    vtk_grid(filename, dh.grid; compress=compress)
end

"""
    reshape_to_nodes(dh::AbstractDofHandler, u::Vector{T}, fieldname::Symbol) where T

Reshape the entries of the dof-vector `u` which correspond to the field `fieldname` in nodal order.
Return a matrix with a column for every node and a row for every dimension of the field.
For superparametric fields only the entries corresponding to nodes of the grid will be returned. Do not use this function for subparametric approximations.
"""
function reshape_to_nodes(dh::DofHandler, u::Vector{T}, fieldname::Symbol) where T
    # make sure the field exists
    fieldname ∈ Ferrite.getfieldnames(dh) || error("Field $fieldname not found.")

    field_idx = findfirst(i->i==fieldname, getfieldnames(dh))
    offset = field_offset(dh, fieldname)
    field_dim = getfielddim(dh, field_idx)

    space_dim = field_dim == 2 ? 3 : field_dim
    data = fill(zero(T), space_dim, getnnodes(dh.grid))

    reshape_field_data!(data, dh, u, offset, field_dim)

    return data
end

function reshape_field_data!(data::Matrix{T}, dh::AbstractDofHandler, u::Vector{T}, field_offset::Int, field_dim::Int, cellset=Set{Int}(1:getncells(dh.grid))) where T

    _celldofs = Vector{Int}(undef, ndofs_per_cell(dh, first(cellset)))
    for cell in CellIterator(dh, collect(cellset))
        celldofs!( _celldofs, cell)
        counter = 1
        for node in getnodes(cell)
            for d in 1:field_dim
                data[d, node] = u[_celldofs[counter + field_offset]]
                @debug println("  exporting $(u[_celldofs[counter + field_offset]]) for dof#$(_celldofs[counter + field_offset])")
                counter += 1
            end
            if field_dim == 2
                # paraview requires 3D-data so pad with zero
                data[3, node] = 0
            end
        end
    end
    return data
end


text:
export lyapunovspectrum!, lyapunovspectrum
export perturbationgrowth!, perturbationgrowth, perturbationevolution

const δqind = SV{Int}(1,2)
const δpind = SV{Int}(3,4)

@inline curvature(::Semicircle{T}) where {T} = -one(T)
@inline curvature(::Disk{T}) where {T} = one(T)

################################################################################
## SPECULAR (LINEAR)
################################################################################

#="""
    specular!(p::AbstractParticle, o::Obstacle, offset::MArray)
Perform specular reflection based on the normal vector of the Obstacle.
The function updates the position and velocity of the particle
together with the components of 4 offset vectors stored in the matrix
`offset` as columns.
"""=#
function specular!(p::Particle{T}, o::Circular{T},
                   offset::Vector{SVector{4, T}}) where {T<:AbstractFloat}
    n = normalvec(o, p.pos)
    ti = SV{T}(-p.vel[2],p.vel[1])

    cosa = -dot(n, p.vel)
    p.vel = p.vel + 2.0*cosa*n

    tf = SV{T}(-p.vel[2], p.vel[1])

    @inbounds for k in eachindex(offset)
        δqprev = offset[k][δqind]
        δpprev = offset[k][δpind]
        # Formulas from Dellago, Posch and Hoover, PRE 53, 2, 1996: 1485-1501 (eq. 27)
        # with norm(p) = 1
        δq  = δqprev - 2*dot(δqprev,n)*n
        δp  = δpprev - 2*dot(δpprev,n)*n - curvature(o)*2/o.r*dot(δqprev,ti)/cosa*tf
        ###
        offset[k] = vcat(δq, δp)
    end
end

function specular!(p::Particle{T}, o::Union{InfiniteWall{T},FiniteWall{T}},
                   offset::Vector{SVector{4, T}}) where {T<:AbstractFloat}

    n = normalvec(o, p.pos)
    specular!(p, o)

    for k in eachindex(offset)
        δqprev = offset[k][δqind]
        δpprev = offset[k][δpind]
        # Formulas from Dellago, Posch and Hoover, PRE 53, 2, 1996: 1485-1501 (eq. 20)
        δq  = δqprev - 2.0*dot(δqprev,n)*n
        δp  = δpprev - 2.0*dot(δpprev,n)*n
        ###
        offset[k] = vcat(δq, δp)
    end

end

################################################################################
## SPECULAR (MAGNETIC)
################################################################################

function specular!(p::MagneticParticle{T}, o::Circular{T},
                   offset::Vector{SVector{4, T}}) where {T<:AbstractFloat}
    n = normalvec(o, p.pos)

    cosa = -dot(n, p.vel)

    Rn = SV{T}(-n[2], n[1])
    ti = SV{T}(-p.vel[2],p.vel[1])

    # magterm = ω*δτ_c*(BR-RB)*p_i / δqprev_normal
    # i.e. everything that can be calculated without using the previous offset vector
    magterm = -2*p.omega*((dot(p.vel, Rn)/(-cosa))*n + Rn)

    #actual specular reflection should not occur before magterm is computed
    p.vel = p.vel + 2*cosa*n
    tf = SV{T}(-p.vel[2], p.vel[1])

    for k in eachindex(offset)
        δqprev = offset[k][δqind]
        δpprev = offset[k][δpind]

        δqprev_normal = dot(δqprev, n)

        # Formulas derived analogously to Dellago et al., PRE 53, 2, 1996: 1485-1501
        δq  = δqprev - 2*δqprev_normal*n
        δp  = δpprev - 2*dot(δpprev,n)*n - curvature(o)* 2/o.r * dot(δqprev,ti)/cosa*tf + magterm*δqprev_normal
        ###
        offset[k] = vcat(δq, δp)
    end

end

function specular!(p::MagneticParticle{T}, o::Union{InfiniteWall{T},FiniteWall{T}},
                   offset::Vector{SVector{4, T}}) where {T<:AbstractFloat}

    n = normalvec(o, p.pos)

    cosa = -dot(n, p.vel)
    Rn = SV{T}(-n[2], n[1])

    # magterm = ω*δτ_c*(BR-RB)*p_i / δqprev_normal
    # i.e. everything that can be calculated without using the previous offset vector
    magterm = -2.0*p.omega*((dot(p.vel, Rn)/(-cosa))*n + Rn)

    #actual specular reflection should not occur before magterm is calculated
    p.vel = p.vel + 2*cosa*n

    for k in eachindex(offset)
        δqprev = offset[k][δqind]
        δpprev = offset[k][δpind]

        δqprev_normal = dot(δqprev, n)

        # Formulas derived analogously to Dellago et al., PRE 53, 2, 1996: 1485-1501
        δq  = δqprev - 2.0*δqprev_normal*n
        δp  = δpprev - 2.0*dot(δpprev,n)*n + magterm*δqprev_normal
        ###
        offset[k] = vcat(δq, δp)
    end
end

################################################################################
## RESOLVECOLLISON
################################################################################

#="""
    resolvecollision!(p::AbstractParticle, o::Union{Disk, InfiniteWall}, offset::MArray)
Resolve the collision between particle `p` and obstacle `o` of type *Circular*,
updating the components of the offset vectors stored in the matrix `offset` as columns.
"""=#
resolvecollision!(p::Particle{T}, o::Obstacle{T},
                  offset::Vector{SVector{4, T}}) where {T} = specular!(p, o, offset)

resolvecollision!(p::Particle{T}, o::PeriodicWall{T},
                  offset::Vector{SVector{4, T}}) where {T} = resolvecollision!(p, o)

function resolvecollision!(p::MagneticParticle{T}, o::Obstacle{T},
    offset::Vector{SVector{4, T}}) where {T}

    specular!(p, o, offset)
    p.center = find_cyclotron(p)
    return
end

# this method only exists to avoid ambiguity for MagneticParticles colliding with
# periodic walls.
resolvecollision!(p::MagneticParticle{T}, o::PeriodicWall{T},
                  offset::Vector{SVector{4, T}}) where {T} = resolvecollision!(p, o)

################################################################################
## OFFSET PROPAGATION
################################################################################

"""
    propagate_offset!(offset::MArray{Tuple{4,4},T}, p::AbstractParticle)
Computes the linearized evolution of the offset vectors during propagation for a
time interval `t`
"""
function propagate_offset!(offset::Vector{SVector{4, T}}, t::T,     #linear case
                           p::Particle{T}) where T
    for k in eachindex(offset)
        δΓ = offset[k]
        offset[k] = SVector{4,T}(δΓ[1] + t*δΓ[3], δΓ[2] + t*δΓ[4], δΓ[3], δΓ[4])
    end
end

#magnetic
function propagate_offset!(offset::Vector{SVector{4, T}}, t::T,
                           p::MagneticParticle{T}) where T
    ω = p.omega
    sω, cω = sincos(ω*t)

    for k in eachindex(offset)
        δΓ = offset[k]
        offset[k] = SVector{4,T}(δΓ[1] + sω/ω*δΓ[3] + (cω - 1)/ω*δΓ[4],
                                 δΓ[2] + (1 - cω)/ω*δΓ[3] + sω/ω*δΓ[4],
                                 cω*δΓ[3] - sω*δΓ[4],
                                 sω*δΓ[3] + cω*δΓ[4])
    end
end

################################################################################
## PROPAGATION & RELOCATION
################################################################################

#="""
    propagate!(p::AbstractParticle{T}, newpos::SV{T}, t::T,
    offset::MArray{Tuple{4,4},T})
Propagate the particle `p` for given time `t`, changing appropriately the the
`p.pos` and `p.vel` fields together with the components of the offset vectors
stored in the `offset` matrix.
"""=#
function propagate!(p::AbstractParticle{T}, newpos::SV{T}, t::T,
                    offset::Vector{SVector{4, T}}) where {T<: AbstractFloat}

    propagate!(p, newpos, t)
    propagate_offset!(offset, t, p)
    return
end

#="""
    relocate(p::AbstractParticle, o::Obstacle, t, cp::SV{T}, offset::MArray)
Propagate the particle's position for time `t` (corrected) and update the components
of the `offset` matrix.
"""=#
function relocate!(p::AbstractParticle{T}, o::Obstacle{T}, tmin, cp::SV{T},
                   offset::Vector{SVector{4, T}}) where {T <: AbstractFloat}

    okay = relocate!(p, o, tmin, cp)
    propagate_offset!(offset, tmin, p)
    return okay
end

################################################################################
## HIGH-LEVEL FUNCTION
################################################################################
function lyapunovspectrum!(p::AbstractParticle{T}, bd::Billiard{T}, t;
                           warning::Bool = false) where {T<:AbstractFloat}

    if t ≤ 0.0
        throw(ArgumentError(
            "`lyapunovspectrum()` cannot evolve backwards in time."))
    end

    # intial offset vectors
    offset = [SVector{4, T}(1,0,0,0), SVector{4, T}(0,1,0,0),
              SVector{4, T}(0,0,1,0), SVector{4, T}(0,0,0,1)]
    λ = zeros(T, 4)

    timecount = zero(T)
    count = zero(t)

    # check for pinning before evolution
    if ispinned(p, bd)
        warning && @warn "Pinned particle!"
        return λ
    end

    ismagnetic = typeof(p) <: MagneticParticle

    while count < t
        # bouncing
        i::Int, tmin::T, cp::SV{T} = next_collision(p, bd)
        relocate!(p, bd[i], tmin, cp, offset)
        resolvecollision!(p, bd[i], offset)

        timecount += tmin
        count += increment_counter(t, tmin)

        # update cyclotron data
        ismagnetic && (p.center = find_cyclotron(p))

        # QR decomposition to get Lyapunov spectrum
        Q, R = qr(hcat(offset[1], offset[2], offset[3], offset[4]))
        offset[1], offset[2], offset[3], offset[4] = Q[:, 1], Q[:, 2], Q[:, 3], Q[:, 4]

        for i ∈ 1:4
            λ[i] += log(abs(R[i,i]))
        end

    end#time loop

    return λ./timecount
end

"""
    lyapunovspectrum([p::AbstractParticle,] bd::Billiard, t)
Returns the finite time lyapunov exponents (averaged over time `t`)
for a given particle in a billiard table using the method outlined in [1].
`t` can be either `Float` or `Int`, meaning total time or total amount
of collisions.

Returns zeros for pinned particles.

If a particle is not given, a random one is picked through [`randominside`](@ref).
See [`parallelize`](@ref) for a parallelized version.

[1] : Ch. Dellago *et al*, [Phys. Rev. E **53** (1996)](http://link.aps.org/doi/10.1103/PhysRevE.53.1485)
"""
lyapunovspectrum(p::AbstractParticle, args...) = lyapunovspectrum!(copy(p), args...)
lyapunovspectrum(bd::Billiard, args...) =
    lyapunovspectrum!(randominside(bd), bd, args...)

################################################################################
## Raw perturbation growth
################################################################################
function perturbationgrowth!(p::AbstractParticle{T}, bd::Billiard{T},
                             t) where {T<:AbstractFloat}

    offset = [SVector{4, T}(1,1,1,1)]

    count = zero(t)
    timecount = zero(T)

    Δ = Vector{SVector{4, T}}() # perturbation vectors
    tim = T[] # sample times
    obst = Int[] # obstacle indices

     # check for pinning before evolution
    if ispinned(p, bd)
        return tim, Δ, obst
    end
    ismagnetic = typeof(p) <: MagneticParticle

    while count < t
        off = offset[1]
        # bounce
        i::Int, tmin::T, cp::SV{T} = next_collision(p, bd)
        relocate!(p, bd[i], tmin, cp, offset)
        timecount += tmin

        # push perturbations before collision
        push!(tim, timecount)
        push!(obst, i)

        # push after propagation evolution
        push!(Δ, offset[1] ./ off)
        off = offset[1]

        resolvecollision!(p, bd[i], offset)
        # push perturbations after collision
        push!(tim, timecount)
        push!(obst, i)
        push!(Δ, offset[1] ./ off)

        ismagnetic && (p.center = find_cyclotron(p))

        # normalize perturbation vector
        offset[1] = offset[1] ./ norm(offset[1])

        count += increment_counter(t, tmin)

    end#time loop

    return tim, Δ, obst
end

"""
    perturbationgrowth([p,] bd, t) -> ts, Rs, is
Calculate the evolution of the perturbation vector `Δ` along the trajectory of `p` in
`bd` for total time `t`. `Δ` is initialised as `[1,1,1,1]`.

If a particle is not given, a random one is picked through [`randominside`](@ref).
Returns empty lists for pinned particles.

## Description
This function *safely* computes the time evolution of a perturbation vector using the
linearized dynamics of the system, as outlined by [1]. Because the dynamics are linear,
we can safely re-normalize the perturbation vector after every collision (otherwise the
perturbations grow to infinity).

Immediately before *and after* every collison, this function computes
* the current time.
* the element-wise ratio of Δ with its previous value
* the obstacle index of the current obstacle
and returns these in three vectors `ts, Rs, is`.

To obtain the *actual* evolution of the perturbation vector you can use
the function `perturbationevolution(Rs)` which simply does
```julia
Δ = Vector{SVector{4,Float64}}(undef, length(R))
Δ[1] = R[1]
for i in 2:length(R)
    Δ[i] = R[i] .* Δ[i-1]
end
```

[1] : Ch. Dellago *et al*, [Phys. Rev. E **53** (1996)](http://link.aps.org/doi/10.1103/PhysRevE.53.1485)
"""
perturbationgrowth(p::AbstractParticle, args...) = perturbationgrowth!(copy(p), args...)
perturbationgrowth(bd::Billiard, args...) =
    perturbationgrowth!(randominside(bd), bd, args...)

function perturbationevolution(R::Vector{SVector{4, T}}) where T
    Δ = Vector{SVector{4,T}}(undef, length(R))
    Δ[1] = R[1]
    @inbounds for i in 2:length(R)
        Δ[i] = R[i] .* Δ[i-1]
    end
    return Δ
end

@deprecate pertubationevolution perturbationevolution

ext: jl


text:
@testset "Constructors and basic properties" begin
    let
        F = FullBinner()
        F2 = FullBinner()
        @test typeof(F) <: BinningAnalysis.AbstractBinner{Float64}
        @test eltype(F) == Float64
        @test ndims(F) == 1
        @test length(F) == 0
        @test size(F) == (0,)
        @test lastindex(F) == 0
        @test axes(F) == (Base.OneTo(0),)
        @test isempty(F)
        @test F == F2
        @test F ≈ F2
    end

    let F = FullBinner(ComplexF64)
        @test typeof(F) <: BinningAnalysis.AbstractBinner{ComplexF64}
        @test eltype(F) == ComplexF64
    end

    let x = [0.84, 0.381169, 0.34174, 0.888868, 0.0973183, 0.722725, 0.0957878, 0.432443, 0.755033, 0.864252]
        F = FullBinner(x)
        F2 = FullBinner(Float64)
        @test length(F) == 10
        @test lastindex(F) == 10
        @test axes(F) == (Base.OneTo(10),)
        @test !isempty(F)
        @test F != F2
        @test !(F ≈ F2)
        append!(F2, x)
        @test F == F2
        @test F ≈ F2
    end

    let x = [rand(rng, 2,3) for _ in 1:100]
        F = FullBinner(x)
        @test length(F) == 100
        @test eltype(F) == Array{Float64,2}
    end
end

@testset "Scalars statistics" begin
    # Real, Complex
    for F in (
            FullBinner(collect(1:10_000)),
            FullBinner(collect((1:10_000) .+ ((10_000:-1:1) .* im)))
        )

        # Unbinned Statistics vs Base
        µ = mean(F.x)
        v = var(F.x, mean = µ)
        Δµ = sqrt(v / 10_000)

        @test mean(F, 1) ≈ µ
        @test var(F, 1) ≈ v
        @test varN(F, 1) ≈ v / 10_000
        @test std_error(F, 1) ≈ Δµ

        # Binned stats (these values are strongly correlated)
        binsize = BinningAnalysis._reliable_level(F)
        @test binsize == div(length(F), 32)
        @test BinningAnalysis._eachlevel(F) == 1:binsize
        @test mean(F) ≈ µ
        @test var(F) ≈ v rtol = 0.05
        @test std_error(F) ≈ Δµ * sqrt(binsize) rtol = 0.05
        @test tau(F) ≈ 0.5(binsize-1) rtol = 0.05
    end
end

@testset "Arrays statistics" begin
    for T in (Float64, ComplexF64)
        StableRNGs.seed!(rng, 123)
        F = FullBinner([rand(rng, T, 2, 3) for _ in 1:1000])

        # Unbinned Statistics vs Base
        µ = mean(F.x)
        v = var(F.x, mean = µ)
        Δµ = sqrt.(v ./ 1000)

        @test mean(F, 1) ≈ µ
        @test var(F, 1) ≈ v
        @test varN(F, 1) ≈ v ./ 1000
        @test std_error(F, 1) ≈ Δµ

        # Binned stats (these values are uncorrelated)
        # Note that the errors can be quite large due to the relatively small
        # sample size
        binsize = BinningAnalysis._reliable_level(F)
        @test binsize == div(length(F), 32)
        @test BinningAnalysis._eachlevel(F) == 1:binsize
        @test mean(F) ≈ µ
        @test var(F) ≈ v/binsize rtol = 0.25
        @test std_error(F) ≈ Δµ rtol = 0.15
        @test all(tau(F) .< 1)

        if T == Float64
            @test var(F) ≈ [0.003208002497613544 0.0029222442318581378 0.002097985378538715; 0.0027042522201602293 0.0035246447808234393 0.003948312674521377]
        else
            @test var(F) ≈ [0.006091716418175418 0.0052399001679973665 0.007313376580285079; 0.0059813002774426445 0.006343840253127798 0.004989114094997521]
        end
    end

    F = FullBinner(Vector{Float64})
    x = Float64[1, 2, 3]
    push!(F, x)
    x .= 5
    push!(F, x)
    @test F.x[1] == Float64[1, 2, 3]
    @test F.x[2] == Float64[5, 5, 5]
end

@testset "Cosmetics (show, print, etc.)" begin
    F = FullBinner();
    # empty binner
    oldstdout = stdout
    (read_pipe, write_pipe) = redirect_stdout()
    println(F) # compact
    show(write_pipe, MIME"text/plain"(), F) # full
    redirect_stdout(oldstdout);
    close(write_pipe);

    # compact
    if VERSION < v"1.7.0"
        @test readline(read_pipe) == "FullBinner{Float64,Array{Float64,1}}()"
        @test readline(read_pipe) == "FullBinner{Float64,Array{Float64,1}}"
    else
        @test readline(read_pipe) == "FullBinner{Float64,Vector{Float64}}()"
        @test readline(read_pipe) == "FullBinner{Float64,Vector{Float64}}"
    end
    # full
    @test readline(read_pipe) == "| Count: 0"
    @test length(readlines(read_pipe)) == 0
    close(read_pipe);

    # filled binner
    StableRNGs.seed!(rng, 123)
    append!(F, rand(rng, 1000))
    (read_pipe, write_pipe) = redirect_stdout()
    show(write_pipe, MIME"text/plain"(), F)
    redirect_stdout(oldstdout);
    close(write_pipe);
    @test if VERSION < v"1.7.0"
        readline(read_pipe) == "FullBinner{Float64,Array{Float64,1}}"
    else
        readline(read_pipe) == "FullBinner{Float64,Vector{Float64}}"
    end
    @test readline(read_pipe) == "| Count: 1000"
    @test readline(read_pipe) == "| Mean: 0.49387"
    @test length(readlines(read_pipe)) == 0
    close(read_pipe);
end

max_stars_repo_name: crstnbr/BinningAnalysis

max_issues_repo_name: crstnbr/BinningAnalysis

max_forks_repo_name: crstnbr/BinningAnalysis


text:
using MagneticResonanceSignals
using AxisArrays
using Test
using Unitful

# UUGH. See the improvements at https://github.com/JuliaArrays/AxisArrays.jl/pull/152
getaxis(samples, name) = AxisArrays.axes(samples, Axis{name}).val

@testset "MagneticResonanceSignals" begin

include("core.jl")
include("mr_load.jl")
include("processing.jl")
include("windows.jl")

include("fixedstring.jl")
include("io_twix.jl")
include("io_rda.jl")
include("io_felix.jl")
include("io_nmrpipe.jl")

end

max_stars_repo_name: TRIImaging/MagneticResonanceSignals

max_issues_repo_name: TRIImaging/MagneticResonanceSignals

max_forks_repo_name: TRIImaging/MagneticResonanceSignals


text:
# cellphone.jl
# Cellphone simulation from PySim

# Parameters

Nchannels = 4
maxN = 1_000
lam = 1.0
mu = 0.6667
meanLifeTime = inv(mu)
Nhours = 10
interv = 60.0
gap = 15.0

# Globals

Nfree = Nchannels
totalBusyVisits = 0
totalBusyTime = 0.0
busyStartTime = 0.0
busyEndTime = 0.0

# Statistics

m = nothing
bn = nothing

@process generator(n::Integer) begin
    iatime = inv(lam)
    for i in 1:n
        @schedule now job(i)
        work(rand(Exponential(iatime)))
    end
end

@process job(i::Integer) begin
    global Nfree, totalBusyVisits, totalBusyTime, busyStartTime, busyEndTime
    if Nfree > 0
        Nfree -= 1
        if Nfree == 0
            # Start busy period
            busyStartTime = current_time()
            totalBusyVisits += 1
            interIdleTime = current_time() - busyEndTime
        end
        work(rand(Exponential(inv(mu))))
        if Nfree == 0
            # End busy period
            busyEndTime = current_time()
            busy = current_time() - busy_start_time
            totalBusyTime += busy
        end
        Nfree += 1
    end
end

@process statistician() begin
    global Nfree, totalBusyVisits, totalBusyTime, busyStartTime, busyEndTime
    for i = 1:Nhours
        # Wait the specified gap time
        work(gap)
        # Initialize
        totalBusyTime = 0.0
        totalBusyVisits = 0
        if Nfree == 0
            busyStartTime = current_time()
        end
        # Wait the specified interval time
        work(interv)
        # Trace busy time and busy visits
        if Nfree == 0
            totalBusyTime += current_time() - busyStartTime
        end
        println("$(current_time()): busy time = $totalBusyTime; busy visits = $totalBusyVisits")
        # Tally statistics
        m.value = totalBusyTime
        bn.value = totalBusyVisits
    end
    # Print final statistics
    println("Busy time:   mean = $(m.mean), variance = $(m.variance)")
    println("Busy number: mean = $(bn.mean), variance = $(bn.variance)")
end

@event stop_sim() begin
    stop_simulation()
end

function run_simulation()
    global Nfree, totalBusyVisits, totalBusyTime, busyStartTime, busyEndTime
    global m, bn
    @simulation begin
        # Print parameters
        println("lambda = $lam")
        println("mu     = $mu")
        println("s      = $s")
        println("NHours = $Nhours")
        println("interv = $interv")
        println("gap    = $gap")
        # Initialize the global variables
        Nfree = Nchannels
        totalBusyVisits = 0
        totalBusyTime = 0.0
        busyEndTime = 0.0
        # Create statistics
        m = Variable{Int64}(data=:tally, history=true)
        bn = Variable{Float64}(data=:tally, history=true)
        # Create initial processes and start the execution
        @schedule at 0.0 generator(maxN)
        @schedule at 0.0 statistician()
        @schedule at 10_000.0 stop_sim()
        start_simulation()
    end
end


text:
# Discriminant AutoRegressive Normal
#

struct DARNormal1D{K}
    arnormals::ARNormal1DSet{K}

    # Prior/posterior over `ξ`, the IIR filter modeling the noise
    ξprior::Normal
    ξposterior::Normal
end

#######################################################################
# Constructors

function DARNormal1D(
    arnormals::ARNormal1DSet{K},
    u₀::Vector{T},
    V₀::Matrix{T};
    u::Union{Vector{T}, Missing} = missing,
    V::Union{Matrix{T}, Missing} = missing
) where {K, T <: AbstractFloat}
    u = ismissing(u) ? copy(u₀) : u
    V = ismissing(V) ? copy(V₀) : V
    DARNormal1D{K}(arnormals, Normal(u₀, V₀), Normal(u, V))
end

#######################################################################
# Expected log-likelihood

function loglikelihood(
    model::DARNormal1D{K},
    x::Vector{T}
) where {K, T<:AbstractFloat}
    r = Regressors1D(K, x)

    η_hs = vcat([gradlognorm(m.hposterior)' for m in model.arnormals]...)
    η_λs = vcat([gradlognorm(m.λposterior)' for m in model.arnormals]...)
    η_ξ = gradlognorm(model.ξposterior)

    ms = [mean(m.hposterior) for m in model.arnormals]
    u = mean(model.ξposterior)

    # Sufficient statistics for h
    s1 = hcat([xₜ * x̂ₜ for (xₜ, x̂ₜ) in zip(x, r)]...)
    s2 = hcat([vec(x̂ₜ* x̂ₜ') for x̂ₜ in r]...)
    stats_hξ = vcat(s1, -.5 * s2)

    iter = zip(eachrow(η_hs), eachrow(η_λs), ms)
    vcat([(hcat(-.5 * x.^2 .+ stats_hξ' * (η_h + η_ξ)  .- s2' *vec(m * u'),
                .5 * ones(T, length(x))) * η_λ)'
          for (η_h, η_λ, m) in iter]...) .- .5 * log(2π)
end

function loglikelihood(
    models::Vector{DARNormal1D{K}},
    X::Matrix{T}
) where {K, T<:AbstractFloat}
    reduce(+, [loglikelihood(models[d], X[d, :]) for d in 1:size(X, 1)])
end

#######################################################################
# Posterior predictive
#
# NOTE: we use the MAP of precision parameters and we integrate over the
#       the filters (h_1, h_2, ...)

function predict(
    model::DARNormal1D{K},
    x::Vector{T}
) where {K, T<:AbstractFloat}

    @warn "erroneous computation"

    r = Regressors1D(K, x)
    λ_map_inv = 1 / mean(model.λposterior)[1]
    m = model.hposterior.μ
    Σ = model.hposterior.Σ

    # Sufficient statistics for h
    s1 = hcat([x̂ₜ for x̂ₜ in r]...)
    s2 = hcat([vec(x̂ₜ* x̂ₜ') for x̂ₜ in r]...)

    ψ = s1' * m
    ϑ_inv = s2' * vec(Σ)

    prec = 1 ./ (λ_map_inv .+ ϑ_inv)

    -.5 .* prec .* (x .- ψ).^2 .+ .5 .* log.(prec) .- .5 * log(2π)
end

function predict(
    models::Vector{DARNormal1D{K}},
    x::Vector{T}
) where {K, T<:AbstractFloat}
    vcat([predict(models[s], x)' for s in 1:length(models)]...)
end

"""
    predict(model, X)

Return the logarithm of the posterior predictive distribution for each
frame of `X`
"""
function predict(
    models::Vector{Vector{DARNormal1D{K}}},
    X::Matrix{T}
) where {K, T<:AbstractFloat}
    reduce(+, [predict(models[d], X[d, :]) for d in 1:size(X, 1)])
end

########################################################################
# Accumulate statistics

function accstats_h(
    model::ARNormal1D{K},
    x::Vector{T},
    γ::Vector{T},
    ξ::Vector{T}
) where {K, T<:AbstractFloat}

    r = Regressors1D(K, x)

    # Sufficient statistics
    λ = mean(model.λposterior)
    s1 = hcat([(xₜ- ξ' * x̂ₜ) * x̂ₜ for (xₜ, x̂ₜ) in zip(x, r)]...)
    s2 = hcat([vec(x̂ₜ* x̂ₜ') for x̂ₜ in r]...)
    stats_h = λ .* vcat(s1, -.5 * s2)

    stats_h * γ
end

function accstats_h(
    model::DARNormal1D{K},
    x::Vector{T},
    γ::Matrix{T}
) where {K, T<:AbstractFloat}
    ξ = mean(model.ξposterior)
    [accstats_h(m, x, γ[d, :], ξ) for (d, m) in enumerate(model.arnormals)]
end

function accstats_ξ(
    model::ARNormal1D{K},
    x::Vector{T},
    γ::Vector{T}
) where {K, T<:AbstractFloat}

    r = Regressors1D(K, x)

    # Sufficient statistics
    h = mean(model.hposterior)
    λ = mean(model.λposterior)[1]
    s1 = hcat([λ * (xₜ- h' * x̂ₜ) * x̂ₜ for (xₜ, x̂ₜ) in zip(x, r)]...)
    s2 = hcat([λ * vec(x̂ₜ* x̂ₜ') for x̂ₜ in r]...)
    stats_ξ =  vcat(s1, -.5 * s2)

    stats_ξ * γ
end

function accstats_ξ(
    model::DARNormal1D{K},
    x::Vector{T},
    γ::Matrix{T}
) where {K, T<:AbstractFloat}
    reduce(+, [accstats_ξ(m, x, γ[s, :]) for (s, m) in enumerate(model.arnormals)])
end

function accstats_λ(
    model::DARNormal1D{K},
    x::Vector{T},
    γ::Matrix{T}
) where {K, T<:AbstractFloat}
    r = Regressors1D(K, x)

    η_hs = vcat([gradlognorm(m.hposterior)' for m in model.arnormals]...)
    η_ξ = gradlognorm(model.ξposterior)

    ms = [mean(m.hposterior) for m in model.arnormals]
    u = mean(model.ξposterior)

    # Sufficient statistics for h
    s1 = hcat([xₜ * x̂ₜ for (xₜ, x̂ₜ) in zip(x, r)]...)
    s2 = hcat([vec(x̂ₜ* x̂ₜ') for x̂ₜ in r]...)
    stats_hξ = vcat(s1, -.5 * s2)

    [hcat(-.5 * x.^2 .+ stats_hξ' * (η_h + η_ξ)  .- s2' *vec(m * u'),
          .5 * ones(T, length(x)))' * γₛ
     for (η_h, m, γₛ) in zip(eachrow(η_hs), ms, eachrow(γ))]
end

function accstats_λ(
    models::Vector{DARNormal1D{K}},
    x::Vector{T},
    γ::Matrix{T}
) where {K, T<:AbstractFloat}
    [accstats_λ(m, x, γ[d, :]) for (d, m) in enumerate(models)]
end

#######################################################################
# Utilities

function save(
    fname::AbstractString,
    model::Vector{DARNormal1D{K}}
) where K
    bson(fname, Dict(:model => model, :order=>K, :type=>typeof(model)))
end

max_stars_repo_name: BUTSpeechFIT/ConvHMM

max_issues_repo_name: BUTSpeechFIT/ConvHMM

max_forks_repo_name: BUTSpeechFIT/ConvHMM


text:
immutable HermitianTensorOperator{S<:HilbertSpace,P,T}<:AbstractTensorMap{S,P,T}
    map::AbstractTensorMap{S,P,T}
    function HermitianTensorOperator(A::AbstractTensorMap{S,P,T})
        domain(A)==dual(conj(codomain(A))) || throw(SpaceError("Not an operator, i.e. domain != codomain"))
        new(A)
    end
end
hermitian{S,P,T}(A::AbstractTensorMap{S,P,T})=HermitianTensorOperator{S,P,T}(A::AbstractTensorMap{S,P,T})

# properties
domain(A::HermitianTensorOperator)=domain(A.map)
codomain(A::HermitianTensorOperator)=codomain(A.map)

Base.isreal(A::HermitianTensorOperator)=isreal(A.map)
Base.issym(A::HermitianTensorOperator)=isreal(A.map)
Base.ishermitian(A::HermitianTensorOperator)=true
Base.isposdef(A::HermitianTensorOperator)=isposdef(A.map)

# comparison
==(A::HermitianTensorOperator,B::HermitianTensorOperator)=A.map==B.map

# multiplication with vector
Base.A_mul_B!{S,P}(y::AbstractTensor{S,P},A::HermitianTensorOperator{S,P},x::AbstractTensor{S,P})=Base.A_mul_B!(y,A.map,x)
Base.Ac_mul_B!{S,P}(y::AbstractTensor{S,P},A::HermitianTensorOperator{S,P},x::AbstractTensor{S,P})=Base.A_mul_B!(y,A.map,x)


text:
using JuMP
using GLPKMathProgInterface

myModel = Model(solver=GLPKSolverLP())
# Name of the model object. All constraints and variables of an optimization problem are associated
# with a particular model object. The name of the model object does not have to be myModel, it can be yourModel too! The argument of Model,
# solver=GLPKsolverLP() means that to solve the optimization problem we will use GLPK solver.

#VARIABLES
#---------

# A variable is modelled using @defVar(name of the model object, variable name and bound, variable type)
# Bound can be lower bound, upper bound or both. If no variable type is defined, then it is treated as
#real. For binary variable write Bin and for integer use Int.

@defVar(myModel, x >= 0) # Models x >=0

# Some possible variations:
# @defVar(myModel, x, Binary) # No bound on x present, but x is a binary variable now
# @defVar(myModel, x <= 10) # This one defines a variable with lower bound x <= 10
# @defVar(myModel, 0 <= x <= 10, Int) # This one has both lower and upper bound, and x is an integer

@defVar(myModel, y >= 0) # Models y >= 0

#OBJECTIVE
#---------

@setObjective(myModel, Min, x + y) # Sets the objective to be minimized. For maximization use Max

#CONSTRAINTS
#-----------

@addConstraint(myModel, x + y <= 1) # Adds the constraint x + y <= 1

#THE MODEL IN A HUMAN-READABLE FORMAT
#------------------------------------
println("The optimization problem to be solved is:")
print(myModel) # Shows the model constructed in a human-readable form

#SOLVE IT AND DISPLAY THE RESULTS
#--------------------------------
status = solve(myModel) # solves the model

println("Objective value: ", getObjectiveValue(myModel)) # getObjectiveValue(model_name) gives the optimum objective value
println("x = ", getValue(x)) # getValue(decision_variable) will give the optimum value of the associated decision variable
println("y = ", getValue(y))


text:
nlower(n::Integer) = (n * (n + 1)) >>> 1
nlower{T}(A::LowerTriangular{T, Matrix{T}}) = nlower(Compat.LinAlg.checksquare(A))

"""
    getθ!{T}(v::AbstractVector{T}, A::LowerTriangular{T, Matrix{T}})

Overwrite `v` with the elements of the lower triangle of `A` (column-major ordering)
"""
function getθ!{T}(v::AbstractVector{T}, A::LowerTriangular{T,Matrix{T}})
    n, Ad = Compat.LinAlg.checksquare(A), A.data
    if length(v) ≠ nlower(n)
        throw(DimensionMismatch("length(v) = $(length(v)) should be $(nlower(n))"))
    end
    k = 0
    for j = 1 : n, i in j : n
        v[k += 1] = Ad[i, j]
    end
    v
end

"""
    getθ(A::LowerTriangular{T, Matrix{T}})

Return a vector of the elements of the lower triangle of `A` (column-major ordering)
"""
getθ{T}(A::LowerTriangular{T, Matrix{T}}) = getθ!(Array(T, nlower(A)), A)

"""
    lowerbd{T}(A::LowerTriangular{T,Matrix{T}})

Return the vector of lower bounds on the parameters, `θ`.

These are the elements in the lower triangle in column-major ordering.
Diagonals have a lower bound of `0`.  Off-diagonals have a lower-bound of `-Inf`.
"""
function lowerbd{T}(A::LowerTriangular{T,Matrix{T}})
    n = Compat.LinAlg.checksquare(A)
    res = fill(convert(T, -Inf), nlower(n))
    k = -n
    for j in n+1:-1:2
        res[k += j] = zero(T)
    end
    res
end

"""
    LT(A)

Create a lower triangular matrix compatible with the blocks of `A`
and initialized to the identity.
"""
LT{T}(A::ScalarReMat{T}) = LowerTriangular(ones(T, (1, 1)))
LT{T}(A::VectorReMat{T}) = LowerTriangular(full(eye(T, size(A.z, 1))))

function setθ!{T}(A::LowerTriangular{T, Matrix{T}}, v::AbstractVector{T})
    Ad = A.data
    n = Compat.LinAlg.checksquare(Ad)
    if length(v) ≠ nlower(n)
        throw(DimensionMismatch("length(v) = $(length(v)) should be $(nlower(n))"))
    end
    offset = 0
    for j in 1 : n, i in j : n
        Ad[i, j] = v[offset += 1]
    end
    A
end

"""
    tscale!(A::LowerTriangular, B::HBlkDiag)
    tscale!(A::LowerTriangular, B::Diagonal)
    tscale!(A::LowerTriangular, B::DenseVecOrMat)
    tscale!(A::LowerTriangular, B::SparseMatrixCSC)
    tscale!(A::HBlkDiag, B::LowerTriangular)
    tscale!(A::Diagonal, B::LowerTriangular)
    tscale!(A::DenseVecOrMat, B::LowerTriangular)
    tscale!(A::SparseMatrixCSC, B::LowerTriangular)

Scale a matrix using the implicit expansion of the lower triangular matrix
to a diagonal or homogeneous block diagonal matrix.

Used in evaluating `Λ'Z'ZΛ` from `Z'Z` without explicitly evaluating the matrix `Λ`.
"""
function tscale!(A::LowerTriangular,B::HBlkDiag)
    Ba = B.arr
    r, s, k = size(Ba)
    n = Compat.LinAlg.checksquare(A)
    if n ≠ r
        throw(DimensionMismatch("size(A,2) ≠ blocksize of B"))
    end
    Ac_mul_B!(A, reshape(Ba,(r,s * k)))
    B
end

function tscale!{T}(A::LowerTriangular{T}, B::Diagonal{T})
    if size(A, 1) ≠ 1
        throw(DimensionMismatch("A must be a 1×1 LowerTriangular"))
    end
    scale!(A.data[1], B.diag)
    B
end

function tscale!{T}(A::LowerTriangular{T}, B::DenseVecOrMat{T})
    if (l = size(A, 1)) == 1
        return scale!(A.data[1], B)
    end
    m, n = size(B, 1), size(B, 2)  # this sets n = 1 when B is a vector
    q, r = divrem(m, l)
    if r ≠ 0
        throw(DimensionMismatch("size(B,1) is not a multiple of size(A,1)"))
    end
    Ac_mul_B!(A, reshape(B, (l, q * n)))
    B
end

function tscale!{T}(A::LowerTriangular{T}, B::SparseMatrixCSC{T})
    if (l = size(A, 1)) == 1
        scale!(A.data[1], B.nzval)
    else
        m, n = size(B)
        q, r = divrem(m, l)
        if r ≠ 0
            throw(DimensionMismatch("size(B, 1) is not a multiple of size(A, 1)"))
        end
        q, r = divrem(nnz(B), l)
        if r ≠ 0
            throw(DimensionMismatch("nnz(B) is not a multiple of size(A, 1)"))
        end
        Ac_mul_B!(A, reshape(B.nzval, (l, q)))
    end
    B
end

function tscale!{T}(A::SparseMatrixCSC{T}, B::LowerTriangular{T})
    if (l = size(B, 1)) == 1
        scale!(A.nzval, B.data[1])
    else
        m, n = size(A)
        q, r = divrem(nnz(A), l)
        if r ≠ 0
            throw(DimensionMismatch("nnz(A) is not a multiple of size(B, 1)"))
        end
        q, r = divrem(n, l)
        if r ≠ 0
            throw(DimensionMismatch("size(A, 2) is not a multiple of size(B, 1)"))
        end
        Ar = rowvals(A)
        Acp = A.colptr
        Anz = nonzeros(A)
        offset = 0
        for k in 1 : q
            nzr1 = nzrange(A, offset + 1)
            Ar1 = Compat.view(Ar, nzr1)
            for j in 2 : l
                if Compat.view(Ar, nzrange(A, offset + j)) ≠ Ar1
                    throw(ArgumentError("A does not have block structure for tscale!"))
                end
            end
            lnzr = length(Ar1)
            Aa = reshape(Compat.view(Anz, nzr1[1] + (0 : (lnzr * l - 1))), (lnzr, l))
            A_mul_B!(Aa, Aa, B)
            offset += l
        end
    end
    A
end

function tscale!{T}(A::Diagonal{T}, B::LowerTriangular{T})
    if (l = Compat.LinAlg.checksquare(B)) ≠ 1
        throw(DimensionMismatch(
        "in tscale!(A::Diagonal,B::LowerTriangular) B must be 1×1"))
    end
    scale!(B.data[1], A.diag)
    A
end

function tscale!{T}(A::HBlkDiag{T}, B::LowerTriangular{T})
    aa = A.arr
    r, s, l = size(aa)
    scr = Array(T, r, s)
    for k in 1 : l
        for j in 1 : s, i in 1 : r
            scr[i, j] = aa[i, j, k]
        end
        A_mul_B!(scr, B)
        for j in 1 : s, i in 1 : r
            aa[i, j, k] = scr[i, j]
        end
    end
    A
end

function tscale!{T}(A::StridedMatrix{T}, B::LowerTriangular{T})
    if (l = size(B,1)) == 1
        return scale!(A, B.data[1])
    end
    m, n = size(A)
    q, r = divrem(n, l)
    if r ≠ 0
        throw(DimensionMismatch("size(A,2) = $n must be a multiple of size(B,1) = $l"))
    end
    for k in 0:(q - 1)
        A_mul_B!(Compat.view(A, : , k * l + (1:l)), B)
    end
    A
end


text:
#This produces a StitchedSeries that incorporates both the stitch transform and slicewise transformations that compensate for
#wiggle in the field of view resulting from fast dynamic piezo operation
#slice_tfms is the set of slicewise transformations that must be applied to img_top (and also to img_bottom after it is aligned to the space of img_top)
#octr is the point around which rotations encoded by slice_tfms are centered.  This is the (x,y) center of the image sequence used for calibration
#This function assumes that sequence was acquired on camera1.  If the (x,y) size of that sequence differes from the (x,y) sizes of the top and bottom
#images then the default octr will not be optimal.  Override it with the true center
function stitch_slicewise(img_top, img_bottom, st_tfm, slice_tfms, out_type=Float64; octr=RegisterQD.center(view(img_top,:,:,1,1)), correct_bias=true, sqrt_tfm=false, flip_y_bottom=true)
    tfms_top, tfms_bottom = recenter_tfms(slice_tfms, st_tfm, img_top, octr)
    za_top = ZWarpedArray(img_top, tfms_top, out_type; correct_bias=false, sqrt_tfm=false)
    return StitchedSeries(za_top, img_bottom, tfms_bottom, Float64; correct_bias=correct_bias, sqrt_tfm=sqrt_tfm, flip_y_bottom=flip_y_bottom)
end

#assume that slice_tfms are already centered on the calibration images, so need to be recentered by the difference between that point and the
#center of camera1's image.
function recenter_tfms(slice_tfms, st_tfm, cam1, octr)
    tfmstop = []
    tfmsbottom = []
    cam1ctr = RegisterQD.center(view(cam1, :,:,1,1)).-octr
    #The below calculation of cam2ctr isn't perfect but is close enough that it's not noticeable by eye
    cam2ctr = cam1ctr + [0.0; -(size(cam1,2)+1)/2] #The center must move negatively in y (assumes top and bottom strips have equal y sz and perfect alignment)
    #Future: maybe for cam2ctr we can just apply the stitching transform to the center of the image? watch for sign errors
    for i =1:length(slice_tfms)
        push!(tfmstop, recenter(slice_tfms[i], cam2ctr))
        push!(tfmsbottom, st_tfm ∘ recenter(slice_tfms[i], cam2ctr))
    end
    return tfmstop, tfmsbottom
end


text:
include("uninformative/out.jl")

include("addition/marginals.jl")
include("addition/out.jl")
include("addition/in1.jl")
include("addition/in2.jl")

include("subtraction/marginals.jl")
include("subtraction/out.jl")
include("subtraction/in1.jl")
include("subtraction/in2.jl")

include("multiplication/marginals.jl")
include("multiplication/out.jl")
include("multiplication/in.jl")

include("gamma/out.jl")
include("gamma/marginals.jl")

include("gamma_shape_rate/out.jl")
include("gamma_shape_rate/marginals.jl")

include("beta/out.jl")
include("beta/marginals.jl")

include("dirichlet/marginals.jl")
include("dirichlet/out.jl")

include("matrix_dirichlet/out.jl")

include("categorical/out.jl")
include("categorical/p.jl")

include("bernoulli/p.jl")
include("bernoulli/out.jl")
include("bernoulli/marginals.jl")

include("gcv/x.jl")
include("gcv/y.jl")
include("gcv/z.jl")
include("gcv/k.jl")
include("gcv/w.jl")
include("gcv/marginals.jl")
include("gcv/gaussian_extension.jl")

include("kernel_gcv/x.jl")
include("kernel_gcv/y.jl")
include("kernel_gcv/z.jl")
include("kernel_gcv/marginals.jl")

include("mv_normal_mean_covariance/out.jl")
include("mv_normal_mean_covariance/mean.jl")
include("mv_normal_mean_covariance/marginals.jl")

include("normal_mean_variance/out.jl")
include("normal_mean_variance/mean.jl")
include("normal_mean_variance/var.jl")
include("normal_mean_variance/marginals.jl")

include("mv_normal_mean_precision/out.jl")
include("mv_normal_mean_precision/mean.jl")
include("mv_normal_mean_precision/precision.jl")
include("mv_normal_mean_precision/marginals.jl")

include("mv_normal_weightedmean_precision/marginals.jl")

include("normal_mean_precision/out.jl")
include("normal_mean_precision/mean.jl")
include("normal_mean_precision/precision.jl")
include("normal_mean_precision/marginals.jl")

include("wishart/out.jl")
include("wishart/marginals.jl")

include("normal_mixture/switch.jl")
include("normal_mixture/m.jl")
include("normal_mixture/p.jl")
include("normal_mixture/out.jl")

include("gamma_mixture/switch.jl")
include("gamma_mixture/a.jl")
include("gamma_mixture/b.jl")
include("gamma_mixture/out.jl")

include("dot_product/marginals.jl")
include("dot_product/out.jl")
include("dot_product/in1.jl")
include("dot_product/in2.jl")

include("transition/marginals.jl")
include("transition/out.jl")
include("transition/in.jl")
include("transition/a.jl")

include("autoregressive/y.jl")
include("autoregressive/x.jl")
include("autoregressive/theta.jl")
include("autoregressive/gamma.jl")
include("autoregressive/marginals.jl")

include("probit/marginals.jl")
include("probit/in.jl")
include("probit/out.jl")

include("flow/marginals.jl")
include("flow/in.jl")
include("flow/out.jl")

include("bifm/marginals.jl")
include("bifm/in.jl")
include("bifm/out.jl")
include("bifm/znext.jl")
include("bifm/zprev.jl")

include("bifm_helper/in.jl")
include("bifm_helper/out.jl")