text:
for rot in (:UnitQuaternion, :MRP, :RodriguesParam)
    @eval rotation_name(::Type{<:$rot}) = $rot
end

struct DiagonalLieCost{n,m,T,nV,nR,Rot} <: QuadraticCostFunction{n,m,T}
    Q::SVector{nV,T}
    R::Diagonal{T,SVector{m,T}}
    q::SVector{nV,T}
    r::SVector{m,T}
    c::T
    w::Vector{T}                     # weights on rotations (1 per rotation)
    vinds::SVector{nV,Int}           # inds of vector states
    qinds::Vector{SVector{nR,Int}}   # inds of rot states
    qrefs::Vector{UnitQuaternion{T}} # reference rotations
    function DiagonalLieCost(s::RD.LieState{Rot,P},
            Q::AbstractVector{<:Real}, R::AbstractVector{<:Real},
            q::AbstractVector{<:Real}, r::AbstractVector{<:Real},
            c::Real, w::AbstractVector,
            qrefs::Vector{<:Rotation}
        ) where {Rot,P}
        n = length(s)
        m = length(R)
        nV = sum(P)
        nR = Rotations.params(Rot)
        num_rots = length(P)-1
        @assert length(Q) == length(q) == nV
        @assert length(r) == m
        @assert length(qrefs) == length(w) == num_rots
        vinds = [RobotDynamics.vec_inds(Rot, P, i) for i = 1:num_rots+1]
        rinds = [RobotDynamics.rot_inds(Rot, P, i) for i = 1:num_rots]
        vinds = SVector{nV}(vcat(vinds...))
        rinds = SVector{nR}.(rinds)
        qrefs = UnitQuaternion.(qrefs)
        T = promote_type(eltype(Q), eltype(R), eltype(q), eltype(r), typeof(c), eltype(w))
        R = Diagonal(SVector{m}(R))
        new{n,m,T,nV,nR,rotation_name(Rot)}(Q, R, q, r, c, w, vinds, rinds, qrefs)
    end
    function DiagonalLieCost{Rot}(Q,R,q,r,c,w,vinds,qinds,qrefs) where Rot
        nV = length(Q)
        num_rots = length(qinds)
        nR = length(qinds[1])
        n = nV + num_rots*nR
        m = size(R,1)
        T = eltype(Q)
        new{n,m,T,nV,nR,rotation_name(Rot)}(Q, R, q, r, c, w, vinds, qinds, qrefs)
    end
end

function Base.copy(c::DiagonalLieCost)
    Rot = RD.rotation_type(c)
    DiagonalLieCost{Rot}(copy(c.Q), c.R, copy(c.q), copy(c.r), c.c, copy(c.w),
        copy(c.vinds), copy(c.qinds), deepcopy(c.qrefs)
    )
end

function DiagonalLieCost(s::RD.LieState{Rot,P},
        Q::Vector{<:AbstractVector{Tq}},
        R::Union{<:Diagonal{Tr},<:AbstractVector{Tr}};
        q = [zeros(Tq,p) for p in P],
        r=zeros(Tr,size(R,1)),
        c=0.0,
        w=ones(RobotDynamics.num_rotations(s)),
        qrefs=[one(UnitQuaternion) for i in 1:RobotDynamics.num_rotations(s)]
    ) where {Rot,P,Tq,Tr}
    @assert length.(Q) == collect(P) == length.(q)
    Q = vcat(Q...)
    q = vcat(q...)
    if R isa Diagonal
        R = R.diag
    end
    if w isa Real
        w = fill(w,length(qrefs))
    end
    DiagonalLieCost(s, Q, R, q, r, c, w, qrefs)
end

function DiagonalLieCost(s::RD.LieState{Rot,P}, Q::Diagonal, R::Diagonal;
        q::AbstractVector=zero(Q.diag), kwargs...
    ) where {Rot,P}
    if length(Q.diag) == sum(P)
        vinds = cumsum(insert(SVector(P), 1, 0))
        vinds = [vinds[i]+1:vinds[i+1] for i = 1:length(P)]
        Qs = [Q.diag[vind] for vind in vinds]
        qs = [q[vind] for vind in vinds]
        DiagonalLieCost(s, Qs, R.diag; q=qs, kwargs...)
    else
        num_rots = length(P) - 1
        vinds = [RobotDynamics.vec_inds(Rot, P, i) for i = 1:num_rots+1]
        rinds = [RobotDynamics.rot_inds(Rot, P, i) for i = 1:num_rots]

        # Split Q and q into vector and rotation components
        Qv = [Q[vind] for vind in vinds]
        Qr = [Q[rind] for rind in rinds]
        qv = [q[vind] for vind in vinds]

        # Use sum of diagonal Q as w
        w = sum.(Qr)
        @show vinds

        DiagonalLieCost(s, Qv, R.diag; q=qv, w=w, kwargs...)
    end
end

function LieLQRCost(s::RD.LieState{Rot,P},
        Q::Diagonal,
        R::Diagonal,
        xf, uf=zeros(size(R,1)); kwargs...) where {Rot,P}
    @assert length(Q.diag) == length(s)
    vinds = [RobotDynamics.vec_inds(Rot, P, i) for i = 1:length(P)]
    Qs = [Q.diag[vind] for vind in vinds]
    LieLQRCost(s, Qs, R, xf, uf; kwargs...)
end

function LieLQRCost(s::RD.LieState{Rot,P},
        Q::Vector{<:AbstractVector},
        R::Union{<:AbstractVector,<:Diagonal},
        xf,
        uf=zeros(size(R,1));
        w=ones(length(Q)-1),
    ) where {Rot,P}
    if R isa AbstractVector
        R = Diagonal(R)
    end
    num_rots = length(P) - 1
    vinds = [RobotDynamics.vec_inds(Rot, P, i) for i = 1:num_rots+1]
    qinds = [RobotDynamics.rot_inds(Rot, P, i) for i = 1:num_rots]
    qrefs = [Rot(xf[qinds[i]]) for i = 1:num_rots]
    q = [-Diagonal(Q[i])*xf[vinds[i]] for i = 1:length(P)]
    r = -R*uf
    c = sum(0.5*xf[vinds[i]]'Diagonal(Q[i])*xf[vinds[i]] for i = 1:length(P))
    c += 0.5*uf'R*uf
    return DiagonalLieCost(s, Q, R; q=q, r=r, c=c, w=w, qrefs=qrefs)
end

RobotDynamics.rotation_type(::DiagonalLieCost{<:Any,<:Any,<:Any, <:Any,<:Any, Rot}) where Rot = Rot
RobotDynamics.state_dim(::DiagonalLieCost{n}) where n = n
RobotDynamics.control_dim(::DiagonalLieCost{<:Any,m}) where m = m
is_blockdiag(::DiagonalLieCost) = true
is_diag(::DiagonalLieCost) = true

function stage_cost(cost::DiagonalLieCost, x::SVector, u::SVector)
    stage_cost(cost, x) + 0.5*u'cost.R*u + cost.r'u
end

function stage_cost(cost::DiagonalLieCost, x::AbstractVector)
    Rot = RobotDynamics.rotation_type(cost)
    Jv = veccost(cost.Q, cost.q, x, cost.vinds)
    Jr = quatcost(Rot, cost.w, x, cost.qinds, cost.qrefs)
    # @show Jv
    # @show Jr
    # @show cost.c
    return Jv + Jr + cost.c
end

function veccost(Q, q, x, vinds)
    xv = x[vinds]
    0.5*xv'Diagonal(Q)*xv + q'xv
end

function quatcost(::Type{Rot}, w, x, qinds, qref) where Rot<:Rotation
    J = zero(eltype(x))
    for i = 1:length(qinds)
        # q = Rotations.params(UnitQuaternion(Rot(x[qinds[i]])))
        q = toquat(Rot, x[qinds[i]])
        qd = Rotations.params(qref[i])
        err = q'qd
        # @show err
        J = w[i]*min(1-err,1+err)
        # @show J
    end
    return J
end

function gradient!(E::QuadraticCostFunction, cost::DiagonalLieCost, x::AbstractVector)
    # Vector states
    Rot = RobotDynamics.rotation_type(cost)
    xv = x[cost.vinds]
    E.q[cost.vinds] .= cost.Q .* xv + cost.q

    # Quaternion states
    for i = 1:length(cost.qinds)
        qind = cost.qinds[i]
        q = toquat(Rot, x[qind])
        qref = Rotations.params(cost.qrefs[i])
        dq = q'qref
        jac = Rotations.jacobian(UnitQuaternion, Rot(q))
        E.q[qind] .= -cost.w[i]*sign(dq)*jac'qref
    end
    return false
end

function hessian!(E::QuadraticCostFunction, cost::DiagonalLieCost{n}, x::AbstractVector) where n
    for (i,j) in enumerate(cost.vinds)
        E.Q[j,j] = cost.Q[i]
    end

    for qind in cost.qinds
        for i in qind
            E.Q[i,i] = 0
        end
    end
    return true
end

toquat(::Type{<:UnitQuaternion}, q::AbstractVector) = q
toquat(::Type{Rot}, q::AbstractVector) where Rot <: Rotation =
    Rotations.params(UnitQuaternion(Rot(q)))

# Jacobians missing from Rotatations
function Rotations.jacobian(::Type{UnitQuaternion}, p::RodriguesParam)
    p = Rotations.params(p)
    np = p'p
    d = 1/sqrt(1 + p'p)
    d3 = d*d*d
    ds = -d3 * p
    pp = -d3*(p*p')
    SA[
        ds[1] ds[2] ds[3];
        pp[1] + d pp[4] pp[7];
        pp[2] pp[5] + d pp[8];
        pp[3] pp[6] pp[9] + d;
    ]
end

function Rotations.jacobian(::Type{RodriguesParam}, q::UnitQuaternion)
    s = 1/q.w
    SA[
        -s*s*q.x s 0 0;
        -s*s*q.y 0 s 0;
        -s*s*q.z 0 0 s
    ]
end

Rotations.jacobian(::Type{R}, q::R) where R = I

############################################################################################
#                        QUADRATIC QUATERNION COST FUNCTION
############################################################################################
"""
    DiagonalQuatCost

Quadratic cost function for states that includes a 3D rotation, that penalizes deviations
    from a provided 3D rotation, represented as a Unit Quaternion.

The cost function penalizes geodesic distance between unit quaternions:

``\\frac{1}{2} \\big( x^T Q x + u^T R u \\big) + q^T x + r^T u + c + w \\min 1 \\pm p_f^T p``

where ``p`` is the quaternion extracted from ``x`` (i.e. `p = x[q_ind]`), and ``p_f``
is the reference quaternion. ``Q`` and ``R`` are assumed to be diagonal.

We've found this perform better than penalizing a quadratic on the quaternion error
state ([`ErrorQuadratic`](@ref)). This cost should still be considered experimental.

# Constructors
* `DiagonalQuatCost(Q::Diagonal, R::Diagonal, q, r, c, w, q_ref, q_ind; terminal)`
* `DiagonalQuatCost(Q::Diagonal, R::Diagonal; q, r, c, w, q_ref, q_ind, terminal)`

where `q_ref` is the reference quaternion (provided as a `SVector{4}`), and
    `q_ind::SVector{4,Int}` provides the indices of the quaternion in the state vector
    (default = `SA[4,5,6,7]`). Note that `Q` and `q` are the size of the full state,
    so `Q.diag[q_ind]` and `q[qind]` should typically be zero.
"""
struct DiagonalQuatCost{N,M,T,N4} <: QuadraticCostFunction{N,M,T}
    Q::Diagonal{T,SVector{N,T}}
    R::Diagonal{T,SVector{M,T}}
    q::SVector{N,T}
    r::SVector{M,T}
    c::T
    w::T
    q_ref::SVector{4,T}
    q_ind::SVector{4,Int}
    Iq::SMatrix{N,4,T,N4}
    terminal::Bool
    function DiagonalQuatCost(Q::Diagonal{T,SVector{N,T}}, R::Diagonal{T,SVector{M,T}},
            q::SVector{N,T}, r::SVector{M,T}, c::T, w::T,
            q_ref::SVector{4,T}, q_ind::SVector{4,Int}; terminal::Bool=false) where {T,N,M}
        Iq = @MMatrix zeros(N,4)
        for i = 1:4
            Iq[q_ind[i],i] = 1
        end
        Iq = SMatrix{N,4}(Iq)
        return new{N,M,T,N*4}(Q, R, q, r, c, w, q_ref, q_ind, Iq, terminal)
    end
end

state_dim(::DiagonalQuatCost{N,M,T}) where {T,N,M} = N
control_dim(::DiagonalQuatCost{N,M,T}) where {T,N,M} = M
is_blockdiag(::DiagonalQuatCost) = true
is_diag(::DiagonalQuatCost) = true

function DiagonalQuatCost(Q::Diagonal{T,SVector{N,T}}, R::Diagonal{T,SVector{M,T}};
        q=(@SVector zeros(N)), r=(@SVector zeros(M)), c=zero(T), w=one(T),
        q_ref=(@SVector [1.0,0,0,0]), q_ind=(@SVector [4,5,6,7])) where {T,N,M}
    DiagonalQuatCost(Q, R, q, r, c, q_ref, q_ind)
end

function stage_cost(cost::DiagonalQuatCost, x::SVector, u::SVector)
    stage_cost(cost, x) + 0.5*u'cost.R*u + cost.r'u
end

function stage_cost(cost::DiagonalQuatCost, x::SVector)
    J = 0.5*x'cost.Q*x + cost.q'x + cost.c
    q = x[cost.q_ind]
    dq = cost.q_ref'q
    J += cost.w*min(1+dq, 1-dq)
end

function gradient!(E::QuadraticCostFunction, cost::DiagonalQuatCost{T,N,M},
        x::AbstractVector) where {T,N,M}
    Qx = cost.Q*x + cost.q
    q = x[cost.q_ind]
    dq = cost.q_ref'q
    if dq < 0
        Qx += cost.w*cost.Iq*cost.q_ref
    else
        Qx -= cost.w*cost.Iq*cost.q_ref
    end
    E.q .= Qx
    return false
end

"""
    QuatLQRCost(Q, R, xf, [uf; w, quat_ind])

Defines a cost function of the form:

``\\frac{1}{2} \\big( (x - x_f)^T Q (x - x_f) + (u - u_f)^T R (u - u_f) \\big) + w \\min 1 \\pm q_f^T q``

where ``Q`` and ``R`` are diagonal, ``x_f`` is the goal state,
``u_f`` is the reference control, and ``q_f``, ``q`` are
the quaternions, extracted from ``x`` using `quat_ind`, i.e. `q = x[quat_ind]`.

The last term is the geodesic distance between quaternions. It's typically recommended that
`Q.diag[quad_ind] == zeros(4)`.

This is just a convenience constructor for [`DiagonalQuatCost`](@ref).

# Example
For a standard rigid body state vector `x = [p; q; v; ω]`, where `q` is a unit quaternion,
we could define a cost function that penalizes the distance to the goal state `xf`.
We can create this cost function as follows:
```julia
Q = Diagonal(SVector(RBState(fill(0.1,3), zeros(4), fill(0.1,3), fill(0.1,3))))
R = Diagonal(@SVector fill(0.01, 6))
xf = RBState([1,2,3], rand(UnitQuaternion), zeros(3), zeros(3))
QuatLQRCost(Q,R,xf)
```
We can add a reference control and change the weight on the rotation error with the optional
arguments:
```julia
QuatLQRCost(Q,R,xf,uf, w=10.0)
```
which is equivalent to
```julia
QuatLQRCost(Q,R,xf,uf, w=10.0, quat_inds=4:7)
```
"""
function QuatLQRCost(Q::Diagonal{T,SVector{N,T}}, R::Diagonal{T,SVector{M,T}}, xf,
        uf=(@SVector zeros(M)); w=one(T), quat_ind=(@SVector [4,5,6,7])) where {T,N,M}
    @assert length(quat_ind) == 4 "quat_ind argument must be of length 4"
    quat_ind = SVector{4}(quat_ind)
    r = -R*uf
    q = -Q*xf
    c = 0.5*xf'Q*xf + 0.5*uf'R*uf
    q_ref = xf[quat_ind]
    return DiagonalQuatCost(Q, R, q, r, c, T(w), q_ref, quat_ind)
end

function change_dimension(cost::DiagonalQuatCost, n, m, ix, iu)
    Qd = zeros(n)
    Rd = zeros(m)
    q = zeros(n)
    r = zeros(m)
    Qd[ix] = diag(cost.Q)
    Rd[iu] = diag(cost.R)
    q[ix] = cost.q
    r[iu] = cost.r
    qind = (1:n)[ix[cost.q_ind]]
    DiagonalQuatCost(Diagonal(SVector{n}(Qd)), Diagonal(SVector{m}(Rd)),
        SVector{n}(q), SVector{m}(r), cost.c, cost.w, cost.q_ref, qind)
end

function (+)(cost1::DiagonalQuatCost, cost2::QuadraticCostFunction)
    @assert state_dim(cost1) == state_dim(cost2)
    @assert control_dim(cost1) == control_dim(cost2)
    is_diag(cost2) || @assert norm(cost2.H) ≈ 0
    DiagonalQuatCost(cost1.Q + cost2.Q, cost1.R + cost2.R,
        cost1.q + cost2.q, cost1.r + cost2.r, cost1.c + cost2.c,
        cost1.w, cost1.q_ref, cost1.q_ind)
end

(+)(cost1::QuadraticCostFunction, cost2::DiagonalQuatCost) = cost2 + cost1

function Base.copy(c::DiagonalQuatCost)
    DiagonalQuatCost(c.Q, c.R, c.q, c.r, c.c, c.w, c.q_ref, c.q_ind)
end

############################################################################################
#                             Error Quadratic
############################################################################################

"""
    ErrorQuadratic{Rot,N,M}

Cost function of the form:

``\\frac{1}{2} (x_k \\ominus x_d)^T Q_k (x_k \\ominus x_d)``

where ``x_k \\ominus x_d`` is the error state, computed using
`RobotDynamics.state_diff`.
This cost function isn't recommended: we've found that `DiagonalQuatCost` usually
    peforms better and is much more computationally efficient.
"""
struct ErrorQuadratic{Rot,N,M} <: CostFunction
    model::RD.RigidBody{Rot}
    Q::Diagonal{Float64,SVector{12,Float64}}
    R::Diagonal{Float64,SVector{M,Float64}}
    r::SVector{M,Float64}
    c::Float64
    x_ref::SVector{N,Float64}
    q_ind::SVector{4,Int}
end
function Base.copy(c::ErrorQuadratic)
    ErrorQuadratic(c.model, c.Q, c.R, c.r, c.c, c.x_ref, c.q_ind)
end

state_dim(::ErrorQuadratic{Rot,N,M}) where {Rot,N,M} = N
control_dim(::ErrorQuadratic{Rot,N,M}) where {Rot,N,M} = M

function ErrorQuadratic(model::RD.RigidBody{Rot}, Q::Diagonal{T,<:SVector{N0}},
        R::Diagonal{T,<:SVector{M}},
        x_ref::SVector{N},
        u_ref=(@SVector zeros(T,M));
        r=(@SVector zeros(T,M)),
        c=zero(T),
        q_ind=(@SVector [4,5,6,7])
    ) where {T,N,N0,M,Rot}
    if Rot <: UnitQuaternion && N0 == N
        Qd = deleteat(Q.diag, 4)
        Q = Diagonal(Qd)
    end
    r += -R*u_ref
    c += 0.5*u_ref'R*u_ref
    return ErrorQuadratic{Rot,N,M}(model, Q, R, r, c, x_ref, q_ind)
end

function stage_cost(cost::ErrorQuadratic, x::SVector)
    dx = RD.state_diff(cost.model, x, cost.x_ref, Rotations.ExponentialMap())
    return 0.5*dx'cost.Q*dx + cost.c
end

function stage_cost(cost::ErrorQuadratic, x::SVector, u::SVector)
    stage_cost(cost, x) + 0.5*u'cost.R*u + cost.r'u
end

function gradient!(E::QuadraticCostFunction, cost::ErrorQuadratic, x)
    f(x) = stage_cost(cost, x)
    ForwardDiff.gradient!(E.q, f, x)
    return false

    model = cost.model
    Q = cost.Q
    q = RD.orientation(model, x)
    q_ref = RD.orientation(model, cost.x_ref)
    dq = Rotations.params(q_ref \ q)
    err = RD.state_diff(model, x, cost.x_ref)
    dx = @SVector [err[1],  err[2],  err[3],
                    dq[1],   dq[2],   dq[3],   dq[4],
                   err[7],  err[8],  err[9],
                   err[10], err[11], err[12]]
    # G = state_diff_jacobian(model, dx) # n × dn

    # Gradient
    dmap = inverse_map_jacobian(model, dx) # dn × n
    # Qx = G'dmap'Q*err
    Qx = dmap'Q*err
    E.q = Qx
    return false
end
function gradient!(E::QuadraticCostFunction, cost::ErrorQuadratic, x, u)
    gradient!(E, cost, x)
    Qu = cost.R*u
    E.r .= Qu
    return false
end

function hessian!(E::QuadraticCostFunction, cost::ErrorQuadratic, x)
    f(x) = stage_cost(cost, x)
    ForwardDiff.hessian!(E.Q, f, x)
    return false

    model = cost.model
    Q = cost.Q
    q = RD.orientation(model, x)
    q_ref = RD.orientation(model, cost.x_ref)
    dq = Rotations.params(q_ref\q)
    err = RD.state_diff(model, x, cost.x_ref)
    dx = @SVector [err[1],  err[2],  err[3],
                    dq[1],   dq[2],   dq[3],   dq[4],
                   err[7],  err[8],  err[9],
                   err[10], err[11], err[12]]
    # G = state_diff_jacobian(model, dx) # n × dn

    # Gradient
    dmap = inverse_map_jacobian(model, dx) # dn × n

    # Hessian
    ∇jac = inverse_map_∇jacobian(model, dx, Q*err)
    # Qxx = G'dmap'Q*dmap*G + G'∇jac*G + ∇²differential(model, x, dmap'Q*err)
    Qxx = dmap'Q*dmap + ∇jac #+ ∇²differential(model, x, dmap'Q*err)
    E.Q = Qxx
    E.H .*= 0
    return false
end

function hessian!(E::QuadraticCostFunction, cost::ErrorQuadratic, x, u)
    hessian!(E, cost, x)
    E.R .= cost.R
    return false
end

function change_dimension(cost::ErrorQuadratic, n, m)
    n0,m0 = state_dim(cost), control_dim(cost)
    Q_diag = diag(cost.Q)
    R_diag = diag(cost.R)
    r = cost.r
    if n0 != n
        dn = n - n0  # assumes n > n0
        pad = @SVector zeros(dn) # assume the new states don't have quaternions
        Q_diag = [Q_diag; pad]
    end
    if m0 != m
        dm = m - m0  # assumes m > m0
        pad = @SVector zeros(dm)
        R_diag = [R_diag; pad]
        r = [r; pad]
    end
    ErrorQuadratic(cost.model, Diagonal(Q_diag), Diagonal(R_diag), r, cost.c,
        cost.x_ref, cost.q_ind)
end

function (+)(cost1::ErrorQuadratic, cost2::QuadraticCost)
    @assert control_dim(cost1) == control_dim(cost2)
    @assert norm(cost2.H) ≈ 0
    @assert norm(cost2.q) ≈ 0
    if state_dim(cost2) == 13
        rm_quat = @SVector [1,2,3,4,5,6,8,9,10,11,12,13]
        Q2 = Diagonal(diag(cost2.Q)[rm_quat])
    else
        Q2 = cost2.Q
    end
    ErrorQuadratic(cost1.model, cost1.Q + Q2, cost1.R + cost2.R,
        cost1.r + cost2.r, cost1.c + cost2.c,
        cost1.x_ref, cost1.q_ind)
end

(+)(cost1::QuadraticCost, cost2::ErrorQuadratic) = cost2 + cost1

@generated function RD.state_diff_jacobian(model::RD.RigidBody{<:UnitQuaternion},
    x0::SVector{N,T}, errmap::D=Rotations.CayleyMap()) where {N,T,D}
    if D <: IdentityMap
        :(I)
    else
        quote
            q0 = RD.orientation(model, x0)
            # G = TrajectoryOptimization.∇differential(q0)
            G = Rotations.∇differential(q0)
            I1 = @SMatrix [1 0 0 0 0 0 0 0 0 0 0 0;
                        0 1 0 0 0 0 0 0 0 0 0 0;
                        0 0 1 0 0 0 0 0 0 0 0 0;
                        0 0 0 G[1] G[5] G[ 9] 0 0 0 0 0 0;
                        0 0 0 G[2] G[6] G[10] 0 0 0 0 0 0;
                        0 0 0 G[3] G[7] G[11] 0 0 0 0 0 0;
                        0 0 0 G[4] G[8] G[12] 0 0 0 0 0 0;
                        0 0 0 0 0 0 1 0 0 0 0 0;
                        0 0 0 0 0 0 0 1 0 0 0 0;
                        0 0 0 0 0 0 0 0 1 0 0 0;
                        0 0 0 0 0 0 0 0 0 1 0 0;
                        0 0 0 0 0 0 0 0 0 0 1 0;
                        0 0 0 0 0 0 0 0 0 0 0 1.]
        end
    end
end
function inverse_map_jacobian(model::RD.RigidBody{<:UnitQuaternion},
    x::SVector, errmap=Rotations.CayleyMap())
    q = RD.orientation(model, x)
    # G = TrajectoryOptimization.inverse_map_jacobian(q)
    G = Rotations.jacobian(inv(errmap), q)
    return @SMatrix [
            1 0 0 0 0 0 0 0 0 0 0 0 0;
            0 1 0 0 0 0 0 0 0 0 0 0 0;
            0 0 1 0 0 0 0 0 0 0 0 0 0;
            0 0 0 G[1] G[4] G[7] G[10] 0 0 0 0 0 0;
            0 0 0 G[2] G[5] G[8] G[11] 0 0 0 0 0 0;
            0 0 0 G[3] G[6] G[9] G[12] 0 0 0 0 0 0;
            0 0 0 0 0 0 0 1 0 0 0 0 0;
            0 0 0 0 0 0 0 0 1 0 0 0 0;
            0 0 0 0 0 0 0 0 0 1 0 0 0;
            0 0 0 0 0 0 0 0 0 0 1 0 0;
            0 0 0 0 0 0 0 0 0 0 0 1 0;
            0 0 0 0 0 0 0 0 0 0 0 0 1;
    ]
end

function inverse_map_∇jacobian(model::RD.RigidBody{<:UnitQuaternion},
    x::SVector, b::SVector, errmap=Rotations.CayleyMap())
    q = RD.orientation(model, x)
    bq = @SVector [b[4], b[5], b[6]]
    # ∇G = TrajectoryOptimization.inverse_map_∇jacobian(q, bq)
    ∇G = Rotations.∇jacobian(inv(errmap), q, bq)
    return @SMatrix [
        0 0 0 0 0 0 0 0 0 0 0 0 0;
        0 0 0 0 0 0 0 0 0 0 0 0 0;
        0 0 0 0 0 0 0 0 0 0 0 0 0;
        0 0 0 ∇G[1] ∇G[5] ∇G[ 9] ∇G[13] 0 0 0 0 0 0;
        0 0 0 ∇G[2] ∇G[6] ∇G[10] ∇G[14] 0 0 0 0 0 0;
        0 0 0 ∇G[3] ∇G[7] ∇G[11] ∇G[15] 0 0 0 0 0 0;
        0 0 0 ∇G[4] ∇G[8] ∇G[12] ∇G[16] 0 0 0 0 0 0;
        0 0 0 0 0 0 0 0 0 0 0 0 0;
        0 0 0 0 0 0 0 0 0 0 0 0 0;
        0 0 0 0 0 0 0 0 0 0 0 0 0;
        0 0 0 0 0 0 0 0 0 0 0 0 0;
        0 0 0 0 0 0 0 0 0 0 0 0 0;
        0 0 0 0 0 0 0 0 0 0 0 0 0;
    ]
end


text:
using Random
using DPMMSubClusters
include("run_experiments.jl")

dir = "/home/vlad/projects/DPMMSubClusters.jl/test_dir/"
SEED = 1234

Random.seed!(SEED)

# model = "pointNet"
model = "setTransformer"

REPEATS=3
ITERS=100
BURNOUT=10
# REPEATS=1
# ITERS=10
# BURNOUT=2
ALPHA=100

imagenet_random = repeat_experiment(:gauss2d,REPEATS, "none", ALPHA, BURNOUT, ITERS)
imagenet_kmeans = repeat_experiment(:gauss2d,REPEATS, "kmeans", ALPHA, BURNOUT, ITERS)
imagenet_splitnet = repeat_experiment(:gauss2d,REPEATS, "splitnet_2d", ALPHA,BURNOUT, ITERS)

# imagenet_random = repeat_experiment(:imagenet_swav_d128,REPEATS, "none", ALPHA, BURNOUT, ITERS)
# imagenet_kmeans = repeat_experiment(:imagenet_swav_d128,REPEATS, "kmeans", ALPHA, BURNOUT, ITERS)
# imagenet_splitnet = repeat_experiment(:imagenet_swav_d128,REPEATS, "splitnet-gpu", ALPHA,BURNOUT, ITERS)

all_results = Dict(
    "vanilla" => imagenet_random,
    "kmeans" => imagenet_kmeans,
    "splitnet" => imagenet_splitnet,
)

title = "ImageNet-Swav-D=128--Burnout-" * string(BURNOUT) * "_iters_"* string(ITERS)* "_alpha_"* string(ALPHA)

f1= plot_nmi_for_iters(all_results, title)
# savefig("./compare_methods_gifs/result_nmi_"* string(title) *".png")

f2 = plot_ll_for_iters(all_results, title)

f3 = plot_K_for_iters(all_results, title)

# f3 = plot(f1, f2, layout = (2, 1), legend = true, size=(600, 800))
f4 = plot(f1, f2, f3, layout = (3, 1), size=(600, 1200))

savefig("../compare_methods_gifs/result_"* string(title) *".png")

ext: jl


text:
W,(WW,EW,FW) = Lar.cuboid([1,1], true)
square = (W,EW)

square = Lar.Struct([  Lar.r(pi/4),  square  ])
squares = Lar.Struct([  square,  Lar.t(√2/2, 0),  square  ])

V,EV = Lar.struct2lar(squares)
VV = [[k] for k=1:size(V,2)]
GL.VIEW(GL.numbering(.4)((V,[VV, EV])));

model = CAGD.Model(V);
CAGD.addModelCells!(model, 1, EV);

# Planar Arrangment

A_model = CAGD.planar_arrangement(model);

A_model.G
Matrix(A_model.T[1])
Matrix(A_model.T[2])

# Display Models

if todisplay
    VV = [[k] for k=1:size(A_model.G, 2)];
    triangulated_faces = CAGD.triangulate2D(A_model);
    FVs = convert(Array{Lar.Cells}, triangulated_faces);

    GL.VIEW(GL.GLExplode(A_model.G, FVs, 1.0, 1.0, 1.0, 99, 1));
    GL.VIEW(GL.GLExplode(A_model.G, FVs, 1.2, 1.2, 1.2, 99, 1));
end

innerpoints = internalpoints2d(A_model.G', A_model.T[1], A_model.T[2][1:end, :])
# associate internal points to 2-cells
listOfModels = Lar.evalStruct(squares)
inputfacenumbers = [length(listOfModels[k][2]) for k=1:length(listOfModels)]
# test input data for containment of reference points
boolmatrix = BitArray(undef, length(innerpoints), length(listOfModels))
containmenttest = Lar.testinternalpoint2d(listOfModels)
for (k,point) in enumerate(innerpoints) # k runs on columns
    cells = containmenttest(point) # contents of columns
    #println(k," ",faces)
    for l in cells
        boolmatrix[k,l] = 1
    end
end

sq1 = boolmatrix[:, 1]
sq2 = boolmatrix[:, 2]

sqXor = sq1 .⊻ sq2

if todisplay
    GL.VIEW(GL.GLExplode(A_model.G, FVs[sparse(sqXor).nzind], 1.0, 1.0, 1.0, 99, 1));
end


text:
function CJS_empirical(F::Matrix,ndim::Int64,nsamples::Int64,nsteps::Int64)

    results = zeros(ndim,ndim)

    for k in 1:ndim
        for l in k:ndim
            lbound = minimum(F[:,k])
            rbound = maximum(F[:,k])

            #print(edf(F[:,k],rbound))

            int=0.0

            for u in 1:nsteps

                x = lbound+(u-1)*(rbound-lbound)/nsteps
                p=edf(F[:,k],x)
                q=edf(F[:,l],x)

                int+=(p*log2(2*p/(p+q))+1/2 * 1/log(2) * (q-p))/nsteps

            end

            results[k,l] = int

        end
    end

    return(results+results')
end

function edf(data::Array{Float64},x::Float64)
    return(length(data[data .<= x])/length(data))
end


text:
"""

In this example, the distance between the C-terminal carbon
and a sodium atom is computed and ploted as a function of time

"""

using M3GTools

println(" Loading Plots... ")
using Plots
ENV["GKSwstype"]="nul" # This supresses the need of a display while plotting

# Initialize simulation data and path to VMD executable
println(" Reading simulation data ... ")
mysim = M3GTools.init(psf="../simulation_files/structure.psf",
                      dcd="../simulation_files/structure.dcd",
                      vmd="vmd")

println(" Defining selections... ")
cterm = M3GTools.select(mysim,"protein and resid 253 and name C")
sod = M3GTools.select(mysim,"index 20455")

# Vector that will contain the distances for each frame
d = Vector{Float64}(undef,mysim.nframes)

println(" Computing distances ... ")
for i in 1:mysim.nframes
  sides, x, y, z = M3GTools.nextframe(mysim)

  # Wrap coordinates of the sodium around cterm

  x_cterm = [ x[cterm[1]], y[cterm[1]], z[cterm[1]] ]
  M3GTools.wrap!(sides,x,y,z;center=x_cterm,sel=sod)

  # Compute distance

  d[i] = sqrt( ( x[cterm[1]] - x[sod[1]] )^2 +
               ( y[cterm[1]] - y[sod[1]] )^2 +
               ( z[cterm[1]] - z[sod[1]] )^2 )

end
M3GTools.closedcd(mysim)

println(" Plotting... ")

x = [ i for i in 1:mysim.nframes ]
plot(x,d,xlabel="Frame",ylabel="Distance")
savefig("distance.pdf")
println(" Created figure: distance.pdf")


text:
struct RandomPQs <: SequentialPQs end

function qs_score(qs::RandomPQs, x::Array{T, 2}, labels::Dict{Symbol, Array{Int, 1}})::Array{Float64, 1} where T <: Real
    return rand(size(x, 2))
end


text:
VGG_CONFIGS = Dict(
    "A" => [
        64, "M",
        128, "M",
        256, 256, "M",
        512, 512, "M",
        512, 512, "M"
    ],
    "B" => [
        64, 64, "M",
        128, 128, "M",
        256, 256, "M",
        512, 512, "M",
        512, 512, "M"
    ],
    "D" => [
        64, 64, "M",
        128, 128, "M",
        256, 256, 256, "M",
        512, 512, 512, "M",
        512, 512, 512, "M"
    ],
    "E" =>  [
        64, 64, "M",
        128, 128, "M",
        256, 256, 256, 256, "M",
        512, 512, 512, 512, "M",
        512, 512, 512, 512, "M"
    ],
)

VGG_URLS = Dict(
    "A" => "https =>//download.pytorch.org/models/vgg11-bbd30ac9.pth",
    "B" => "https =>//download.pytorch.org/models/vgg13-c768596a.pth",
    # "D" => "https =>//download.pytorch.org/models/vgg16-6487cb99.mat",
    "D" => "http://www.vlfeat.org/matconvnet/models/imagenet-vgg-verydeep-16.mat",
    "E" => "https =>//download.pytorch.org/models/vgg19-dcbb9e9d.pth",
)

struct VGG
    features::Sequential
    classifier::Sequential
end

function VGG(; features::Sequential, num_classes::Int=1000)
    classifier = Sequential([
        Linear(in_features=512*7*7, out_features=4096),
        Relu(),
        Linear(in_features=4096, out_features=4096),
        Relu(),
        Linear(in_features=4096, out_features=num_classes),
    ])

    return VGG(features, classifier)
end

function (vgg::VGG)(x)
    x = vgg.features(x)
    x = mat(x)
    x = vgg.classifier(x)
end

function topk(scores, k)
    scores = vec(Array(scores))
    ranked = sortperm(scores, rev=true)[1:k]
end

function make_vgg_layers(config)
    layers = []
    in_channels = 3
    for out_channels in config
        if out_channels == "M"
            push!(layers, MaxPool())
        else
            conv = Conv(
                in_channels=in_channels,
                out_channels=out_channels,
                kernel_size=3,
                padding=1,
                mode=1)
            append!(layers, (conv, Relu()))
            in_channels = out_channels
        end
    end
    return Sequential(layers)
end

function _vgg(config; pretrained::Bool=false)
    vgg = VGG(features=make_vgg_layers(VGG_CONFIGS[config]))
    if pretrained
        url = VGG_URLS[config]
        filepath = dir("..", "data", split(url, "/")[end])
        !isfile(filepath) && download(url, filepath)
        _load_weights!(vgg, matread(filepath))
    end
    return vgg
end

function _load_weights!(vgg::VGG, matdata)
    _layers = matdata["layers"]
    features, classifier = _layers[1:end-6], _layers[end-5:end-1]
    for (layer, pretrained) in zip(vgg.features.layers, features)
        length(pretrained["weights"]) != 2 && continue
        weight, bias = pretrained["weights"]
        layer.weight.value .= Knet.atype(weight)
        layer.bias.value .= reshape(Knet.atype(bias), size(layer.bias))
    end
    for (layer, pretrained) in zip(vgg.classifier.layers, classifier)
        length(pretrained["weights"]) != 2 && continue
        weight, bias = pretrained["weights"]
        layer.weight.value .= Knet.atype(transpose(mat(weight)))
        layer.bias.value .= reshape(Knet.atype(bias), size(layer.bias))
    end
end

function _get_imagenet_metadata()
    url = VGG_URLS["D"]
    filepath = dir("..", "data", split(url, "/")[end])
    !isfile(filepath) && download(url, filepath)
    meta = matread(filepath)["meta"]
    description = meta["classes"]["description"]
    average_image = meta["normalization"]["averageImage"]
    return average_image, description
end

vgg11(; pretrained::Bool=false) = _vgg("A"; pretrained=pretrained)
vgg13(; pretrained::Bool=false) = _vgg("B"; pretrained=pretrained)
vgg16(; pretrained::Bool=false) = _vgg("D"; pretrained=pretrained)
vgg19(; pretrained::Bool=false) = _vgg("E"; pretrained=pretrained)


text:
using Gaston
g = Gaston;
set_terminal("x11");

c = g.CurveConf();
c.legend = "Random";
c.plotstyle = "errorbars";

y = exp(-(1:.1:4.9));
ylow = y - 0.05*rand(40);
yhigh = y + 0.05*rand(40);
g.addcoords(1:40,y,c);
g.adderror(0.1*rand(40));

a = g.AxesConf();
a.title = "Error bars (ydelta)";
g.addconf(a);
g.llplot();


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

# Various Unicode functionality from the utf8proc library
module UTF8proc

import Base: show, ==, hash, string, Symbol, isless, length, eltype, start, next, done, convert, isvalid, lowercase, uppercase, titlecase

export isgraphemebreak, category_code, category_abbrev, category_string

# also exported by Base:
export normalize_string, graphemes, is_assigned_char, charwidth, isvalid,
   islower, isupper, isalpha, isdigit, isnumber, isalnum,
   iscntrl, ispunct, isspace, isprint, isgraph

# whether codepoints are valid Unicode scalar values, i.e. 0-0xd7ff, 0xe000-0x10ffff
isvalid(::Type{Char}, ch::Unsigned) = !((ch - 0xd800 < 0x800) | (ch > 0x10ffff))
isvalid(::Type{Char}, ch::Integer) = isvalid(Char, Unsigned(ch))
isvalid(::Type{Char}, ch::Char) = isvalid(Char, UInt32(ch))

isvalid(ch::Char) = isvalid(Char, ch)

# utf8 category constants
const UTF8PROC_CATEGORY_CN = 0
const UTF8PROC_CATEGORY_LU = 1
const UTF8PROC_CATEGORY_LL = 2
const UTF8PROC_CATEGORY_LT = 3
const UTF8PROC_CATEGORY_LM = 4
const UTF8PROC_CATEGORY_LO = 5
const UTF8PROC_CATEGORY_MN = 6
const UTF8PROC_CATEGORY_MC = 7
const UTF8PROC_CATEGORY_ME = 8
const UTF8PROC_CATEGORY_ND = 9
const UTF8PROC_CATEGORY_NL = 10
const UTF8PROC_CATEGORY_NO = 11
const UTF8PROC_CATEGORY_PC = 12
const UTF8PROC_CATEGORY_PD = 13
const UTF8PROC_CATEGORY_PS = 14
const UTF8PROC_CATEGORY_PE = 15
const UTF8PROC_CATEGORY_PI = 16
const UTF8PROC_CATEGORY_PF = 17
const UTF8PROC_CATEGORY_PO = 18
const UTF8PROC_CATEGORY_SM = 19
const UTF8PROC_CATEGORY_SC = 20
const UTF8PROC_CATEGORY_SK = 21
const UTF8PROC_CATEGORY_SO = 22
const UTF8PROC_CATEGORY_ZS = 23
const UTF8PROC_CATEGORY_ZL = 24
const UTF8PROC_CATEGORY_ZP = 25
const UTF8PROC_CATEGORY_CC = 26
const UTF8PROC_CATEGORY_CF = 27
const UTF8PROC_CATEGORY_CS = 28
const UTF8PROC_CATEGORY_CO = 29

# strings corresponding to the category constants
const category_strings = [
    "Other, not assigned",
    "Letter, uppercase",
    "Letter, lowercase",
    "Letter, titlecase",
    "Letter, modifier",
    "Letter, other",
    "Mark, nonspacing",
    "Mark, spacing combining",
    "Mark, enclosing",
    "Number, decimal digit",
    "Number, letter",
    "Number, other",
    "Punctuation, connector",
    "Punctuation, dash",
    "Punctuation, open",
    "Punctuation, close",
    "Punctuation, initial quote",
    "Punctuation, final quote",
    "Punctuation, other",
    "Symbol, math",
    "Symbol, currency",
    "Symbol, modifier",
    "Symbol, other",
    "Separator, space",
    "Separator, line",
    "Separator, paragraph",
    "Other, control",
    "Other, format",
    "Other, surrogate",
    "Other, private use"
]

const UTF8PROC_STABLE    = (1<<1)
const UTF8PROC_COMPAT    = (1<<2)
const UTF8PROC_COMPOSE   = (1<<3)
const UTF8PROC_DECOMPOSE = (1<<4)
const UTF8PROC_IGNORE    = (1<<5)
const UTF8PROC_REJECTNA  = (1<<6)
const UTF8PROC_NLF2LS    = (1<<7)
const UTF8PROC_NLF2PS    = (1<<8)
const UTF8PROC_NLF2LF    = (UTF8PROC_NLF2LS | UTF8PROC_NLF2PS)
const UTF8PROC_STRIPCC   = (1<<9)
const UTF8PROC_CASEFOLD  = (1<<10)
const UTF8PROC_CHARBOUND = (1<<11)
const UTF8PROC_LUMP      = (1<<12)
const UTF8PROC_STRIPMARK = (1<<13)

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

utf8proc_error(result) = error(unsafe_string(ccall(:utf8proc_errmsg, Cstring, (Cssize_t,), result)))

function utf8proc_map(str::String, options::Integer)
    nwords = ccall(:utf8proc_decompose, Int, (Ptr{UInt8}, Int, Ptr{UInt8}, Int, Cint),
                   str, sizeof(str), C_NULL, 0, options)
    nwords < 0 && utf8proc_error(nwords)
    buffer = Base.StringVector(nwords*4)
    nwords = ccall(:utf8proc_decompose, Int, (Ptr{UInt8}, Int, Ptr{UInt8}, Int, Cint),
                   str, sizeof(str), buffer, nwords, options)
    nwords < 0 && utf8proc_error(nwords)
    nbytes = ccall(:utf8proc_reencode, Int, (Ptr{UInt8}, Int, Cint), buffer, nwords, options)
    nbytes < 0 && utf8proc_error(nbytes)
    return String(resize!(buffer, nbytes))
end

utf8proc_map(s::AbstractString, flags::Integer) = utf8proc_map(String(s), flags)

function normalize_string(s::AbstractString; stable::Bool=false, compat::Bool=false, compose::Bool=true, decompose::Bool=false, stripignore::Bool=false, rejectna::Bool=false, newline2ls::Bool=false, newline2ps::Bool=false, newline2lf::Bool=false, stripcc::Bool=false, casefold::Bool=false, lump::Bool=false, stripmark::Bool=false)
    flags = 0
    stable && (flags = flags | UTF8PROC_STABLE)
    compat && (flags = flags | UTF8PROC_COMPAT)
    if decompose
        flags = flags | UTF8PROC_DECOMPOSE
    elseif compose
        flags = flags | UTF8PROC_COMPOSE
    elseif compat || stripmark
        throw(ArgumentError("compat=true or stripmark=true require compose=true or decompose=true"))
    end
    stripignore && (flags = flags | UTF8PROC_IGNORE)
    rejectna && (flags = flags | UTF8PROC_REJECTNA)
    newline2ls + newline2ps + newline2lf > 1 && throw(ArgumentError("only one newline conversion may be specified"))
    newline2ls && (flags = flags | UTF8PROC_NLF2LS)
    newline2ps && (flags = flags | UTF8PROC_NLF2PS)
    newline2lf && (flags = flags | UTF8PROC_NLF2LF)
    stripcc && (flags = flags | UTF8PROC_STRIPCC)
    casefold && (flags = flags | UTF8PROC_CASEFOLD)
    lump && (flags = flags | UTF8PROC_LUMP)
    stripmark && (flags = flags | UTF8PROC_STRIPMARK)
    utf8proc_map(s, flags)
end

"""
    normalize_string(s::AbstractString, normalform::Symbol)

Normalize the string `s` according to one of the four "normal forms" of the Unicode
standard: `normalform` can be `:NFC`, `:NFD`, `:NFKC`, or `:NFKD`.  Normal forms C
(canonical composition) and D (canonical decomposition) convert different visually identical
representations of the same abstract string into a single canonical form, with form C being
more compact.  Normal forms KC and KD additionally canonicalize "compatibility equivalents":
they convert characters that are abstractly similar but visually distinct into a single
canonical choice (e.g. they expand ligatures into the individual characters), with form KC
being more compact.

Alternatively, finer control and additional transformations may be be obtained by calling
`normalize_string(s; keywords...)`, where any number of the following boolean keywords
options (which all default to `false` except for `compose`) are specified:

* `compose=false`: do not perform canonical composition
* `decompose=true`: do canonical decomposition instead of canonical composition (`compose=true`
  is ignored if present)
* `compat=true`: compatibility equivalents are canonicalized
* `casefold=true`: perform Unicode case folding, e.g. for case-insensitive string comparison
* `newline2lf=true`, `newline2ls=true`, or `newline2ps=true`: convert various newline sequences
  (LF, CRLF, CR, NEL) into a linefeed (LF), line-separation (LS), or paragraph-separation (PS)
  character, respectively
* `stripmark=true`: strip diacritical marks (e.g. accents)
* `stripignore=true`: strip Unicode's "default ignorable" characters (e.g. the soft hyphen
  or the left-to-right marker)
* `stripcc=true`: strip control characters; horizontal tabs and form feeds are converted to
  spaces; newlines are also converted to spaces unless a newline-conversion flag was specified
* `rejectna=true`: throw an error if unassigned code points are found
* `stable=true`: enforce Unicode Versioning Stability

For example, NFKC corresponds to the options `compose=true, compat=true, stable=true`.
"""
function normalize_string(s::AbstractString, nf::Symbol)
    utf8proc_map(s, nf == :NFC ? (UTF8PROC_STABLE | UTF8PROC_COMPOSE) :
                    nf == :NFD ? (UTF8PROC_STABLE | UTF8PROC_DECOMPOSE) :
                    nf == :NFKC ? (UTF8PROC_STABLE | UTF8PROC_COMPOSE
                                   | UTF8PROC_COMPAT) :
                    nf == :NFKD ? (UTF8PROC_STABLE | UTF8PROC_DECOMPOSE
                                   | UTF8PROC_COMPAT) :
                    throw(ArgumentError(":$nf is not one of :NFC, :NFD, :NFKC, :NFKD")))
end

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

"""
    charwidth(c)

Gives the number of columns needed to print a character.
"""
charwidth(c::Char) = Int(ccall(:utf8proc_charwidth, Cint, (UInt32,), c))

lowercase(c::Char) = isascii(c) ? ('A' <= c <= 'Z' ? c + 0x20 : c) : Char(ccall(:utf8proc_tolower, UInt32, (UInt32,), c))
uppercase(c::Char) = isascii(c) ? ('a' <= c <= 'z' ? c - 0x20 : c) : Char(ccall(:utf8proc_toupper, UInt32, (UInt32,), c))
titlecase(c::Char) = isascii(c) ? ('a' <= c <= 'z' ? c - 0x20 : c) : Char(ccall(:utf8proc_totitle, UInt32, (UInt32,), c))

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

# returns UTF8PROC_CATEGORY code in 0:30 giving Unicode category
category_code(c) = ccall(:utf8proc_category, Cint, (UInt32,), c)

# more human-readable representations of the category code
category_abbrev(c) = unsafe_string(ccall(:utf8proc_category_string, Cstring, (UInt32,), c))
category_string(c) = category_strings[category_code(c)+1]

"""
    is_assigned_char(c) -> Bool

Returns `true` if the given char or integer is an assigned Unicode code point.
"""
is_assigned_char(c) = category_code(c) != UTF8PROC_CATEGORY_CN

## libc character class predicates ##

"""
    islower(c::Char) -> Bool

Tests whether a character is a lowercase letter.
A character is classified as lowercase if it belongs to Unicode category Ll,
Letter: Lowercase.
"""
islower(c::Char) = (category_code(c) == UTF8PROC_CATEGORY_LL)

# true for Unicode upper and mixed case

"""
    isupper(c::Char) -> Bool

Tests whether a character is an uppercase letter.
A character is classified as uppercase if it belongs to Unicode category Lu,
Letter: Uppercase, or Lt, Letter: Titlecase.
"""
function isupper(c::Char)
    ccode = category_code(c)
    return ccode == UTF8PROC_CATEGORY_LU || ccode == UTF8PROC_CATEGORY_LT
end

"""
    isdigit(c::Char) -> Bool

Tests whether a character is a numeric digit (0-9).
"""
isdigit(c::Char)  = ('0' <= c <= '9')

"""
    isalpha(c::Char) -> Bool

Tests whether a character is alphabetic.
A character is classified as alphabetic if it belongs to the Unicode general
category Letter, i.e. a character whose category code begins with 'L'.
"""
isalpha(c::Char)  = (UTF8PROC_CATEGORY_LU <= category_code(c) <= UTF8PROC_CATEGORY_LO)

"""
    isnumber(c::Char) -> Bool

Tests whether a character is numeric.
A character is classified as numeric if it belongs to the Unicode general category Number,
i.e. a character whose category code begins with 'N'.
"""
isnumber(c::Char) = (UTF8PROC_CATEGORY_ND <= category_code(c) <= UTF8PROC_CATEGORY_NO)

"""
    isalnum(c::Char) -> Bool

Tests whether a character is alphanumeric.
A character is classified as alphabetic if it belongs to the Unicode general
category Letter or Number, i.e. a character whose category code begins with 'L' or 'N'.
"""
function isalnum(c::Char)
    ccode = category_code(c)
    return (UTF8PROC_CATEGORY_LU <= ccode <= UTF8PROC_CATEGORY_LO) ||
           (UTF8PROC_CATEGORY_ND <= ccode <= UTF8PROC_CATEGORY_NO)
end

# following C++ only control characters from the Latin-1 subset return true

"""
    iscntrl(c::Char) -> Bool

Tests whether a character is a control character.
Control characters are the non-printing characters of the Latin-1 subset of Unicode.
"""
iscntrl(c::Char) = (c <= Char(0x1f) || Char(0x7f) <= c <= Char(0x9f))

"""
    ispunct(c::Char) -> Bool

Tests whether a character belongs to the Unicode general category Punctuation, i.e. a
character whose category code begins with 'P'.
"""
ispunct(c::Char) = (UTF8PROC_CATEGORY_PC <= category_code(c) <= UTF8PROC_CATEGORY_PO)

# \u85 is the Unicode Next Line (NEL) character

"""
    isspace(c::Char) -> Bool

Tests whether a character is any whitespace character. Includes ASCII characters '\\t',
'\\n', '\\v', '\\f', '\\r', and ' ', Latin-1 character U+0085, and characters in Unicode
category Zs.
"""
@inline isspace(c::Char) = c == ' ' || '\t' <= c <='\r' || c == '\u85' || '\ua0' <= c && category_code(c) == UTF8PROC_CATEGORY_ZS

"""
    isprint(c::Char) -> Bool

Tests whether a character is printable, including spaces, but not a control character.
"""
isprint(c::Char) = (UTF8PROC_CATEGORY_LU <= category_code(c) <= UTF8PROC_CATEGORY_ZS)

# true in principal if a printer would use ink

"""
    isgraph(c::Char) -> Bool

Tests whether a character is printable, and not a space.
Any character that would cause a printer to use ink should be
classified with `isgraph(c)==true`.
"""
isgraph(c::Char) = (UTF8PROC_CATEGORY_LU <= category_code(c) <= UTF8PROC_CATEGORY_SO)

############################################################################
# iterators for grapheme segmentation

isgraphemebreak(c1::Char, c2::Char) =
    ccall(:utf8proc_grapheme_break, Bool, (UInt32, UInt32), c1, c2)

# Stateful grapheme break required by Unicode-9 rules: the string
# must be processed in sequence, with state initialized to Ref{Int32}(0).
# Requires utf8proc v2.0 or later.
isgraphemebreak!(state::Ref{Int32}, c1::Char, c2::Char) =
    ccall(:utf8proc_grapheme_break_stateful, Bool, (UInt32, UInt32, Ref{Int32}), c1, c2, state)

struct GraphemeIterator{S<:AbstractString}
    s::S # original string (for generation of SubStrings)
end

"""
    graphemes(s::AbstractString) -> GraphemeIterator

Returns an iterator over substrings of `s` that correspond to the extended graphemes in the
string, as defined by Unicode UAX #29. (Roughly, these are what users would perceive as
single characters, even though they may contain more than one codepoint; for example a
letter combined with an accent mark is a single grapheme.)
"""
graphemes(s::AbstractString) = GraphemeIterator{typeof(s)}(s)

eltype{S}(::Type{GraphemeIterator{S}}) = SubString{S}

function length(g::GraphemeIterator)
    c0 = Char(0x00ad) # soft hyphen (grapheme break always allowed after this)
    n = 0
    state = Ref{Int32}(0)
    for c in g.s
        n += isgraphemebreak!(state, c0, c)
        c0 = c
    end
    return n
end

start(g::GraphemeIterator) = (start(g.s), Ref{Int32}(0))
done(g::GraphemeIterator, i) = done(g.s, i[1])

function next(g::GraphemeIterator, i_)
    s = g.s
    i, state = i_
    j = i
    c0, k = next(s, i)
    while !done(s, k) # loop until next grapheme is s[i:j]
        c, ℓ = next(s, k)
        isgraphemebreak!(state, c0, c) && break
        j = k
        k = ℓ
        c0 = c
    end
    return (SubString(s, i, j), (k, state))
end

==(g1::GraphemeIterator, g2::GraphemeIterator) = g1.s == g2.s
hash(g::GraphemeIterator, h::UInt) = hash(g.s, h)
isless(g1::GraphemeIterator, g2::GraphemeIterator) = isless(g1.s, g2.s)

convert{S<:AbstractString}(::Type{S}, g::GraphemeIterator) = convert(S, g.s)

show{S}(io::IO, g::GraphemeIterator{S}) = print(io, "length-$(length(g)) GraphemeIterator{$S} for \"$(g.s)\"")

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

end # module


text:
# logical functions

"""
Can use matrix mul to create constraints:
julia> @constraint(model, con, A * x .== b)
2-element Array{ConstraintRef{Model,MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64},MathOptInterface.EqualTo{Float64}},ScalarShape},1}:
 x[1] + 2 x[2] == 5.0
 3 x[1] + 4 x[2] == 6.0
"""
mutable struct Constraint
    # e.g. 5*x + 3*theta < 3.2
    coeffs::Array{T} where T <: Real # 5, 3
    relation::Symbol # <
    scalar::Real # 3.2
end

# mutable struct Property
#     temporal # e.g. "always_true" or "eventually_true"
#     type # e.g. "safe" or "avoid"
#     constraint # e.g. hyperrectangle, for now, or list of Constraints
# end


text:
module Hilbert

"""
    hilbert(x)
    hilbert(x, n)
Analytic signal, computed using the Hilbert transform.
`y = hilbert(x)` gives the analytic signal `y = x + i*xi` where `xi` is the
Hilbert transform of vector `x`. If `x` is complex, only the real part is used.
If `x` is a matrix, then `hilbert` operates along columns. If n is provided, the
n-point Hilbert transform is computed.

**Example**
```julia
x = randn(10,10)
y = hilbert(x)
```
"""
function hilbert{T<:Real}(x::AbstractArray{T,2}, n::Int64=0)

  x_ = copy(x)

  if(eltype(x_) <: Complex)
    warn("Using real part, ignoring complex part")
    x_ = real(x_)
  end

  # work along columns
  size(x_,1)==1 ? x_ = permutedims(x_,[2 1]) : nothing

  if(n>0 && n<size(x_,1))
    x_ = x_[1:n,:]
  elseif(n>0 && n>size(x_,1))
    x_ = cat(1,x_,zeros(n-size(x_,1),size(x_,2)))
  else
    n = size(x_,1)
  end

  xf = fft(x_,1)
  h = zeros(Int64,n)
  if n>0 && n % 2 == 0
    #even, nonempty
    h[1:div(n,2)+1] = 1
    h[2:div(n,2)] = 2
  elseif n>0
    #odd, nonempty
    h[1] = 1
    h[2:div(n + 1,2)] = 2
  end
  x_ = ifft(xf .* h[:,ones(Int64,size(xf,2))],1)

  # restore to original shape if necessary
  size(x,1)==1 ? x_ = permutedims(x_,[2 1]) : nothing

  return x_

end

function hilbert{T<:Complex}(x::AbstractArray{T}, n::Int64=0)
  warn("Using real part, ignoring complex part")
  Hilbert.hilbert(real(x))
end

end # module


text:
g = EvolvingGraph{Node{String},String}()

add_edge!(g, "a", "b", "t1")
add_edge!(g, "b", "c", "t1")
add_edge!(g, "c", "d", "t2")
add_edge!(g, "a", "b", "t2")

egwrite(g, "test1.csv")
g = egread("test1.csv")
@test num_nodes(g) == 4
@test num_timestamps(g) == 2

rm("test1.csv")

# twitter data
twitterdata = joinpath(Pkg.dir("EvolvingGraphs"), "data", "manunited_cont.csv")
twitter = egread(twitterdata)
ns = nodes(twitter)
es = edges(twitter)
ts = timestamps(twitter)


text:
__precompile__()

module Intervals

using Compat
using Compat.Dates
using Compat.Printf
using TimeZones

using Compat.Dates: value, coarserperiod
using Compat: AbstractDateTime

import Base: ⊆, ⊇, ⊈, ⊉, union, union!, merge

abstract type AbstractInterval{T} end

Base.eltype(::AbstractInterval{T}) where {T} = T

if VERSION >= v"0.7"
    Base.broadcastable(x::AbstractInterval) = Ref(x)
end

include("inclusivity.jl")
include("endpoint.jl")
include("interval.jl")
include("anchoredinterval.jl")
include("description.jl")

export AbstractInterval,
       Interval,
       AnchoredInterval,
       HourEnding,
       HourBeginning,
       HE,
       HB,
       Inclusivity,
       inclusivity,
       first,
       last,
       isclosed,
       anchor,
       span,
       merge,
       union,
       union!,
       less_than_disjoint,
       greater_than_disjoint,
       superset,
       .., ≪, ≫, ⊆, ⊇, ⊈, ⊉
end


text:
using Images, ImageView

PROJECT_IMG_FOLDER = joinpath("data", "book-covers")
DATASETS = ["test", "train"]

preview_count = 5
shape = (100, 78)
preview_img = zeros(typeof(RGB4{N0f8}(0.,0.,0.)), (100, 1));

original_shapes = []

for dataset_name in DATASETS

    current_folder = joinpath(PROJECT_IMG_FOLDER, dataset_name)
    files = filter(x -> contains(x, ".jpg"), readdir(current_folder))

    shuffle!(files)

    for i=1:preview_count

        img_original = load(joinpath(current_folder, files[i]))
        img_resized = imresize(img_original, shape)

        preview_img = hcat(preview_img, img_resized)
        push!(original_shapes, size(img_original))
    end
end

imshow(preview_img)
println(original_shapes)


text:
# Simulation submodules
include("integrators.jl")
include("propagators.jl")


text:
module QuDemonCooling
using QuBase
import QuBase.AbstractQuVector
export cooling!

# for general use
function cooling!(state::AbstractQuVector,U::AbstractMatrix,t::Real,gamma::Real)
    dice = rand()
    if dice <= 0.5*(1+sin(gamma))
        heat!(state,U,t,gamma)
        return 0.5*(1+sin(gamma))
    else
        cool!(state,U,t,gamma)
        return 0.5*(1-sin(gamma))
    end
end

function heat!(state::AbstractQuVector,U::AbstractMatrix,t::Real,gamma::Real)
    next_state = 0.5*( coeffs(state)+im*exp(im*gamma)*U*coeffs(state) ) |> normalize!
    CQST = QuBase.similar_type(state)
    return CQST(next_state,bases(state))
end

function cool!(state::AbstractQuVector,U::AbstractMatrix,t::Real,gamma::Real)
    next_state = 0.5*( coeffs(state)-im*exp(im*gamma)*U*coeffs(state) ) |> normalize!
    CQST = QuBase.similar_type(state)
    return CQST(next_state,bases(state))
end

end # module


text:
using Distances
using Test

# I copied these tests :x

@testset "Distances tests" begin
    @test GetDistance(Euclidean::DistanceType, [2,2], [2,4]) == 2
    @test GetDistance(Cityblock::DistanceType, [2,2], [2,4]) == 2
    @test GetDistance(TotalVariation::DistanceType, [2,2], [2,4]) == 1
    @test GetDistance(Chebyshev::DistanceType, [2,2], [2,4]) == 2
    @test GetDistance(Jaccard::DistanceType, [2,2], [2,4]) == 1 - 2/3
    @test GetDistance(BrayCurtis::DistanceType, [2,2], [2,4]) == 1/5
    @test GetDistance(SpanNormDist::DistanceType, [2,2], [2,4]) == 2
end

@testset "Errors test" begin
    @test_throws DimensionMismatch GetDistance(Cityblock::DistanceType, [2, 2, 2], [2, 4])
    @test_throws DimensionMismatch GetDistance(Cityblock::DistanceType, [2], [2, 4])
    @test_throws MethodError GetDistance(Cityblock::DistanceType, [2, "hello"], [2, 4])
    @test_throws MethodError GetDistance(Cityblock::DistanceType, [], [2, 4])
    @test_throws MethodError GetDistance(Cityblock::DistanceType, [2], [])
end