text:
# invarianttensor.jl
#
# ITensor provides a dense implementation of an AbstractTensor type living
# in an InvariantSpace, i.e. the invariant subspace of the tensor product
# of its index spaces.
# Currenlty only defined for abelian sectors.

#+++++++++++++++++++++++
# InvariantTensor type:
#+++++++++++++++++++++++
# Type definition and constructors:
#-----------------------------------
immutable InvariantTensor{G,S,T,N} <: AbstractTensor{S,InvariantSpace,T,N}
    data::Vector{T}
    space::InvariantSpace{G,S,N}
    _datasectors::Dict{NTuple{N,G},Array{T,N}}
    function InvariantTensor(data::Array{T},space::InvariantSpace{G,S,N})
        if !(length(data)==dim(space) || ndims(data) == dim(space) == 0)
            throw(DimensionMismatch("data not of right size"))
        end
        if promote_type(T,eltype(S))!=eltype(S)
            warn("For a tensor in $(space), the entries should not be of type $(T)")
        end
        _datasectors=Dict{NTuple{N,G},Array{T,N}}()
        ind=1
        for s in sectors(space)
            dims=ntuple(n->dim(space[n],s[n]),N)
            _datasectors[s]=pointer_to_array(pointer(data,ind),dims)
            ind+=prod(dims)
        end
        return new(vec(data),space,_datasectors)
    end
end

# Show method:
#-------------
function Base.show{G,S,T,N}(io::IO,t::InvariantTensor{G,S,T,N})
    print(io," InvariantTensor ∈ $T")
    print(io,"[")
    showcompact(io,t.space.dims)
    println(io,"]:")
    sec = sectors(t)
    if length(sec) == 0
        println(io,t.data[1])
    else
        for s in sectors(t)
            println(io,"$s:")
            Base.showarray(io,t[s];header=false)
            println(io,"")
        end
    end
end

# Basic methods for characterising a tensor:
#--------------------------------------------
space(t::InvariantTensor,ind::Int)=t.space[ind]
space(t::InvariantTensor)=t.space

# General constructors
#---------------------
# with data
tensor{G,S,T,N}(data::Array{T},P::InvariantSpace{G,S,N})=InvariantTensor{G,S,T,N}(data,P)

vdim(dim::Int) = dim == 0 ? () : (dim)

# without data
tensor{T}(::Type{T},P::InvariantSpace)=tensor(Array(T,vdim(dim(P))),P)
tensor(P::InvariantSpace)=tensor(Float64,V)

Base.similar{G,S,T,N}(t::InvariantTensor{G,S},::Type{T},P::InvariantSpace{G,S,N}=space(t))=tensor(similar(t.data,T,vdim(dim(P))),P)
Base.similar{G,S,T,N}(t::InvariantTensor{G,S},::Type{T},P::ProductSpace{S,N})=similar(t,T,invariant(P))
Base.similar{G,S,T}(t::InvariantTensor{G,S},::Type{T},V::S)=similar(t,T,invariant(V))
Base.similar{G,S,N}(t::InvariantTensor{G,S},P::InvariantSpace{G,S,N}=space(t))=similar(t,eltype(t),P)
Base.similar{G,S,N}(t::InvariantTensor{G,S},P::ProductSpace{S,N})=similar(t,eltype(t),P)
Base.similar{G,S}(t::InvariantTensor{G,S},V::S)=similar(t,eltype(t),V)

Base.zero(t::InvariantTensor)=tensor(zero(t.data),space(t))

Base.zeros{T}(::Type{T},P::InvariantSpace)=tensor(zeros(T,vdim(dim(P))),P)
Base.zeros(P::InvariantSpace)=zeros(Float64,P)

Base.rand{T}(::Type{T},P::InvariantSpace)=tensor(rand(T,vdim(dim(P))),P)
Base.rand(P::InvariantSpace)=rand(Float64,P)

function Base.eye{S<:UnitaryRepresentationSpace,T}(::Type{T},::Type{InvariantSpace},V::S)
    t=zeros(T,invariant(V⊗dual(V)))
    for s in sectors(V)
        for n=1:dim(V,s)
            t[(s,conj(s))][n,n] = one(T)
        end
    end
    return t
end
Base.eye{S<:UnitaryRepresentationSpace}(::Type{InvariantSpace},V::S)=eye(Float64,InvariantSpace,V)

Base.eye{G<:Sector,S<:UnitaryRepresentationSpace,T}(::Type{T},P::InvariantSpace{G,S,2})=(P[1]==dual(P[2]) ? eye(T,InvariantSpace,P[1]) : throw(SpaceError("Cannot construct eye-tensor when second space is not the dual of the first space")))
Base.eye{G<:Sector,S<:UnitaryRepresentationSpace}(P::InvariantSpace{G,S,2})=eye(Float64,P)

# TO BE DONE
# # tensors from concatenation
# function tensorcat{S}(catind, X::Tensor{S}...)
#     catind = collect(catind)
#     isempty(catind) && error("catind should not be empty")
#     # length(unique(catdims)) != length(catdims) && error("every dimension should appear only once")

#     nargs = length(X)
#     numindX = map(numind, X)

#     all(n->(n == numindX[1]), numindX) || throw(SpaceError("all tensors should have the same number of indices for concatenation"))

#     numindC = numindX[1]
#     ncatind = setdiff(1:numindC,catind)
#     spaceCvec = Array(S, numindC)
#     for n = 1:numindC
#         spaceCvec[n] = space(X[1],n)
#     end
#     for i = 2:nargs
#         for n in catind
#             spaceCvec[n] = directsum(spaceCvec[n], space(X[i],n))
#         end
#         for n in ncatind
#             spaceCvec[n] == space(X[i],n) || throw(SpaceError("space mismatch for index $n"))
#         end
#     end
#     spaceC = ⊗(spaceCvec...)
#     typeC = mapreduce(eltype, promote_type, X)
#     dataC = zeros(typeC, map(dim,spaceC))

#     offset = zeros(Int,numindC)
#     for i=1:nargs
#         currentdims=ntuple(n->dim(space(X[i],n)),numindC)
#         currentrange=[offset[n]+(1:currentdims[n]) for n=1:numindC]
#         dataC[currentrange...] = X[i].data
#         for n in catind
#             offset[n]+=currentdims[n]
#         end
#     end
#     return tensor(dataC,spaceC)
# end

# Copy and fill tensors:
#------------------------
function Base.copy!(tdest::InvariantTensor,tsource::InvariantTensor)
    # Copies data of tensor tsource to tensor tdest if compatible
    space(tdest)==space(tsource) || throw(SpaceError())
    for s in sectors(tdest)
        copy!(tdest[s],tsource[s])
    end
    return tdest
end
Base.fill!{G,S,T}(tdest::InvariantTensor{G,S,T},value::Number)=fill!(tdest.data,convert(T,value))

# Vectorization:
#----------------
Base.vec(t::InvariantTensor)=t.data
# Convert the non-trivial degrees of freedom in a tensor to a vector to be passed to eigensolvers etc.

@generated function Base.full{G,S,T,N}(t::InvariantTensor{G,S,T,N})
    quote
        dims=@ntuple $N n->dim(space(t,n))
        a=zeros(T,dims)
        for s in sectors(t)
            @nexprs $N n->(r_{n}=to_range(s[n],space(t,n)))
            (@nref $N a r) = t[s]
        end
        return a
    end
end

# Conversion and promotion:
#---------------------------
Base.promote_rule{G,S,T1,T2,N}(::Type{InvariantTensor{G,S,T1,N}},::Type{InvariantTensor{G,S,T2,N}})=InvariantTensor{G,S,promote_type(T1,T2),N}
Base.promote_rule{G,S,T1,T2,N1,N2}(::Type{InvariantTensor{G,S,T1,N1}},::Type{InvariantTensor{G,S,T2,N2}})=InvariantTensor{G,S,promote_type(T1,T2)}
Base.promote_rule{G,S,T1,T2}(::Type{InvariantTensor{G,S,T1}},::Type{InvariantTensor{G,S,T2}})=InvariantTensor{G,S,promote_type(T1,T2)}

Base.promote_rule{G,S,T1,T2,N}(::Type{AbstractTensor{S,InvariantSpace,T1,N}},::Type{InvariantTensor{G,S,T2,N}})=AbstractTensor{S,InvariantSpace,promote_type(T1,T2),N}
Base.promote_rule{G,S,T1,T2,N1,N2}(::Type{AbstractTensor{S,InvariantSpace,T1,N1}},::Type{InvariantTensor{G,S,T2,N2}})=AbstractTensor{S,InveriantSpace,promote_type(T1,T2)}
Base.promote_rule{G,S,T1,T2}(::Type{AbstractTensor{S,InvariantSpace,T1}},::Type{InvariantTensor{G,S,T2}})=AbstractTensor{S,InvariantSpace,promote_type(T1,T2)}

Base.convert{G,S,T,N}(::Type{InvariantTensor{G,S,T,N}},t::InvariantTensor{G,S,T,N})=t
Base.convert{G,S,T1,T2,N}(::Type{InvariantTensor{G,S,T1,N}},t::InvariantTensor{G,S,T2,N})=copy!(similar(t,T1),t)
Base.convert{G,S,T}(::Type{InvariantTensor{G,S,T}},t::InvariantTensor{G,S,T})=t
Base.convert{G,S,T1,T2}(::Type{InvariantTensor{G,S,T1}},t::InvariantTensor{G,S,T2})=copy!(similar(t,T1),t)

Base.float{G,S,T<:AbstractFloat}(t::InvariantTensor{G,S,T})=t
Base.float(t::InvariantTensor)=tensor(float(t.data),space(t))

Base.real{G,S,T<:Real}(t::InvariantTensor{G,S,T})=t
Base.real(t::InvariantTensor)=tensor(real(t.data),space(t))
Base.complex{G,S,T<:Complex}(t::InvariantTensor{G,S,T})=t
Base.complex(t::InvariantTensor)=tensor(complex(t.data),space(t))

for (f,T) in ((:float32,    Float32),
              (:float64,    Float64),
              (:complex64,  Complex64),
              (:complex128, Complex128))
    @eval (Base.$f){G,S}(t::InvariantTensor{G,S,$T}) = t
    @eval (Base.$f)(t::InvariantTensor) = tensor(($f)(t.data),space(t))
end

invariant(t::InvariantTensor)=t

function invariant{G<:Abelian,T,N}(t::Tensor{AbelianSpace{G},T,N})
    for s in sectors(t)
        prod(s)==one(G) || all(t[s].==0) || throw(InexactError())
    end
    tdest=tensor(T,invariant(space(t)))
    for s in sectors(tdest)
        copy!(tdest[s],t[s])
    end
    return tdest
end

# Basic algebra:
#----------------
function Base.conj!(t1::InvariantTensor,t2::InvariantTensor)
    space(t1)==conj(space(t2)) || throw(SpaceError())
    for s in sectors(t2)
        copy!(t1[conj(s)],t2[s])
    end
    conj!(t1.data)
    return t1
end

# transpose inverts order of indices, compatible with graphical notation
function Base.transpose!(tdest::InvariantTensor,tsource::InvariantTensor)
    space(tdest)==space(tsource).' || throw(SpaceError())
    N=numind(tsource)
    for s in sectors(tsource)
        TensorOperations.tensorcopy!(tsource[s],1:N,tdest[reverse(s)],reverse(1:N))
    end
    return tdest
end
function Base.ctranspose!(tdest::InvariantTensor,tsource::InvariantTensor)
    space(tdest)==space(tsource)' || throw(SpaceError())
    N=numind(tsource)
    for s in sectors(tsource)
        TensorOperations.tensorcopy!(tsource[s],1:N,tdest[reverse(map(conj,s))],reverse(1:N))
    end
    conj!(tdest.data)
    return tdest
end

Base.scale!{G,S,T,N}(t1::InvariantTensor{G,S,T,N},t2::InvariantTensor{G,S,T,N},a::Number)=(space(t1)==space(t2) ? scale!(t1.data,t2.data,a) : throw(SpaceError());return t1)
Base.scale!{G,S,T,N}(t1::InvariantTensor{G,S,T,N},a::Number,t2::InvariantTensor{G,S,T,N})=(space(t1)==space(t2) ? scale!(t1.data,a,t2.data) : throw(SpaceError());return t1)

Base.LinAlg.axpy!(a::Number,x::InvariantTensor,y::InvariantTensor)=(space(x)==space(y) ? Base.LinAlg.axpy!(a,x.data,y.data) : throw(SpaceError()); return y)

-(t::InvariantTensor)=tensor(-t.data,space(t))
+(t1::InvariantTensor,t2::InvariantTensor)= space(t1)==space(t2) ? tensor(t1.data+t2.data,space(t1)) : throw(SpaceError())
-(t1::InvariantTensor,t2::InvariantTensor)= space(t1)==space(t2) ? tensor(t1.data-t2.data,space(t1)) : throw(SpaceError())

# Scalar product and norm: only implemented for AbelianSpace
Base.dot{G<:Abelian}(t1::InvariantTensor{G,AbelianSpace{G}},t2::InvariantTensor{G,AbelianSpace{G}})= (space(t1)==space(t2) ? dot(vec(t1),vec(t2)) : throw(SpaceError()))
Base.vecnorm{G<:Abelian}(t::InvariantTensor{G,AbelianSpace{G}})=vecnorm(t.data) # frobenius norm

# Indexing
#----------
# using sectors
Base.getindex{G,S,T,N}(t::InvariantTensor{G,S,T,N},s::NTuple{N,G})=t._datasectors[s]
Base.setindex!{G,S,T,N}(t::InvariantTensor{G,S,T,N},v::Array,s::NTuple{N,G})=(length(v)==length(t[s]) ? copy!(t[s],v) : throw(DimensionMismatch()))
Base.setindex!{G,S,T,N}(t::InvariantTensor{G,S,T,N},v::Number,s::NTuple{N,G})=fill!(t[s],v)

# Tensor Operations
#-------------------
TensorOperations.scalar(t::InvariantTensor)=iscnumber(space(t)) ? t.data[1] : throw(SpaceError("Not a scalar"))

function TensorOperations.add!{CA}(α, A::InvariantTensor, ::Type{Val{CA}}, β, C::InvariantTensor, indCinA)
    # Implements C = β*C + α*permute(op(A))
    NA = numind(A)

    spaceA = CA == :C ? conj(space(A)) : space(A)

    for i = 1:NA
        spaceA[indCinA[i]] == space(C,i) || throw(SpaceError("incompatible index spaces of tensors"))
    end

    if NA == 0 #scalars
        C.data[1] = β*C.data[1] + α*A.data[1]
    elseif indCinA == collect(1:NA) #trivial permutation
        if β == 0
            scale!(copy!(C,A), α)
        else
            Base.LinAlg.axpy!(α, A, scale!(C, β))
        end
    else
        for s in sectors(A)
            TensorOperations.add!(α, A[s], Val{CA}, β, C[s[indCinA]], indCinA)
        end
    end

    return C
end

function TensorOperations.trace!{CA}(α, A::InvariantTensor, ::Type{Val{CA}}, β, C::InvariantTensor, indCinA, cindA1, cindA2)
    NA = numind(A)
    NC = numind(C)

    spaceA = CA == :C ? conj(space(A)) : space(A)

    for i = 1:NC
        spaceA[indCinA[i]] == space(C,i) || throw(SpaceError("space mismatch"))
    end

    for i = 1:div(NA-NC, 2)
        spaceA[cindA1[i]] == dual(spaceA[cindA2[i]]) || throw(SpaceError("space mismatch"))
    end

    if length(indCinA) == 0
        Cscal = pointer_to_array(pointer(C.data, 1), ())
        for sA in sectors(A)
            if sA[cindA1] != conj(sA[cindA2])
                continue
            end
            TensorOperations.trace!(α, A[sA], Val{CA}, β, Cscal, indCinA, cindA1, cindA2)
            β = one(β)
        end
    else
        βdict = [s=>β for s in sectors(C)]
        for sA in sectors(A)
            if sA[cindA1] != conj(sA[cindA2])
                continue
            end
            sC = sA[indCinA]
            TensorOperations.trace!(α, A[sA], Val{CA}, βdict[sC], C[sC], indCinA, cindA1, cindA2)
            βdict[sC] = one(β)
        end
    end

    return C
end

function TensorOperations.contract!{CA,CB,ME}(α, A::InvariantTensor, ::Type{Val{CA}}, B::InvariantTensor, ::Type{Val{CB}}, β, C::InvariantTensor, oindA, cindA, oindB, cindB, indCinoAB, ::Type{Val{ME}}=Val{:BLAS})
    # check size compatibility
    spaceA = CA == :C ? conj(space(A)) : space(A)
    spaceB = CB == :C ? conj(space(B)) : space(B)
    spaceC = space(C)

    cspaceA = spaceA[cindA]
    cspaceB = spaceB[cindB]

    ospaceA = spaceA[oindA]
    ospaceB = spaceB[oindB]

    ospaceAB = ospaceA ⊗ ospaceB

    for i = 1:length(cspaceA)
        cspaceA[i] == dual(cspaceB[i]) || throw(SpaceError("A and B have incompatible index spaces"))
    end

    for i in 1:length(indCinoAB)
        spaceC[i] == ospaceAB[indCinoAB[i]] || throw(SpaceError("C has incompatible index space"))
    end

    if length(indCinoAB) == 0
        Cscal = pointer_to_array(pointer(C.data, 1), ())
        for sA in sectors(spaceA), sB in sectors(spaceB)
            if sA[cindA] == conj(sB)[cindB]
                sA_ = (CA==:C?conj(sA):sA)
                sB_ = (CB==:C?conj(sB):sB)
                TensorOperations.contract!(α, A[sA_], Val{CA}, B[sB_], Val{CB}, β, Cscal, oindA, cindA, oindB, cindB, indCinoAB, Val{ME})
                β = one(β)
            end
        end
    else
        βdict = [s=>β for s in sectors(C)]
        for sA in sectors(spaceA), sB in sectors(spaceB)
            if sA[cindA] == conj(sB)[cindB]
                sC = tuple(sA[oindA]..., sB[oindB]...)[indCinoAB]
                sA_ = (CA == :C ? conj(sA) : sA) #undo conjugation for indexing of original tensors
                sB_ = (CB == :C ? conj(sB) : sB)
                TensorOperations.contract!(α, A[sA_], Val{CA}, B[sB_], Val{CB}, βdict[sC], C[sC], oindA, cindA, oindB, cindB, indCinoAB, Val{ME})
                βdict[sC] = one(β)
            end
        end
    end

    return C
end

function TensorOperations.similar_from_indices{T,CA}(::Type{T}, indices, A::InvariantTensor, ::Type{Val{CA}}=Val{:N})
    spaceA = (CA == :C ? conj(space(A)) : space(A))
    return similar(A, T, spaceA[indices])
end

function TensorOperations.similar_from_indices{T,CA,CB}(::Type{T}, indices, A::InvariantTensor, B::InvariantTensor, ::Type{Val{CA}}=Val{:N}, ::Type{Val{CB}}=Val{:N})
    spaceA = (CA == :C ? conj(space(A)) : space(A))
    spaceB = (CB == :C ? conj(space(B)) : space(B))
    spaceAB = spaceA ⊗ spaceB
    return similar(A, T, spaceAB[indices])
end

# Index methods
#---------------
@eval function insertind{G,S}(t::InvariantTensor{G,S},ind::Int,V::S)
    N=numind(t)
    0<=ind<=N || throw(IndexError("index out of range"))
    iscnumber(V) || throw(SpaceError("can only insert index with c-number index space"))
    spacet=space(t)
    newspace=spacet[1:ind] ⊗ V ⊗ spacet[ind+1:N]
    tdest=similar(t,newspace)
    for s in sectors(tdest)
        st=tuple(s[1:ind]...,s[ind+2:N+1]...)
        copy!(tdest[s],t[st])
    end
    return tdest
end
@eval function deleteind(t::InvariantTensor,ind::Int)
    1<=ind<=numind(t) || throw(IndexError("index out of range"))
    iscnumber(space(t,ind)) || throw(SpaceError("can only delete index with c-number index space"))
    spacet=space(t)
    newspace=spacet[1:ind-1] ⊗ spacet[ind+1:N]
    tdest=similar(t,newspace)
    for s in sectors(tsource)
        st=tuple(s[1:ind-1]...,s[ind+1:N]...)
        copy!(tdest[st],t[s])
    end
    return tensor(t.data,newspace)
end

# Methods below are only implemented for AbelianTensor:
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++
typealias AbelianTensor{G<:Abelian,T,N} InvariantTensor{G,AbelianSpace{G},T,N}

# for (S,TT) in ((CartesianSpace,CartesianTensor),(ComplexSpace,ComplexTensor))
#     @eval function fuseind(t::$TT,ind1::Int,ind2::Int,V::$S)
#         N=numind(t)
#         ind2==ind1+1 || throw(IndexError("only neighbouring indices can be fused"))
#         1<=ind1<=N-1 || throw(IndexError("index out of range"))
#         fuse(space(t,ind1),space(t,ind2),V) || throw(SpaceError("index spaces $(space(t,ind1)) and $(space(t,ind2)) cannot be fused to $V"))
#         spacet=space(t)
#         newspace=spacet[1:ind1-1]*V*spacet[ind2+1:N]
#         return tensor(t.data,newspace)
#     end
#     @eval function splitind(t::$TT,ind::Int,V1::$S,V2::$S)
#         1<=ind<=numind(t) || throw(IndexError("index out of range"))
#         fuse(V1,V2,space(t,ind)) || throw(SpaceError("index space $(space(t,ind)) cannot be split into $V1 and $V2"))
#         spacet=space(t)
#         newspace=spacet[1:ind-1]*V1*V2*spacet[ind+1:N]
#         return tensor(t.data,newspace)
#     end
# end

# Factorizations:
#-----------------
function _reorderdata!{G<:Abelian,T,N}(t::AbelianTensor{G,T,N},tmp::Vector{T},n::Int)
    # Reorders data in t, using a temporary array of size dim(t), as preparation for
    # a factorization according to a cut between index n and index n+1. The resulting
    # tensor t is no longer safe to use, as t will no longer be equal to tensor(vec(t),space(t)).

    spacet=space(t)
    leftspace=spacet[1:n]
    rightspace=spacet[n+1:N]

    length(tmp)==dim(spacet) || throw(DimensionMismatch())

    t2=tensor(tmp,space(t))
    copy!(t2,t)

    leftdims=Dict{NTuple{n,G},Int}()
    rightdims=Dict{NTuple{N-n,G},Int}()
    for s in sectors(t2)
        ls=s[1:n]
        rs=s[n+1:N]
        leftdims[ls]=dim(leftspace,ls)
        rightdims[rs]=dim(rightspace,rs)
    end

    # Determine sizes of blocks of constant fusing charge
    blockleftdims=Dict{G,Int}()
    blockrightdims=Dict{G,Int}()
    leftoffsets=Dict{NTuple{n,G},Int}()
    rightoffsets=Dict{NTuple{N-n,G},Int}()
    for s in keys(leftdims)
        c=prod(s)
        leftoffsets[s]=get(blockleftdims,c,0)
        blockleftdims[c]=get(blockleftdims,c,0)+leftdims[s]
    end
    for s in keys(rightdims)
        c=conj(prod(s))
        rightoffsets[s]=get(blockrightdims,c,0)
        blockrightdims[c]=get(blockrightdims,c,0)+rightdims[s]
    end

    # Build new blocks in t
    blocks=Dict{G,Array{T,2}}()
    offset=0
    for c in keys(blockleftdims)
        blocks[c]=pointer_to_array(pointer(t.data,offset+1),(blockleftdims[c],blockrightdims[c]))
        offset+=length(blocks[c])
    end
    for s in sectors(t2)
        ls=s[1:n]
        rs=s[n+1:N]
        c=prod(ls)
        lo=leftoffsets[ls]
        ro=rightoffsets[rs]
        ld=leftdims[ls]
        rd=rightdims[rs]

        src=reshape(t2[s],(ld,rd))
        copy!(blocks[c],lo+(1:ld),ro+(1:rd),src,1:ld,1:rd)
    end
    return blocks, blockleftdims, blockrightdims, leftoffsets, rightoffsets, leftdims, rightdims
end

function svd!{G<:Abelian,T,N}(t::AbelianTensor{G,T,N},n::Int,::NoTruncation)
    # Perform singular value decomposition corresponding to bipartion of the
    # tensor indices into the left indices 1:n and remaining right indices,
    # thereby destroying the original tensor.

    spacet=space(t)
    leftspace=spacet[1:n]
    rightspace=spacet[n+1:N]
    truncerr=abs(zero(T))
    ST=typeof(truncerr)

    # prepare data
    tmp=similar(vec(t))
    blocks, blockleftdims, blockrightdims, leftoffsets, rightoffsets, leftdims, rightdims = _reorderdata!(t,tmp,n)

    # perform singular value decomposition
    dims=Dict{G,Int}()
    facts=Dict{G,Base.LinAlg.SVD{T,ST}}()
    for c in keys(blocks)
        F=svdfact!(blocks[c])
        facts[c]=F
        dims[c]=length(F[:S])
    end
    newspace=AbelianSpace(dims)

    # bring output in tensor form
    spaceU=invariant(leftspace*newspace')
    spaceV=invariant(newspace*rightspace)
    spaceS=invariant(newspace*newspace')
    if abs(dim(spaceU)-length(tmp)) < abs(dim(spaceV)-length(tmp))
        dataU=resize!(tmp,dim(spaceU))
        dataV=Array(T,dim(spaceV))
    else
        dataU=Array(T,dim(spaceU))
        dataV=resize!(tmp,dim(spaceV))
    end
    dataS=zeros(ST,dim(spaceS))
    U=tensor(dataU,spaceU)
    V=tensor(dataV,spaceV)
    S=tensor(dataS,spaceS)

    for c in sectors(newspace)
        Sblock=S[tuple(c,conj(c))]
        F=facts[c]
        for i=1:dims[c]
            Sblock[i,i]=F[:S][i]
        end
    end
    for s in sectors(t)
        ls=s[1:n]
        rs=s[n+1:N]
        c=prod(ls)
        lo=leftoffsets[ls]
        ro=rightoffsets[rs]
        ld=leftdims[ls]
        rd=rightdims[rs]
        d=dims[c]

        if d>0
            Ublock=reshape(U[tuple(ls...,conj(c))],(ld,d))
            Vblock=reshape(V[tuple(c,rs...)],(d,rd))

            F=facts[c]
            copy!(Ublock,1:ld,1:d,F[:U],lo+(1:ld),1:d)
            copy!(Vblock,1:d,1:rd,F[:Vt],1:d,ro+(1:rd))
        end
    end

    return U,S,V,truncerr
end

function svd!{G<:Abelian,T,N}(t::AbelianTensor{G,T,N},n::Int,trunc::TruncationDimension)
    # Perform singular value decomposition corresponding to bipartion of the
    # tensor indices into the left indices 1:n and remaining right indices,
    # thereby destroying the original tensor.
    spacet=space(t)
    leftspace=spacet[1:n]
    rightspace=spacet[n+1:N]
    truncerr=abs(zero(T))
    ST=typeof(truncerr)

    tmp=similar(vec(t))
    blocks, blockleftdims, blockrightdims, leftoffsets, rightoffsets, leftdims, rightdims = _reorderdata!(t,tmp,n)

    # perform singular value decomposition
    dims=Dict{G,Int}()
    facts=Dict{G,Base.LinAlg.SVD{T,ST}}()
    for c in keys(blocks)
        F=svdfact!(blocks[c])
        facts[c]=F
        dims[c]=length(F[:S])
    end
    newspace=AbelianSpace(dims)

    # perform truncation
    sing=Array(ST,sum(values(dims)))
    ind=1
    for c in keys(facts)
        F=facts[c]
        for sigma in F[:S]
            sing[ind]=sigma
            ind+=1
        end
    end
    sing=sort!(sing,rev=true)
    numsing=length(sing)
    normsing=vecnorm(sing)
    if eps(trunc)!=0
        truncdim=0
        while truncdim<length(sing) && vecnorm(slice(sing,(truncdim+1):numsing))>eps(trunc)*normsing
            truncdim+=1
        end
        truncdim=min(truncdim,dim(trunc))
    else
        truncdim=dim(trunc)
    end
    if truncdim<length(sing)
        truncerr=vecnorm(slice(sing,(truncdim+1):numsing))/normsing
        trunctol=sqrt(sing[truncdim]*sing[truncdim+1])
    else
        trunctol=zero(ST)
    end
    for c in keys(facts)
        F=facts[c]
        dims[c]=sum(F[:S].>trunctol)
    end
    newspace=AbelianSpace(dims)

    # bring output in tensor form
    spaceU=invariant(leftspace*newspace')
    spaceV=invariant(newspace*rightspace)
    spaceS=invariant(newspace*newspace')
    if abs(dim(spaceU)-length(tmp)) < abs(dim(spaceV)-length(tmp))
        dataU=resize!(tmp,dim(spaceU))
        dataV=Array(T,dim(spaceV))
    else
        dataU=Array(T,dim(spaceU))
        dataV=resize!(tmp,dim(spaceV))
    end
    dataS=zeros(ST,dim(spaceS))
    U=tensor(dataU,spaceU)
    V=tensor(dataV,spaceV)
    S=tensor(dataS,spaceS)

    for c in sectors(newspace)
        Sblock=S[tuple(c,conj(c))]
        F=facts[c]
        for i=1:dims[c]
            Sblock[i,i]=F[:S][i]
        end
    end
    for s in sectors(t)
        ls=s[1:n]
        rs=s[n+1:N]
        c=prod(ls)
        lo=leftoffsets[ls]
        ro=rightoffsets[rs]
        ld=leftdims[ls]
        rd=rightdims[rs]
        d=dims[c]

        if d>0
            Ublock=reshape(U[tuple(ls...,conj(c))],(ld,d))
            Vblock=reshape(V[tuple(c,rs...)],(d,rd))

            F=facts[c]
            copy!(Ublock,1:ld,1:d,F[:U],lo+(1:ld),1:d)
            copy!(Vblock,1:d,1:rd,F[:Vt],1:d,ro+(1:rd))
        end
    end

    return U,S,V,truncerr
end

function svd!{G<:Abelian,T,N}(t::AbelianTensor{G,T,N},n::Int,trunc::TruncationSpace)
    # Perform singular value decomposition corresponding to bipartion of the
    # tensor indices into the left indices 1:n and remaining right indices,
    # thereby destroying the original tensor.
    spacet=space(t)
    leftspace=spacet[1:n]
    rightspace=spacet[n+1:N]
    truncerr=abs(zero(T))
    ST=typeof(truncerr)

    tmp=similar(vec(t))
    blocks, blockleftdims, blockrightdims, leftoffsets, rightoffsets, leftdims, rightdims = _reorderdata!(t,tmp,n)

    # perform singular value decomposition
    dims=Dict{G,Int}()
    facts=Dict{G,Base.LinAlg.SVD{T,ST}}()
    for c in keys(blocks)
        F=svdfact!(blocks[c])
        facts[c]=F
        dims[c]=length(F[:S])
    end
    newspace=AbelianSpace(dims)

    # perform truncation
    sing=Array(ST,sum(values(dims)))
    ind=1
    for c in keys(facts)
        F=facts[c]
        for sigma in F[:S]
            sing[ind]=sigma
            ind+=1
        end
    end
    sing=sort!(sing,rev=true)
    numsing=length(sing)
    normsing=vecnorm(sing)
    if eps(trunc)!=0
        truncdim=0
        while truncdim<length(sing) && vecnorm(slice(sing,(truncdim+1):numsing))>eps(trunc)*normsing
            truncdim+=1
        end
        if truncdim<length(sing)
            trunctol=sqrt(sing[truncdim]*sing[truncdim+1])
        else
            trunctol=zero(ST)
        end
        for c in keys(facts)
            F=facts[c]
            dims[c]=min(sum(F[:S].>=trunctol),dim(space(trunc),c))
        end
    else
        for c in keys(facts)
            F=facts[c]
            dims[c]=min(dims[c],dim(space(trunc),c))
        end
    end
    newspace=AbelianSpace(dims)
    fill!(sing,0)
    ind=1
    for c in keys(facts)
        F=facts[c]
        for i=1:dims[c]
            sing[ind]=F[:S][i]
            ind+=1
        end
    end
    truncerr=vecnorm(sing)/normsing

    # bring output in tensor form
    spaceU=invariant(leftspace*newspace')
    spaceV=invariant(newspace*rightspace)
    spaceS=invariant(newspace*newspace')
    if abs(dim(spaceU)-length(tmp)) < abs(dim(spaceV)-length(tmp))
        dataU=resize!(tmp,dim(spaceU))
        dataV=Array(T,dim(spaceV))
    else
        dataU=Array(T,dim(spaceU))
        dataV=resize!(tmp,dim(spaceV))
    end
    dataS=zeros(ST,dim(spaceS))
    U=tensor(dataU,spaceU)
    V=tensor(dataV,spaceV)
    S=tensor(dataS,spaceS)

    for c in sectors(newspace)
        Sblock=S[tuple(c,conj(c))]
        F=facts[c]
        for i=1:dims[c]
            Sblock[i,i]=F[:S][i]
        end
    end
    for s in sectors(t)
        ls=s[1:n]
        rs=s[n+1:N]
        c=prod(ls)
        lo=leftoffsets[ls]
        ro=rightoffsets[rs]
        ld=leftdims[ls]
        rd=rightdims[rs]
        d=dims[c]

        if d>0
            Ublock=reshape(U[tuple(ls...,conj(c))],(ld,d))
            Vblock=reshape(V[tuple(c,rs...)],(d,rd))

            F=facts[c]
            copy!(Ublock,1:ld,1:d,F[:U],lo+(1:ld),1:d)
            copy!(Vblock,1:d,1:rd,F[:Vt],1:d,ro+(1:rd))
        end
    end

    return U,S,V,truncerr
end

function svd!{G<:Abelian,T,N}(t::AbelianTensor{G,T,N},n::Int,trunc::TruncationError)
    # Perform singular value decomposition corresponding to bipartion of the
    # tensor indices into the left indices 1:n and remaining right indices,
    # thereby destroying the original tensor.
    spacet=space(t)
    leftspace=spacet[1:n]
    rightspace=spacet[n+1:N]
    truncerr=abs(zero(T))
    ST=typeof(truncerr)

    tmp=similar(vec(t))
    blocks, blockleftdims, blockrightdims, leftoffsets, rightoffsets, leftdims, rightdims = _reorderdata!(t,tmp,n)

    # perform singular value decomposition
    dims=Dict{G,Int}()
    facts=Dict{G,Base.LinAlg.SVD{T,ST}}()
    for c in keys(blocks)
        F=svdfact!(blocks[c])
        facts[c]=F
        dims[c]=length(F[:S])
    end
    newspace=AbelianSpace(dims)

    # perform truncation
    sing=Array(ST,sum(values(dims)))
    ind=1
    for c in keys(facts)
        F=facts[c]
        for sigma in F[:S]
            sing[ind]=sigma
            ind+=1
        end
    end
    sing=sort!(sing,rev=true)
    normsing=vecnorm(sing)
    numsing=length(sing)
    truncdim=0
    while truncdim<length(sing) && vecnorm(slice(sing,(truncdim+1):numsing))>eps(trunc)*normsing
        truncdim+=1
    end
    if truncdim<length(sing)
        truncerr=vecnorm(slice(sing,(truncdim+1):numsing))/normsing
        trunctol=sqrt(sing[truncdim]*sing[truncdim+1])
    else
        trunctol=zero(ST)
    end
    for c in keys(facts)
        F=facts[c]
        dims[c]=sum(F[:S].>=trunctol)
    end
    newspace=AbelianSpace(dims)

    # bring output in tensor form
    spaceU=invariant(leftspace*newspace')
    spaceV=invariant(newspace*rightspace)
    spaceS=invariant(newspace*newspace')
    if abs(dim(spaceU)-length(tmp)) < abs(dim(spaceV)-length(tmp))
        dataU=resize!(tmp,dim(spaceU))
        dataV=Array(T,dim(spaceV))
    else
        dataU=Array(T,dim(spaceU))
        dataV=resize!(tmp,dim(spaceV))
    end
    dataS=zeros(ST,dim(spaceS))
    U=tensor(dataU,spaceU)
    V=tensor(dataV,spaceV)
    S=tensor(dataS,spaceS)

    for c in sectors(newspace)
        Sblock=S[tuple(c,conj(c))]
        F=facts[c]
        for i=1:dims[c]
            Sblock[i,i]=F[:S][i]
        end
    end
    for s in sectors(t)
        ls=s[1:n]
        rs=s[n+1:N]
        c=prod(ls)
        lo=leftoffsets[ls]
        ro=rightoffsets[rs]
        ld=leftdims[ls]
        rd=rightdims[rs]
        d=dims[c]

        if d>0
            Ublock=reshape(U[tuple(ls...,conj(c))],(ld,d))
            Vblock=reshape(V[tuple(c,rs...)],(d,rd))

            F=facts[c]
            copy!(Ublock,1:ld,1:d,F[:U],lo+(1:ld),1:d)
            copy!(Vblock,1:d,1:rd,F[:Vt],1:d,ro+(1:rd))
        end
    end

    return U,S,V,truncerr
end

function leftorth!{G<:Abelian,T,N}(t::AbelianTensor{G,T,N},n::Int,::NoTruncation)
    # Perform singular value decomposition corresponding to bipartion of the
    # tensor indices into the left indices 1:n and remaining right indices,
    # thereby destroying the original tensor.
    spacet=space(t)
    leftspace=spacet[1:n]
    rightspace=spacet[n+1:N]
    truncerr=abs(zero(T))
    ST=typeof(truncerr)

    # prepare data
    tmp=similar(vec(t))
    blocks, blockleftdims, blockrightdims, leftoffsets, rightoffsets, leftdims, rightdims = _reorderdata!(t,tmp,n)

    # perform singular value decomposition
    dims=Dict{G,Int}()
    Us=Dict{G,Array{T,2}}()
    Cs=Dict{G,Array{T,2}}()
    for c in keys(blocks)
        if blockleftdims[c]>blockrightdims[c]
            F=qrfact!(blocks[c])
            Us[c]=full(F[:Q])
            Cs[c]=full(F[:R])
        else
            Us[c]=eye(T,blockleftdims[c])
            Cs[c]=blocks[c]
        end
        dims[c]=min(blockleftdims[c],blockrightdims[c])
    end
    newspace=AbelianSpace(dims)

    # bring output in tensor form
    spaceU=invariant(leftspace*newspace')
    spaceC=invariant(newspace*rightspace)
    dataU=resize!(tmp,dim(spaceU))
    dataC=Array(T,dim(spaceC))
    U=tensor(dataU,spaceU)
    C=tensor(dataC,spaceC)

    for s in sectors(t)
        ls=s[1:n]
        rs=s[n+1:N]
        c=prod(ls)
        lo=leftoffsets[ls]
        ro=rightoffsets[rs]
        ld=leftdims[ls]
        rd=rightdims[rs]
        d=dims[c]

        if d>0
            Ublock=reshape(U[tuple(ls...,conj(c))],(ld,d))
            Cblock=reshape(C[tuple(c,rs...)],(d,rd))

            copy!(Ublock,1:ld,1:d,Us[c],lo+(1:ld),1:d)
            copy!(Cblock,1:d,1:rd,Cs[c],1:d,ro+(1:rd))
        end
    end

    return U,C,truncerr
end

function leftorth!{G<:Abelian,T,N}(t::AbelianTensor{G,T,N},n::Int,trunc::TruncationDimension)
    # Perform singular value decomposition corresponding to bipartion of the
    # tensor indices into the left indices 1:n and remaining right indices,
    # thereby destroying the original tensor.
    spacet=space(t)
    leftspace=spacet[1:n]
    rightspace=spacet[n+1:N]
    truncerr=abs(zero(T))
    ST=typeof(truncerr)

    tmp=similar(vec(t))
    blocks, blockleftdims, blockrightdims, leftoffsets, rightoffsets, leftdims, rightdims = _reorderdata!(t,tmp,n)

    # perform singular value decomposition
    dims=Dict{G,Int}()
    facts=Dict{G,Base.LinAlg.SVD{T,ST}}()
    Us=Dict{G,Array{T,2}}()
    Cs=Dict{G,Array{T,2}}()
    for c in keys(blocks)
        F=svdfact!(blocks[c])
        facts[c]=F
        dims[c]=length(F[:S])
        Us[c]=F[:U]
        Cs[c]=scale!(F[:S],F[:Vt])
    end
    newspace=AbelianSpace(dims)

    # perform truncation
    sing=Array(ST,sum(values(dims)))
    ind=1
    for c in keys(facts)
        F=facts[c]
        for sigma in F[:S]
            sing[ind]=sigma
            ind+=1
        end
    end
    sing=sort!(sing,rev=true)
    numsing=length(sing)
    normsing=vecnorm(sing)
    if eps(trunc)!=0
        truncdim=0
        while truncdim<length(sing) && vecnorm(slice(sing,(truncdim+1):numsing))>eps(trunc)*normsing
            truncdim+=1
        end
        truncdim=min(truncdim,dim(trunc))
    else
        truncdim=dim(trunc)
    end
    if truncdim<length(sing)
        truncerr=vecnorm(slice(sing,(truncdim+1):numsing))/normsing
        trunctol=sqrt(sing[truncdim]*sing[truncdim+1])
    else
        trunctol=zero(ST)
    end
    for c in keys(facts)
        F=facts[c]
        dims[c]=sum(F[:S].>trunctol)
    end
    newspace=AbelianSpace(dims)

    # bring output in tensor form
    spaceU=invariant(leftspace*newspace')
    spaceC=invariant(newspace*rightspace)
    dataU=resize!(tmp,dim(spaceU))
    dataC=Array(T,dim(spaceC))
    U=tensor(dataU,spaceU)
    C=tensor(dataC,spaceC)

    for s in sectors(t)
        ls=s[1:n]
        rs=s[n+1:N]
        c=prod(ls)
        lo=leftoffsets[ls]
        ro=rightoffsets[rs]
        ld=leftdims[ls]
        rd=rightdims[rs]
        d=dims[c]

        if d>0
            Ublock=reshape(U[tuple(ls...,conj(c))],(ld,d))
            Cblock=reshape(C[tuple(c,rs...)],(d,rd))

            copy!(Ublock,1:ld,1:d,Us[c],lo+(1:ld),1:d)
            copy!(Cblock,1:d,1:rd,Cs[c],1:d,ro+(1:rd))
        end
    end

    return U,C,truncerr
end

function leftorth!{G<:Abelian,T,N}(t::AbelianTensor{G,T,N},n::Int,trunc::TruncationSpace)
    # Perform singular value decomposition corresponding to bipartion of the
    # tensor indices into the left indices 1:n and remaining right indices,
    # thereby destroying the original tensor.
    spacet=space(t)
    leftspace=spacet[1:n]
    rightspace=spacet[n+1:N]
    truncerr=abs(zero(T))
    ST=typeof(truncerr)

    tmp=similar(vec(t))
    blocks, blockleftdims, blockrightdims, leftoffsets, rightoffsets, leftdims, rightdims = _reorderdata!(t,tmp,n)

    # perform singular value decomposition
    dims=Dict{G,Int}()
    facts=Dict{G,Base.LinAlg.SVD{T,ST}}()
    Us=Dict{G,Array{T,2}}()
    Cs=Dict{G,Array{T,2}}()
    for c in keys(blocks)
        F=svdfact!(blocks[c])
        facts[c]=F
        dims[c]=length(F[:S])
        Us[c]=F[:U]
        Cs[c]=scale!(F[:S],F[:Vt])
    end
    newspace=AbelianSpace(dims)

    # perform truncation
    sing=Array(ST,sum(values(dims)))
    ind=1
    for c in keys(facts)
        F=facts[c]
        for sigma in F[:S]
            sing[ind]=sigma
            ind+=1
        end
    end
    sing=sort!(sing,rev=true)
    numsing=length(sing)
    normsing=vecnorm(sing)
    if eps(trunc)!=0
        truncdim=0
        while truncdim<length(sing) && vecnorm(slice(sing,(truncdim+1):numsing))>eps(trunc)*normsing
            truncdim+=1
        end
        if truncdim<length(sing)
            trunctol=sqrt(sing[truncdim]*sing[truncdim+1])
        else
            trunctol=zero(ST)
        end
        for c in keys(facts)
            F=facts[c]
            dims[c]=min(sum(F[:S].>=trunctol),dim(space(trunc),c))
        end
    else
        for c in keys(facts)
            F=facts[c]
            dims[c]=min(dims[c],dim(space(trunc),c))
        end
    end
    newspace=AbelianSpace(dims)
    fill!(sing,0)
    ind=1
    for c in keys(facts)
        F=facts[c]
        for i=1:dims[c]
            sing[ind]=F[:S][i]
            ind+=1
        end
    end
    truncerr=vecnorm(sing)/normsing

    # bring output in tensor form
    spaceU=invariant(leftspace*newspace')
    spaceC=invariant(newspace*rightspace)
    dataU=resize!(tmp,dim(spaceU))
    dataC=Array(T,dim(spaceC))
    U=tensor(dataU,spaceU)
    C=tensor(dataC,spaceC)

    for s in sectors(t)
        ls=s[1:n]
        rs=s[n+1:N]
        c=prod(ls)
        lo=leftoffsets[ls]
        ro=rightoffsets[rs]
        ld=leftdims[ls]
        rd=rightdims[rs]
        d=dims[c]

        if d>0
            Ublock=reshape(U[tuple(ls...,conj(c))],(ld,d))
            Cblock=reshape(C[tuple(c,rs...)],(d,rd))

            copy!(Ublock,1:ld,1:d,Us[c],lo+(1:ld),1:d)
            copy!(Cblock,1:d,1:rd,Cs[c],1:d,ro+(1:rd))
        end
    end

    return U,C,truncerr
end

function leftorth!{G<:Abelian,T,N}(t::AbelianTensor{G,T,N},n::Int,trunc::TruncationError)
    # Perform singular value decomposition corresponding to bipartion of the
    # tensor indices into the left indices 1:n and remaining right indices,
    # thereby destroying the original tensor.
    spacet=space(t)
    leftspace=spacet[1:n]
    rightspace=spacet[n+1:N]
    truncerr=abs(zero(T))
    ST=typeof(truncerr)

    tmp=similar(vec(t))
    blocks, blockleftdims, blockrightdims, leftoffsets, rightoffsets, leftdims, rightdims = _reorderdata!(t,tmp,n)

    # perform singular value decomposition
    dims=Dict{G,Int}()
    facts=Dict{G,Base.LinAlg.SVD{T,ST}}()
    Us=Dict{G,Array{T,2}}()
    Cs=Dict{G,Array{T,2}}()
    for c in keys(blocks)
        F=svdfact!(blocks[c])
        facts[c]=F
        dims[c]=length(F[:S])
        Us[c]=F[:U]
        Cs[c]=scale!(F[:S],F[:Vt])
    end
    newspace=AbelianSpace(dims)

    # perform truncation
    sing=Array(ST,sum(values(dims)))
    ind=1
    for c in keys(facts)
        F=facts[c]
        for sigma in F[:S]
            sing[ind]=sigma
            ind+=1
        end
    end
    sing=sort!(sing,rev=true)
    numsing=length(sing)
    normsing=vecnorm(sing)
    truncdim=0
    while truncdim<length(sing) && vecnorm(slice(sing,(truncdim+1):numsing))>eps(trunc)*normsing
        truncdim+=1
    end
    if truncdim<length(sing)
        truncerr=vecnorm(slice(sing,(truncdim+1):numsing))/normsing
        trunctol=sqrt(sing[truncdim]*sing[truncdim+1])
    else
        trunctol=zero(ST)
    end
    for c in keys(facts)
        F=facts[c]
        dims[c]=sum(F[:S].>=trunctol)
    end
    newspace=AbelianSpace(dims)

    # bring output in tensor form
    spaceU=invariant(leftspace*newspace')
    spaceC=invariant(newspace*rightspace)
    dataU=resize!(tmp,dim(spaceU))
    dataC=Array(T,dim(spaceC))
    U=tensor(dataU,spaceU)
    C=tensor(dataC,spaceC)

    for s in sectors(t)
        ls=s[1:n]
        rs=s[n+1:N]
        c=prod(ls)
        lo=leftoffsets[ls]
        ro=rightoffsets[rs]
        ld=leftdims[ls]
        rd=rightdims[rs]
        d=dims[c]

        if d>0
            Ublock=reshape(U[tuple(ls...,conj(c))],(ld,d))
            Cblock=reshape(C[tuple(c,rs...)],(d,rd))

            copy!(Ublock,1:ld,1:d,Us[c],lo+(1:ld),1:d)
            copy!(Cblock,1:d,1:rd,Cs[c],1:d,ro+(1:rd))
        end
    end

    return U,C,truncerr
end

function rightorth!{G<:Abelian,T,N}(t::AbelianTensor{G,T,N},n::Int,::NoTruncation)
    # Perform singular value decomposition corresponding to bipartion of the
    # tensor indices into the left indices 1:n and remaining right indices,
    # thereby destroying the original tensor.
    spacet=space(t)
    leftspace=spacet[1:n]
    rightspace=spacet[n+1:N]
    truncerr=abs(zero(T))
    ST=typeof(truncerr)

    # prepare data
    tmp=similar(vec(t))
    t2=tensor(tmp,space(t))
    copy!(t2,t)
    blocks, blockleftdims, blockrightdims, leftoffsets, rightoffsets, leftdims, rightdims = _reorderdata!(t2,t.data,n)
    offset=0
    tblocks=Dict{G,Array{T,2}}()
    for c in keys(blockleftdims)
        tblocks[c]=pointer_to_array(pointer(t.data,offset+1),(blockrightdims[c],blockleftdims[c]))
        Base.transpose!(tblocks[c],blocks[c])
        offset+=length(tblocks[c])
    end

    # perform singular value decomposition
    dims=Dict{G,Int}()
    Cs=Dict{G,Array{T,2}}()
    Us=Dict{G,Array{T,2}}()
    for c in keys(blocks)
        if blockleftdims[c]<blockrightdims[c]
            F=qrfact!(tblocks[c])
            Us[c]=full(F[:Q])
            Cs[c]=full(F[:R])
        else
            Us[c]=eye(T,blockrightdims[c])
            Cs[c]=tblocks[c]
        end
        dims[c]=min(blockleftdims[c],blockrightdims[c])
    end
    newspace=AbelianSpace(dims)

    # bring output in tensor form
    spaceC=invariant(leftspace*newspace')
    spaceU=invariant(newspace*rightspace)
    dataC=Array(T,dim(spaceC))
    dataU=resize!(tmp,dim(spaceU))
    C=tensor(dataC,spaceC)
    U=tensor(dataU,spaceU)

    for s in sectors(t)
        ls=s[1:n]
        rs=s[n+1:N]
        c=prod(ls)
        lo=leftoffsets[ls]
        ro=rightoffsets[rs]
        ld=leftdims[ls]
        rd=rightdims[rs]
        d=dims[c]

        if d>0
            Cblock=reshape(C[tuple(ls...,conj(c))],(ld,d))
            Ublock=reshape(U[tuple(c,rs...)],(d,rd))

            Base.transpose!(Cblock,slice(Cs[c],1:d,lo+(1:ld)))
            Base.transpose!(Ublock,slice(Us[c],ro+(1:rd),1:d))
        end
    end

    return C,U,truncerr
end

function rightorth!{G<:Abelian,T,N}(t::AbelianTensor{G,T,N},n::Int,trunc::TruncationDimension)
    # Perform singular value decomposition corresponding to bipartion of the
    # tensor indices into the left indices 1:n and remaining right indices,
    # thereby destroying the original tensor.
    spacet=space(t)
    leftspace=spacet[1:n]
    rightspace=spacet[n+1:N]
    truncerr=abs(zero(T))
    ST=typeof(truncerr)

    tmp=similar(vec(t))
    blocks, blockleftdims, blockrightdims, leftoffsets, rightoffsets, leftdims, rightdims = _reorderdata!(t,tmp,n)

    # perform singular value decomposition
    dims=Dict{G,Int}()
    facts=Dict{G,Base.LinAlg.SVD{T,ST}}()
    Cs=Dict{G,Array{T,2}}()
    Us=Dict{G,Array{T,2}}()
    for c in keys(blocks)
        F=svdfact!(blocks[c])
        facts[c]=F
        dims[c]=length(F[:S])
        Cs[c]=scale!(F[:U],F[:S])
        Us[c]=F[:Vt]
    end
    newspace=AbelianSpace(dims)

    # perform truncation
    sing=Array(ST,sum(values(dims)))
    ind=1
    for c in keys(facts)
        F=facts[c]
        for sigma in F[:S]
            sing[ind]=sigma
            ind+=1
        end
    end
    sing=sort!(sing,rev=true)
    numsing=length(sing)
    normsing=vecnorm(sing)
    if eps(trunc)!=0
        truncdim=0
        while truncdim<length(sing) && vecnorm(slice(sing,(truncdim+1):numsing))>eps(trunc)*normsing
            truncdim+=1
        end
        truncdim=min(truncdim,dim(trunc))
    else
        truncdim=dim(trunc)
    end
    if truncdim<length(sing)
        truncerr=vecnorm(slice(sing,(truncdim+1):numsing))/normsing
        trunctol=sqrt(sing[truncdim]*sing[truncdim+1])
    else
        trunctol=zero(ST)
    end
    for c in keys(facts)
        F=facts[c]
        dims[c]=sum(F[:S].>trunctol)
    end
    newspace=AbelianSpace(dims)

    # bring output in tensor form
    spaceC=invariant(leftspace*newspace')
    spaceU=invariant(newspace*rightspace)
    dataC=Array(T,dim(spaceC))
    dataU=resize!(tmp,dim(spaceU))
    C=tensor(dataC,spaceC)
    U=tensor(dataU,spaceU)

    for s in sectors(t)
        ls=s[1:n]
        rs=s[n+1:N]
        c=prod(ls)
        lo=leftoffsets[ls]
        ro=rightoffsets[rs]
        ld=leftdims[ls]
        rd=rightdims[rs]
        d=dims[c]

        if d>0
            Cblock=reshape(C[tuple(ls...,conj(c))],(ld,d))
            Ublock=reshape(U[tuple(c,rs...)],(d,rd))

            copy!(Cblock,1:ld,1:d,Cs[c],lo+(1:ld),1:d)
            copy!(Ublock,1:d,1:rd,Us[c],1:d,ro+(1:rd))
        end
    end

    return C,U,truncerr
end

function rightorth!{G<:Abelian,T,N}(t::AbelianTensor{G,T,N},n::Int,trunc::TruncationSpace)
    # Perform singular value decomposition corresponding to bipartion of the
    # tensor indices into the left indices 1:n and remaining right indices,
    # thereby destroying the original tensor.
    spacet=space(t)
    leftspace=spacet[1:n]
    rightspace=spacet[n+1:N]
    truncerr=abs(zero(T))
    ST=typeof(truncerr)

    tmp=similar(vec(t))
    blocks, blockleftdims, blockrightdims, leftoffsets, rightoffsets, leftdims, rightdims = _reorderdata!(t,tmp,n)

    # perform singular value decomposition
    dims=Dict{G,Int}()
    facts=Dict{G,Base.LinAlg.SVD{T,ST}}()
    Cs=Dict{G,Array{T,2}}()
    Us=Dict{G,Array{T,2}}()
    for c in keys(blocks)
        F=svdfact!(blocks[c])
        facts[c]=F
        dims[c]=length(F[:S])
        Cs[c]=scale!(F[:U],F[:S])
        Us[c]=F[:Vt]
    end
    newspace=AbelianSpace(dims)

    # perform truncation
    sing=Array(ST,sum(values(dims)))
    ind=1
    for c in keys(facts)
        F=facts[c]
        for sigma in F[:S]
            sing[ind]=sigma
            ind+=1
        end
    end
    sing=sort!(sing,rev=true)
    numsing=length(sing)
    normsing=vecnorm(sing)
    if eps(trunc)!=0
        truncdim=0
        while truncdim<length(sing) && vecnorm(slice(sing,(truncdim+1):numsing))>eps(trunc)*normsing
            truncdim+=1
        end
        if truncdim<length(sing)
            trunctol=sqrt(sing[truncdim]*sing[truncdim+1])
        else
            trunctol=zero(ST)
        end
        for c in keys(facts)
            F=facts[c]
            dims[c]=min(sum(F[:S].>=trunctol),dim(space(trunc),c))
        end
    else
        for c in keys(facts)
            F=facts[c]
            dims[c]=min(dims[c],dim(space(trunc),c))
        end
    end
    newspace=AbelianSpace(dims)
    fill!(sing,0)
    ind=1
    for c in keys(facts)
        F=facts[c]
        for i=1:dims[c]
            sing[ind]=F[:S][i]
            ind+=1
        end
    end
    truncerr=vecnorm(sing)/normsing

    # bring output in tensor form
    spaceC=invariant(leftspace*newspace')
    spaceU=invariant(newspace*rightspace)
    dataC=Array(T,dim(spaceC))
    dataU=resize!(tmp,dim(spaceU))
    C=tensor(dataC,spaceC)
    U=tensor(dataU,spaceU)

    for s in sectors(t)
        ls=s[1:n]
        rs=s[n+1:N]
        c=prod(ls)
        lo=leftoffsets[ls]
        ro=rightoffsets[rs]
        ld=leftdims[ls]
        rd=rightdims[rs]
        d=dims[c]

        if d>0
            Cblock=reshape(C[tuple(ls...,conj(c))],(ld,d))
            Ublock=reshape(U[tuple(c,rs...)],(d,rd))

            copy!(Cblock,1:ld,1:d,Cs[c],lo+(1:ld),1:d)
            copy!(Ublock,1:d,1:rd,Us[c],1:d,ro+(1:rd))
        end
    end

    return C,U,truncerr
end

function rightorth!{G<:Abelian,T,N}(t::AbelianTensor{G,T,N},n::Int,trunc::TruncationError)
    # Perform singular value decomposition corresponding to bipartion of the
    # tensor indices into the left indices 1:n and remaining right indices,
    # thereby destroying the original tensor.
    spacet=space(t)
    leftspace=spacet[1:n]
    rightspace=spacet[n+1:N]
    truncerr=abs(zero(T))
    ST=typeof(truncerr)

    tmp=similar(vec(t))
    blocks, blockleftdims, blockrightdims, leftoffsets, rightoffsets, leftdims, rightdims = _reorderdata!(t,tmp,n)

    # perform singular value decomposition
    dims=Dict{G,Int}()
    facts=Dict{G,Base.LinAlg.SVD{T,ST}}()
    Cs=Dict{G,Array{T,2}}()
    Us=Dict{G,Array{T,2}}()
    for c in keys(blocks)
        F=svdfact!(blocks[c])
        facts[c]=F
        dims[c]=length(F[:S])
        Cs[c]=scale!(F[:U],F[:S])
        Us[c]=F[:Vt]
    end
    newspace=AbelianSpace(dims)

    # perform truncation
    sing=Array(ST,sum(values(dims)))
    ind=1
    for c in keys(facts)
        F=facts[c]
        for sigma in F[:S]
            sing[ind]=sigma
            ind+=1
        end
    end
    sing=sort!(sing,rev=true)
    numsing=length(sing)
    normsing=vecnorm(sing)
    truncdim=0
    while truncdim<length(sing) && vecnorm(slice(sing,(truncdim+1):numsing))>eps(trunc)*normsing
        truncdim+=1
    end
    if truncdim<length(sing)
        truncerr=vecnorm(slice(sing,(truncdim+1):numsing))/normsing
        trunctol=sqrt(sing[truncdim]*sing[truncdim+1])
    else
        trunctol=zero(ST)
    end
    for c in keys(facts)
        F=facts[c]
        dims[c]=sum(F[:S].>=trunctol)
    end
    newspace=AbelianSpace(dims)

    # bring output in tensor form
    spaceC=invariant(leftspace*newspace')
    spaceU=invariant(newspace*rightspace)
    dataC=Array(T,dim(spaceC))
    dataU=resize!(tmp,dim(spaceU))
    C=tensor(dataC,spaceC)
    U=tensor(dataU,spaceU)

    for s in sectors(t)
        ls=s[1:n]
        rs=s[n+1:N]
        c=prod(ls)
        lo=leftoffsets[ls]
        ro=rightoffsets[rs]
        ld=leftdims[ls]
        rd=rightdims[rs]
        d=dims[c]

        if d>0
            Cblock=reshape(C[tuple(ls...,conj(c))],(ld,d))
            Ublock=reshape(U[tuple(c,rs...)],(d,rd))

            copy!(Cblock,1:ld,1:d,Cs[c],lo+(1:ld),1:d)
            copy!(Ublock,1:d,1:rd,Us[c],1:d,ro+(1:rd))
        end
    end

    return C,U,truncerr
end

# # Methods below are only implemented for CartesianMatrix or ComplexMatrix:
# #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# typealias ComplexMatrix{T} ComplexTensor{T,2}
# typealias CartesianMatrix{T} CartesianTensor{T,2}
#
# function Base.pinv(t::Union{ComplexMatrix,CartesianMatrix})
#     # Compute pseudo-inverse
#     spacet=space(t)
#     data=copy(t.data)
#     leftdim=dim(spacet[1])
#     rightdim=dim(spacet[2])
#
#     F=svdfact!(data)
#     Sinv=F[:S]
#     for k=1:length(Sinv)
#         if Sinv[k]>eps(Sinv[1])*max(leftdim,rightdim)
#             Sinv[k]=one(Sinv[k])/Sinv[k]
#         end
#     end
#     data=F[:V]*scale(F[:S],F[:U]')
#     return tensor(data,spacet')
# end
#
# function Base.eig(t::Union{ComplexMatrix,CartesianMatrix})
#     # Compute eigenvalue decomposition.
#     spacet=space(t)
#     spacet[1] == spacet[2]' || throw(SpaceError("eigenvalue factorization only exists if left and right index space are dual"))
#     data=copy(t.data)
#
#     F=eigfact!(data)
#
#     Lambda=tensor(diagm(F[:values]),spacet)
#     V=tensor(F[:vectors],spacet)
#     return Lambda, V
# end
#
# function Base.inv(t::Union{ComplexMatrix,CartesianMatrix})
#     # Compute inverse.
#     spacet=space(t)
#     spacet[1] == spacet[2]' || throw(SpaceError("inverse only exists if left and right index space are dual"))
#
#     return tensor(inv(t.data),spacet)
# end


text:
# Defines structures and functions relating to display.

export render, ScreenCoords, ScreenPixelCoords, ScreenPixelPos, ScreenPixelDims,
        UIPixelCoords, UIPixelPos, UIPixelDims, Camera, dpiScale, worldScale,
        screenCenter, screenCenterX, screenCenterY, screenOffsetFromCenter,
        winWidth, winHeight, winWidth_highDPI, winHeight_highDPI,
        toScreenPos, toScreenPixelDims, toUIPixelPos, toUIPixelDims,
        SetRenderDrawColor, blendAlphaColors, topLeftPos, rectOrigin, render,
        renderRectCentered, renderRectFromOrigin, renderText, renderTextSurface,
        createText, sizeText, hcat_render_text, loadFont

"""
    render(obj, pos::AbstractPos, cam::Camera, renderer;
           size = nothing, src_rect = C_NULL)

Anything that can be drawn to the screen in the game should support the `render` method,
which should draw the `obj` onto the `renderer`. This will be called once each frame to
update the screen.

Note that SDL2 will handle things like drawing outside the current screen, so you don't need
to worry about checking if the object is visible inside the `cam`, etc. Simply draw the
`obj` onto the `renderer` using any available functionality.
"""
function render end

""" Coordinate system to represent the screen, not the Game World.
Importantly, these are "upside-down" from WorldCoords: (0,0) is top-left, and
higher numbers go down and to the right.
"""
abstract type ScreenCoords end

""" ScreenPixelCoords are absolute space on screen, in actual pixels. """
struct ScreenPixelCoords <: ScreenCoords end

"""
    ScreenPixelPos(1200, 1150)
Absolute position on screen, in pixels. (0,0) is top-left of screen.
"""
struct ScreenPixelPos <: AbstractPos{ScreenPixelCoords}  # 0,0 == top-left
    x::Int
    y::Int
end
ScreenPixelPos(x::Number, y::Number) = ScreenPixelPos(convert.(Int, floor.((x,y)))...)

# TODO: CONSIDER SWITCHING TO DOWN-SCALING INSTEAD of up-scaling, so you can get
# more precise sizes (1.5 "pixels" -> 3 pixels). Or consider using Floats for UIPixelCoords.
""" UIPixelCoords are space on screen, in un-dpi-scaled "pixels". """
struct UIPixelCoords <: ScreenCoords end

"""
    struct UIPixelPos <: AbstractPos{UIPixelCoords}

Position on screen, in un-dpi-scaled "pixels". (0,0) is top-left of screen.

These should be used whenever placing anything on the screen, since they are
indepdent of resolution-scaling.
"""
struct UIPixelPos <: AbstractPos{UIPixelCoords}  # 0,0 == top-left
    x::Int
    y::Int
end
UIPixelPos(x::Number, y::Number) = UIPixelPos(convert.(Int, floor.((x,y)))...)

+(a::UIPixelPos, b::UIPixelPos) = UIPixelPos(a.x+b.x, a.y+b.y)
# TODO: How to define .+ like this in Julia 1.0+?
#.+(a::UIPixelPos, x::Number) = UIPixelPos(a.x+x, a.y+x)
+(a::UIPixelPos, x::Number) = UIPixelPos(a.x+x, a.y+x)

"""
    struct ScreenPixelDims <: AbstractDims{ScreenPixelCoords}

Absolute size on screen, in pixels. Use with ScreenPixelPos.
"""
struct ScreenPixelDims <: AbstractDims{ScreenPixelCoords}
    w::Int32  # Int32 to match SDL
    h::Int32
end
"""
    struct UIPixelDims <: AbstractDims{UIPixelCoords}

Size on screen, in un-dpi-scaled "pixels". Use with UIPixelPos.
"""
struct UIPixelDims <: AbstractDims{UIPixelCoords}
    w::Int
    h::Int
end

mutable struct Camera
    pos::WorldPos
    # Note, these are in WorldPos size.
    w::Threads.Atomic{Float32}   # Note: These are Atomics, since they can be modified by the
    h::Threads.Atomic{Float32}   # windowEventWatcher callback, which can run in another thread!
end
Camera() = Camera(WorldPos(0,0),Threads.Atomic{Float32}(100),Threads.Atomic{Float32}(100))

# Note: These are all Atomics, since they can be modified by the
# windowEventWatcher callback, which can run in another thread!
const winWidth, winHeight = Threads.Atomic{Int32}(800), Threads.Atomic{Int32}(600)
const winWidth_highDPI, winHeight_highDPI = Threads.Atomic{Int32}(800), Threads.Atomic{Int32}(600)

screenCenter() = UIPixelPos(winWidth[]/2, winHeight[]/2)
screenCenterX() = winWidth[]/2
screenCenterY() = winHeight[]/2
screenOffsetFromCenter(x::Int,y::Int) = UIPixelPos(screenCenterX()+x,screenCenterY()+y)

dpiScale() = winWidth_highDPI[] / winWidth[];
worldScale(c::Camera) = dpiScale() * (winWidth[] / c.w[]);
function toScreenPos(p::WorldPos, c::Camera)
    scale = worldScale(c)
    ScreenPixelPos(
        round(winWidth_highDPI[]/2. + scale*(p.x-c.pos.x)), round(winHeight_highDPI[]/2. - scale*(p.y-c.pos.y)))
end
function toScreenPos(p::UIPixelPos, c::Camera)
    scale = dpiScale()
    ScreenPixelPos(round(scale*p.x), round(scale*p.y))
end
toScreenPos(p::ScreenPixelPos, c::Camera) = p
function toWorldPos(p::ScreenPixelPos, c::Camera)
    scale = worldScale(c)
    WorldPos((p.x - winWidth_highDPI[]/2.)/scale + c.pos.x, -(p.y - winHeight_highDPI[]/2.)/scale + c.pos.y)
end
function toWorldPos(p::UIPixelPos, c::Camera)
    toWorldPos(toScreenPos(p, c), c)
end
toWorldPos(p::WorldPos, c::Camera) = p
function toUIPixelPos(p::ScreenPixelPos, c::Camera)
    scale = dpiScale()
    UIPixelPos(round(p.x/scale), round(p.y/scale))
end
function toUIPixelPos(p::WorldPos, c::Camera)
    toUIPixelPos(toScreenPos(p, c), c)
end
toUIPixelPos(p::UIPixelPos, c::Camera) = p
function toScreenPixelDims(dims::UIPixelDims,c::Camera)
    scale = dpiScale()
    ScreenPixelDims(round(scale*dims.w), round(scale*dims.h))
end
function toScreenPixelDims(dims::WorldDims,c::Camera)
    scale = worldScale(c)
    ScreenPixelDims(round(scale*dims.w), round(scale*dims.h))
end
toScreenPixelDims(dims::ScreenPixelDims,c::Camera) = dims
function toUIPixelDims(dims::ScreenPixelDims,c::Camera)
    scale = dpiScale()
    UIPixelDims(round(dims.w/scale), round(dims.h/scale))
end
function toUIPixelDims(d::WorldDims, c::Camera)
    toUIPixelDims(toScreenPixelDims(d, c), c)
end
toUIPixelDims(dims::UIPixelDims,c::Camera) = dims
function toWorldDims(dims::ScreenPixelDims,c::Camera)
    scale = worldScale(c)
    WorldDims(dims.w/scale, dims.h/scale)
end
function toWorldDims(d::UIPixelDims, c::Camera)
    toWorldDims(toScreenPixelDims(d, c), c)
end
toWorldDims(dims::WorldDims,c::Camera) = dims

SetRenderDrawColor(renderer::Ptr{SDL2.Renderer}, c::SDL2.Color) = SDL2.SetRenderDrawColor(
    renderer, Int64(c.r), Int64(c.g), Int64(c.b), Int64(c.a))

# Convenience functions to allow `WorldPos(x,y)...` to become `x,y`
Base.iterate(p::AbstractPos, i=1) = if i==1 return (p.x,2) elseif i==2 return (p.y,3) else nothing end
Base.iterate(p::AbstractDims, i=1) = if i==1 return (p.w,2) elseif i==2 return (p.h,3) else nothing end
# To allow .+ (but doesn't work)
#Base.length(::Union{AbstractPos,AbstractDims}) = 2

topLeftPos(center::P, dims::D) where P<:AbstractPos{Coord} where D<:AbstractDims{Coord} where Coord<:ScreenCoords = P(center.x - dims.w/2., center.y - dims.h/2.)  # positive is down
topLeftPos(center::P, dims::D) where P<:AbstractPos{Coord} where D<:AbstractDims{Coord} where Coord = P(center.x - dims.w/2., center.y + dims.h/2.)  # positive is up
rectOrigin(center::P, dims::D) where P<:AbstractPos{C} where D<:AbstractDims{C} where {C} = P(center.x - dims.w/2., center.y - dims.h/2.)  # always minus

#makeSDL2RectFromOrigin(p::ScreenPixelPos, d::ScreenPixelDims) = SDL2.Rect(p..., d...)
#makeSDL2RectFromCenter(p::ScreenPixelPos, d::ScreenPixelDims) = SDL2.Rect(rectOrigin(p,d)..., d...)
#makeSDL2RectFromOrigin(p,d) = makeSDL2RectFromCenter(toScreenPos(p), toScreenPixelDims(d))
#makeSDL2RectFromCenter(p,d) =
function renderRectCentered(cam, renderer, center::AbstractPos{C}, dims::AbstractDims{C}, color; outlineColor=nothing) where C
    origin = topLeftPos(center, dims)
    renderRectFromOrigin(cam, renderer, origin, dims, color; outlineColor=outlineColor)
end
function renderRectFromOrigin(cam, renderer, origin::AbstractPos{C}, dims::AbstractDims{C}, color; outlineColor=nothing) where C
    screenPos = toScreenPos(origin, cam)
    rect = SDL2.Rect(screenPos.x, screenPos.y, toScreenPixelDims(dims, cam)...)
    if color != nothing
        SetRenderDrawColor(renderer, color)
        SDL2.RenderFillRect(renderer, Ref(rect) )
    end
    if outlineColor != nothing
        SetRenderDrawColor(renderer, outlineColor)
        SDL2.RenderDrawRect(renderer, Ref(rect) )
    end
end

function blendAlphaColors(x::SDL2.Color, y::SDL2.Color)
    xAlphaPercent = x.a / 255
    yAlphaPercent = y.a / 255
    r = round(x.r*xAlphaPercent + (y.r - Int32(x.r) * xAlphaPercent) * yAlphaPercent)
    g = round(x.g*xAlphaPercent + (y.g - Int32(x.g) * xAlphaPercent) * yAlphaPercent)
    b = round(x.b*xAlphaPercent + (y.b - Int32(x.b) * xAlphaPercent) * yAlphaPercent)
    a = round(x.a + (1-xAlphaPercent)*yAlphaPercent * 255)
    SDL2.Color(r,g,b,a)
end

# pointwise subtraction with bounds checking (floors to 0)
-(a::SDL2.Color, b::Int) = SDL2.Color(a.r-min(b,a.r), a.g-min(b,a.g), a.b-min(b,a.b), a.a-min(b,a.a))
SDL2.Color(1,5,1,1) - 2 == SDL2.Color(0,3,0,0)

# ---- Text Rendering ----
function set_default_font(fontname, fontsize)
    global defaultFontName = fontname
    global defaultFontSize = fontsize
end

fonts_cache = Dict()
txt_cache = Dict()
function sizeText(txt, fontName, fontSize)
    sizeText(txt, loadFont(dpiScale(), fontName, fontSize))
end
function sizeText(txt, font::Ptr{SDL2.TTF_Font})
   fw,fh = Cint[1], Cint[1]
   SDL2.TTF_SizeText(font, txt, pointer(fw), pointer(fh))
   return ScreenPixelDims(fw[1],fh[1])
end
function loadFont(scale, fontName, fontSize)
   fontSize = scale*fontSize
   fontKey = (fontName, Cint(round(fontSize)))
   if haskey(fonts_cache, fontKey)
       font = fonts_cache[fontKey]
   else
       font = SDL2.TTF_OpenFont(fontKey...)
       font == C_NULL && throw(ErrorException("Failed to load font '$fontKey'"))
       fonts_cache[fontKey] = font
   end
   return font
end
function createText(renderer, cam, txt, fontName, fontSize)
   font = loadFont(dpiScale(), fontName, fontSize)
   txtKey = (font, txt)
   if haskey(txt_cache, txtKey)
       tex = txt_cache[txtKey]
   else
       text = SDL2.TTF_RenderText_Blended(font, txt, SDL2.Color(20,20,20,255))
       tex = SDL2.CreateTextureFromSurface(renderer,text)
       #SDL2.FreeSurface(text)
       txt_cache[txtKey] = tex
   end

   fw,fh = Cint[1], Cint[1]
   SDL2.TTF_SizeText(font, txt, pointer(fw), pointer(fh))
   fw,fh = fw[1],fh[1]

   return tex, ScreenPixelDims(fw, fh)
end
@enum TextAlign centered leftJustified rightJustified
function renderText(renderer, cam::Camera, txt::String, pos::UIPixelPos
                    ; fontName = defaultFontName,
                     fontSize=defaultFontSize, align::TextAlign = centered)
   tex, fDims = createText(renderer, cam, txt, fontName, fontSize)
   renderTextSurface(renderer, cam, pos, tex, toUIPixelDims(fDims,cam), align)
end

function renderTextSurface(renderer, cam::Camera, pos::AbstractPos{C},
                           tex::Ptr{SDL2.Texture}, dims::AbstractDims{C}, align::TextAlign) where C
   screenPos = toScreenPos(pos, cam)
   screenDims = toScreenPixelDims(dims, cam)
   x,y, fw, fh = screenPos.x, screenPos.y, screenDims.w, screenDims.h
   renderPos = SDL2.Rect(0,0,0,0)
   if align == centered
       renderPos = SDL2.Rect(Int(floor(x-fw/2.)), Int(floor(y-fh/2.)), fw,fh)
   elseif align == leftJustified
       renderPos = SDL2.Rect(Int(floor(x)), Int(floor(y-fh/2.)), fw,fh)
   else # align == rightJustified
       renderPos = SDL2.Rect(Int(floor(x-fw)), Int(floor(y-fh/2.)), fw,fh)
   end
   SDL2.RenderCopy(renderer, tex, C_NULL, pointer_from_objref(renderPos))
   #SDL2.DestroyTexture(tex)
end

function hcat_render_text(lines, renderer, cam, gap, pos::UIPixelPos;
         fixedWidth=nothing, fontName=defaultFontName, fontSize=defaultFontSize)
    numLines = size(lines)[1]
    if fixedWidth != nothing
        widths = fill(fixedWidth,size(lines))
    else
        widths = [toUIPixelDims(sizeText(line, fontName, fontSize), cam).w for line in lines]
    end
    totalWidth = sum(widths) + gap*(numLines-1)
    runningWidth = 0
    leftMostPos = pos.x - totalWidth/2.0
    text_centers = []
    for i in 1:numLines
        linePos = leftMostPos + runningWidth
        leftPos = UIPixelPos(linePos, pos.y)
        renderText(renderer, cam, lines[i], leftPos; fontName=fontName, fontSize=fontSize, align=leftJustified)
        runningWidth += widths[i] + gap
        push!(text_centers, UIPixelPos(leftPos.x + widths[i]÷2, leftPos.y))
    end
    return text_centers
end

#  ------- Image rendering ---------

function render(t::Ptr{SDL2.Texture}, pos::AbstractPos{C}, cam::Camera, renderer;
                size::Union{Cvoid, AbstractDims{C}, Number} = nothing,
                src_rect::Union{Ptr{Cvoid}, SDL2.Rect} = C_NULL,
                ) where {C}
    if (t == C_NULL) return end
    pos = toScreenPos(pos, cam)
    if size isa AbstractDims
        size = toScreenPixelDims(size, cam)
        w = size.w
        h = size.h
    else
        w,h,access = Cint[1], Cint[1], Cint[1]
        format = Cuint[1]
        SDL2.QueryTexture( t, format, access, w, h );
        w,h = w[], h[]

        if size isa Number
            w,h = Cint(round(w*size)), Cint(round(h*size))
        end
    end
    rect = SDL2.Rect(pos.x - w÷2,pos.y - h÷2,w,h)
    SDL2.RenderCopy(renderer, t, src_rect, pointer_from_objref(rect))
end

# --------- Display navigation ----------

ext: jl


text:
using ADIOS2
using Base.Filesystem
using MPI
using Printf
using Test

################################################################################

# Initialize MPI
const mpi_initialized = MPI.Initialized()
if !mpi_initialized
    MPI.Init()
end
const comm = MPI.COMM_WORLD
const comm_rank = MPI.Comm_rank(comm)
const comm_size = MPI.Comm_size(comm)
const comm_root = 0
const use_mpi = comm_size > 1

################################################################################

"""
Convert an object to a string as the REPL would
"""
function showmime(obj)
    buf = IOBuffer()
    show(buf, MIME"text/plain"(), obj)
    return String(take!(buf))
end

################################################################################

include("internal.jl")
include("basic.jl")
include("highlevel.jl")

################################################################################

# Finalize MPI
const mpi_finalized = MPI.Finalized()
if mpi_initialized && !mpi_finalized
    MPI.Finalize()
end


text:
function get_frequency(obj::ObservationType, arg0::SatelliteSystem)
    return jcall(obj, "getFrequency", Frequency, (SatelliteSystem,), arg0)
end

function get_measurement_type(obj::ObservationType)
    return jcall(obj, "getMeasurementType", MeasurementType, ())
end

function get_signal_code(obj::ObservationType)
    return jcall(obj, "getSignalCode", SignalCode, ())
end

function value_of(::Type{ObservationType}, arg0::JString)
    return jcall(ObservationType, "valueOf", ObservationType, (JString,), arg0)
end

function values(::Type{ObservationType})
    return jcall(ObservationType, "values", Vector{ObservationType}, ())
end


text:
export PropensityGradient, TimeInvariantPropensityGradient, TimeVaryingPropensityGradient, statefactorgrad, tfactorgrad

"""
Base type for all propensity gradients.
"""
abstract type PropensityGradient end
abstract type TimeInvariantPropensityGradient <: PropensityGradient end
abstract type TimeVaryingPropensityGradient <: PropensityGradient end

"""
Store the partial derivatives of a `StandardTimeInvariantPropensity` instance.

# Fields
- `pardiffs`: Vector of callables, whose length is the number of model parameters. Each of these callables has syntax `pardiffs[i](x,p)::Real` where `x` is the state vector and `p` the parameter vector.
"""
struct StandardTimeInvariantPropensityGradient <: TimeInvariantPropensityGradient
    pardiffs::Vector{Any}
end

"""
Store the partial derivatives of a `SeparableTimeVaryingPropensity` instance.

# Fields
- `tfactor`: Callable in the form `tfactor(t,p)` where `t` is time, `p` is parameter vector.
- `statefactor`: Callable in the form `statefactor(x,p)` where `x` is CME state, `p` the parameter vector.
- `tfactor_pardiffs`: Vector of callables. Each of these callables has syntax `tfactor_pardiffs[i](t,p)::Real` where `t` is time variable and `p` the parameter vector.
- `statefactor_pardiffs`: Vector of callables, whose length is the number of model parameters. Each of these callables has syntax `statefactor_pardiffs[i](x,p)::Real` where `x` is the state vector and `p` the parameter vector.
"""
struct SeparableTimeVaryingPropensityGradient <: TimeVaryingPropensityGradient
    tfactor::Any
    statefactor::Any
    tfactor_pardiffs::Vector{Any}
    statefactor_pardiffs::Vector{Any}
end

"""
Store the partial derivatives of a `JointTimeVaryingPropensity` instance.

# Fields
- `pardiffs`: Vector of callables, whose length is the number of model parameters. Each of these callables has syntax `pardiffs[i](t,x,p)::Real` where `t` is time, `x` is the state vector and `p` the parameter vector.
"""
struct JointTimeVaryingPropensityGradient <: TimeVaryingPropensityGradient
    pardiffs::Vector{Any}
end

function propensitygrad(∇f::Vector{Any})
    return StandardTimeInvariantPropensityGradient(∇f)
end
function propensitygrad_timevarying(∇f::Vector{Any})
    return JointTimeVaryingPropensityGradient(∇f)
end
function propensitygrad_timevarying(tfactor::Any, statefactor::Any, ∇tfactor::Vector{Any}, ∇statefactor::Vector{Any})
    return SeparableTimeVaryingPropensityGradient(tfactor, statefactor, ∇tfactor, ∇statefactor)
end

function (g::StandardTimeInvariantPropensityGradient)(x::AbstractVector, p::AbstractVector)
    return [g.pardiffs[i](x, p) for i in 1:length(g.pardiffs)]
end
function (g::JointTimeVaryingPropensityGradient)(t::AbstractFloat, x::AbstractVector, p::AbstractVector)
    return [g.pardiffs[i](t, x, p) for i in 1:length(g.pardiffs)]
end
function (g::SeparableTimeVaryingPropensityGradient)(t::AbstractFloat, x::AbstractVector, p::AbstractVector)
    return [
        g.statefactor(x, p) * g.tfactor_pardiffs[i](t, p) + g.tfactor(t, p) * g.statefactor_pardiffs[i](x, p)
        for i in 1:length(g.tfactor_pardiffs)
    ]
end

get_pardiffs(g::Union{StandardTimeInvariantPropensityGradient,JointTimeVaryingPropensityGradient}) = g.pardiffs
get_single_pardiff(g::Union{StandardTimeInvariantPropensityGradient, JointTimeVaryingPropensityGradient}, ip) = g.pardiffs[ip]
get_tfactor_pardiffs(g::SeparableTimeVaryingPropensityGradient) = g.tfactor_pardiffs
get_statefactor_pardiffs(g::SeparableTimeVaryingPropensityGradient) = g.statefactor_pardiffs
get_single_tfactor_pardiff(g::SeparableTimeVaryingPropensityGradient, ip) = g.tfactor_pardiffs[ip]
get_single_statefactor_pardiff(g::SeparableTimeVaryingPropensityGradient, ip) = g.statefactor_pardiffs[ip]

include("senstools/sparsity_pattern.jl")
include("senstools/forwarddiffpropensity.jl")


text:
"""
    simpledown

Compute the mesh topology structures that correspond to just the downward
adjacency between the tetrahedra and the vertices.
Memory is used only to store the incidence relation and the locations of the
vertices.
"""
module simpledown
using StaticArrays
using MeshCore: nvertices, nshapes, manifdim, attribute
using MeshCore: ir_skeleton, ir_bbyfacets, ir_transpose
using MeshSteward: vtkwrite, T4block
using BenchmarkTools
using Test

include("usedbytes.jl")

function test()
    n = 3
    membytes = 0; summembytes = 0
    @info "Initial (3, 0)"
    @time connectivity = T4block(1.0, 2.0, 3.0, n*7, n*9, n*10, :a; intbytes = 4)
    ir30 = connectivity
    @show "($(manifdim(ir30.left)), $(manifdim(ir30.right)))"
    @show (nshapes(ir30.left), nshapes(ir30.right))
    @show membytes = usedbytes(ir30._v)
    summembytes += membytes
    geom = attribute(ir30.right, "geom")
    @show membytes = usedbytes(geom.v)
    summembytes += membytes

    # Print the total number of megabytes used to store the database
    @show summembytes/2^20

    # vtkwrite("speedtest1", connectivity)
    true
end
end
using .simpledown
# using BenchmarkTools
# @btime simpledown.test()
simpledown.test()


text:
## hclust.jl (c) 2014, 2017 David A. van Leeuwen
## Hierarchical clustering, similar to R's hclust()

## Algorithms are based upon C. F. Olson, Parallel Computing 21 (1995) 1313--1325.

"""
The output of [`hclust`](@ref), hierarchical clustering of data points.

Provides the bottom-up definition of the dendrogram as the sequence of
merges of the two lower subtrees into a higher level subtree.

This type mostly follows R's `hclust` class.

# Fields
 - `merges::Matrix{Int}`: ``N×2`` matrix encoding subtree merges:
   * each row specifies the left and right subtrees (referenced by their ``id``s)
     that are merged
   * negative subtree ``id`` denotes the leaf node and corresponds to the data
     point at position ``-id``
   * positive ``id`` denotes nontrivial subtree (the row `merges[id, :]`
     specifies its left and right subtrees)
 - `linkage::Symbol`: the name of *cluster linkage* function used to construct
   the hierarchy (see [`hclust`](@ref))
 - `heights::Vector{T}`: subtree heights, i.e. the distances between the left
    and right branches of each subtree calculated using the specified `linkage`
 - `order::Vector{Int}`: the data point indices ordered so that there are no
    intersecting branches on the *dendrogram* plot. This ordering also puts
    the points of the same cluster close together.

See also: [`hclust`](@ref).
"""
struct Hclust{T<:Real}
    merges::Matrix{Int} # the tree merge sequence. 1st column: left subtree, 2nd column: right subtree
    heights::Vector{T}  # subtrees heights (aggregated distance between its elements)
    order::Vector{Int}  # the order of datapoint (leaf node) indices in the final tree
    linkage::Symbol     # subtree distance type (cluster linkage)
end

nmerges(h::Hclust) = length(h.heights)  # number of tree merges
nnodes(h::Hclust) = length(h.order)     # number of datapoints (leaf nodes)
height(h::Hclust) = isempty(h.heights) ? typemin(eltype(h.heights)) : last(h.heights)

# FIXME remove after deprecation period for merge/height/method
Base.propertynames(hclu::Hclust, private::Bool = false) =
    (:merges, :heights, :order, :linkage,
     #= deprecated =# :height, :labels, :merge, :method)

# FIXME remove after deprecation period for height/labels/merge/method
function Base.getproperty(hclu::Hclust, prop::Symbol)
    if prop == :height
        Base.depwarn("Hclust::height is deprecated, use Hclust::heights", Symbol("Hclust::height"))
        return getfield(hclu, :heights)
    elseif prop == :labels
        Base.depwarn("Hclust::labels is deprecated and will be removed in future versions", Symbol("Hclust::labels"))
        return 1:nnodes(hclu)
    elseif prop == :merge
        Base.depwarn("Hclust::merge is deprecated, use Hclust::merges", Symbol("Hclust::merge"))
        return getfield(hclu, :merges)
    elseif prop == :method
        Base.depwarn("Hclust::method is deprecated, use Hclust::linkage", Symbol("Hclust::method"))
        return getfield(hclu, :linkage)
    else
        return getfield(hclu, prop)
    end
end

function assertdistancematrix(d::AbstractMatrix)
    nr, nc = size(d)
    nr == nc || throw(DimensionMismatch("Distance matrix should be square."))
    issymmetric(d) || throw(ArgumentError("Distance matrix should be symmetric."))
end

## R's order of trees
_isrordered(i::Integer, j::Integer) =
    i < 0 && j < 0 && i > j ||  # leaves (datapoints) are sorted in ascending order
    i > 0 && j > 0 && i < j ||  # if i-th tree was created before j-th one, it goes first
    i < 0 && j > 0              # leaves go before trees

# the sequence of tree merges
struct HclustMerges{T<:Real}
    nnodes::Int         # number of datapoints (leaf nodes)
    heights::Vector{T}  # tree height
    mleft::Vector{Int}  # ID of the left subtree merged
    mright::Vector{Int} # ID of the right subtree merged

    function HclustMerges{T}(nnodes::Integer) where {T<:Real}
        ntrees = max(nnodes-1, 0)
        new{T}(nnodes, sizehint!(T[], ntrees),
               sizehint!(Int[], ntrees), sizehint!(Int[], ntrees))
    end
end

nmerges(hmer::HclustMerges) = length(hmer.heights)
nnodes(hmer::HclustMerges) = hmer.nnodes

# merges i-th and j-th subtrees into a new tree of height h and returns its index
function push_merge!(hmer::HclustMerges{T}, i::Integer, j::Integer, h::T) where T<:Real
    if !_isrordered(i, j)
        i, j = j, i
    end
    push!(hmer.mleft, i)
    push!(hmer.mright, j)
    push!(hmer.heights, h)
    return nmerges(hmer)
end

#= utilities for working with the vector of clusters =#

cluster_size(cl::AbstractVector{Vector{Int}}, i::Integer) =
    return i > 0 ? length(cl[i]) : #= leaf node =# 1

# indices of nodes assigned to i-th cluster
# if i-th cluster is a leaf node (i < 0), return leafcluster setting its contents to [-i]
function cluster_elems(clusters::AbstractVector{Vector{Int}}, i::Integer,
                       leafcluster::AbstractVector{Int})
    if i > 0
        return clusters[i]
    else # i < 0 means it's a leaf node
        @assert length(leafcluster) == 1
        @inbounds leafcluster[1] = -i
        return leafcluster
    end
end

# merges i-th and j-th clusters and adds the result to the end of the `cl` list;
# i-th and j-th clusters are deactivated (emptied or replaced by `noels` vector)
# if either i or j are negative, the corresponding cluster is a leaf node (-i or -j, resp.)
function merge_clusters!(cl::AbstractVector{Vector{Int}},
                         i::Integer, j::Integer,
                         noels::Vector{Int} = Int[])
    if j < 0 # negative == cluster is a leaf node (-j)
        newclu = i < 0 ? [-i, -j] : push!(cl[i], -j)
    else
        clj = cl[j]
        if i < 0
            newclu = pushfirst!(clj, -i)
            cl[j] = noels
        else
            newclu = append!(cl[i], clj)
            empty!(clj) # not used anymore
        end
    end
    if i > 0
        cl[i] = noels # not used anymore
    end
    return push!(cl, newclu)
end

# compute resulting leaves (original datapoints) permutation
# given a sequence of tree nodes merges
function hclust_perm(hmer::HclustMerges)
    n = nmerges(hmer)
    perm = fill(1, nnodes(hmer)) # resulting permutation
    clusters = Vector{Int}[]     # clusters elements
    onel = [0]                   # placeholder for the elements of a leaf node
    noels = Int[]                # placeholder for empty decativated trees
    for i in 1:n
        ml = hmer.mleft[i]
        mr = hmer.mright[i]
        # elements in the right subtree are moved length(ml) positions from the start
        nl = cluster_size(clusters, ml)
        @inbounds for i in cluster_elems(clusters, mr, onel)
            perm[i] += nl
        end
        merge_clusters!(clusters, ml, mr, noels)
    end
    return perm
end

# convert HclustMerges to Hclust
function Hclust(hmer::HclustMerges, method::Symbol)
    Hclust(hcat(hmer.mleft, hmer.mright), hmer.heights,
           invperm(hclust_perm(hmer)), method)
end

# active trees of hclust algorithm
struct HclustTrees{T<:Real}
    merges::HclustMerges{T}     # history of tree merges
    id::Vector{Int}             # IDs of active trees
    cl::Vector{Vector{Int}}     # elements in the non-trivial trees
    noels::Vector{Int}          # empty placeholder for elements of deactivated trees

    HclustTrees{T}(n::Integer) where T<:Real =
        new{T}(HclustMerges{T}(n),
               collect(-(1:n)), # init with all leaves
               sizehint!(Vector{Int}[], n),
               Vector{Int}())
end

nmerges(htre::HclustTrees) = nmerges(htre.merges)
ntrees(htre::HclustTrees) = length(htre.id)
nnodes(htre::HclustTrees) = nnodes(htre.merges)

tree_size(htre::HclustTrees, i::Integer) = cluster_size(htre.cl, htre.id[i])

# ids of elements assigned to the tree with i-th index
# if i-th the is a leaf node, return leafcluster setting its contents to the id of that node
tree_elems(htre::HclustTrees, i::Integer,
           leafcluster::AbstractVector{Int}) =
    cluster_elems(htre.cl, htre.id[i], leafcluster)

# merges the trees referenced by indices i and j in htre.id into a new tree of height h;
# the i-th and j-th trees are deactived, their containers are emptied or replaced by `noels` placeholder
function merge_trees!(htre::HclustTrees, i::Integer, j::Integer, h::Real)
    # get tree ids
    ci = htre.id[i]
    cj = htre.id[j]
    cnew = push_merge!(htre.merges, ci, cj, h)
    # in the tree list, replace ci by cnew and cj by cnew the last tree
    htre.id[i] = cnew
    htre.id[j] = htre.id[end]
    pop!(htre.id)
    merge_clusters!(htre.cl, ci, cj, htre.noels)
    return htre
end

## This seems to work like R's implementation, but it is extremely inefficient
## This probably scales O(n^3) or worse. We can use it to check correctness
function hclust_n3(d::AbstractMatrix, linkage::Function)
    assertdistancematrix(d)
    T = eltype(linkage(d, 1:0, 1:0))
    htre = HclustTrees{T}(size(d, 1))
    onecol = [0]
    onerow = [0]
    while ntrees(htre) > 1
        # find the closest pair of trees mi and mj, mj < mi
        NNmindist = typemax(T)
        NNi = NNj = 0           # indices of nearest neighbors clusters
        for j in 1:ntrees(htre)
            cols = tree_elems(htre, j, onecol)
            for i in (j+1):ntrees(htre)
                rows = tree_elems(htre, i, onerow)
                dist = linkage(d, rows, cols) # very expensive
                if (NNi == 0) || (dist < NNmindist)
                    NNmindist = dist
                    NNi = i
                    NNj = j
                end
            end
        end
        merge_trees!(htre, NNj, NNi, NNmindist)
    end
    return htre.merges
end

"""
Base type for _reducible_ Lance–Williams cluster metrics.

The metric `d` is called _reducible_ if for any clusters `A`, `B` and `C` and
some `ρ > 0` s.t.
```
d(A, B) < ρ, d(A, C) > ρ, d(B, C) > ρ
```
it follows that
```
d(A∪B, C) > ρ
```

If the cluster metrics belongs to Lance-Williams family, there is an efficient
formula that defines `d(A∪B, C)` using `d(A, C)`, `d(B, C)` and `d(A, B)`.
"""
abstract type ReducibleMetric{T <: Real} end

"""
Distance between the clusters is the minimal distance between any pair of their
points.
"""
struct MinimalDistance{T} <: ReducibleMetric{T}
    MinimalDistance(d::AbstractMatrix{T}) where T<:Real = new{T}()
end

# update `metric` distance between `k`-th cluster and `i`-th cluster
# (`d[k, i]`, `k < i`) after `j`-th cluster was merged into `i`-th cluster
@inline update!(metric::MinimalDistance{T}, d::AbstractMatrix{T},
    k::Integer, i::Integer, d_ij::T, d_kj::T,
    ni::Integer, nj::Integer, nk::Integer
) where T =
    (d[k, i] > d_kj) && (d[k, i] = d_kj)

"""
Ward distance between the two clusters `A` and `B` is the amount by
which merging the two clusters into a single larger cluster `A∪B` would increase
the average squared distance of a point to its cluster centroid.
"""
struct WardDistance{T} <: ReducibleMetric{T}
    WardDistance(d::AbstractMatrix{T}) where T<:Real = new{typeof(one(T)/2)}()
end

# update `metric` distance between `k`-th cluster and `i`-th cluster
# (`d[k, i]`, `k < i`) after `j`-th cluster was merged into `i`-th cluster
@inline function update!(metric::WardDistance{T}, d::AbstractMatrix{T},
    k::Integer, i::Integer, d_ij::T, d_kj::T,
    ni::Integer, nj::Integer, nk::Integer
) where T
    nall = ni + nj + nk
    d[k, i] = ((ni + nk) * d[k, i] + (nj + nk) * d_kj - nk * d_ij) / nall
end

"""
Average distance between a pair of points from each clusters.
"""
struct AverageDistance{T} <: ReducibleMetric{T}
    AverageDistance(d::AbstractMatrix{T}) where T<:Real = new{typeof(one(T)/2)}()
end

# update `metric` distance between `k`-th cluster and `i`-th cluster
# (`d[k, i]`, `k < i`) after `j`-th cluster was merged into `i`-th cluster
@inline function update!(metric::AverageDistance{T}, d::AbstractMatrix{T},
    k::Integer, i::Integer, d_ij::T, d_kj::T,
    ni::Integer, nj::Integer, nk::Integer
) where T
    nij = ni + nj
    d[k, i] = (ni * d[k, i] + nj * d_kj) / nij
end

"""
Maximum distance between a pair of point from each clusters.
"""
struct MaximumDistance{T} <: ReducibleMetric{T}
    MaximumDistance(d::AbstractMatrix{T}) where T<:Real = new{T}()
end

# update `metric` distance between `k`-th cluster and `i`-th cluster
# (`d[k, i]`, `k < i`) after `j`-th cluster was merged into `i`-th cluster
@inline update!(metric::MaximumDistance{T}, d::AbstractMatrix{T},
    k::Integer, i::Integer, d_ij::T, d_kj::T,
    ni::Integer, nj::Integer, nk::Integer
) where T =
    (d[k, i] < d_kj) && (d[k, i] = d_kj)

# Update upper-triangular matrix `d` of cluster-cluster `metric`-based distances
# after merging cluster `j` into cluster `i` and
# moving the last cluster (`N`) into the `j`-th slot
function update_distance_after_merge!(
    d::AbstractMatrix{T},
    metric::ReducibleMetric{T},
    clu_size::Function,
    i::Integer, j::Integer, N::Integer
) where {T <: Real}
    @assert 1 <= i < j <= N <= size(d, 1) "1 ≤ i=$i < j=$j <= N=$N ≤ $(size(d, 1))"
    @inbounds d_ij = d[i, j]
    @inbounds ni = clu_size(i)
    @inbounds nj = clu_size(j)
    ## update d, split in ranges k<i, i<k<j, j<k≤newj
    for k in 1:i          # k ≤ i
        @inbounds update!(metric, d, k, i, d_ij, d[k,j], ni, nj, clu_size(k))
    end
    for k in (i+1):(j-1)  # i < k < j
        @inbounds update!(metric, d, i, k, d_ij, d[k,j], ni, nj, clu_size(k))
    end
    for k in (j+1):N      # j < k ≤ N
        @inbounds update!(metric, d, i, k, d_ij, d[j,k], ni, nj, clu_size(k))
    end
    ## move N-th row/col into j
    if j < N
        @inbounds d[j, j] = d[N, N]
        for k in 1:(j-1)         # k < j < N
            @inbounds d[k,j] = d[k,N]
        end
        for k in (j+1):(N-1)     # j < k < N
            @inbounds d[j,k] = d[k,N]
        end
    end
    return d
end

# nearest neighbor to i-th node given symmetric distance matrix d;
# returns 0 if no nearest neighbor (1×1 matrix)
function nearest_neighbor(d::AbstractMatrix, i::Integer, N::Integer=size(d, 1))
    (N <= 1) && return 0, NaN
    # initialize with the first non-i node
    @inbounds if i > 1
        NNi = 1
        NNdist = d[NNi, i]
    else
        NNi = 2
        NNdist = d[i, NNi]
    end
    @inbounds for j in (NNi+1):(i-1)
        if NNdist > d[j, i]
            NNi = j
            NNdist = d[j, i]
        end
    end
    @inbounds for j in (i+1):N
        if NNdist > d[i, j]
            NNi = j
            NNdist = d[i, j]
        end
    end
    return NNi, NNdist
end

## Efficient single link algorithm, according to Olson, O(n^2), fig 2.
## Verified against R's implementation, correct, and about 2.5 x faster
## For each i < j compute D(i,j) (this is already given)
## For each 0 < i ≤ n compute Nearest Neighbor NN(i)
## Repeat n-1 times
##   find i,j that minimize D(i,j)
##   merge clusters i and j
##   update D(i,j) and NN(i) accordingly
function hclust_minimum(ds::AbstractMatrix{T}) where T<:Real
    d = Matrix(ds)      # active trees distances, only upper (i < j) is used
    mindist = MinimalDistance(d)
    hmer = HclustMerges{T}(size(d, 1))
    n = nnodes(hmer)
    ## For each 0 < i ≤ n compute Nearest Neighbor NN[i]
    NN = [nearest_neighbor(d, i, n)[1] for i in 1:n]
    ## the main loop
    trees = collect(-(1:n))  # indices of active trees, initialized to all leaves
    while length(trees) > 1  # O(n)
        # find a pair of nearest trees, i and j
        i = 1
        NNmindist = i < NN[i] ? d[i, NN[i]] : d[NN[i], i]
        for k in 2:length(trees) # O(n)
            @inbounds dist = k < NN[k] ? d[k,NN[k]] : d[NN[k],k]
            if dist < NNmindist
                NNmindist = dist
                i = k
            end
        end
        j = NN[i]
        if i > j
            i, j = j, i     # make sure i < j
        end
        last_tree = length(trees)
        update_distance_after_merge!(d, mindist, i -> 0, i, j, last_tree)
        trees[i] = push_merge!(hmer, trees[i], trees[j], NNmindist)
        # reassign the last tree to position j
        trees[j] = trees[last_tree]
        NN[j] = NN[last_tree]
        pop!(NN)
        pop!(trees)
        ## update NN[k]
        for k in eachindex(NN)
            if NN[k] == j        # j is merged into i (only valid for the min!)
                NN[k] = i
            elseif NN[k] == last_tree # last_tree is moved into j
                NN[k] = j
            end
        end
        ## finally we need to update NN[i], because it was nearest to j
        NNmindist = typemax(T)
        NNi = 0
        for k in 1:(i-1)
            @inbounds if (NNi == 0) || (d[k,i] < NNmindist)
                NNmindist = d[k,i]
                NNi = k
            end
        end
        for k in (i+1):length(trees)
            @inbounds if (NNi == 0) || (d[i,k] < NNmindist)
                NNmindist = d[i,k]
                NNi = k
            end
        end
        NN[i] = NNi
#        for n in NN[1:nc] print(n, " ") end; println()
    end
    return hmer
end

## functions to compute maximum, minimum, mean for just a slice of an array
## FIXME: method(view(d, cl1, cl2)) would be much more generic, but it leads to extra allocations

function slicemaximum(d::AbstractMatrix, cl1::AbstractVector{Int}, cl2::AbstractVector{Int})
    maxdist = typemin(eltype(d))
    @inbounds for j in cl2, i in cl1
        if d[i,j] > maxdist
            maxdist = d[i,j]
        end
    end
    maxdist
end

function sliceminimum(d::AbstractMatrix, cl1::AbstractVector{Int}, cl2::AbstractVector{Int})
    mindist = typemax(eltype(d))
    @inbounds for j in cl2, i in cl1
        if d[i,j] < mindist
            mindist = d[i,j]
        end
    end
    mindist
end

function slicemean(d::AbstractMatrix, cl1::AbstractVector{Int}, cl2::AbstractVector{Int})
    s = zero(eltype(d))
    @inbounds for j in cl2
        sj = zero(eltype(d))
        for i in cl1
            sj += d[i,j]
        end
        s += sj
    end
    s / (length(cl1)*length(cl2))
end

## reorders the tree merges by the height of resulting trees
## (to be compatible with R's hclust())
function rorder!(hmer::HclustMerges)
    o = sortperm(hmer.heights)
    io = invperm(o)
    ml = hmer.mleft
    mr = hmer.mright
    for i in eachindex(ml)
        if ml[i] > 0
            ml[i] = io[ml[i]]
        end
        if mr[i] > 0
            mr[i] = io[mr[i]]
        end
        if !_isrordered(ml[i], mr[i])
            ml[i], mr[i] = mr[i], ml[i]
        end
    end
    permute!(ml, o)
    permute!(mr, o)
    Base.permute!!(hmer.heights, o)
    return hmer
end

## Another nearest neighbor algorithm, for reducible metrics
## From C. F. Olson, Parallel Computing 21 (1995) 1313--1325, fig 5
## Verfied against R implementation for mean and maximum, correct but ~ 5x slower
## Pick c1: 0 ≤ c1 ≤ n random
## i <- 1
## repeat n-1 times
##   repeat
##     i++
##     c[i] = nearest neighbor c[i-1]
##   until c[i] = c[i-2] ## nearest of nearest is cluster itself
##   merge c[i] and nearest neighbor c[i]
##   if i>3 i -= 3 else i <- 1
function hclust_nn(d::AbstractMatrix, linkage::Function)
    T = eltype(linkage(d, 1:0, 1:0))
    htre = HclustTrees{T}(size(d, 1))
    onerow = [0]  # placeholder for a leaf node of cl_i
    onecol = [0]  # placeholder for a leaf node of cl_j
    NN = [1]      # nearest neighbors chain of tree indices, init by random tree index
    while ntrees(htre) > 1
        # search for a pair of closest clusters,
        # they would be mutual nearest neighbors on top of the NN stack
        NNmindist = typemax(T)
        while true
            NNtop = NN[end]
            els_top = tree_elems(htre, NNtop, onecol)
            ## find NNnext: the nearest neighbor of NNtop and the next stack top
            NNnext = NNtop > 1 ? 1 : 2
            NNmindist = linkage(d, els_top, tree_elems(htre, NNnext, onerow))
            for k in (NNnext+1):ntrees(htre) if k != NNtop
                dist = linkage(d, tree_elems(htre, k, onerow), els_top)
                if dist < NNmindist
                    NNmindist = dist
                    NNnext = k
                end
            end end
            if length(NN) > 1 && NNnext == NN[end-1] # NNnext==NN[end-1] and NNtop=NN[end] are mutual n.neighbors
                break
            else
                push!(NN, NNnext)   # grow the chain
            end
        end
        ## merge NN[end] and its nearest neighbor, i.e., NN[end-1]
        NNlo = pop!(NN)
        NNhi = pop!(NN)
        if NNlo > NNhi
             NNlo, NNhi = NNhi, NNlo
        end
        last_tree = ntrees(htre)
        merge_trees!(htre, NNlo, NNhi, NNmindist)
        ## replace any nearest neighbor referring to the last_tree with NNhi
        if NNhi < last_tree
            for k in eachindex(NN)
                if NN[k] == last_tree
                    NN[k] = NNhi
                end
            end
        end
        isempty(NN) && push!(NN, 1) # restart NN chain
    end
    return rorder!(htre.merges)
end

## Nearest neighbor chain algorithm for reducible Lance-Williams metrics.
## In comparison to hclust_nn() maintains the upper-triangular matrix
## of cluster-cluster distances, so it requires O(N²) memory, but it's faster,
## because distance calculation is more efficient.
function hclust_nn_lw(d::AbstractMatrix, metric::ReducibleMetric{T}) where {T<:Real}
    dd = copyto!(Matrix{T}(undef, size(d)...), d)
    htre = HclustTrees{T}(size(d, 1))
    NN = [1]      # nearest neighbors chain of tree indices, init by random tree index
    while ntrees(htre) > 1
        # search for a pair of closest clusters,
        # they would be mutual nearest neighbors on top of the NN stack
        NNmindist = typemax(T)
        while true
            ## find NNnext: nearest neighbor of NN[end] (and the next stack top)
            NNnext, NNmindist = nearest_neighbor(dd, NN[end], ntrees(htre))
            @assert NNnext > 0
            if length(NN) > 1 && NNnext == NN[end-1] # NNnext==NN[end-1] and NN[end] are mutual n.neighbors
                break
            else
                push!(NN, NNnext)
            end
        end
        ## merge NN[end] and its nearest neighbor, i.e., NN[end-1]
        NNlo = pop!(NN)
        NNhi = pop!(NN)
        if NNlo > NNhi
            NNlo, NNhi = NNhi, NNlo
        end
        last_tree = ntrees(htre)
        ## update the distance matrix (while the trees are not merged yet)
        update_distance_after_merge!(dd, metric, i -> tree_size(htre, i), NNlo, NNhi, last_tree)
        merge_trees!(htre, NNlo, NNhi, NNmindist)
        ## replace any nearest neighbor referring to the last cluster with NNhi
        for k in eachindex(NN)
            if NN[k] == last_tree
                NN[k] = NNhi
            end
        end
        isempty(NN) && push!(NN, 1) # restart NN chain
    end
    return rorder!(htre.merges)
end

"""
    hclust(d::AbstractMatrix; [linkage], [uplo])

Perform hierarchical clustering using the distance matrix `d` and
the cluster `linkage` function.

Returns the dendrogram as an object of type [`Hclust`](@ref).

# Arguments
 - `d::AbstractMatrix`: the pairwise distance matrix. ``d_{ij}`` is the distance
    between ``i``-th and ``j``-th points.
 - `linkage::Symbol`: *cluster linkage* function to use. `linkage` defines how
   the distances between the data points are aggregated into the distances between
   the clusters. Naturally, it affects what clusters are merged on each
   iteration. The valid choices are:
   * `:single` (the default): use the minimum distance between any of the
     cluster members
   * `:average`: use the mean distance between any of the cluster members
   * `:complete`: use the maximum distance between any of the members
   * `:ward`: the distance is the increase of the average squared distance of
     a point to its cluster centroid after merging the two clusters
   * `:ward_presquared`: same as `:ward`, but assumes that the distances
     in `d` are already squared.
 - `uplo::Symbol` (optional): specifies whether the upper (`:U`) or the
   lower (`:L`) triangle of `d` should be used to get the distances.
   If not specified, the method expects `d` to be symmetric.
"""
function hclust(d::AbstractMatrix; linkage::Symbol = :single,
                uplo::Union{Symbol, Nothing} = nothing)
    if uplo !== nothing
        sd = Symmetric(d, uplo) # use upper/lower part of d
    else
        assertdistancematrix(d)
        sd = d
    end
    if linkage == :single
        hmer = hclust_minimum(sd)
    elseif linkage == :complete
        hmer = hclust_nn_lw(sd, MaximumDistance(sd))
    elseif linkage == :average
        hmer = hclust_nn_lw(sd, AverageDistance(sd))
    elseif linkage == :ward_presquared
        hmer = hclust_nn_lw(sd, WardDistance(sd))
    elseif linkage == :ward
        if sd === d
            sd = abs2.(sd)
        else
            sd .= abs2.(sd)
        end
        hmer = hclust_nn_lw(sd, WardDistance(sd))
        hmer.heights .= sqrt.(hmer.heights)
    else
        throw(ArgumentError("Unsupported cluster linkage $linkage"))
    end

    Hclust(hmer, linkage)
end

@deprecate hclust(d, method::Symbol, uplo::Union{Symbol, Nothing} = nothing) hclust(d, linkage=method, uplo=uplo)

"""
    cutree(hclu::Hclust; [k], [h])

Cuts the `hclu` dendrogram to produce clusters at the specified level of
granularity.

Returns the cluster assignments vector ``z`` (``z_i`` is the index of the
cluster for the ``i``-th data point).

# Arguments
 - `k::Integer` (optional) the number of desired clusters.
 - `h::Real` (optional) the height at which the tree is cut.

If both `k` and `h` are specified, it's guaranteed that the
number of clusters is not less than `k` and their height is not above `h`.

See also: [`hclust`](@ref)
"""
function cutree(hclu::Hclust;
                k::Union{Integer, Nothing} = nothing,
                h::Union{Real, Nothing} = nothing)
    # check k and h
    (k !== nothing || h !== nothing) ||
        throw(ArgumentError("Either `k` or `h` must be specified"))
    n = nnodes(hclu)
    m = nmerges(hclu)
    # use k and h to calculate how many merges to do before cutting
    if k !== nothing
        k >= min(1, n) || throw(ArgumentError("`k` should be greater or equal $(min(1,n))"))
        cutm = n - k
    else
        cutm = m
    end
    if h !== nothing
        # adjust cutm w.r.t h
        hix = findlast(hh -> hh ≤ h, hclu.heights)
        if hix !== nothing && hix < cutm
            cutm = hix
        elseif nmerges(hclu) >= 1 && first(hclu.heights) > h
            # corner case, the requested h smaller that the smallest nontrivial subtree
            cutm = 0
        end
    end
    clusters = Vector{Int}[]
    unmerged = fill(true, n) # if a node is not merged to a cluster
    noels = Int[]            # placeholder for empty deactivated trees
    i = 1
    while i ≤ cutm
        c1 = hclu.merges[i, 1]
        c2 = hclu.merges[i, 2]
        (c1 < 0) && (unmerged[-c1] = false)
        (c2 < 0) && (unmerged[-c2] = false)
        merge_clusters!(clusters, c1, c2, noels)
        i += 1
    end
    ## build an array of cluster indices (R's order)
    res = fill(0, n)
    # sort non-empty clusters by the minimal element index
    filter!(!isempty, clusters)
    permute!(clusters, sortperm(minimum.(clusters)))
    i = findfirst(unmerged)
    next = 1
    for clu in clusters
        cl1 = minimum(clu)
        while (i !== nothing) && (i < cl1)
            res[i] = next
            next += 1
            i = findnext(unmerged, i+1)
        end
        res[clu] .= next
        next += 1
    end
    while i !== nothing
        res[i] = next
        next += 1
        i = findnext(unmerged, i+1)
    end
    return res
end

## some diagnostic functions, not exported
function printupper(d::Matrix)
    n = size(d,1)
    for i in 1:(n-1)
        print(" " ^ ((i-1) * 6))
        for j in (i+1):n
            print(@sprintf("%5.2f ", d[i,j]))
        end
        println()
    end
end


text:
const DATA_DIR = joinpath(dirname(dirname(pathof(CDMrdata))),"data")

function load_data(dataset_name)
    rda = dataset_name * ".rda"
    try
        data = load(joinpath(DATA_DIR,rda))[ "data." * dataset_name]
        return data
    catch
        data = load(joinpath(DATA_DIR,rda))[dataset_name]
        return data
    end
end

list_datasets() =show(stdout, "text/plain", readdir(DATA_DIR))