text:
using LinearAlgebra
import Base: isless
@enum EdgeType CircEdge = 1 InsideEdge = 2

"""
    AbstractInterpolation{T, edge, order, nd}

Abstract supertype for all interpolation type

# Type parameters
- `T` : type of number on witch interpolation works
- `edge` : type of edge treatment
- `order` : order of interpolation

# Implementation constraint
- `tabfct::Vector` : this attribut must be on the implementation, it is a table of function of size order+1
"""
abstract type AbstractInterpolation{T,edge,order} end

"""
    get_order(_::AbstractInterpolation{T, edge, order}) where{T, edge, order}
Return the order of interpolation implementation
"""
get_order(_::AbstractInterpolation{T,edge,order}) where {T,edge,order} = order

"""
    sol(_::AbstractInterpolation, line::AbstractVector)

Interface method to transform the treated line, by default this method does nothing

# Arguments :
- `_::AbstractInterpolation` : interpolation implementation
- `line::AbstractVector` : line to transform

# Return :
The transformed line

"""
sol(_::AbstractInterpolation, b::AbstractArray) = b

issolidentity(_::AbstractInterpolation) = true

isbspline(_::AbstractInterpolation) = false

Base.show(io::IO, interp::AbstractInterpolation) = print(io, typeof(interp))

function sol!(
    Y::AbstractArray{T2,N},
    interp_t::AbstractVector{I},
    b::AbstractArray{T2,N},
) where {T,N,I<:AbstractInterpolation{T},T2<:Union{T,OpTuple{N,T}}}
    sz = size(Y)
    p = circshift(1:N, -1)
    perm = p
    bufout = b
    for i = 1:N
        interp = interp_t[i]
        if issolidentity(interp)
            if i == 1
                bufout = b
            else
                perm = perm[p]
            end
        else
            bufin = (i == 1) ? b : permutedims(bufout, perm)
            bufout = zeros(T2, sz)
            for c in CartesianIndices(sz[2:end])
                bo = view(bufout, :, c)
                bi = view(bufin, :, c)
                sol!(bo, interp, bi)
            end
            perm = p
        end
        sz = sz[p]
    end
    permutedims!(Y, bufout, perm)
    return Y
end
function sol(
    interp_t::AbstractVector{I},
    b::Union{AbstractArray{T,N},AbstractArray{Complex{T},N},AbstractArray{OpTuple{N,T},N}},
) where {T<:Real,N,I<:AbstractInterpolation{T}}
    if all(issolidentity.(interp_t))
        return b
    else
        return sol!(zeros(eltype(b), size(b)), interp_t, b)
    end
end

@inline function getprecal(interp::AbstractInterpolation{T}, decf::T) where {T}
    return [T(fct(decf)) for fct in interp.tabfct]
end

@inline function getprecal!(v::Vector{T}, bsp::AbstractInterpolation{T}, decf::T) where {T}
    v .= getprecal(bsp, decf)
end
@inline function getprecal!(
    v::Vector{T},
    bsp::Vector{I},
    decf::T,
) where {T,I<:AbstractInterpolation{T}}
    getprecal!(v, bsp[1], decf)
end
@inline function getprecal!(
    v::Array{T,N},
    bsp::Vector{I},
    decf::NTuple{N,T},
) where {T,N,I<:AbstractInterpolation{T}}
    v .= dotprod(ntuple(x -> getprecal(bsp[x], decf[x]), N))
end

function get_allprecal(
    interp::AbstractInterpolation{T,InsideEdge,order},
    decint::Int,
    decfloat::T,
) where {T,order}
    origin = -div(order, 2)
    indbeg = origin + decint
    indend = indbeg + order
    return [getprecal(interp, decfloat + i) for i = indbeg:indend]
end
struct ValInv{T}
    a::Int
    b::T
    val::Complex{T}
    ind::CartesianIndex{2}
    function ValInv(ind::CartesianIndex{2}, val::Complex{T}) where {T}
        return new{T}(round(Int, real(val)), imag(val), val, ind)
    end
end
function Base.isless(a::ValInv{T}, b::ValInv{T}) where {T}
    return a.a == b.a ? Base.isless(a.b, b.b) : Base.isless(a.a, b.a)
end

function corind(indref::CartesianIndex{2}, indorigin::CartesianIndex{2}, sz::Tuple{Int,Int})
    return CartesianIndex(modone.(indref.I .- indorigin.I, sz))
end

"""
    interpolate!( fp::AbstractVector{T},
        fi::AbstractVector{T},
        decint::Int,
        precal::Vector{T},
        interp::AbstractInterpolation{T, CircEdge, order},
        self::AdvectionData{T},
        tabmod=gettabmod(size(fi)) ) where {T, order}

apply an offset to the function fi interpolate by interp struct, the result is in fp vector,
decint and precal are precompute with getprecal method, the TypeEdge is CircEdge

# Arguments
- `fp::AbstractVector` : output vector
- `fi::AbstractVector` : input vector
- `decint` : offset in units of dx
- `precal::Vector` : vector of length order+1 precompute with getprecal(interp, dec) (dec is the offset)
- `interp::AbstractInterpolation{T, CircEdge, order}` : interpolation implementation, note that TypeEdge is CircEdge
- `tabmod=gettabmod(length(fi))` : precompute for "begin at one" modulo

# Returns :
- No return
"""
function interpolate!(
    fp::AbstractVector{T},
    fi::AbstractVector{T},
    decint::Int,
    precal::Vector{T},
    interp::AbstractInterpolation{T,CircEdge,order},
    tabmod = gettabmod(length(fi))
) where {T,order}
    res = sol(interp, fi)
    # @show typeof(res)
    # @show typeof(tabmod)
    origin = -div(order, 2)
    lg = length(fi)
    decal = (origin + decint + 5lg) % lg
    @inbounds for i = 1:lg
        indbeg = i + decal
        indend = indbeg + order
        fp[i] = sum(res[tabmod[indbeg:indend]] .* precal)
    end
    missing
end
@inline function interpolate!(
    fp::AbstractVector{T},
    fi::AbstractVector{T},
    decint::Int,
    precal::Vector{T},
    tinterp::Vector{I},
    tabmod = gettabmod(length(fi))
) where {T,I<:AbstractInterpolation{T,CircEdge}}
    return interpolate!(fp, fi, decint, precal, tinterp[1], tabmod[1])
end

function interpolate!(
    fp::AbstractArray{T,N},
    fi::AbstractArray{T,N},
    decint::NTuple{N,Int},
    precal::Array{T,N},
    interp::Vector{I},
    tabmod = gettabmod.(size(fi)),
) where {T, N, I <: AbstractInterpolation{T,CircEdge}}
    res = sol(interp, fi)
    sz = size(fi)
    order = get_order.(interp)
    origin = -div.(order, (2,))
    decall = (5 .* sz .+ origin) .% sz .+ sz
    for ind in CartesianIndices(fi)
        deb_i = decint .+ decall .+ ind.I
        end_i = deb_i .+ order
        fp[ind] = sum(res[ntuple(x -> tabmod[x][deb_i[x]:end_i[x]], N)...] .* precal)
    end
    missing
end

"""
    interpolate!(
    fp::AbstractVector{T}, fi::AbstractVector{T}, decint::Int,
    allprecal::Vector{Vector{T}},
    interp::AbstractInterpolation{T, InsideEdge, order},
    tabmod=gettabmod(length(fi))
    ) where {T, order}

apply an offset to the function fi interpolate by interp struct, the result is in fp vector,
decint and precal are precompute with getprecal method, the TypeEdge is InsideEdge, it is a marginal case

# Arguments
- `fp::AbstractVector` : output vector
- `fi::AbstractVector` : input vector
- `decint` : offset in units of dx
- `allprecal::Vector{Vector{T}}` : vector of vector of length order+1 precompute with getprecal(interp, dec) (dec is the offset)
- `interp::AbstractInterpolation{T, InsideEdge, order}` : interpolation implementation, note that TypeEdge is CircEdge
- `tabmod=gettabmod(length(fi))` : precompute for "begin at one" modulo

# Returns :
- No return
"""
function interpolate!(
    fp::AbstractVector{T},
    fi::AbstractVector{T},
    decint::Int,
    allprecal::Vector{Vector{T}},
    interp::AbstractInterpolation{T,InsideEdge,order},
    tabmod = gettabmod(length(fi));
) where {T,order}
    res = sol(interp, fi)
    origin = -div(order, 2)
    lg = length(fi)
    lgp = length(allprecal)
    borne1 = -decint - origin
    borne2 = lg - decint + origin - 1
    for i = 1:borne1
        indbeg = 1
        indend = order + 1
        ind = i
        #        @show 1, i, ind, indbeg, indend, decint, lgp
        fp[i] = sum(res[indbeg:indend] .* allprecal[ind])
    end
    for i = borne1+1:borne2
        indbeg = i - borne1
        indend = indbeg + order
        ind = borne1 + 1
        #        @show 2, i, ind, indbeg, indend, decint, lgp
        fp[i] = sum(res[indbeg:indend] .* allprecal[ind])
    end
    for i = borne2+1:lg
        indbeg = lg - order
        indend = lg
        ind = lgp - (lg - i)
        #        @show 3, i, ind, indbeg, indend, decint, lgp
        fp[i] = sum(res[indbeg:indend] .* allprecal[ind])
    end
    missing
end

"""
    interpolate!( fp, fi, dec, interp)

apply the offset dec to the function fi interpolate by interp struct, the result is in fp Vector

# Arguments
- `fp` : output vector of length n
- `fi` : input vector of length n
- `dec` : offset in units of dx
- `interp::AbstractInterpolation` : interpolation implementation

# Returns :
- No return
"""
function interpolate!(
    fp,
    fi,
    dec,
    interp::AbstractInterpolation{T,edge,order},
) where {T,edge,order}
    decint = convert(Int, floor(dec))
    decfloat = dec - decint
    if edge == CircEdge
        return interpolate!(fp, fi, decint, getprecal(interp, decfloat), interp)
    else
        return interpolate!(fp, fi, decint, get_allprecal(interp, decint, decfloat), interp)
    end
end

mutable struct CachePrecal{T,N,I}
    interps::Vector{I}
    cache_alpha::Union{NTuple{N,T},T}
    cache_int::Union{NTuple{N,Int},Int}
    precal::Array{T,N}
    function CachePrecal(
        interps::Vector{I},
        ::T,
    ) where {T<:Real,I<:AbstractInterpolation{T}}
        N = length(interps)
        cache_alpha = (N == 1) ? zero(T) : ntuple(x -> zero(T), N)
        cache_int = (N == 1) ? 0 : ntuple(x -> 0, N)
        sz = totuple(get_order.(interps) .+ 1)
        precal = zeros(T, sz)
        getprecal!(precal, interps, cache_alpha)
        return new{T,N,I}(interps, cache_alpha, cache_int, precal)
    end
    CachePrecal(
        interps::Vector{I},
        x::Complex{T},
    ) where {T<:Real,I<:AbstractInterpolation{T}} = CachePrecal(interps, real(x))
    CachePrecal(
        interps::Vector{I},
        x::OpTuple{N,T},
    ) where {N,T<:Real,I<:AbstractInterpolation{T}} = CachePrecal(interps, x.v[1])
end
@inline function getprecal(self::CachePrecal{T,N}, alpha::NTuple{N,T}) where {T,N}
    if alpha != self.cache_alpha
        self.cache_alpha = alpha
        self.cache_int = Int.(floor.(alpha))
        decfloat = alpha .- self.cache_int
        getprecal!(self.precal, self.interps, decfloat)
    end
    return self.cache_int, self.precal
end
@inline getprecal(self::CachePrecal{T,2}, alpha::Complex{T}) where {T<:Real} =
    getprecal(self, reim(alpha))
@inline getprecal(self::CachePrecal{T,N}, alpha::OpTuple{N,T}) where {N,T<:Real} =
    getprecal(self, alpha.v)
@inline function getprecal(self::CachePrecal{T,1}, alpha::T) where {T}
    if alpha != self.cache_alpha
        self.cache_alpha = alpha
        self.cache_int = Int(floor(alpha))
        decfloat = alpha - self.cache_int
        getprecal!(self.precal, self.interps[1], decfloat)
    end
    return self.cache_int, self.precal
end
@inline getprecal(self::CachePrecal{T,1}, alpha::Tuple{T}) where {T} =
    getprecal(self, alpha[1])

function interpolate!(
    fp::Union{AbstractArray{T,N},AbstractArray{Complex{T},N}},
    fi::Union{AbstractArray{T,N},AbstractArray{Complex{T},N}},
    dec::Function,
    interp_t::AbstractVector{I};
    tabmod::NTuple{N,Vector{Int}} = gettabmod.(size(fi)),
    cache::CachePrecal{T,N} = CachePrecal(interp_t, one(eltype(fp))),
) where {T,N,I<:AbstractInterpolation{T}}

    N == length(interp_t) || thrown(
        ArgumentError(
            "The number of Interpolation $(length(interp_t)) is different of N=$N",
        ),
    )
    sz = size(fp)
    #    @show "interpolate!", sz
    res = sol(interp_t, fi)

    order = get_order.(interp_t)
    origin = -div.(order, (2,))
    decall = (5 .* sz .+ origin) .% sz .+ sz

    for ind in CartesianIndices(sz)
        dint, tab = getprecal(cache, dec(ind))
        deb_i = dint .+ decall .+ ind.I
        end_i = deb_i + order
        fp[ind] = sum(res[ntuple(x -> tabmod[x][deb_i[x]:end_i[x]], N)...] .* tab)
    end
    #    @show "std2d", diffmax, decminmin, decmaxmax
end

function interpolate!(
    fp::Union{AbstractArray{T,N},AbstractArray{OpTuple{N,T},N},AbstractArray{Complex{T},N}},
    fi::Union{AbstractArray{T,N},AbstractArray{OpTuple{N,T},N},AbstractArray{Complex{T},N}},
    bufdec::Union{AbstractArray{OpTuple{N,T},N},AbstractArray{Complex{T},N}} ,
    interp_t::AbstractVector{I};
    tabmod::NTuple{N,Vector{Int}} = gettabmod.(size(fi)),
    mpid::Union{MPIData,Missing} = missing,
    t_split::Union{Tuple,Missing} = missing,
    cachethreads::Union{Vector{CachePrecal{T}},Missing} = missing,
    #    itr::AbstractArray=CartesianIndices(fi)
) where {T,N,I<:AbstractInterpolation{T}}
    N == length(interp_t) || thrown(
        ArgumentError(
            "The number of Interpolation $(length(interp_t)) is different of N=$N",
        ),
    )

    isthreads = !ismissing(cachethreads)
    sz = size(fp)
    #    @show "interpolate!", sz
    res = sol(interp_t, fi)

    order = get_order.(interp_t)
    origin = -div.(order, (2,))
    decall = (5 .* sz .+ origin) .% sz .+ sz

    function fct(ind::CartesianIndex{N2}, cache::CachePrecal{T2,N2,I2}) where {N2,T2,I2 <: AbstractInterpolation{T2}}
        dint, tab = getprecal(cache, bufdec[ind])
        deb_i = dint .+ decall .+ ind.I
        end_i = deb_i + order
        fp[ind] = sum(res[ntuple(x -> tabmod[x][deb_i[x]:end_i[x]], N2)...] .* tab)

    end
    ci = CartesianIndices(sz)
    if isthreads
        if ismissing(t_split)
#            println("SimpleThreads")
            @threads for ind in ci
                local cache = cachethreads[Threads.threadid()]
                fct(ind, cache)
            end
        else
#            println("SplitThreads")
            @threads for it in t_split
                local cache = cachethreads[Threads.threadid()]
                for ind in ci[it]
                    fct(ind, cache)
                end
            end
        end
    else
        local cache = CachePrecal(interp_t, zero(T))
        itr = ismissing(mpid) ? ci : ci[t_split[mpid.ind]]
        for ind in itr
            fct(ind, cache)
        end
    end

    if !ismissing(mpid)
        mpibroadcast(mpid, t_split, fp)
    end
    true
    #    @show "std2d", diffmax, decminmin, decmaxmax
end
function interpolatemod!(
    fp::AbstractArray{T,N},
    fi::AbstractArray{T,N},
    bufc::Tuple{Array{T,N},Array{T,N}},
    interp_t::AbstractVector{I},
    lgmesh::Union{T,UniformMesh{T}},
) where {T,N,I<:AbstractInterpolation{T,CircEdge}}

    (N == 2 && N == length(interp_t)) || thrown(
        ArgumentError(
            "The number of Interpolation $(length(interp_t)) is different of N=$N",
        ),
    )

    sz = size(fp)

    res = sol(interp_t, fi)

    extr = extrema.(bufc)
    decbegin = ntuple(x -> -Int(floor(extr[x][1])), N)
    decend = ntuple(x -> Int(ceil(extr[x][2])), N)

    order = totuple(get_order.(interp_t))
    origin = 0 .- div.(order, (2,))

    dec_b = decbegin .- origin
    dec_e = decend .+ order .+ origin
    fiext = getextarray(res, dec_b, dec_e)

    decall = decbegin

    traitmodbegin!(lgmesh, fiext)

    cache::CachePrecal{T,N} = CachePrecal(interp_t, one(eltype(fp)))
    diffmax = 0
    decminmin = ntuple(x -> Inf, N)
    decmaxmax = ntuple(x -> -Inf, N)
    decabs = ntuple(x -> 0, N)

    # @show size(fi), decbegin, decend
    # @show order, origin, dec_b, dec_e, size(fiext), decall
    for ind in CartesianIndices(sz)
        dint, tab = getprecal(cache, (bufc[1][ind], bufc[2][ind]))
        decabs = max.(decabs, abs.((bufc[1][ind], bufc[2][ind])))
        decmin, decmax = extrema(tab)
        decminmin = min.(decminmin, dint .+ decmin)
        decmaxmax = max.(decmaxmax, dint .+ decmax)

        deb_i = dint .+ decall .+ ind.I
        end_i = deb_i .+ order
        vmin, vmax = extrema(fiext[ntuple(x -> deb_i[x]:end_i[x], N)...])
        diff = vmax - vmin
        diffmax = max(diff, diffmax)
        fp[ind] = sum(fiext[ntuple(x -> deb_i[x]:end_i[x], N)...] .* tab)
        # if ind.I[1] <= 10 && ind.I[2] <= 10
        #     indcor = ind +CartesianIndex(decall .- origin)
        #     @show decall, dint, origin
        #     @show ind, fp[ind], indcor, fiext[indcor]
        # end
    end
    # @show diffmax, decminmin, decmaxmax, decabs
    # @show fp[1:10,1:10]
    traitmodend!(lgmesh, fp)
    # @show fp[1:10,1:10]
    missing
end

function getinverse(
    dec::NTuple{N,Array{T,N}},
    interp::Vector{I},
) where {T,N,I<:AbstractInterpolation{T,CircEdge}}
    sz = size(dec[1])
    res = ntuple(x -> [T(ind.I[x] - 1) for ind in CartesianIndices(sz)], N)
    decinv = ntuple(x -> [T(ind.I[x] - 1) for ind in CartesianIndices(sz)], N)
    buf = zeros(T, sz)
    res2 = ntuple(x -> zeros(T, sz), N)
    res3 = ntuple(x -> zeros(T, sz), N)
    res4 = ntuple(x -> zeros(T, sz), N)
    decfmr = ntuple(x -> zeros(T, sz), N)

    sz_2 = div.(sz, 2)

    for i = 1:N
        interpolatemod!(res2[i], res[i], dec, interp, T(sz[i]))
        interpolatemod!(decfmr[i], -dec[i], dec, interp, T(sz[i]))
        decfmr[i] .= mod.(decfmr[i] .+ sz_2[i], sz[i]) .- sz_2[i]

        res4[i] .= res2[i]

        res3[i] .= mod.(res2[i] .- res[i] .- dec[i] .+ sz_2[i], sz[i]) .- sz_2[i]

    end

    @show norm(res3)

    borne = log(eps(T) * prod(sz))

    for z = 1:2
        ind = 1
        note = Inf
        while note > borne
            for i = 1:N
                interpolatemod!(res3[i], res2[i], decfmr, interp, T(sz[i]))
                interpolatemod!(buf, decinv[i], decfmr, interp, T(sz[i]))
                decinv[i] .= buf
            end
            for i = 1:N
                decfmr[i] .= mod.(res[i] .- res3[i] .+ sz_2[i], sz[i]) .- sz_2[i]
                res2[i] .= res3[i]
            end
            note = log(2, norm(decfmr))
            if note > 1
                for i = 1:N
                    decfmr[i] .= decfmr[i] / note
                end
            end
            @show ind, note, borne
            ind += 1
        end
        for i = 1:N
            decinv[i] .= mod.(decinv[i] .- res[i] .+ sz_2[i], sz[i]) .- sz_2[i]
            if z == 1
                decfmr[i] .= decinv[i]
                decinv[i] .= res[i]
                res2[i] .= res4[i]
            end
        end

    end
    return decinv
end

function autointerp!(
    to::Array{OpTuple{N,T},N},
    from::Array{OpTuple{N,T},N},
    nb::Int,
    interp_t::AbstractVector{I};
    mpid::Union{MPIData,Missing} = missing,
    t_split::Union{Tuple,Missing} = missing,
    cachethreads::Union{Vector{CachePrecal{T}},Missing} = missing,
) where {N,T,I<:AbstractInterpolation}
    if nb < 1
        to .= from
    end
    fmr = copy(from)
    for i = 1:nb
        interpolate!(
            to,
            from,
            fmr,
            interp_t;
            mpid = mpid,
            t_split = t_split,
            cachethreads = cachethreads,
        )
        #        @show "autointerp!", i, nb, norm(fmr-to)
        if i != nb
            fmr .= to
        end
    end
end
function interpbufc!(
    t_buf::Vector{Array{OpTuple{N,T},N}},
    bufdec::Array{OpTuple{N,T},N},
    interp_t::AbstractVector{I},
    nb::Int = length(t_buf);
    mpid::Union{MPIData,Missing} = missing,
    t_split::Union{Tuple,Missing} = missing,
    cachethreads::Union{Vector{CachePrecal{T}},Missing} = missing,
) where {N,T,I<:AbstractInterpolation{T}}

    for i = 0:nb-1
        buf = t_buf[end-i]
        interpolate!(
            buf,
            copy(buf),
            bufdec,
            interp_t,
            mpid = mpid,
            t_split = t_split,
            cachethreads = cachethreads,
        )
    end
end


text:
function nk_glfw3_char_callback(win::GLFW.Window, codepoint::Char)
    if glfw.text_len < NK_GLFW_TEXT_MAX
        glfw.text_len += 1
        glfw.text[glfw.text_len] = codepoint
    end
    return nothing
end

function nk_gflw3_scroll_callback(win::GLFW.Window, xoff::Cdouble, yoff::Cdouble)
    new_x = glfw.scroll.x + Cfloat(xoff)
    new_y = glfw.scroll.y + Cfloat(yoff)
    glfw.scroll = nk_vec2(new_x, new_y)
    return nothing
end

function nk_glfw3_mouse_button_callback(window::GLFW.Window, button::GLFW.MouseButton, action::GLFW.Action, mods::Cint)
    button != GLFW.MOUSE_BUTTON_LEFT && return nothing
    x, y = GLFW.GetCursorPos(window)
    if action == GLFW.PRESS
        dt = time_ns()/1e9 - glfw.last_button_click
        if NK_GLFW_DOUBLE_CLICK_LO < dt < NK_GLFW_DOUBLE_CLICK_HI
            glfw.is_double_click_down = true
            glfw.double_click_pos = nk_vec2(Cfloat(x), Cfloat(y))
        end
        glfw.last_button_click = time_ns()/1e9
    else
        glfw.is_double_click_down = false
    end
    return nothing
end

function nk_glfw3_clipboard_paste(usr::nk_handle, edit::Ptr{nk_text_edit})::Cvoid
    text = clipboard()
    !isempty(text) && nk_textedit_paste(edit, text, nk_strlen(text))
    return nothing
end

function nk_glfw3_clipboard_copy(usr::nk_handle, text::Ptr{UInt8}, len::Cint)::Cvoid
    len != 0 && clipboard(unsafe_string(text, len))
    return nothing
end

ext: jl


text:
export
ConvexCombMeasure, ConvexCombScore,
S3Measure, S3Score, DS3Measure, DS3Score,
WECMeasure, WECScore, DWECMeasure, DWECScore,
NetalMeasure, NetalScore, NullMeasure, NullScore,
LCCSMeasure, LCCSScore, CIQMeasure, CIQScore,
GOCMeasure, GOCScore,
GhostMeasure, GhostScore

"""
Null measure
Always returns 0.0 score
"""
immutable NullMeasure <: NetalignMeasure
    G
    function NullMeasure(G::SparseMatrixCSC...)
        if length(G) < 2 error("Need 2 of more networks") end
        for i = 2:length(G)
            if size(G[i-1],1) > size(G[i],1) error("Network size order") end
        end
        new(G)
    end
end

immutable NullScore <: NetalignScore
    score :: Float64
end

measure(m::NullMeasure,fs::Tuple{Vector{<:Integer}}) = NullScore(0.0)
measure(m::NullMeasure,fs::Vector{<:Integer}...) = NullScore(0.0)
dim(m::NullMeasure,d::Int) = size(m.G[end],1)
dim(m::NullMeasure) = size(m.G[end],length(m.G))

"""
    ConvexCombMeasure(S::A,T::B,alpha::Float64) where {A<:NetalignMeasure,B<:NetalignMeasure}

Convex combination of two `NetalignMeasure`s. `0 <= alpha <= 1`.
"""
immutable ConvexCombMeasure{A<:NetalignMeasure,B<:NetalignMeasure} <: NetalignMeasure
    S :: A
    T :: B
    alpha :: Float64 # alpha S + (1-alpha) T
    function ConvexCombMeasure{A,B}(S::A, T::B,alpha::Float64) where
        {A<:NetalignMeasure,B<:NetalignMeasure}
        if dim(S)!=dim(T) error("Sizes of measures A and B do not match: $(dim(S)) vs $(dim(T))") end
        new(S,T,alpha)
    end
end
ConvexCombMeasure(S::A,T::B,alpha::Float64) where {A<:NetalignMeasure,B<:NetalignMeasure} =
    ConvexCombMeasure{A,B}(S,T,alpha)

immutable ConvexCombScore{A<:NetalignScore,B<:NetalignScore} <: NetalignScore
    s :: A
    t :: B
    score :: Float64
end
ConvexCombScore(s::A,t::B) where {A<:NetalignScore,B<:NetalignScore} =
    ConvexCombScore{A,B}(s,t)

function measure(m::ConvexCombMeasure,f::Vector{Int})
    a = measure(m.S,f)
    b = measure(m.T,f)
    s = m.alpha * score(a) + (1.0-m.alpha) * score(b)
    ConvexCombScore(a,b,s)
end

dim(m::ConvexCombMeasure,d::Int) = dim(m.S,d)
dim(m::ConvexCombMeasure) = dim(m,2)

"""
    S3Measure(G1::SparseMatrixCSC,G2::SparseMatrixCSC)

S3 edge conservation measure.
"""
immutable S3Measure <: NetalignMeasure
    G1 :: SparseMatrixCSC{Int,Int}
    G2 :: SparseMatrixCSC{Int,Int}
    function S3Measure(G1::SparseMatrixCSC,G2::SparseMatrixCSC)
        if size(G1,1) > size(G2,1) error("Network size order") end
        new(G1,G2)
    end
end

immutable S3Score <: NetalignScore
    Nc :: Int
    Nn :: Int
    score :: Float64 # Nc/(Nc+Nn)
end

function measure(meas::S3Measure,f::Vector{Int})
    # Assumes that all edge weights are 1
    h = view(f,1:size(meas.G1,1))
    w = nonzeros(meas.G1 + meas.G2[h,h])
    Nc = count(x->x==2,w)
    Nn = length(w) - Nc
    Nc,cr = divrem(Nc,2)
    Nn,nr = divrem(Nn,2) # need this & above so it works well with sim. ann. code
    if cr!=0 || nr!=0; error("G1 and G2 need to be symmetric"); end
    score = Nc/(Nc + Nn)
    S3Score(Nc,Nn,score)
end

dim(m::S3Measure,d::Int) = d==1 ? size(m.G1,1) : d==2 ? size(m.G2,1) : error("d")
dim(m::S3Measure) = dim(m,2)

"""
    DS3Measure(G1::SparseMatrixCSC,G2::SparseMatrixCSC)
DS3 event conservation measure.
"""
immutable DS3Measure <: NetalignMeasure
    G1 :: SparseMatrixCSC{Events,Int}
    G2 :: SparseMatrixCSC{Events,Int}
    function DS3Measure(G1::SparseMatrixCSC,G2::SparseMatrixCSC)
        if size(G1,1) > size(G2,1) error("Bad args") end
        new(G1,G2)
    end
end

immutable DS3Score <: NetalignScore
    Tc :: Float64
    Tn :: Float64
    score :: Float64 # Tc/(Tc+Tn)
end

function measure(m::DS3Measure,f::Vector{Int})
    G1 = m.G1
    G2 = m.G2
    h = view(f,1:size(G1,1))
    # Important: assumes edges (of G1 and G2) are symmetric,
    I1,J1,V1 = findnz(G1)
    I2,J2,V2 = findnz(G2[h,h])
    S = sparse(vcat(I1,I2),vcat(J1,J2),vcat(V1,V2),
               size(G1,1),size(G1,1),
               mergesorted)
    Tc,Tn = 0.0,0.0
    for tis in nonzeros(S)
        Tcp,Tnp = cet_ncet(tis)
        Tc += Tcp
        Tn += Tnp
    end
    Tc /= 2
    Tn /= 2 # this & above required for sim. ann. code
    score = Tc / (Tc + Tn)
    DS3Score(Tc,Tn,score)
end

dim(m::DS3Measure,d::Int) =
    d==1 ? size(m.G1,1) : d==2 ? size(m.G2,1) : error("d")
dim(m::DS3Measure) = dim(m,2)

"""
    WECMeasure(G1::SparseMatrixCSC,G2::SparseMatrixCSC,S::AbstractMatrix)

Weighted edge conservation from WAVE.
"""
immutable WECMeasure{T<:AbstractMatrix} <: NetalignMeasure
    G1 :: SparseMatrixCSC{Int,Int}
    G2 :: SparseMatrixCSC{Int,Int}
    S :: T # matrix type
    function WECMeasure{T}(G1::SparseMatrixCSC,
                           G2::SparseMatrixCSC,S::T) where {T<:AbstractMatrix}
        if size(G1,1) > size(G2,1) || size(S,1)!=size(G1,1) ||
             size(S,2)!=size(G2,1) error("Network/matrix dims") end
        new(G1,G2,S)
    end
end
WECMeasure(G1::SparseMatrixCSC,G2::SparseMatrixCSC,S::AbstractMatrix) =
    WECMeasure{typeof(S)}(G1,G2,S)

immutable WECScore <: NetalignScore
    score :: Float64
end

function measure(m::WECMeasure,f::Vector{Int})
    # Assumes that all edge weights are 1
    h = view(f,1:size(m.G1,1))
    I,J,w = findnz(m.G1 + m.G2[h,h])
    wix = w.==2  # conserved edges
    I = I[wix]
    KI = sub2ind(size(m.S), I, h[I])
    score = sum(m.S[KI]) / min(nnz(m.G1),nnz(m.G2))
    WECScore(score)
end

dim(m::WECMeasure,d::Int) = size(m.S,d)
dim(m::WECMeasure) = dim(m,2)

"""
    DWECMeasure(G1::SparseMatrixCSC{Events},G2::SparseMatrixCSC{Events},S::AbstractMatrix)

Dynamic weighted edge conservation from DynaWAVE.
"""
immutable DWECMeasure{T} <: NetalignMeasure
    # IMP: symmetric "adjacency" matrices
    G1 :: SparseMatrixCSC{Events,Int}
    G2 :: SparseMatrixCSC{Events,Int}
    S :: T # matrix type
    activitysum1 :: Float64
    activitysum2 :: Float64
    function DWECMeasure{T}(G1::SparseMatrixCSC{Events},
                            G2::SparseMatrixCSC{Events},S::T) where {T}
        if size(G1,1) > size(G2,1) || size(S,1)!=size(G1,1) ||
             size(S,2)!=size(G2,1) error("Bad args") end
        new(G1, G2, S, networkactivity(G1), networkactivity(G2))
    end
end
DWECMeasure(G1::SparseMatrixCSC{Events}, G2::SparseMatrixCSC{Events},S::AbstractMatrix) =
    DWECMeasure{typeof(S)}(G1,G2,S)

immutable DWECScore <: NetalignScore
    score :: Float64
end

function measure(m::DWECMeasure,f::Vector{Int})
    G1 = m.G1
    G2 = m.G2
    h = view(f,1:size(G1,1))
    I1,J1,V1 = findnz(G1)
    I2,J2,V2 = findnz(G2[h,h])
    I,J,W = findnz(sparse(vcat(I1,I2),vcat(J1,J2),vcat(V1,V2),
                          size(G1,1),size(G1,1), mergesorted))
    w = map(es -> cet_ncet(es)[1], W)
    KI = sub2ind(size(m.S), I, h[I])
    score = dot(w,m.S[KI]) / min(m.activitysum1,m.activitysum2)
    DWECScore(score)
end

dim(m::DWECMeasure,d::Int) = size(m.S,d)
dim(m::DWECMeasure) = size(m.S,2)

"""
    NetalMeasure(G1::SparseMatrixCSC,G2::SparseMatrixCSC)

Measure from NETAL paper.
"""
immutable NetalMeasure <: NetalignMeasure
    G1 :: SparseMatrixCSC{Int,Int}
    G2 :: SparseMatrixCSC{Int,Int}
    dep1 :: Vector{Float64}
    dep2 :: Vector{Float64}
    expci1 :: Vector{Float64}
    expci2 :: Vector{Float64}
    maxdeg1 :: Int
    maxdeg2 :: Int
    function NetalMeasure(G1::SparseMatrixCSC,G2::SparseMatrixCSC)
        if size(G1,1) > size(G2,1) error("Bad args") end
        deg1 = degree(G1)
        deg2 = degree(G2)
        dep1 = dependency(G1,deg1)
        dep2 = dependency(G2,deg2)
        expci1 = expci(G1,dep1)
        expci2 = expci(G2,dep2)
        new(G1,G2, dep1,dep2, expci1,expci2, maximum(deg1),maximum(deg2))
    end
end

immutable NetalScore <: NetalignScore
    score :: Float64
end

function measure(meas::NetalMeasure,f::AbstractVector{Int})
    # Assumes that all edge weights are 1
    h = view(f,1:size(meas.G1,1))
    I,J,w = findnz(meas.G1 + meas.G2[h,h])
    score = count(x -> x==2, w) / min(nnz(meas.G1),nnz(meas.G2))
    NetalScore(score)
end

dim(m::NetalMeasure,d::Int) =
    d==1 ? size(m.G1,1) : d==2 ? size(m.G2,1) : error("d")
dim(m::NetalMeasure) = dim(m,2)

"""
   LCCSMeasure(G::SparseMatrixCSC...)

Largest common conserved sub-graph, from multiMAGNA++ paper. Works for >2 networks.
"""
immutable LCCSMeasure <: NetalignMeasure
    G
    function LCCSMeasure(G::SparseMatrixCSC...)
        if length(G) < 2 error("Need 2 of more networks") end
        new(G)
    end
end

immutable LCCSScore <: NetalignScore
    edgefrac :: Float64
    nodefrac :: Float64
    edgecov :: Int
    nodecov :: Int
    score :: Float64
end

function compositegraph(G::Tuple,
                        fs::Tuple{Vector{<:Integer}})
    k = length(G)
    ncs = size(G[end],1)
    I,J,_ = findnz(G[end])
    M = Vector(k)
    M[k] = sparse(I,J,1,ncs,ncs)
    for l = 1:k-1
        f = fs[l]
        Ic,Jc,Vc = findnz(G[l])
        fi = f[Ic]
        fj = f[Jc]
        M[l] = sparse(fi,fj,1,ncs,ncs)
        append!(I, fi)
        append!(J, fj)
    end
    M,sparse(I,J,1,ncs,ncs)
end

# Only works for 1-1 MNA
# fs[i] is permutation from G[i] to G[end]
function measure(meas::LCCSMeasure,fs::Tuple{Vector{<:Integer}})
    for i = 2:length(meas.G)
        if size(meas.G[i-1],1) > size(meas.G[i],1) error("Network size order") end
    end
    k = length(meas.G)
    length(fs) == k-1 || error("number of permutations")
    ncs = size(meas.G[end],1)
    M,R = compositegraph(meas.G, fs) # wrt G[end]

    I,J,V = findnz(R)
    ix = V.==k
    I = I[ix]
    J = J[ix]
    A = sparse(I,J,1,size(R,1),size(R,2))

    if nnz(A)==0
        efrac,nfrac,nc,ec = (0.0,0.0,0.0,0,0)
    end

    Acc,p = largest_component(A)
    nc = sum(p)
    ec = nnz(Acc)

    # find the maximum possible number of edges for each graph
    # want the inverse map from cluster to each graph
    # see magna 2014 bioinf paper
    emin = minimum([nnz(M[l][p,p]) for l in 1:k])
    nmin = minimum(map(G->size(G,1), meas.G))
    if nmin!=0 && emin!=0
        efrac = ec/emin
        nfrac = nc/nmin
    else
        nfrac = 0.0
        efrac = 0.0
    end
    return LCCSScore(efrac,nfrac,ec,nc,sqrt(efrac*nfrac))
end

immutable CIQMeasure <: NetalignMeasure
    G
    function CIQMeasure(G::SparseMatrixCSC...)
        if length(G) < 2 error("Need 2 of more networks") end
        new(G)
    end
end
immutable CIQScore <: NetalignScore
    ciqn :: Float64
    ciqd :: Float64
    score :: Float64
end
# Only works for 1-1 MNA
# fs[i] is permutation from G[i] to G[end]
function measure(meas::CIQMeasure,fs::Tuple{Vector{<:Integer}})
    for i = 2:length(meas.G)
        if size(meas.G[i-1],1) > size(meas.G[i],1) error("Network size order") end
    end

    k = length(meas.G)
    length(fs) == k-1 || error("number of permutations")
    ncs = size(meas.G[end],1)
    M,R = compositegraph(meas.G, fs)

    if nnz(R)==0
        ciqn,ciqd = (0.0,0.0)
    end

    scl = ones(Int,ncs,k)
    for l = 1:k-1
        f = fs[l]
        for i = 1:size(meas.G[l],1)
            scl[f[i],l] += 1
        end
    end

    u,v,r = findnz(triu(R,1))
    scl = scl .> 0
    # >0 if there is a node in each cluster that both come from the same network
    s = vec(sum(scl[u,:] .& scl[v,:],2))

    cs = r./s
    cs[(s.==0) .| (r.<=1)] = 0.0

    ciqd = sum(r)
    ciqn = dot(r,cs)

    ciq = ciqd!=0 ? ciqn/ciqd : 0.0
    CIQScore(ciqn,ciqd,ciq)
end

function gofilter(known)
    counts = Dict{String,Int}()
    for node in known, goterm in node
        counts[goterm] = get(counts,goterm,0) + 1
    end
    goodterms = keys(filter((k,v) -> v>1, counts))
    map(goterms -> intersect(goterms,goodterms), known)
end

immutable GOCMeasure <: NetalignMeasure
    known1 :: Vector # known[i] contains annotated go for node i
    known2 :: Vector
    k :: Int
    function GOCMeasure(known1::Vector,known2::Vector,k::Int=1)
        if length(known1) > length(known2) error("Network size order") end
        new(gofilter(known1),gofilter(known2),k)
    end
end
immutable GOCScore <: NetalignScore
    num :: Int
    den :: Int
    score :: Float64
end
function measure(meas::GOCMeasure,f::Vector{<:Integer})
    s = 0
    c = 0
    for (x,y) in zip(meas.known1,meas.known2[view(f,1:length(meas.known1))])
        ixy = length(intersect(x,y))
        mxy = min(length(x),length(y))
        if mxy >= meas.k
            s += Int(ixy > 0) #/mxy
            c += 1
        end
    end
    GOCScore(s, c, c==0 ? 0.0 : s/c)
end

include("ghost.jl")

"""
GHOST signature and node similarity
"""
immutable GhostMeasure <: NetalignMeasure
    G1 :: SparseMatrixCSC{Int,Int}
    G2 :: SparseMatrixCSC{Int,Int}
    evs1::Vector{Vector{Vector{Float64}}}
    evs2::Vector{Vector{Vector{Float64}}}
    k :: Int
    function GhostMeasure(G1::SparseMatrixCSC,G2::SparseMatrixCSC,k=4;verbose=true)
        if size(G1,1) > size(G2,1) error("Network size order") end
        verbose && println("Calculating GHOST signatures")
        m = size(G1,1)
        n = size(G2,1)
        evs1 = ghost_signature(G1,k)
        evs2 = ghost_signature(G2,k)
        new(G1,G2,evs1,evs2,k)
    end
end

immutable GhostScore <: NetalignScore
    score :: Float64
end

function measure(meas::GhostMeasure,f::Vector{Int})
    s = 0.0
    n = length(f)
    for i = 1:n
        s += score(meas, i, f[i])
    end
    GhostScore(s/n)
end

function score(meas::GhostMeasure, u::Int, v::Int)
    ev1 = meas.evs1[u]
    ev2 = meas.evs2[v]
    d = 0.0
    for r = 1:meas.k
        d += sqrt(js_divergence(ev1[r], ev2[r]))
    end
    1.0 - d/meas.k
end

dim(m::GhostMeasure,d::Int) = d==1 ? size(m.G1,1) : d==2 ? size(m.G2,1) : error("d")
dim(m::GhostMeasure) = dim(m,2)


text:
struct IHTOptions
 epsTol::Float64
 maxIter::Int64
 checkEvery::Int64
end

IHTOptions(;epsTol::Float64=1e-2, maxIter::Int64=500, checkEvery::Int64=10) =
  IHTOptions(epsTol, maxIter, checkEvery)

  function covselLatentIHT_grad_S!(gS, Σ, S, L)
    T = S + L
    T = inv(T)
    gS .= Σ - T
    gS
  end

  function covselLatentIHT_grad_U!(gU, gS, Σ, S, L, U)
    covselLatentIHT_grad_S!(gS, Σ, S, L)
    scale!(gS, 2.)
    A_mul_B!(gU, gS, U)
  end

  function covselLatentIHT!(Σ::StridedMatrix, ηS, ηU, s, r;
    epsTol=1e-2, maxIter=100, checkEvery=5, callback=nothing)

    p = size(Σ, 1)

    gS = zeros(p, p)
    gU = zeros(p, r)

    S = zeros(p, p)
    L = zeros(p, p)

    # init
    tmp = full(inv(Symmetric(Σ)))
    HardThreshold!(S, tmp, s, false)

    tmp .= tmp - S
    U, d, V = svd(tmp)
    U = U[:,1:r] .* sqrt.(d[1:r])'
    L = U * U'

    fvals = []
    # gradient descent
    #
    fv = covSel_objective(Σ, S, L)
    push!(fvals, fv)
    for iter=1:maxIter
      # update S
      covselLatentIHT_grad_S!(gS, Σ, S, L)
      @. S -= ηS * gS
      HardThreshold!(S, S, s, false)
      max_gS = maximum( abs.(gS) )

      # update U and L
      covselLatentIHT_grad_U!(gU, gS, Σ, S, L, U)
      @. U -= ηU * gU
      A_mul_Bt!(L, U, U)
      max_gU = maximum( abs.(gU) )
      if callback != nothing
        callback(Σ, S, L)
      end

      done = max(max_gS, max_gU) < epsTol
      # check for convergence
      if mod(iter, checkEvery) == 0
        fv_new = covSel_objective(Σ, S, L)
        push!(fvals, fv_new)
        done = abs(fv_new - fv) <= epsTol
        fv = fv_new
      end
      done && break
    end

    (S, L, U, fvals, iter)
  end


text:
using ExprTools
using Base.Test

macro splitcombine(fundef) # should be a no-op
    dict = splitdef(fundef)
    esc(ExprTools.combinedef(dict))
end

let
    # Ideally we'd compare the result against :(function f(x)::Int 10 end),
    # but it fails because of :line and :block differences
    @test ExprTools.longdef1(:(f(x)::Int = 10)).head == :function
    # That test would belong in MacroTools.jl
    #@test shortdef(:(function f(x)::Int 10 end)).head != :function
    @test map(splitarg, (:(f(a=2, x::Int=nothing, y, args...))).args[2:end]) ==
        [(:a, :Any, false, 2), (:x, :Int, false, :nothing),
         (:y, :Any, false, nothing), (:args, :Any, true, nothing)]
    @test splitarg(:(::Int)) == (nothing, :Int, false, nothing)

    @splitcombine foo(x) = x+2
    @test foo(10) == 12
    @splitcombine add(a, b=2; c=3, d=4)::Float64 = a+b+c+d
    @test add(1; d=10) === 16.0
    @splitcombine fparam{T}(a::T) = T
    @test fparam([]) == Vector{Any}
    immutable Orange end
    @splitcombine (::Orange)(x) = x+2
    @test Orange()(10) == 12
    if VERSION >= v"0.6.0"
        include_string("""
        @splitcombine fwhere(a::T) where T = T
        @test fwhere(10) == Int
        """)
    end
end


text:
using CreditMetrics
using Statistics

# ========================================================
# Using CreditMetrics for simulating a single counterparty
# ========================================================

# Probabilities of being in each rating class after a period
# Note that ProbabilityVector is very unforgiving when checking if probabilities sum to one
pv = ProbabilityVector([0.05, 0.15, 0.25, 0.25, 1.0 - (0.05 + 0.15 + 0.25 + 0.25)])

# Value of the exposure to the counterparty in each rating class after a period
vals = [-1.0, -0.2, -0.1, 0.0, 0.05]

# ---------------------------------------------
# Random outcomes without macroeconomic effects
# ---------------------------------------------

# Compute a vector of values in a number of scenarios
my_base_outcomes = outcome_value(pv, vals, rand(1_000_000))

# Analytically predicted mean value from probabilities and values
meanvalue(pv, vals)

# Observed mean value is the same
mean(my_base_outcomes)

# Observed standard deviation
std(my_base_outcomes)

# --------------------------------------------------------------
# Single counterparty, but in a -2 sigma macroeconomic scenarios
# --------------------------------------------------------------

# Correlation to macroeconomic environment
rho = 0.3

# Convert to z (Gauss) based thresholds
link = NormalThresholdLink{Float64}()
thr = RealLineThresholds(UnitIntervalThresholds(pv), link)

# Deform thresholds according to macroeconomic environment
dthr = deform(thr, 2.0 * rho, sqrt(1 - rho^2))

# At this time we can either convert back to unit interval thresholds or stay in transformed coordinates
# Normally, it will be more efficient to just procede in transformed coordinates
my_stressed_outcomes = outcome_value(dthr, vals, randn(1_000_000))

# Observed mean value differs from original mean due to macroeconomic stress
mean(my_stressed_outcomes)

# Observed standard deviation
std(my_stressed_outcomes)

# Simulating over macroeconomic scenarios, but taking mean wrt the idiosyncratic effect

# For each scenario, compute value after deformation
my_macro_outcomes = deformation_value(thr, vals, link, -rho * randn(1_000_000), sqrt(1 - rho^2))

# Observed mean corresponds to unstressed mean
mean(my_macro_outcomes)

# Observed unexpected loss at 0.05%
quantile(my_macro_outcomes, 0.0005) - mean(my_macro_outcomes)


text:
#######################
# AbstractInstruction #
#######################

abstract type AbstractInstruction end

const InstructionTape = Vector{AbstractInstruction}

function record!(tp::InstructionTape, ::Type{InstructionType}, args...) where InstructionType
    tp !== NULL_TAPE && push!(tp, InstructionType(args...))
    return nothing
end

function Base.:(==)(a::AbstractInstruction, b::AbstractInstruction)
    return (a.func == b.func &&
            a.input == b.input &&
            a.output == b.output &&
            a.cache == b.cache)
end

# Ensure that the external state is "captured" so that external
# reference-breaking (e.g. destructive assignment) doesn't break
# internal instruction state. By default, `capture` is a no-op.
@inline capture(state) = state
@inline capture(state::Tuple) = map(capture, state)

# ScalarInstruction #
#-------------------#

struct ScalarInstruction{F,I,O,C} <: AbstractInstruction
    func::F
    input::I
    output::O
    cache::C
    # disable default outer constructor
    function ScalarInstruction{F,I,O,C}(func, input, output, cache) where {F,I,O,C}
        return new{F,I,O,C}(func, input, output, cache)
    end
end

@inline function _ScalarInstruction(func::F, input::I, output::O, cache::C) where {F,I,O,C}
    return ScalarInstruction{F,I,O,C}(func, input, output, cache)
end

function ScalarInstruction(func, input, output, cache = nothing)
    return _ScalarInstruction(func, capture(input), capture(output), cache)
end

# SpecialInstruction #
#--------------------#

struct SpecialInstruction{F,I,O,C} <: AbstractInstruction
    func::F
    input::I
    output::O
    cache::C
    # disable default outer constructor
    function SpecialInstruction{F,I,O,C}(func, input, output, cache) where {F,I,O,C}
        return new{F,I,O,C}(func, input, output, cache)
    end
end

@inline function _SpecialInstruction(func::F, input::I, output::O, cache::C) where {F,I,O,C}
    return SpecialInstruction{F,I,O,C}(func, input, output, cache)
end

function SpecialInstruction(func, input, output, cache = nothing)
    return _SpecialInstruction(func, capture(input), capture(output), cache)
end

##########
# passes #
##########

function forward_pass!(tape::InstructionTape)
    for instruction in tape
        forward_exec!(instruction)
    end
    return nothing
end

@noinline forward_exec!(instruction::ScalarInstruction) = scalar_forward_exec!(instruction)
@noinline forward_exec!(instruction::SpecialInstruction) = special_forward_exec!(instruction)

function reverse_pass!(tape::InstructionTape)
    for i in length(tape):-1:1
        reverse_exec!(tape[i])
    end
    return nothing
end

@noinline reverse_exec!(instruction::ScalarInstruction) = scalar_reverse_exec!(instruction)
@noinline reverse_exec!(instruction::SpecialInstruction) = special_reverse_exec!(instruction)

###################
# Pretty Printing #
###################

# extra spaces here accomodates padding in show(::IO, ::AbstractInstruction)
compactrepr(x::Tuple) = "("*join(map(compactrepr, x), ",\n           ")*")"
compactrepr(x::AbstractArray) = length(x) < 5 ? match(r"\[.*?\]", repr(x)).match : summary(x)
compactrepr(x) = repr(x)

function Base.show(io::IO, instruction::AbstractInstruction, pad = "")
    name = isa(instruction, ScalarInstruction) ? "ScalarInstruction" : "SpecialInstruction"
    println(io, pad, "$(name)($(instruction.func)):")
    println(io, pad, "  input:  ", compactrepr(instruction.input))
    println(io, pad, "  output: ", compactrepr(instruction.output))
    print(io,   pad, "  cache:  ", compactrepr(instruction.cache))
end

function Base.show(io::IO, tp::InstructionTape)
    println("$(length(tp))-element InstructionTape:")
    i = 1
    for instruction in tp
        print(io, "$i => ")
        show(io, instruction)
        println(io)
        i += 1
    end
end


text:
module GPT

using  KnetLayers, WordTokenizers, BytePairEncoding, JSON, NPZ
import KnetLayers: IndexedDict, AbstractTransformer, Layer, arrtype, dir, Datasets
using KnetLayers.Datasets
import KnetLayers.Datasets: Container

struct Gpt <: AbstractTransformer
    ts::Chain
    drop::Dropout
    loss::CrossEntropyLoss
end

function Gpt(size::Int, head::Int, ps::Int, layer::Int;
             max_len::Int = 512, trainable = true, act = ELU(), pdrop = 0.1)
    rem(size, head) != 0 && error("size not divisible by head")
    Gpt(size, head, div(size, head), ps, layer; max_len=max_len, trainable=trainable, act=act, pdrop=pdrop)
end

function Gpt(size::Int, head::Int, hs::Int, ps::Int, layer::Int;
             max_len::Int = 512, trainable = true, act = ELU(), pdrop = 0.1)
    Gpt(Chain([Transformer(size, head, hs, ps; future=false, act=act, pdrop=pdrop) for i = 1:layer]...),
        Dropout(pdrop), CrossEntropyLoss())
end

function (gpt::Gpt)(x, mask=nothing; all::Bool=false)
    e = gpt.drop(x)
    t = gpt.ts(e)
    if mask != nothing
        mask = convert(arrtype,mask)
        t = t .* reshape(mask,1,size(mask)...)
    end
    return t
end
"""
    lmloss(embed, onehot, encoding, mask)
compute the language modeling loss for Gpt, onehot is the onehot array of the origin
input sentence. encoding the output of Gpt, mask is the mask between input sentences.
"""
lmloss(embed::Embed, o, t, mask) = lmloss(embed, o, t, mask)
function lmloss(embed::Embed, et, t, mask)
    N = ndims(t)
    if N == 3
        t = t[:, 1:end-1, :]
        s = size(t)
        #FIXME: transpose
        out = Embed(embed.weight')
        return KnetLayers.nllmask(out(t), et[2:end,:].*mask[2:end,:]) #(vocab, seq_len*batch)
    elseif N == 2
        t = t[:, 1:end-1]
        s = size(t)
        out = Embed(embed.weight')
        return KnetLayers.nllmask(out(t), et[2:end,:] .*mask[2:end]) #(vocab, seq_len*batch)
    end
end

function Base.show(io::IO, gpt::Gpt)
    hs = div(size(gpt.ts[1].mh.iqproj)[1], gpt.ts[1].mh.head)
    h, ps = size(gpt.ts[1].pw.layers[end])
    print(io, "Gpt(")
    print(io, "layers=$(length(gpt.ts.layers)), ")
    print(io, "head=$(gpt.ts[1].mh.head), ")
    print(io, "head_size=$(hs), ")
    print(io, "pwffn_size=$(ps), ")
    print(io, "size=$(h))")
end

function load_gpt_pretrain_params()
    shapes  = JSON.parsefile(dir("examples/gpt/pretrain/params_shapes.json"))
    offsets = accumulate(+, prod.(shapes))
    shapes  = map(s -> length(s) > 1 ? (s[end], s[end-1]) : s, shapes)
    params  = cat([npzread(joinpath(dirname(@__FILE__), "pretrain/params_$(i).npy")) for i = 0:9]..., dims=1)
    params  = [collect(reshape(selectdim(params, 1, a+1:b), s...)) for (a, b, s) in zip([0;offsets[1:end-1]], offsets, shapes)]
    return params
end

gpt_tokenizer(x) = toktok_tokenize(text_standardize(x))

function load_gpt_pretrain(n::Int=12;
                           startsym="_start_",
                           delisym="_delimiter_",
                           clfsym="_classify_",
                           unksym="<unk>")
    n > 12 && error("pretrain maximum layer: 12")
    #set_tokenizer((x)->nltk_word_tokenize(text_standardize(x)))
    emp = JSON.parsefile(dir("examples/gpt/pretrain/encoder_bpe_40000.json"))
    vocab = map(first, sort!(collect(emp), by=(x)->x.second))
    push!(vocab, startsym)
    push!(vocab, delisym)
    push!(vocab, clfsym)
    if unksym ∉ vocab
         pushfirst!(vocab,unk)
    end
    bpe = Bpe(joinpath(dir("examples/gpt/pretrain/vocab_40000.bpe")))
    vocab = IndexedDict(vocab)
    embed = Embed(output=768, input=length(vocab))
    pe = PositionEmbedding(768, 512; trainable=true)
    ce = CompositeEmbedding(embed, pe, unksym)
    gpt = Gpt(768, 12, 768*4, 12; act=GeLU(), pdrop=0.1)
    pms = load_gpt_pretrain_params()
    loadparams!(embed, [hcat(pms[2],randn(Float32,768, 3) .* 0.02)])
    loadparams!(pe, [pms[1]])
    for i = 1:n
        mhW = pms[12(i-1) + 3]
        mhb = pms[12(i-1) + 4]
        loadparams!(gpt.ts[i].mh.iqproj,[selectdim(mhW, 1, 1:768),
                                         selectdim(mhb, 1, 1:768)])
        loadparams!(gpt.ts[i].mh.ikproj,[selectdim(mhW, 1, 768+1:2*768),
                                         selectdim(mhb, 1, 768+1:2*768)])
        loadparams!(gpt.ts[i].mh.ivproj,[selectdim(mhW, 1, 2*768+1:3*768),
                                         selectdim(mhb, 1, 2*768+1:3*768)])
        loadparams!(gpt.ts[i].mh.oproj,[pms[12(i-1) + 5],
                                        pms[12(i-1) + 6]])
        loadparams!(gpt.ts[i].mhn,[pms[12(i-1) + 7],
                                   pms[12(i-1) + 8]])
        loadparams!(gpt.ts[i].pw.layers[1],[pms[12(i-1) + 9],
                                      pms[12(i-1) + 10]])
        loadparams!(gpt.ts[i].pw.layers[2],[pms[12(i-1) + 11],
                                       pms[12(i-1) + 12]])
        loadparams!(gpt.ts[i].pwn,[pms[12(i-1) + 13],
                                   pms[12(i-1) + 14]])
    end
    #gpt, embed, bpe, vocab
    TransformerModel(ce, gpt), bpe, vocab, gpt_tokenizer
end

struct CompositeEmbedding <: Layer
    tok::Embed
    pos::PositionEmbedding
    unksym::String
end

(m::CompositeEmbedding)(x) =  m.pos(m.tok(x))

function loadparams!(m::Layer, ws)
    for (w,wo) in zip(params(m),ws)
        copyto!(w, convert(arrtype,wo))
    end
end

"""
    getmask(ls::Container{<:Container})
get the mask for batched data.
"""
function getmask(ls::Container{<:Container})
    lens = map(length, ls)
    m = zeros(Int, maximum(lens), length(lens))

    for (i, l) ∈ enumerate(ls)
        m[1:length(l),i] .= 1
        #selectdim(selectdim(m, 2, i), 1, 1:length(l)) .= 1
    end
    m
end

function text_standardize(text)
    text = lowercase(text)
    text = replace(text, r"([a-z])(,|\.)"=>s"\1 \2")
    text = replace(text, "—"=>"-")
    text = replace(text, "–"=>"-")
    text = replace(text, "―"=>"-")
    text = replace(text, "…"=>"...")
    text = replace(text, "´"=>"'")
    text = replace(text, r"""(-+|~+|!+|"+|;+|\?+|\++|,+|\)+|\(+|\\+|\/+|\*+|\[+|\]+|}+|{+|\|+|_+)"""=>s" \1 ")
    text = replace(text, r"\s*\n\s*"=>" \n ")
    text = replace(text, r"[^\S\n]+"=>" ")
    strip(text)
end

end


text:
println("=========== testing multiobjectiveGraphs package =============")

#=
import multiobjectiveGraphs
using multiobjectiveGraphs
=#

greet() ### encore une fois ###

#=
import Pkg
Pkg.add("Test")
import Test
=#

using Test

#=

@testset "priority lists" begin
	PL =  priorityList{UInt8}()
	add!(PL,UInt8(5)) ; pretty(PL)
	add!(PL,UInt8(12)) ; pretty(PL)
	add!(PL,UInt8(7)) ; pretty(PL)
	add!(PL,UInt8(2)) ; pretty(PL)
	add!(PL,UInt8(20)) ; pretty(PL)
	add!(PL,UInt8(9)) ; pretty(PL)
	@test size(PL) == 6
	@test first(PL) == 2
	@test removefirst!(PL) == 2
	@test removefirst!(PL) == 5
	@test size(PL) == 4
	for x::UInt8 in [(3*i) % 100 for i in 20:80] add!(PL,x) end
	@test removefirst!(PL) == 1
	@test first(PL) == 2
	pretty(PL) ; println("taille : ",length(PL))
end

=#

@testset "objectives" begin
	anObj = weightMinSum{Int32}()
	@test typeof(objTypes(anObj)[1]) == typeof(		weightMinSum{Int32}() 		) 	# --> 1
  	@test typeof(objTypes(anObj)[2]) == typeof(		weightMinProduct{Int32}() 	)  	# --> 2
  	@test typeof(objTypes(anObj)[3]) == typeof(		weightMinMax{Int32}() 		)  	# --> 3
  	@test typeof(objTypes(anObj)[4]) == typeof(		weightMaxSum{Int32}() 		)	# --> 4
  	@test typeof(objTypes(anObj)[5]) == typeof(		weightMaxProduct{Int32}() 	)	# --> 5
  	@test typeof(objTypes(anObj)[6]) == typeof(		weightMaxMin{Int32}() 		) 	# --> 6
end

@testset "multiObjectives" begin
    multObj = multiobj(2,1,2, typeofvalues=Float32)

    @test string(typeof(multObj)) == "multiobj"
    @test string(multObj) ==  "multiobj(5, genericWeightCategory[weightMinSum{Float32}(0.0f0), weightMinSum{Float32}(0.0f0), weightMinProduct{Float32}(1.0f0), weightMinMax{Float32}(0.0f0), weightMinMax{Float32}(0.0f0)])"
    @test descr(multiobj())=="[ min∑(0.0)]"
    @test descr(multiobj(multiobj()))=="[ min∑(0.0)]"
    @test descr(multiobj(5))=="[ min∑(0.0),min∑(0.0),min∑(0.0),min∑(0.0),min∑(0.0)]"
    @test descr(multiobj(3,typeofobjectives=typeof(weightMinSum{UInt8}(5))))=="[ min∑(0),min∑(0),min∑(0)]"
end

#=

@testset "MOgraph" begin
    GG = MOgraph{ Int64, Float64 }(5, 3) # graph with 5 vertices, 3 objectives
    ar = edge{  Int64, Float64 }( 2, multiobj{Float64}(3), 5 ) # edge from 2 to 5 with default tri-objective
	addEdge!( GG, ar )
	moType = [weightMinSum{Float64}() , weightMinProduct{Float64}(), weightMinMax{Float64}() ] # default specific tri-objective
	for i in 2:nv(GG)
			mo = mimic( moType, [ i*i, i, 99-i ] ) # same of multiobjective but with these 3 values
			ar=mimic(ar,i-1, mo, i)						# same type of edge but from i-1 to i with this objective
			addEdge!( GG, ar )
			mo = mimic( moType, [ 99-i, i*i, i ] )
			ar = mimic(ar, i, mo, nv(GG)+1-i )
			addEdge!( GG, ar )
	end

	summarize(GG)
	println(typeof(GG))
	t=plot(GG)
	GG=WheelGraph(10)
	println(typeof(GG))
    t=plot(GG)
    save(SVG("TikzDessin.svg"), t)

    @test nv(GG) == 5
    @test ne(GG) == 9
end

=#

max_stars_repo_name: DRobbes/multiobjectiveGraphs

max_issues_repo_name: DRobbes/multiobjectiveGraphs

max_forks_repo_name: DRobbes/multiobjectiveGraphs


text:
module KeywordDispatch

export @kwdispatch, @kwmethod

function ntsort(nt::NamedTuple{N}) where {N}
    if @generated
        names = tuple(sort(collect(N))...)
        types = Tuple{Any[ fieldtype(nt, n) for n in names ]...}
        vals = Any[ :(getfield(nt, $(QuoteNode(n)))) for n in names ]
        :( NamedTuple{$names,$types}(($(vals...),)) )
    else
        names = tuple(sort(collect(N))...)
        types = Tuple{Any[ fieldtype(typeof(nt), n) for n in names ]...}
        vals = map(n->getfield(nt, n), names)
        NamedTuple{names,types}(vals)
    end
end

# given an argument, output a form appropriate for a method definition in pass-through methods
# currently:
#  - replaces underscores with new symbols
#    g(_) => g(newsym)
#  - replaces tuples with a single symbol
#    g((a,b)) => g(newsym)
#  - inserts symbols for 1-arg ::,
#    g(::T)   => g(newsym::T)
argmeth(x::Symbol) = x == :_ ? gensym() : x
function argmeth(x::Expr)
    x.head == :tuple && return gensym()
    x.head == :(::)  && return :($(length(x.args)==1 ? gensym() : argmeth(x.args[1]))::$(x.args[end]))
    x.head == :...   && return Expr(:..., argmeth.(x.args)...)
    return x
end

# given an argument, output a form appropriate for a method call in pass-through methods
# currently:
#  - removes the :: part
#    g(x::T)   => g(x)
argsym(x::Symbol) = x
function argsym(x::Expr)
    x.head == :(::) && return length(x.args) > 1 ? x.args[1] : :_
    x.head == :...  && return Expr(:..., argsym.(x.args)...)
    return x
end

# given an argument, output an argument type
# currently:
#  - g(x::T)  => T
#  - g(x)     => Any
argtype(x::Symbol) = Any
argtype(x::Expr) = x.head == :(::) ? x.args[end] : error("unexpected expression $x")

function unwrap_where(expr)
    stack = Any[]
    while expr isa Expr && expr.head == :where
        push!(stack, expr.args[2])
        expr = expr.args[1]
    end
    expr, stack
end

function wrap_where(expr, stack)
    for w in Iterators.reverse(stack)
        expr = Expr(:where, expr, esc(w))
    end
    expr
end

struct KeywordMethodError <: Exception
    f
    args
    kwargs
    KeywordMethodError(@nospecialize(f), @nospecialize(args), @nospecialize(kwargs)) = new(f, args, kwargs)
end

function Base.showerror(io::IO, err::KeywordMethodError)
    print(io, "KeywordMethodError: no keyword method matching ")
    if err.f isa Function || err.f isa Type
        print(io, err.f)
    else
        print(io, "(::", typeof(err.f), ")")
    end
    print(io, "(")
    for (i, arg) in enumerate(err.args)
        print(io, "::", typeof(arg))
        if i != length(err.args)
            print(io, ", ")
        end
    end
    print(io, "; ")
    for (i, (kw,kwarg)) in enumerate(pairs(err.kwargs))
        print(io, kw, "::", typeof(kwarg))
        if i != length(err.kwargs)
            print(io, ", ")
        end
    end
    print(io, ")")
end

function kwcall(nt::NamedTuple, f, args...)
    throw(KeywordMethodError(f, args, nt))
end

"""
    @kwdispatch sig [methods]

Designate a function signature `sig` that should dispatch on keyword arguments. A
function can also be provided, in which case _all_ calls to that function are will
dispatch to keyword methods.

It is possible to specify keyword aliases by specifying `from => to` pairs in the keyword position.

The optional `methods` argument allows a block of keyword methods specified as anonymous
functions. To define additional keyword methods, use the [`@kwmethod`](@ref) macro.

# Examples

```julia
@kwdispatch f(_)
@kwdispatch f(_::Real)

@kwdispacth f # equivalent to @kwdispatch f(_...)

@kwdispatch f(x) begin
    (a) -> x+a
    (b) -> x-b
end
# equivalent to
#  @kwdispatch f(x)
#  @kwmethod f(x;a) = x+a
#  @kwmethod f(x;b) = x-b

@kwdispatch f(; alpha=>α) # specifies alpha as an alias to α
```
"""
macro kwdispatch(fexpr,methods=nothing)
    fcall, wherestack = unwrap_where(fexpr)

    # handle: `fun`, `Mod.fun`, `(a::B)`
    if fcall isa Symbol || fcall isa Expr && (fcall.head in (:., :(::), :curly))
        fcall = :($fcall(_...))
    end

    @assert fcall isa Expr && fcall.head == :call

    rename_expr = :(kw)

    f = fcall.args[1]
    posargs = fcall.args[2:end]
    if length(posargs) >= 1 && posargs[1] isa Expr && posargs[1].head == :parameters
        parameters = popfirst!(posargs)
        for p in parameters.args
            if p isa Expr && p.head == :call && p.args[1] == :(=>)
                rename_expr = :(kw == $(QuoteNode(p.args[2])) ? $(QuoteNode(p.args[3])) : $rename_expr)
            else
                error("Only renames (from => to) are allowed in keyword position of `@kwdispatch`")
            end
        end
    end
    f = argmeth(f)

    if f isa Expr && f.head == :(::)
        ftype = esc(f.args[end])
    else
        ftype = :(typeof($(esc(f))))
    end

    posargs_method = argmeth.(posargs)
    ff = esc(argsym(f))

    quote
        $(wrap_where(:($(esc(f))($(esc.(posargs_method)...); kwargs...)), wherestack)) = begin
            N = map(kw -> $rename_expr, propertynames(kwargs.data))
            nt = NamedTuple{N}(Tuple(kwargs.data))
            KeywordDispatch.kwcall(ntsort(nt), $ff, $(esc.(argsym.(posargs_method))...))
        end
        $(generate_kwmethods(methods, f, posargs, wherestack))
    end
end

generate_kwmethods(other, f, posargs, wherestack) = other
function generate_kwmethods(expr::Expr, f, posargs, wherestack)
    if expr.head == :block
        for (i, ex) in enumerate(expr.args)
            expr.args[i] = generate_kwmethods(ex, f, posargs, wherestack)
        end
        return expr
    elseif expr.head in (:->, :function)
        if expr.args[1] isa Symbol
            return kwmethod_expr(f, posargs, [expr.args[1]], wherestack, esc(expr.args[2]))
        elseif expr.args[1] isa Expr && expr.args[1].head == :tuple
            return kwmethod_expr(f, posargs, expr.args[1].args, wherestack, esc(expr.args[2]))
        end
    end
    error("Invalid keyword definition $expr.")
end

"""
    @kwmethod expr

Define a keyword method for dispatch. `expr` should be a standard function definition
(either block or inline) with a keyword argument block (which may be empty).

The positional signature should first be designated by the [`@kwdispatch`](@ref) macro.

# Examples

```julia
@kwdispatch f() # designate the no positional argument form

@kwmethod f(;) = nothing # no keyword arguments
@kwmethod f(;a) = a
@kwmethod f(;b) = 2*b
@kwmethod f(;a,b) = a+2*b
```
"""
macro kwmethod(fexpr)
    @assert fexpr isa Expr && fexpr.head in (:function, :(=))
    body = esc(fexpr.args[2])

    fcall, wherestack = unwrap_where(fexpr.args[1])

    @assert fcall isa Expr && fcall.head == :call

    f = fcall.args[1]

    length(fcall.args) >= 2 && fcall.args[2] isa Expr && fcall.args[2].head == :parameters ||
        error("@kwmethod requires functions specify a keyword block.\nUse @kwmethod `f(args...;)` to specify no keywords.")

    kwargs = fcall.args[2].args
    posargs = fcall.args[3:end]
    kwmethod_expr(f, posargs, kwargs, wherestack, body)
end

function kwmethod_expr(f, posargs, kwargs, wherestack, body)
    sort!(kwargs, by=argsym)

    kwsyms = argsym.(kwargs)
    kwtypes = argtype.(kwargs)

    if f isa Expr && f.head == :(::)
        F = esc(f)
    else
        F = :(::($(esc(f)) isa Type ? Type{$(esc(f))} : typeof($(esc(f)))))
    end

    quote
        $(wrap_where(:(KeywordDispatch.kwcall(($(esc.(kwsyms)...),)::NamedTuple{($(QuoteNode.(kwsyms)...),),T},
                                              $F,
                                              $(esc.(posargs)...)) where {T<:Tuple{$(esc.(kwtypes)...)}}), wherestack)) =
                                                  $body
    end
end

end # module


text:
# Many of these are redundant. They should be replaced in the rest of the code by `isa`

is_SSJSym(x) = isa(x,SSJSym)
is_SJSym(x) = isa(x,SJSym)

is_Mxpr(x) = isa(x,Mxpr)

## I think we should prefer isa(Mxpr{:Sym}) here
is_Mxpr(mx::Mxpr{T},s::Symbol) where {T} = T == s
is_Mxpr(x,s::Symbol) = false

## Ruleq or ruleq ? need lc somewhere to signify it is lower level
ruleq(x::Mxpr{:Rule}) = true
ruleq(x) = false

# we could generate some of these with a loop and eval
Base.isempty(mx::Mxpr) = isempty(margs(mx))

"""
    mxpr_head_freeq(mx::Mxpr)

return false if any element at level 1 of mx is an Mxpr with head `head`.
"""
mxpr_head_freeq(mx::Mxpr, head) = ! any(t -> isa(t,Mxpr{head}), margs(mx))

## probably should not use these. isa signals that we are interested in Julia type rather
## than a paraphyletic class of numbers.
is_Number(x) = isa(x,Number)
is_Real(x) = isa(x,Real)
is_Complex(x) = isa(x,Complex)
is_Float(x) = isa(x,AbstractFloat)

## TODO julia has some function like this. We should use, if possible:
## isinteger
## isapprox
## isassigned, for arrays
## isequal
## isempty (instead of length == 0)
is_imaginary_integer(z::Complex{T}) where {T<:Integer} = real(z) == 0
is_imaginary_integer(x) = false

atomq(x) = ! isa(x,Mxpr)

is_Indeterminate(x::Symbol) = x == Indeterminate
is_Indeterminate(x) = false

is_Infintity(x) = false
is_Infintity(mx::Mxpr{:DirectedInfinity}) = ! isempty(mx) && mx[1] == 1 ? true : false

is_ComplexInfinity(x) = false
is_ComplexInfinity(mx::Mxpr{:DirectedInfinity}) =  isempty(mx) ? true : false

# use hasConstant
# is_Constant(x::Symbol) = get_attribute(x, :Constant)
# is_Constant(x) = false

# use hasProtected instead
#is_protected(sj::SJSym) = get_attribute(sj,:Protected)

is_blankxxx(x) = isa(x,BlankXXX)

## question: shold we do isxxx, is_xxx, or xxxq ?
## Julia does isxxx (or is_xxx in a few cases). Mma uses XxxQ uniformly (almost)
## Answer:
## Base.isinteger(1.0) -> true
## So we should do xxxq consistently because we want a clear distinction from Julia
## Julia's isinteger makes perfect sense within Julia. But, the name would collide
## isinteger for us.

## These are equivalent to a single isa, but they can be used directly as predicate functions to
## supply as arguments
listq(x::Mxpr{:List}) = true
listq(x) = false
integerq(x::Integer) = true
integerq(x) = false

listofpredq(x::Mxpr{:List},pred) = isempty(x) ? false : all(pred,x)
listofpredq(x,pred) = false
listoflistsq(x::Mxpr{:List}) = listofpredq(x,listq)

function listoflistsofsamelengthq(x::Mxpr{:List})
    listoflistsq(x) || return false
    n = symlength(x[1])
    all( a -> symlength(a) == n, x)
end
listoflistsofsamelengthq(x) = false

listofintegersq(x::Mxpr{:List}) = listofpredq(x,integerq)  ## maybe we should not define this

listoflistsofpredq(x::Mxpr{:List},pred) = all(a -> listofpredq(a,pred), x)
listoflistsofpredq(x) = false # maybe not necessary, but maybe the compiler needs help ?

"""
    listoflistsofsamelengthpredq(x,pred)

true if x is List of Lists of same length with predq giving true on all elements.
"""
function listoflistsofsamelengthpredq(x::Mxpr{:List},pred)
    listoflistsofsamelengthq(x) && listoflistsofpredq(x,pred)
end
listoflistsofsamelengthpredq(x) = false

####  Symata Predicates

### ConstantQ

@mkapprule ConstantQ :nargs => 1

@sjdoc ConstantQ """
    ConstantQ(x)

return `True` if `x` is a numerical constant.
"""
do_ConstantQ(mx::Mxpr{:ConstantQ}, x) = isConstant(x)

### AtomQ

@mkapprule AtomQ  :nargs => 1
@sjdoc AtomQ """
    AtomQ(expr)

return true if `expr` has no parts accessible with `Part`.

However, currently, Julia `Array`s can be accessed with `Part`, and return `True` under `AtomQ`.
"""
@doap AtomQ(x) = atomq(x)

### OddQ

@mkapprule OddQ  :nargs => 1
@doap OddQ(x) = isa(x, Integer) && ! iseven(x)
@sjdoc OddQ """
    OddQ(expr)

return `True` if `expr` is an odd integer.
"""

### EvenQ

@mkapprule EvenQ  :nargs => 1
@doap EvenQ(x) = isa(x, Integer) && iseven(x)
@sjdoc EvenQ """
    EvenQ(expr)

return `True` if `expr` is an even integer.
"""

@sjseealso_group(AtomQ, EvenQ, OddQ)

### DirtyQ

@mkapprule DirtyQ :nargs => 1
@sjdoc DirtyQ """
    DirtyQ(m)

return `True `if the timestamp of any symbol that `m` depends on
is more recent than the timestamp of `m`. This is for diagnostics.
`DirtyQ` evaluates `m` once.
"""
@doap DirtyQ(expr) = expr |> symval |> checkdirtysyms

### NumericQ

@mkapprule NumericQ  :nargs => 1

@sjdoc NumericQ """
    NumericQ(expr)

return true if `N(expr)` would return a number.
"""
do_NumericQ(mx::Mxpr{:NumericQ}, x) = is_Numeric(x)
is_Numeric(x) = false
is_Numeric(x::Number) = true
is_Numeric(x::Symbol) = isConstant(x)
is_Numeric(x::Mxpr) = get_attribute(x,:NumericFunction)  && all( t -> is_Numeric(t), margs(x))

### NumberQ

@mkapprule NumberQ :nargs => 1

@sjdoc NumberQ """
    NumberQ(x)

return true if `x` is an explicit number. i.e. it is a subtype of Julia type `Number`.
"""
@doap NumberQ(x) = isa(x,Number)

### MachineNumberQ

@mkapprule MachineNumberQ :nargs => 1
# Should we check for smaller floats ?

@sjdoc MachineNumberQ """
    MachineNumberQ(x)

return `True` if `x` is a machine-precision floating point number.
"""
@doap MachineNumberQ(x::Float64) = true
@doap MachineNumberQ(x::Complex{Float64}) = true
@doap MachineNumberQ(x) = false

### InexactNumberQ

@mkapprule InexactNumberQ :nargs => 1

@sjdoc InexactNumberQ """
    InexactNumberQ(x)

return `True` if `x` is an inexact number. Inexact numbers are floating
point numbers and floating point complex numbers.
"""
@doap InexactNumberQ(x::AbstractFloat) = true
@doap InexactNumberQ(x::Complex{<:AbstractFloat}) = true
@doap InexactNumberQ(x) = false

### IntegerQ

@mkapprule IntegerQ :nargs => 1
@sjdoc IntegerQ """
    IntegerQ(x)

return true if `x` is an `Integer`.
"""
@doap IntegerQ(x) = isa(x,Integer)

### ListQ

@mkapprule ListQ :nargs => 1
@sjdoc ListQ """
    ListQ(x)

return `True` if the `Head` of `x` is `List`.
"""
@doap ListQ(x::Mxpr{:List}) = true
@doap ListQ(x) = false

### Positive

@mkapprule Positive :nargs => 1
@sjdoc Positive """
    Positive(x)

return `True` if `x` is a positive number. return `Null` if `x` is a `String`.
Otherwise, return the input expression.
"""
@doap Positive(x::Real) = x > 0
@doap Positive(x::Number) = false
@doap Positive(x::String) = nothing

### PermuationQ

@mkapprule PermutationQ :nargs => 1
@sjdoc PermutationQ """
    PermutationQ(list)

return `True` if and only if `list` is a permuation of the integers from `1` through `Length(list)`.
"""
do_PermutationQ(mx::Mxpr{:PermutationQ}, lst::Mxpr{:List}) =  all(t -> typeof(t) <: Union{Integer,AbstractFloat}, margs(lst)) && isperm(margs(lst))

### VectorQ

@mkapprule VectorQ :nargs => 1:2

@sjdoc VectorQ """
    VectorQ(x)

return `True` if `x` is a `List`, none of whose elements is a `List`.

    VectorQ(x,test)

return `True` if `x` is a vector, all of whose elements satisfy the predicate `test`.
For example, `VectorQ(expr, IntegerQ)` returns `True` if `x` is a vector of `Integers`.
"""
vectorq(x::Mxpr{:List}) = all( t -> (! isa(t,Mxpr{:List})), margs(x))
vectorq(x::Mxpr{:List}, test) = isa(test,Function) ? all(t -> test(t), margs(x)) : all(t ->doeval(mxpr(test,t)), margs(x))
@doap VectorQ(x) = false
@doap VectorQ(x::Mxpr{:List}) = vectorq(x)
@doap VectorQ(x::Mxpr{:List}, test) = vectorq(x,test)

### MatrixQ

@mkapprule MatrixQ :nargs => 1:2  :nodefault => true
@sjdoc MatrixQ """
    MatrixQ(m)

return `True` if `m` has the form of a matrix.
`m` has the form of a matrix if it is a `List`, each
of whose elements is a `List` of the same length,
none of whose elements is a `List`.

    MatrixQ(m,test)

return `True` if `test` applied to each element of the matrix `m`.
"""
@doap MatrixQ(args...) = matrixq(args...)

function matrixq(x::Mxpr{:List})
    a = margs(x)
    length(a) < 1 && return false
    len = symlength(a[1])
    all(t -> ( symlength(t) == len && vectorq(t) ), a)
end

function matrixq(x::Mxpr{:List},test)
    a = margs(x)
    length(a) < 1 && return false
    len = symlength(a[1])
    all(t -> ( symlength(t) == len && vectorq(t,test) ), a)
end
matrixq(x) = false
matrixq(x,test) = false

### TrueQ

@mkapprule TrueQ :nargs => 1
@doap TrueQ(x) = trueq(x)

trueq(x::Bool) = x
trueq(x) = false

### Boole

@mkapprule Boole
@doap Boole(x::Bool) = x ? 1 : 0

### BooleanQ

@mkapprule BooleanQ :nodefault => true
@doap BooleanQ(x::Bool) = true
@doap BooleanQ(x...) = false

### FreeQ

## FIXME: Need Heads => True

# don't know how to put this in the macro string
# const levelstring = """
# - `n`       levels `0` through `n`.
# - `[n]`     level `n` only.
# - `[n1,n2]` levels `n1` through `n2`
# """

@mkapprule FreeQ :nargs => 1:3
@doap FreeQ(expr, pattern) = freeq(LevelSpecAll(),expr,pattern)
@doap FreeQ(expr, pattern, inlevelspec)  = freeq(make_level_specification(expr, inlevelspec), expr, pattern)

@curry_second FreeQ
@sjdoc FreeQ """
    FreeQ(expr, pattern)

return `True` if no subexpression of `expr` matches `pattern`.

    FreeQ(expr, pattern, levelspec)

return `True` if `expr` does not match `pattern` on levels specified by `levelspec`.

- `n`       levels `1` through `n`.
- `[n]`     level `n` only.
- `[n1,n2]` levels `n1` through `n2`

    FreeQ(pattern)

operator form.
"""

mutable struct FreeQData
    pattern
    gotmatch::Bool
end

freeq(expr, pat) = freeq(LevelSpecAll(), expr, pat)

function freeq(levelspec::LevelSpec, expr, pat)
    data = FreeQData(pat,false)
    action = LevelAction(data,
                         function (data, expr)
                             (gotmatch,cap) = match_and_capt(expr,patterntoBlank(data.pattern))
                             if gotmatch
                               data.gotmatch = true
                               action.levelbreak = true
                             end
                         end)
    if has_level_zero(levelspec)  # Do level zero separately
        (gotmatch,cap) = match_and_capt(expr,patterntoBlank(data.pattern))
        gotmatch && return false
    end
    traverse_levels!(action,levelspec,expr)
    data.gotmatch && return false
    return true
end

### Element

## We may want to do this with a dictionary. As it is the Julia compiler
## may have a chance to reason.

@sjdoc Element """
    Element(x,dom),  x ∈ dom

asserts that `x` is an element of the domain `dom`.
"""
@mkapprule Element :nargs => 2

@doap function Element(x::Integer,sym::Symbol)
    if sym == :Integers || sym == :Reals || sym == :Rationals || sym == :Complexes ||
        sym == :Algebraics
        return true
    end
    sym == :Primes && return isprime(x)
    sym == :Booleans && return false
    mx
end

@doap function Element(x::AbstractFloat,sym::Symbol)
    sym == :Integers && return round(x) == x ? mx : false
    sym == :Rationals && return mx
    sym == :Complexes && return true
    sym == :Booleans && return false
    sym == :Reals && return true
    sym == :Algebraics && return mx
    if sym == :Primes
        if round(x) == x
            return isprime(convert(Integer,x))
        else
            return false
        end
    end
    mx
end

@doap function Element(x::Complex{T},sym::Symbol) where T
    sym == :Integers && return false
    sym == :Rationals && return false
    sym == :Reals && return false
    sym == :Booleans && return false
    if sym == :Algebraics
        T <: Integer && return true
        return mx
    end
    sym == :Complexes && return true
    mx
end

@doap function Element(x::Bool,sym::Symbol)
    sym == :Booleans && return true
    sym == :Integers && return false
    sym == :Rationals && return false
    sym == :Reals && return false
    sym == :Complexes && return false
    sym == :Algebraics && return false
    mx
end

@doap function Element(x::Rational,sym::Symbol)
    sym == :Booleans && return false
    sym == :Integers && return false
    sym == :Rationals && return true
    sym == :Reals && return true
    sym == :Complexes && return true
    sym == :Algebraics && return true
    mx
end

@doap function Element(x::Symbol,sym::Symbol)
    if x == :Pi || x == :E
        (sym == :Integers || sym == :Rationals || sym == :Booleans || sym == :Algebraics || sym == :Primes) && return false
        (sym == :Reals || sym == :Complexes)  && return true
    end
    if x == :GoldenRatio
        sym == :Algebraics && return true
        sym == :Reals && return true
        sym == :Rationals && return false
        sym == :Integers && return false
        sym == :Complexes && return true
    end
    if x == :Catalan
        sym == :Reals && return true
        sym == :Integers && return false
        sym == :Complexes && return true
    end
    mx
end

## This will return unevaluated some inputs for which the answer is known, eg, powers involving Pi
@doap Element(x::PowerT, sym::Symbol) = _element_power(mx,x,margs(x)...,sym)
function _element_power(mx,x,b::Number,expt::Number,sym)
    (sym == :Integers || sym == :Rationals || sym == :Booleans || sym == :Primes) && return false
    (sym == :Reals || sym == :Complexes || sym == :Algebraics)  && return true
    x
end
_element_power(mx,args...) = mx

## This will probably give the wrong answer for some inputs
@doap function Element(x::TimesT, sym::Symbol)
    if is_Numeric(x) && freeq(LevelSpecAll(),x, mxpr(:Blank,:Symbol))
        (sym == :Integers || sym == :Rationals || sym == :Booleans || sym == :Primes) && return false
        (sym == :Reals || sym == :Complexes || sym == :Algebraics)  && return true
    end
    mx
end


text:
# Methods for the computation of Gram matrices and continuous projections in general

# By convention Gram functionality is only implemented for dictionaries that are
# associated with a measure.
hasmeasure(Φ::Dictionary) = false

# Determine a measure to use when two dictionaries are given
defaultmeasure(Φ1::Dictionary, Φ2::Dictionary) =
    _defaultmeasure(Φ1, Φ2, measure(Φ1), measure(Φ2))

function _defaultmeasure(Φ1, Φ2, m1, m2)
    if iscompatible(m1, m2)
        m1
    else
        if iscompatible(support(Φ1),support(Φ2))
            lebesguemeasure(support(Φ1))
        else
            error("Please specify which measure to use for the combination of $(Φ1) and $(Φ2).")
        end
    end
end

# Shortcut: Dictionaries of the same type have just one measure
defaultmeasure(Φ1::D, Φ2::D) where {D <: Dictionary} = measure(Φ1)

innerproduct(Φ1::Dictionary, i, Φ2::Dictionary, j; options...) =
    innerproduct(Φ1, i, Φ2, j, defaultmeasure(Φ1, Φ2); options...)

# Convert linear indices to native indices, then call innerproduct_native
innerproduct(Φ1::Dictionary, i::Int, Φ2::Dictionary, j::Int, measure; options...) =
    innerproduct_native(Φ1, native_index(Φ1, i), Φ2, native_index(Φ2, j), measure; options...)
innerproduct(Φ1::Dictionary, i, Φ2::Dictionary, j::Int, measure; options...) =
    innerproduct_native(Φ1, i, Φ2, native_index(Φ2, j), measure; options...)
innerproduct(Φ1::Dictionary, i::Int, Φ2::Dictionary, j, measure; options...) =
    innerproduct_native(Φ1, native_index(Φ1, i), Φ2, j, measure; options...)
innerproduct(Φ1::Dictionary, i, Φ2::Dictionary, j, measure; options...) =
    innerproduct_native(Φ1, i, Φ2, j, measure; options...)

# - innerproduct_native: if not specialized, called innerproduct1
innerproduct_native(Φ1::Dictionary, i, Φ2::Dictionary, j, measure; options...) =
    innerproduct1(Φ1, i,  Φ2, j, measure; options...)
# - innerproduct1: possibility to dispatch on the first dictionary without ambiguity.
#                  If not specialized, we call innerproduct2
innerproduct1(Φ1::Dictionary, i, Φ2, j, measure; options...) =
    innerproduct2(Φ1, i, Φ2, j, measure; options...)
# - innerproduct2: possibility to dispatch on the second dictionary without ambiguity
innerproduct2(Φ1, i, Φ2::Dictionary, j, measure; options...) =
    default_dict_innerproduct(Φ1, i, Φ2, j, measure; options...)

# Make a quadrature strategy using user-supplied atol and rtol if they were given
function quadstrategy(::Type{T}; atol = DomainIntegrals.default_atol(T), rtol = DomainIntegrals.default_rtol(T), maxevals=DomainIntegrals.default_maxevals() , options...) where {T}
	QuadAdaptive{T}(atol, rtol, maxevals)
end

# We make this a separate routine so that it can also be called directly, in
# order to compare to the value reported by a dictionary overriding innerproduct
function default_dict_innerproduct(Φ1::Dictionary, i, Φ2::Dictionary, j, measure::Measure; options...)
	@boundscheck checkbounds(Φ1, i)
	@boundscheck checkbounds(Φ2, j)
    domain1 = support(Φ1, i)
    domain2 = support(Φ2, j)
    domain3 = support(measure)

	# This code is delicate: it is necessary when using basis functions with compact
	# support, because numerical integration over the whole domain might miss the
	# part where the integrand is non-zero. However, the intersection of three domains
	# below might fail.
	try
    	domain = intersectdomain(domain1, domain2, domain3)
		qs = quadstrategy(prectype(Φ1, Φ2); options...)
    	unsafe_default_dict_innerproduct1(Φ1, i, Φ2, j, measure, domain1, domain2, domain3, domain, qs)
	catch e
		# Intersection or integration failed, use safe evaluation
		@warn "Error thrown, trying to recover"
		# qs = quadstrategy(prectype(Φ1, Φ2); options...)
		safe_default_dict_innerproduct(Φ1, i, Φ2, j, measure, qs)
	end
end

# Routine below is safe because it uses eval_element, and not unsafe_eval_element
safe_default_dict_innerproduct(Φ1::Dictionary, i, Φ2::Dictionary, j, measure::Measure, qs) =
	integral(qs, x->conj(eval_element(Φ1, i, x))*eval_element(Φ2, j, x), measure)

function default_dict_innerproduct(Φ1::Dictionary, i, Φ2::Dictionary, j, measure::DiscreteWeight; options...)
	@boundscheck checkbounds(Φ1, i)
	@boundscheck checkbounds(Φ2, j)
    domain1 = support(Φ1, i)
    domain2 = support(Φ2, j)
    unsafe_default_dict_innerproduct2(Φ1, i, Φ2, j, measure, domain1 ∩ domain2)
end

# unsafe for indexing
function unsafe_default_dict_innerproduct1(Φ1::Dictionary, i, Φ2::Dictionary, j, measure::Weight, d1, d2, d3, domain, qs)
    if d1 == d2 == d3
        # -> domains are the same, don't convert the measure
        unsafe_default_dict_innerproduct2(Φ1, i, Φ2, j, measure, qs)
    else
        # -> do compute on the smaller domain and convert the measure
        integral(qs, x->conj(unsafe_eval_element(Φ1, i, x))*unsafe_eval_element(Φ2, j, x)*unsafe_weightfun(measure,x), domain)
    end
end

unsafe_default_dict_innerproduct1(Φ1::Dictionary, i, Φ2::Dictionary, j, measure::Measure, d1, d2, d3, domain::IntersectDomain, qs) =
    # -> disregard the intersection domain, but use the safe eval instead to guarantee correctness
    integral(qs, x->conj(eval_element(Φ1, i, x))*eval_element(Φ2, j, x), measure)
unsafe_default_dict_innerproduct1(Φ1::Dictionary, i, Φ2::Dictionary, j, measure::Weight, d1, d2, d3, domain::IntersectDomain, qs) =
    # -> disregard the intersection domain, but use the safe eval instead to guarantee correctness
    integral(qs, x->conj(eval_element(Φ1, i, x))*eval_element(Φ2, j, x), measure)

# unsafe for indexing and for support of integral
unsafe_default_dict_innerproduct2(Φ1::Dictionary, i, Φ2::Dictionary, j, measure, qs) =
    integral(qs, x->conj(unsafe_eval_element(Φ1, i, x))*unsafe_eval_element(Φ2, j, x), measure)

function unsafe_default_dict_innerproduct2(Φ1::Dictionary, i, Φ2::Dictionary, j, measure::DiscreteWeight, domain)
	# Note: using "unsafe_eval" here is still slightly dangerous: we rely on "integral"
	# not evaluating the integrand outside the given domain for discrete weights
    integral(x->conj(unsafe_eval_element(Φ1, i, x))*unsafe_eval_element(Φ2, j, x), domain, measure)
end

gramelement(dict::Dictionary, i, j, μ = measure(dict); options...) =
    innerproduct(dict, i, dict, j, μ; options...)

# Call this routine in order to evaluate the Gram matrix entry numerically
default_gramelement(dict::Dictionary, i, j, μ = measure(dict); options...) =
    default_dict_innerproduct(dict, i, dict, j, μ; options...)

function grammatrix(Φ::Dictionary, μ = measure(Φ), T = codomaintype(Φ); options...)
    G = zeros(T, length(Φ), length(Φ))
    grammatrix!(G, Φ, μ; options...)
end

function grammatrix!(G, Φ::Dictionary, μ = measure(Φ); options...)
    n = length(Φ)
    for i in 1:n
        for j in 1:i-1
            G[i,j] = gramelement(Φ, i, j, μ; options...)
            G[j,i] = conj(G[i,j])
        end
        G[i,i] = gramelement(Φ, i, i, μ; options...)
    end
    G
end

gram(Φ::Dictionary, args...; options...) = gram(operatoreltype(Φ), Φ, args...; options...)

gram(::Type{T}, Φ::Dictionary; options...) where {T} = gram(T, Φ, measure(Φ); options...)

gram(::Type{T}, Φ::Dictionary, μ::Measure; options...) where {T} =
	gram1(T, Φ, μ; options...)

gram1(T, Φ::Dictionary, μ; options...) =
	gram2(T, Φ, μ; options...)

gram2(T, Φ, μ::Measure; options...) =
	default_gram(T, Φ, μ; options...)

gram2(T, Φ, μ::DiscreteWeight; options...) =
    gram(T, Φ, μ, points(μ), weights(μ); options...)

gram(::Type{T}, Φ, μ::DiscreteWeight, grid::AbstractGrid, weights; options...) where {T} =
    default_mixedgram_discretemeasure(T, Φ, Φ, μ, grid, weights; options...)

default_gram(Φ::Dictionary, args...; options...) =
	default_gram(operatoreltype(Φ), Φ, args...; options...)

function default_gram(::Type{T}, Φ::Dictionary, μ::Measure = measure(Φ); warnslow = true, options...) where {T}
    # warnslow && @debug "Slow computation of Gram matrix entrywise of $Φ with measure $μ)."
    A = grammatrix(Φ, μ, T; options...)
    R = ArrayOperator(A, Φ, Φ)
end

default_diagonal_gram(Φ::Dictionary, μ::Measure = measure(Φ); options...) =
	default_diagonal_gram(operatoreltype(Φ), μ; options...)

function default_diagonal_gram(::Type{T}, Φ::Dictionary, μ::Measure; options...) where {T}
    @assert isorthogonal(Φ, μ)
	n = length(Φ)
	diag = zeros(T, n)
	for i in 1:n
		diag[i] = innerproduct(Φ, i, Φ, i, μ; options...)
	end
	DiagonalOperator(Φ, diag)
end

function default_mixedgram_discretemeasure(::Type{T}, Φ1::Dictionary, Φ2::Dictionary,
			μ::DiscreteWeight, grid::AbstractGrid, weights; options...) where {T}
    E1 = evaluation(T, Φ1, grid; options...)
    E2 = evaluation(T, Φ2, grid; options...)
    W = DiagonalOperator{T}(dest(E2),dest(E1), weights)
    E1'*W*E2
end

"""
Project the function onto the space spanned by the given dictionary.
"""
project(Φ::Dictionary, f, m = measure(Φ); T = coefficienttype(Φ), options...) =
    project!(zeros(T,Φ), Φ, f, m; options...)

function project!(result, Φ, f, μ; options...)
    for i in eachindex(result)
        result[i] = innerproduct(Φ[i], f, μ; options...)
    end
    result
end

########################
# Mixed gram operators
########################

mixedgram(Φ1::Dictionary, Φ2::Dictionary, args...; options...) =
	mixedgram(operatoreltype(Φ1,Φ2), Φ1, Φ2, args...; options...)

mixedgram(::Type{T}, Φ1::Dictionary, Φ2::Dictionary; options...) where {T} =
    mixedgram(T, Φ1, Φ2, defaultmeasure(Φ1, Φ2); options...)

"""
Compute the mixed Gram matrix corresponding to two dictionaries. The matrix
has elements given by the inner products between the elements of the dictionaries,
relative to the given measure.
"""
mixedgram(::Type{T}, Φ1::Dictionary, Φ2::Dictionary, μ::Measure; options...) where {T} =
    mixedgram1(T, Φ1, Φ2, μ; options...)

# The routine mixedgram1 can be specialized by concrete subtypes of the
# first dictionary, while mixedgram2 can be specialized on the second dictionary.
# mixedgram3 can be specialized on the measure
mixedgram1(T, Φ1::Dictionary, Φ2, μ; options...) =
    mixedgram2(T, Φ1, Φ2, μ; options...)

mixedgram2(T, Φ1, Φ2::Dictionary, μ; options...) =
    mixedgram3(T, Φ1, Φ2, μ; options...)

mixedgram3(T, Φ1, Φ2, μ::Measure; options...) =
    default_mixedgram(T, Φ1, Φ2, μ; options...)

mixedgram3(T, Φ1, Φ2, μ::DiscreteWeight; options...) =
    mixedgram(T, Φ1, Φ2, μ, points(μ), weights(μ); options...)

mixedgram(T, Φ1, Φ2, μ::DiscreteWeight, grid::AbstractGrid, weights; options...) =
    default_mixedgram_discretemeasure(T, Φ1, Φ2, μ, grid, weights; options...)

function mixedgram(::Type{T}, Φ1::D, Φ2::D, μ::Measure; options...) where {T,D<:Dictionary}
    if Φ1 == Φ2
        gram(T, Φ1, μ; options...)
    else
        mixedgram1(T, Φ1, Φ2, μ; options...)
    end
end

default_mixedgram(Φ1::Dictionary, Φ2::Dictionary, args...; options...) =
	default_mixedgram(operatoreltype(Φ1, Φ2), Φ1, Φ2, args...; options...)

function default_mixedgram(::Type{T}, Φ1::Dictionary, Φ2::Dictionary, μ; warnslow = true, options...) where {T}
    # warnslow && @debug "Slow computation of mixed Gram matrix entrywise."
    A = mixedgrammatrix(Φ1, Φ2, μ, T; options...)
    ArrayOperator(A, Φ2, Φ1)
end

function mixedgrammatrix(Φ1::Dictionary, Φ2::Dictionary, μ, T = operatoreltype(Φ1,Φ2); options...)
    G = zeros(T, length(Φ1), length(Φ2))
    mixedgrammatrix!(G, Φ1, Φ2, μ; options...)
end

function mixedgrammatrix!(G, Φ1::Dictionary, Φ2::Dictionary, μ; options...)
    m = length(Φ1)
    n = length(Φ2)
    for i in 1:m
        for j in 1:n
            G[i,j] = innerproduct(Φ1, i, Φ2, j, μ; options...)
        end
    end
    G
end


text:
function random_unit_vector3(rng::MersenneTwister, T::DataType)
    normalize(Vec{3, T}(randn(rng, T, 3)))
end

function supersample_crisp(image::AbstractArray{C, 2}, upsize) where {C<:Color}
    bigimage = zeros(C, (size(image) .* upsize)...)
    height, width = size(image)
    for i in 1:height, j in 1:width
        bigimage[(i-1) * upsize + 1:i * upsize, (j-1) * upsize + 1:j * upsize] .= image[i, j]
    end
    bigimage
end