text:
@noinline not_planar_err() = throw(ArgumentError("not a planar diagram expression"))
@noinline not_planar_err(ex) = throw(ArgumentError("not a planar diagram expression: $ex"))

@nospecialize

macro planar(ex::Expr)
    return esc(planar_parser(ex))
end

function planar_parser(ex::Expr)
    parser = TO.TensorParser()

    pop!(parser.preprocessors) # remove TO.extracttensorobjects
    push!(parser.preprocessors, _conj_to_adjoint)
    treebuilder = parser.contractiontreebuilder
    treesorter = parser.contractiontreesorter
    push!(parser.preprocessors, ex->TO.processcontractions(ex, treebuilder, treesorter))
    push!(parser.preprocessors, ex->_check_planarity(ex))
    push!(parser.preprocessors, _extract_tensormap_objects)
    temporaries = Vector{Symbol}()
    push!(parser.preprocessors, ex->_decompose_planar_contractions(ex, temporaries))

    pop!(parser.postprocessors) # remove TO.addtensoroperations
    push!(parser.postprocessors, ex->_update_temporaries(ex, temporaries))
    push!(parser.postprocessors, ex->_annotate_temporaries(ex, temporaries))
    push!(parser.postprocessors, _add_modules)
    return parser(ex)
end

macro planar2(ex::Expr)
    return esc(planar2_parser(ex))
end

function planar2_parser(ex::Expr)
    parser = TO.TensorParser()

    braidingtensors = Vector{Any}()

    pop!(parser.preprocessors) # remove TO.extracttensorobjects
    push!(parser.preprocessors, _conj_to_adjoint)
    push!(parser.preprocessors, _extract_tensormap_objects2)
    push!(parser.preprocessors, _construct_braidingtensors)
    treebuilder = parser.contractiontreebuilder
    treesorter = parser.contractiontreesorter
    push!(parser.preprocessors, ex->TO.processcontractions(ex, treebuilder, treesorter))
    push!(parser.preprocessors, ex->_check_planarity(ex))
    temporaries = Vector{Symbol}()
    push!(parser.preprocessors, ex->_decompose_planar_contractions(ex, temporaries))

    pop!(parser.postprocessors) # remove TO.addtensoroperations
    push!(parser.postprocessors, ex->_update_temporaries(ex, temporaries))
    push!(parser.postprocessors, ex->_annotate_temporaries(ex, temporaries))
    push!(parser.postprocessors, _add_modules)
    return parser(ex)
end

function _conj_to_adjoint(ex::Expr)
    if ex.head == :call && ex.args[1] == :conj && TO.istensor(ex.args[2])
        obj, leftind, rightind = TO.decomposetensor(ex.args[2])
        return Expr(:typed_vcat, Expr(TO.prime, obj),
                        Expr(:tuple, rightind...), Expr(:tuple, leftind...))
    else
        return Expr(ex.head, [_conj_to_adjoint(a) for a in ex.args]...)
    end
end
_conj_to_adjoint(ex) = ex

function get_possible_planar_indices(ex)
    if !TO.istensorexpr(ex)
        return [[]]
    elseif TO.isgeneraltensor(ex)
        _,leftind,rightind = TO.decomposegeneraltensor(ex)
        ind = planar_unique2(vcat(leftind, reverse(rightind)))
        return length(ind) == length(unique(ind)) ? Any[ind] : Any[]
    elseif ex.head == :call && (ex.args[1] == :+ || ex.args[1] == :-)
        inds = get_possible_planar_indices(ex.args[2])
        keep = fill(true, length(inds))
        for i = 3:length(ex.args)
            inds′ = get_possible_planar_indices(ex.args[i])
            keepᵢ = fill(false, length(inds))
            for (j, ind) in enumerate(inds), ind′ in inds′
                if iscyclicpermutation(ind′, ind)
                    keepᵢ[j] = true
                end
            end
            keep .&= keepᵢ
            any(keep) || break # give up early if keep is all false
        end
        return inds[keep]
    elseif ex.head == :call && ex.args[1] == :*
        @assert length(ex.args) == 3
        inds1 = get_possible_planar_indices(ex.args[2])
        inds2 = get_possible_planar_indices(ex.args[3])
        inds = Any[]
        for ind1 in inds1, ind2 in inds2
            for (oind1, oind2, cind1, cind2) in possible_planar_complements(ind1, ind2)
                push!(inds, vcat(oind1, oind2))
            end
        end
        return inds
    else
        return Any[]
    end
end

# remove double indices (trace indices) from cyclic set
function planar_unique2(allind)
    oind = collect(allind)
    removing = true
    while removing
        removing = false
        i = 1
        while i <= length(oind) && length(oind) > 1
            j = mod1(i+1, length(oind))
            if oind[i] == oind[j]
                deleteat!(oind, i)
                deleteat!(oind, mod1(i, length(oind)))
                removing = true
            else
                i += 1
            end
        end
    end
    return oind
end

# remove intersection (contraction indices) from two cyclic sets
function possible_planar_complements(ind1, ind2)
    # quick return path
    (isempty(ind1) || isempty(ind2)) && return Any[(ind1, ind2, Any[], Any[])]
    # general case:
    j1 = findfirst(in(ind2), ind1)
    if j1 === nothing # disconnected diagrams, can be made planar in various ways
        return Any[(circshift(ind1, i-1), circshift(ind2, j-1), Any[], Any[])
                    for i ∈ eachindex(ind1), j ∈ eachindex(ind2)]
    else # genuine contraction
        N1, N2 = length(ind1), length(ind2)
        j2 = findfirst(==(ind1[j1]), ind2)
        jmax1 = j1
        jmin2 = j2
        while jmax1 < N1 && ind1[jmax1+1] == ind2[mod1(jmin2-1, N2)]
            jmax1 += 1
            jmin2 -= 1
        end
        jmin1 = j1
        jmax2 = j2
        if j1 == 1 && jmax1 < N1
            while ind1[mod1(jmin1-1, N1)] == ind2[mod1(jmax2 + 1, N2)]
                jmin1 -= 1
                jmax2 += 1
            end
        end
        if jmax2 > N2
            jmax2 -= N2
            jmin2 -= N2
        end
        indo1 = jmin1 < 1 ? ind1[(jmax1+1):mod1(jmin1-1, N1)] :
                    vcat(ind1[(jmax1+1):N1], ind1[1:(jmin1-1)])
        cind1 = jmin1 < 1 ? vcat(ind1[mod1(jmin1, N1):N1], ind1[1:jmax1]) : ind1[jmin1:jmax1]
        indo2 = jmin2 < 1 ? ind2[(jmax2+1):mod1(jmin2-1, N2)] :
                    vcat(ind2[(jmax2+1):N2], ind2[1:(jmin2-1)])
        cind2 = reverse(cind1)
        return isempty(intersect(indo1, indo2)) ? Any[(indo1, indo2, cind1, cind2)] : Any[]
    end
end

function _check_planarity(ex::Expr)
    if ex.head == :macrocall && ex.args[1] == Symbol("@notensor")
    elseif TO.isassignment(ex) || TO.isdefinition(ex)
        lhs, rhs = TO.getlhs(ex), TO.getrhs(ex)
        if TO.istensorexpr(rhs)
            if TO.istensorexpr(lhs)
                @assert TO.istensor(lhs)
                indlhs = only(get_possible_planar_indices(lhs)) # should have only one element
            else
                indlhs = Any[]
            end
            indsrhs = get_possible_planar_indices(rhs)
            isempty(indsrhs) && not_planar_err(rhs)
            i = findfirst(ind -> iscyclicpermutation(ind, indlhs), indsrhs)
            i === nothing && not_planar_err(ex)
        end
    else
        foreach(ex.args) do a
            _check_planarity(a)
        end
    end
    return ex
end
_check_planarity(ex) = ex

# decompose contraction trees in order to fix index order of temporaries
# to ensure that planarity is guaranteed
_decompose_planar_contractions(ex, temporaries) = ex
function _decompose_planar_contractions(ex::Expr, temporaries)
    if ex.head == :macrocall && ex.args[1] == Symbol("@notensor")
        return ex
    end
    if TO.isassignment(ex) || TO.isdefinition(ex)
        lhs, rhs = TO.getlhs(ex), TO.getrhs(ex)
        if TO.istensorexpr(rhs)
            pre = Vector{Any}()
            rhs = _extract_contraction_pairs(rhs, lhs, pre, temporaries)
            return Expr(:block, pre..., Expr(ex.head, lhs, rhs))
        else
            return ex
        end
    end
    if TO.istensorexpr(ex)
        pre = Vector{Any}()
        rhs = _extract_contraction_pairs(ex, (Any[], Any[]), pre, temporaries)
        return Expr(:block, pre..., rhs)
    end
    if ex.head == :block
        return Expr(ex.head,
                    [_decompose_planar_contractions(a, temporaries) for a in ex.args]...)
    end
    if ex.head == :for || ex.head == :function
        return Expr(ex.head, ex.args[1],
                        _decompose_planar_contractions(ex.args[2], temporaries))
    end
    return ex
end

# decompose a contraction into elementary binary contractions of tensors without inner traces
# if lhs is an expression, it contains the existing lhs and thus the index order
# if lhs is a tuple, the result is a temporary object and the tuple (lind, rind) gives a suggestion for the preferred index order
function _extract_contraction_pairs(rhs, lhs, pre, temporaries)
    if TO.isscalarexpr(rhs)
        return rhs
    elseif TO.isgeneraltensor(rhs)
        if TO.hastraceindices(rhs) && lhs isa Tuple
            s = gensym()
            newlhs = Expr(:typed_vcat, s, Expr(:tuple, lhs[1]...), Expr(:tuple, lhs[2]...))
            push!(temporaries, s)
            push!(pre, Expr(:(:=), newlhs, rhs))
            return newlhs
        else
            return rhs
        end
    elseif rhs.head == :call && rhs.args[1] == :*
        @assert length(rhs.args) == 3

        if lhs isa Expr
            _, leftind, rightind = TO.decomposetensor(lhs)
        else
            leftind, rightind = lhs
        end
        lhs_ind = vcat(leftind, reverse(rightind))

        # find possible planar order
        rhs_inds = Any[]
        for ind1 in get_possible_planar_indices(rhs.args[2])
            for ind2 in get_possible_planar_indices(rhs.args[3])
                for (oind1, oind2, cind1, cind2) in possible_planar_complements(ind1, ind2)
                    if iscyclicpermutation(vcat(oind1, oind2), lhs_ind)
                        push!(rhs_inds, (ind1, ind2, oind1, oind2, cind1, cind2))
                    end
                    isempty(rhs_inds) || break
                end
                isempty(rhs_inds) || break
            end
            isempty(rhs_inds) || break
        end
        ind1, ind2, oind1, oind2, cind1, cind2 = only(rhs_inds) # inds_rhs should hold exactly one match
        if all(in(leftind), oind2) && all(in(rightind), oind1) # reverse order
            a1 = _extract_contraction_pairs(rhs.args[3], (oind2, reverse(cind2)), pre, temporaries)
            a2 = _extract_contraction_pairs(rhs.args[2], (cind1, reverse(oind1)), pre, temporaries)
            oind1, oind2 = oind2, oind1
            cind1, cind2 = cind2, cind1
        else
            a1 = _extract_contraction_pairs(rhs.args[2], (oind1, reverse(cind1)), pre, temporaries)
            a2 = _extract_contraction_pairs(rhs.args[3], (cind2, reverse(oind2)), pre, temporaries)
        end

        if TO.isscalarexpr(a1) || TO.isscalarexpr(a2)
            rhs = Expr(:call, :*, a1, a2)
            s = gensym()
            newlhs = Expr(:typed_vcat, s, Expr(:tuple, oind1...),
                                        Expr(:tuple, reverse(oind2)...))
            push!(temporaries, s)
            push!(pre, Expr(:(:=), newlhs, rhs))
            return newlhs
        end

        # note that index order in _extract... is only a suggestion, now we have actual index order
        _, l1, r1, = TO.decomposegeneraltensor(a1)
        _, l2, r2, = TO.decomposegeneraltensor(a2)
        if all(in(r1), oind1) && all(in(l2), oind2) # reverse order
            a1, a2 = a2, a1
            ind1, ind2 = ind2, ind1
            oind1, oind2 = oind2, oind1
        end
        if lhs isa Tuple
            rhs = Expr(:call, :*, a1, a2)
            s = gensym()
            newlhs = Expr(:typed_vcat, s, Expr(:tuple, oind1...),
                                        Expr(:tuple, reverse(oind2)...))
            push!(temporaries, s)
            push!(pre, Expr(:(:=), newlhs, rhs))
            return newlhs
        else
            if leftind == oind1 && rightind == reverse(oind2)
                rhs = Expr(:call, :*, a1, a2)
                return rhs
            elseif leftind == oind2 && rightind == reverse(oind1) # probably this can not happen anymore
                rhs = Expr(:call, :*, a2, a1)
                return rhs
            else
                rhs = Expr(:call, :*, a1, a2)
                s = gensym()
                newlhs = Expr(:typed_vcat, s, Expr(:tuple, oind1...),
                                            Expr(:tuple, reverse(oind2)...))
                push!(temporaries, s)
                push!(pre, Expr(:(:=), newlhs, rhs))
                return newlhs
            end
        end
    elseif rhs.head == :call && rhs.args[1] ∈ (:+, :-)
        args = [_extract_contraction_pairs(a, lhs, pre, temporaries) for
                    a in rhs.args[2:end]]
        return Expr(rhs.head, rhs.args[1], args...)
    else
        throw(ArgumentError("unknown tensor expression"))
    end
end

# replacement of TensorOperations functionality:
# adds checks for matching number of domain and codomain indices
# special cases adjoints so that t and t' are considered the same object
function _extract_tensormap_objects(ex)
    inputtensors = _remove_adjoint.(TO.getinputtensorobjects(ex))
    outputtensors = _remove_adjoint.(TO.getoutputtensorobjects(ex))
    newtensors = TO.getnewtensorobjects(ex)
    @assert !any(_is_adjoint, newtensors)
    existingtensors = unique!(vcat(inputtensors, outputtensors))
    alltensors = unique!(vcat(existingtensors, newtensors))
    tensordict = Dict{Any,Any}(a => gensym() for a in alltensors)
    pre = Expr(:block, [Expr(:(=), tensordict[a], a) for a in existingtensors]...)
    pre2 = Expr(:block)
    ex = TO.replacetensorobjects(ex) do obj, leftind, rightind
        _is_adj = _is_adjoint(obj)
        if _is_adj
            leftind, rightind = rightind, leftind
            obj = _remove_adjoint(obj)
        end
        newobj = get(tensordict, obj, obj)
        if (obj in existingtensors)
            nl = length(leftind)
            nr = length(rightind)
            nlsym = gensym()
            nrsym = gensym()
            objstr = string(obj)
            errorstr1 = "incorrect number of input-output indices: ($nl, $nr) instead of "
            errorstr2 = " for $objstr."
            checksize = quote
                $nlsym = numout($newobj)
                $nrsym = numin($newobj)
                ($nlsym == $nl && $nrsym == $nr) ||
                    throw(IndexError($errorstr1 * string(($nlsym, $nrsym)) * $errorstr2))
            end
            push!(pre2.args, checksize)
        end
        return _is_adj ? _restore_adjoint(newobj) : newobj
    end
    post = Expr(:block, [Expr(:(=), a, tensordict[a]) for a in newtensors]...)
    pre = Expr(:macrocall, Symbol("@notensor"), LineNumberNode(@__LINE__, Symbol(@__FILE__)), pre)
    pre2 = Expr(:macrocall, Symbol("@notensor"), LineNumberNode(@__LINE__, Symbol(@__FILE__)), pre2)
    post = Expr(:macrocall, Symbol("@notensor"), LineNumberNode(@__LINE__, Symbol(@__FILE__)), post)
    return Expr(:block, pre, pre2, ex, post)
end
_is_adjoint(ex) = isa(ex, Expr) && ex.head == TO.prime
_remove_adjoint(ex) = _is_adjoint(ex) ? ex.args[1] : ex
_restore_adjoint(ex) = Expr(TO.prime, ex)

function _extract_tensormap_objects2(ex)
    inputtensors = collect(filter(!=(:τ), _remove_adjoint.(TO.getinputtensorobjects(ex))))
    outputtensors = _remove_adjoint.(TO.getoutputtensorobjects(ex))
    newtensors = TO.getnewtensorobjects(ex)
    (any(==(:τ), newtensors) || any(==(:τ), outputtensors)) &&
        throw(ArgumentError("The name τ is reserved for the braiding, and should not be assigned to."))
    @assert !any(_is_adjoint, newtensors)
    existingtensors = unique!(vcat(inputtensors, outputtensors))
    alltensors = unique!(vcat(existingtensors, newtensors))
    tensordict = Dict{Any,Any}(a => gensym() for a in alltensors)
    pre = Expr(:block, [Expr(:(=), tensordict[a], a) for a in existingtensors]...)
    pre2 = Expr(:block)
    ex = TO.replacetensorobjects(ex) do obj, leftind, rightind
        _is_adj = _is_adjoint(obj)
        if _is_adj
            leftind, rightind = rightind, leftind
            obj = _remove_adjoint(obj)
        end
        newobj = get(tensordict, obj, obj)
        if (obj in existingtensors)
            nl = length(leftind)
            nr = length(rightind)
            nlsym = gensym()
            nrsym = gensym()
            objstr = string(obj)
            errorstr1 = "incorrect number of input-output indices: ($nl, $nr) instead of "
            errorstr2 = " for $objstr."
            checksize = quote
                $nlsym, $nrsym = numout($newobj), numin($newobj)
                ($nlsym == $nl && $nrsym == $nr) ||
                    throw(IndexError($errorstr1 * string(($nlsym, $nrsym)) * $errorstr2))
            end
            push!(pre2.args, checksize)
        end
        return _is_adj ? _restore_adjoint(newobj) : newobj
    end
    post = Expr(:block, [Expr(:(=), a, tensordict[a]) for a in newtensors]...)
    pre = Expr(:macrocall, Symbol("@notensor"), LineNumberNode(@__LINE__, Symbol(@__FILE__)), pre)
    pre2 = Expr(:macrocall, Symbol("@notensor"), LineNumberNode(@__LINE__, Symbol(@__FILE__)), pre2)
    post = Expr(:macrocall, Symbol("@notensor"), LineNumberNode(@__LINE__, Symbol(@__FILE__)), post)
    return Expr(:block, pre, pre2, ex, post)
end

function _construct_braidingtensors(ex::Expr)
    if TO.isdefinition(ex) || TO.isassignment(ex)
        lhs, rhs = TO.getlhs(ex), TO.getrhs(ex)
        if TO.istensorexpr(rhs)
            list = TO.gettensors(_conj_to_adjoint(rhs))
            if TO.isassignment(ex) && TO.istensor(lhs)
                obj, l, r = TO.decomposetensor(lhs)
                lhs_as_rhs = Expr(:typed_vcat, Expr(TO.prime, obj), Expr(:tuple, r...), Expr(:tuple, l...))
                push!(list, lhs_as_rhs)
            end
        else
            return ex
        end
    elseif TO.istensorexpr(ex)
        list = TO.gettensors(_conj_to_adjoint(ex))
    else
        return Expr(ex.head, map(_construct_braidingtensors, ex.args)...)
    end

    i = 1
    translatebraidings = Dict{Any,Any}()
    while i <= length(list)
        t = list[i]
        if _remove_adjoint(TO.gettensorobject(t)) == :τ
            translatebraidings[t] = Expr(:call, GlobalRef(TensorKit, :BraidingTensor))
            deleteat!(list, i)
        else
            i += 1
        end
    end

    unresolved = Any[] # list of indices that we couldn't yet figure out
    indexmaps = Dict{Any,Any}() # contains the expression to resolve the space at index indexmaps[i]
    for (t, construct_expr) in translatebraidings
        obj, leftind, rightind = TO.decomposetensor(t)
        length(leftind) == length(rightind) == 2 ||
            throw(ArgumentError("The name τ is reserved for the braiding, and should have two input and two output indices."))
        if _is_adjoint(obj)
            i1b, i2b, = leftind
            i2a, i1a, = rightind
        else
            i2b, i1b, = leftind
            i1a, i2a, = rightind
        end

        obj_and_pos1a = _findindex(i1a, list)
        obj_and_pos2a = _findindex(i2a, list)
        obj_and_pos1b = _findindex(i1b, list)
        obj_and_pos2b = _findindex(i2b, list)

        if !isnothing(obj_and_pos1a)
            indexmaps[i1b] = Expr(:call, :space, obj_and_pos1a...)
            indexmaps[i1a] = Expr(:call, :space, obj_and_pos1a...)
        elseif !isnothing(obj_and_pos1b)
            indexmaps[i1b] = Expr(TO.prime, Expr(:call, :space, obj_and_pos1b...))
            indexmaps[i1a] = Expr(TO.prime, Expr(:call, :space, obj_and_pos1b...))
        else
            push!(unresolved,(i1a,i1b))
        end

        if !isnothing(obj_and_pos2a)
            indexmaps[i2b] = Expr(:call, :space, obj_and_pos2a...)
            indexmaps[i2a] = Expr(:call, :space, obj_and_pos2a...)
        elseif !isnothing(obj_and_pos2b)
            indexmaps[i2b] = Expr(TO.prime, Expr(:call, :space, obj_and_pos2b...))
            indexmaps[i2a] = Expr(TO.prime, Expr(:call, :space, obj_and_pos2b...))
        else
            push!(unresolved,(i2a,i2b))
        end
    end
    # simple loop that tries to resolve as many indices as possible
    changed = true
    while changed == true
        changed = false
        i = 1
        while i <= length(unresolved)
            (i1, i2) = unresolved[i]
            if i1 in keys(indexmaps)
                changed = true
                indexmaps[i2] = indexmaps[i1]
                deleteat!(unresolved, i)
            elseif i2 in keys(indexmaps)
                changed = true
                indexmaps[i1] = indexmaps[i2]
                deleteat!(unresolved, i)
            else
                i += 1
            end
        end
    end
    !isempty(unresolved) &&
        throw(ArgumentError("cannot determine the spaces of indices " *
                                string(tuple(unresolved...)) *
                                    "for the braiding tensors in $(ex)"))

    pre = Expr(:block)
    for (t, construct_expr) in translatebraidings
        obj, leftind, rightind = TO.decomposetensor(t)
        if _is_adjoint(obj)
            i1b, i2b, = leftind
            i2a, i1a, = rightind
        else
            i2b, i1b, = leftind
            i1a, i2a, = rightind
        end
        push!(construct_expr.args, indexmaps[i1b])
        push!(construct_expr.args, indexmaps[i2b])
        s = gensym()
        push!(pre.args, Expr(:(=), s, construct_expr))
        ex = TO.replacetensorobjects(ex) do o, l, r
            if o == obj && l == leftind && r == rightind
                return obj  == :τ ? s : Expr(TO.prime, s)
            else
                return o
            end
        end
    end
    return Expr(:block, pre, ex)
end
_construct_braidingtensors(x) = x

function _findindex(i, list) # finds an index i in a list of tensor expressions
    for t in list
        obj, l, r = TO.decomposetensor(t)
        pos = findfirst(==(i), l)
        isnothing(pos) || return (obj, pos)
        pos = findfirst(==(i), r)
        isnothing(pos) || return (obj, pos + length(l))
    end
    return nothing
end

# since temporaries were taken out in preprocessing, they are not identified by the parsing
# step of TensorOperations, and we have to manually fix this
# Step 1: we have to find the new name that TO.tensorify assigned to these temporaries
# since it parses `tmp[] := a[] * b[]` as `newtmp = similar...; tmp = contract!(.... , newtmp, ...)`
function _update_temporaries(ex, temporaries)
    if ex isa Expr && ex.head == :(=)
        lhs = ex.args[1]
        i = findfirst(==(lhs), temporaries)
        if i !== nothing
            rhs = ex.args[2]
            if rhs isa Expr && rhs.head == :call && rhs.args[1] == :add!
                newname = rhs.args[6]
                temporaries[i] = newname
            elseif rhs isa Expr && rhs.head == :call && rhs.args[1] == :contract!
                newname = rhs.args[8]
                temporaries[i] = newname
            elseif rhs isa Expr && rhs.head == :call && rhs.args[1] == :trace!
                newname = rhs.args[6]
                temporaries[i] = newname
            else
                @error "lhs = $lhs , rhs = $rhs"
            end
        end
    elseif ex isa Expr
        for a in ex.args
            _update_temporaries(a, temporaries)
        end
    end
    return ex
end
# Step 2: we find `newtmp = similar_from_...` and replace with `newtmp = cached_similar_from...`
function _annotate_temporaries(ex, temporaries)
    if ex isa Expr && ex.head == :(=)
        lhs = ex.args[1]
        i = findfirst(==(lhs), temporaries)
        if i !== nothing
            rhs = ex.args[2]
            if !(rhs isa Expr && rhs.head == :call && rhs.args[1] == :similar_from_indices)
                @error "lhs = $lhs , rhs = $rhs"
            end
            newrhs = Expr(:call, :cached_similar_from_indices,
                            QuoteNode(lhs), rhs.args[2:end]...)
            return Expr(:(=), lhs, newrhs)
        end
    elseif ex isa Expr
        return Expr(ex.head, [_annotate_temporaries(a, temporaries) for a in ex.args]...)
    end
    return ex
end

const _TOFUNCTIONS = (:similar_from_indices, :cached_similar_from_indices,
                        :scalar, :IndexError)

function _add_modules(ex::Expr)
    if ex.head == :call && ex.args[1] in _TOFUNCTIONS
        return Expr(ex.head, GlobalRef(TensorOperations, ex.args[1]),
                        (_add_modules(ex.args[i]) for i in 2:length(ex.args))...)
    elseif ex.head == :call && ex.args[1] == :add!
        @assert ex.args[4] == :(:N)
        argind = [2,3,5,6,7,8]
        return Expr(ex.head, GlobalRef(TensorKit, Symbol(:planar_add!)),
                        (_add_modules(ex.args[i]) for i in argind)...)
    elseif ex.head == :call && ex.args[1] == :trace!
        @assert ex.args[4] == :(:N)
        argind = [2,3,5,6,7,8,9,10]
        return Expr(ex.head, GlobalRef(TensorKit, Symbol(:planar_trace!)),
                        (_add_modules(ex.args[i]) for i in argind)...)
    elseif ex.head == :call && ex.args[1] == :contract!
        @assert ex.args[4] == :(:N) && ex.args[6] == :(:N)
        argind = vcat([2,3,5], 7:length(ex.args))
        return Expr(ex.head, GlobalRef(TensorKit, Symbol(:planar_contract!)),
                        (_add_modules(ex.args[i]) for i in argind)...)
    else
        return Expr(ex.head, (_add_modules(e) for e in ex.args)...)
    end
end
_add_modules(ex) = ex

@specialize

planar_add!(α, tsrc::AbstractTensorMap{S},
            β, tdst::AbstractTensorMap{S, N₁, N₂},
            p1::IndexTuple{N₁}, p2::IndexTuple{N₂}) where {S, N₁, N₂} =
    add_transpose!(α, tsrc, β, tdst, p1, p2)

function planar_trace!(α, tsrc::AbstractTensorMap{S},
                       β, tdst::AbstractTensorMap{S, N₁, N₂},
                       p1::IndexTuple{N₁}, p2::IndexTuple{N₂},
                       q1::IndexTuple{N₃}, q2::IndexTuple{N₃}) where {S, N₁, N₂, N₃}
    if BraidingStyle(sectortype(S)) == Bosonic()
        return trace!(α, tsrc, β, tdst, p1, p2, q1, q2)
    end

    @boundscheck begin
        all(i->space(tsrc, p1[i]) == space(tdst, i), 1:N₁) ||
            throw(SpaceMismatch("trace: tsrc = $(codomain(tsrc))←$(domain(tsrc)),
                    tdst = $(codomain(tdst))←$(domain(tdst)), p1 = $(p1), p2 = $(p2)"))
        all(i->space(tsrc, p2[i]) == space(tdst, N₁+i), 1:N₂) ||
            throw(SpaceMismatch("trace: tsrc = $(codomain(tsrc))←$(domain(tsrc)),
                    tdst = $(codomain(tdst))←$(domain(tdst)), p1 = $(p1), p2 = $(p2)"))
        all(i->space(tsrc, q1[i]) == dual(space(tsrc, q2[i])), 1:N₃) ||
            throw(SpaceMismatch("trace: tsrc = $(codomain(tsrc))←$(domain(tsrc)),
                    q1 = $(q1), q2 = $(q2)"))
    end

    if iszero(β)
        fill!(tdst, β)
    elseif β != 1
        rmul!(tdst, β)
    end
    pdata = (p1..., p2...)
    for (f1, f2) in fusiontrees(tsrc)
        for ((f1′, f2′), coeff) in planar_trace(f1, f2, p1, p2, q1, q2)
            TO._trace!(α*coeff, tsrc[f1, f2], true, tdst[f1′, f2′], pdata, q1, q2)
        end
    end
    return tdst
end

_cyclicpermute(t::Tuple) = (Base.tail(t)..., t[1])
_cyclicpermute(t::Tuple{}) = ()
function reorder_indices(codA, domA, codB, domB, oindA, oindB, p1, p2)
    N₁ = length(oindA)
    N₂ = length(oindB)
    @assert length(p1) == N₁ && all(in(p1), 1:N₁)
    @assert length(p2) == N₂ && all(in(p2), N₁ .+ (1:N₂))
    oindA2 = TupleTools.getindices(oindA, p1)
    oindB2 = TupleTools.getindices(oindB, p2 .- N₁)
    indA = (codA..., reverse(domA)...)
    indB = (codB..., reverse(domB)...)
    # cycle indA to be of the form (oindA2..., reverse(cindA2)...)
    while length(oindA2) > 0 && indA[1] != oindA2[1]
        indA = _cyclicpermute(indA)
    end
    # cycle indB to be of the form (cindB2..., reverse(oindB2)...)
    while length(oindB2) > 0 && indB[end] != oindB2[1]
        indB = _cyclicpermute(indB)
    end
    for i = 2:N₁
        @assert indA[i] == oindA2[i]
    end
    for j = 2:N₂
        @assert indB[end + 1 - j] == oindB2[j]
    end
    Nc = length(indA) - N₁
    @assert Nc == length(indB) - N₂
    pc = ntuple(identity, Nc)
    cindA2 = reverse(TupleTools.getindices(indA, N₁ .+ pc))
    cindB2 = TupleTools.getindices(indB, pc)
    return oindA2, cindA2, oindB2, cindB2
end
function reorder_indices(codA, domA, codB, domB, oindA, cindA, oindB, cindB, p1, p2)
    oindA2, cindA2, oindB2, cindB2 =
        reorder_indices(codA, domA, codB, domB, oindA, oindB, p1, p2)

    #if oindA or oindB are empty, then reorder indices can only order it correctly up to a cyclic permutation!
    if isempty(oindA2) && !isempty(cindA)
         # isempty(cindA) is a cornercase which I'm not sure if we can encounter
        hit = cindA[findfirst(==(first(cindB2)), cindB)]
        while hit != first(cindA2)
            cindA2 = _cyclicpermute(cindA2)
        end
    end
    if isempty(oindB2) && !isempty(cindB)
        hit = cindB[findfirst(==(first(cindA2)), cindA)]
        while hit != first(cindB2)
            cindB2 = _cyclicpermute(cindB2)
        end
    end
    @assert TupleTools.sort(cindA) == TupleTools.sort(cindA2)
    @assert TupleTools.sort(tuple.(cindA2, cindB2)) == TupleTools.sort(tuple.(cindA, cindB))
    return oindA2, cindA2, oindB2, cindB2
end

function planar_contract!(α, A::AbstractTensorMap{S}, B::AbstractTensorMap{S},
                            β, C::AbstractTensorMap{S},
                            oindA::IndexTuple, cindA::IndexTuple,
                            oindB::IndexTuple, cindB::IndexTuple,
                            p1::IndexTuple, p2::IndexTuple,
                            syms::Union{Nothing, NTuple{3, Symbol}} = nothing) where {S}

    codA, domA = codomainind(A), domainind(A)
    codB, domB = codomainind(B), domainind(B)
    oindA, cindA, oindB, cindB =
        reorder_indices(codA, domA, codB, domB, oindA, cindA, oindB, cindB, p1, p2)

    if oindA == codA && cindA == domA
        A′ = A
    else
        if isnothing(syms)
            A′ = TO.similar_from_indices(eltype(A), oindA, cindA, A, :N)
        else
            A′ = TO.cached_similar_from_indices(syms[1], eltype(A), oindA, cindA, A, :N)
        end
        add_transpose!(true, A, false, A′, oindA, cindA)
    end
    if cindB == codB && oindB == domB
        B′ = B
    else
        if isnothing(syms)
            B′ = TO.similar_from_indices(eltype(B), cindB, oindB, B, :N)
        else
            B′ = TO.cached_similar_from_indices(syms[2], eltype(B), cindB, oindB, B, :N)
        end
        add_transpose!(true, B, false, B′, cindB, oindB)
    end
    mul!(C, A′, B′, α, β)
    return C
end

macro plansor(ex::Expr)
    return esc(plansor_parser(ex))
end

function plansor_parser(ex)
    in_tensors = TO.getinputtensorobjects(ex);
    new_tensors = TO.getnewtensorobjects(ex);

    # want the first non-braiding tensor to determine the braidingstyle
    target_tensor = in_tensors[findfirst(x->x!=:τ,in_tensors)];
    temp_t = gensym();

    ex = TO.replacetensorobjects(ex) do obj,leftind,rightind
        obj == target_tensor ? temp_t : obj
    end

    defaultparser = TO.TensorParser();
    insert!(defaultparser.preprocessors,3,_remove_braidingtensors)
    default = defaultparser(ex);
    planar = planar2_parser(ex);

    ex = quote
        $(temp_t) = $(target_tensor)

        if BraidingStyle(sectortype($temp_t)) isa Bosonic
            $(default)
        else
            $(planar)
        end
    end

    if target_tensor in new_tensors
        ex = quote
            $(ex)

            $(target_tensor) = $(temp_t)
            $(new_tensors[end]) # this is - afaik - the curent behaviour for @tensor blocks
        end
    end

    ex
end

function _remove_braidingtensors(ex::Expr)
    outgoing = [];

    if TO.isdefinition(ex) || TO.isassignment(ex)
        lhs, rhs = TO.getlhs(ex), TO.getrhs(ex)
        if TO.istensorexpr(rhs)
            list = TO.gettensors(_conj_to_adjoint(rhs))
            if TO.istensor(lhs)
                obj, l, r = TO.decomposetensor(lhs)
                outgoing = [l;r];
            end
        else
            return ex
        end
    elseif TO.istensorexpr(ex)
        list = TO.gettensors(_conj_to_adjoint(ex))
    else
        return Expr(ex.head, map(_remove_braidingtensors, ex.args)...)
    end

    τs = Any[];
    i = 1;
    while i <= length(list)
        t = list[i]
        if _remove_adjoint(TO.gettensorobject(t)) == :τ
            push!(τs,t)
            deleteat!(list, i)
        else
            i += 1
        end
    end

    indicemaps = Dict{Any,Any}() # to remove the braidingtensors, we need to map certain indices to other indices
    for t in τs
        obj, leftind, rightind = TO.decomposetensor(t)
        length(leftind) == length(rightind) == 2 ||
            throw(ArgumentError("The name τ is reserved for the braiding, and should have two input and two output indices."))
        if _is_adjoint(obj)
            i1b, i2b, = leftind
            i2a, i1a, = rightind
        else
            i2b, i1b, = leftind
            i1a, i2a, = rightind
        end

        i1a = get(indicemaps,i1a,i1a);
        i1b = get(indicemaps,i1b,i1b);
        i2a = get(indicemaps,i2a,i2a);
        i2b = get(indicemaps,i2b,i2b);

        obj_and_pos1a = _findindex(i1a, list)
        obj_and_pos2a = _findindex(i2a, list)
        obj_and_pos1b = _findindex(i1b, list)
        obj_and_pos2b = _findindex(i2b, list)

        if i1a in outgoing
            indicemaps[i1a] = i1a;
            indicemaps[i1b] = i1a;
        elseif i1b in outgoing
            indicemaps[i1a] = i1b;
            indicemaps[i1b] = i1b;
        else
            if i1a isa Int && i1b isa Int
                indicemaps[i1a] = max(i1a,i1b);
                indicemaps[i1b] = max(i1a,i1b);
            else
                indicemaps[i1a] = i1a;
                indicemaps[i1b] = i1a;
            end
        end

        if i2a in outgoing
            indicemaps[i2a] = i2a;
            indicemaps[i2b] = i2a;
        elseif i2b in outgoing
            indicemaps[i2a] = i2b;
            indicemaps[i2b] = i2b;
        else
            if i2a isa Int && i2b isa Int
                indicemaps[i2a] = max(i2a,i2b);
                indicemaps[i2b] = max(i2a,i2b);
            else
                indicemaps[i2a] = i2a;
                indicemaps[i2b] = i2a;
            end
        end
    end

    # simple loop that tries to simplify the indicemaps (a=>b,b=>c -> a=>c,b=>c)
    changed = true
    while changed == true
        changed = false
        i = 1
        for (k,v) in indicemaps
            if v in keys(indicemaps) && indicemaps[v] != v
                changed = true;
                indicemaps[k] = indicemaps[v]
            end
        end
    end

    ex = TO.replaceindices(i -> get(indicemaps,i,i),ex)
    ex = _purge_braidingtensors!(ex)
    return ex
end
_remove_braidingtensors(x) = x

function _purge_braidingtensors!(ex::Expr) # actually remove the braidingtensors
    filter!(ex.args) do a
        !(a isa Expr) || (
            !(a.args[1] == :τ) &&
            !(a.head == :call && a.args[1] == :conj && a.args[2] isa Expr &&
                a.args[2].args[1] == :τ && a.args[2].head in (:typed_vcat,:typed_hcat))
            )
    end

    if ex.head == :call && ex.args[1] == :* && length(ex.args) == 2 # multiplication with only a single argument is (rightfully) seen as invalid syntax
        _purge_braidingtensors!(ex.args[2])
    else
        Expr(ex.head,map(_purge_braidingtensors!,ex.args)...)
    end
end
_purge_braidingtensors!(x) = x


text:
function convective_adjust!(x)
    # remove negative gradients from temperature profile
    for i in length(x)-3:-1:2
        if x[i] > x[i+1]
            if x[i-1] > x[i]; x[i] = x[i+1]
            else; x[i] = (x[i-1]+x[i+1])/2
            end
        end
    end
end

ext: jl

max_stars_repo_name: adelinehillier/LearnConvection

max_issues_repo_name: adelinehillier/LearnConvection

max_forks_repo_name: adelinehillier/LearnConvection


text:
import DrWatson: quickactivate, savename
quickactivate(@__DIR__, "Chemostat_Heerden2013")

import SparseArrays
import Base.Threads: @threads, threadid, SpinLock
using Serialization

# -------------------------------------------------------------------
# Run add https://github.com/josePereiro/Chemostat_Heerden2013.jl in the Julia Pkg REPL to install the
# package, then you must activate the package enviroment (see README)
@time begin
    import Chemostat_Heerden2013
    const Hd  = Chemostat_Heerden2013.HeerdenData;
    const BD  = Chemostat_Heerden2013.BegData;
    const iJR = Chemostat_Heerden2013.iJR904
end

# -------------------------------------------------------------------
# run add "https://github.com/josePereiro/Chemostat" in the
# julia Pkg REPL for installing the package
@time begin
    import Chemostat
    import Chemostat.LP.MathProgBase
    const ChU = Chemostat.Utils
    const ChSS = Chemostat.SteadyState
    const ChLP = Chemostat.LP
    const ChEP = Chemostat.MaxEntEP
    const ChSU = Chemostat.SimulationUtils
end

## ----------------------------------------------------------------------
# Tooling
function err_str(err; max_len = 500)
    s = sprint(showerror, err, catch_backtrace())
    return length(s) > max_len ? s[1:max_len] * "\n[...]" : s
end

function partial_test(model, title  = "PARTIAL TEST")
    withcost = iJR.COST_IDER in model.rxns
    iders = withcost ? [iJR.BIOMASS_IDER, iJR.COST_IDER] : [iJR.BIOMASS_IDER]
    fbaout = ChLP.fba(model, iders...);
    ChU.tagprintln_inmw(title,
        "\nsize:             ", size(model),
        "\nobj_ider:         ", iJR.BIOMASS_IDER,
        "\nfba obj_val:      ", ChU.av(model, fbaout, iJR.BIOMASS_IDER),
        "\ncost_ider:        ", withcost ? iJR.COST_IDER : "not found",
        "\nfba cost_val:     ", withcost ? ChU.av(model, fbaout, iJR.COST_IDER) : "not found",
    )
end

## ----------------------------------------------------------------------
function prepare_model(exp; verbose = true)
    # prepare model
    model = ChU.MetNet(ChU.load_data(iJR.BASE_MODELS_FILE; verbose))
    model = ChU.uncompressed_model(model)

    exp_growth = Hd.val("D", exp)
    expξ = Hd.val("xi", exp)
    intake_info = deepcopy(iJR.base_intake_info)
    intake_info["EX_glc_LPAREN_e_RPAREN_"]["c"] = Hd.val("cGLC", exp)
    ChSS.apply_bound!(model, expξ, intake_info)
    return model
end

## ----------------------------------------------------------------------
let
    # setup
    exp = 1
    model = prepare_model(exp)
    M, N = size(model)
    obj_ider = iJR.BIOMASS_IDER
    obj_idx = ChU.rxnindex(model, obj_ider)
    alpha = Inf
    maxiter = Int(1e5)
    minvar, maxvar = 1e-45, 1e45
    damp = 0.99
    fbaout = ChLP.fba(model, obj_idx)
    partial_test(model); println()

    # seed
    seed_file = joinpath(iJR.MODEL_PROCESSED_DATA_DIR, "ept0_epout_seed.bson")
    if !isfile(seed_file)
        seed = ChEP.maxent_ep(model; alpha, epsconv = 1e4, damp, maxiter = Int(1e7))
        ChU.save_data(seed_file, seed)
    else
        seed = ChU.load_data(seed_file)
    end
    println()

    # Simulations
    write_lock = ReentrantLock()

    # @threads
    i_epsconvs = collect(enumerate(10.0.^-(4:7)))
    Nepsconvs = length(i_epsconvs)
    @threads for (epsconvi, epsconv) in i_epsconvs

        # cache
        fname = savename("exploration_", (;epsconv = string(epsconv)), "jld")
        fpath = joinpath(iJR.MODEL_PROCESSED_DATA_DIR, fname)
        if isfile(fpath)
            lock(write_lock) do
                @info "Cache found SKIPING" threadid() epsconvi Nepsconvs epsconv fname
                println()
            end
            continue
        end

        beta_range = [0.0; 10.0.^(-1:0.1:30)]
        Nbeta_range = length(beta_range)

        lmodel = deepcopy(model)

        # MaxEnt EP
        beta_vec = zeros(N)
        epoutT0 = deepcopy(seed)

        D = Dict()
        D["epouts"] = Dict()
        D["error"] = ""
        for (betai, beta) in enumerate(beta_range)
            beta_vec[obj_idx] = beta

            lock(write_lock) do
                @info("Doing",
                    threadid(), epsconvi,
                    Nepsconvs, epsconv, betai,
                    Nbeta_range, beta,
                    length(D["epouts"]),
                    epoutT0.iter
                ); println()
            end

            try
                epoutT0 = ChEP.maxent_ep(
                    lmodel;
                    alpha, beta_vec,
                    epsconv, minvar, maxvar,
                    solution = epoutT0,
                    maxiter, damp,
                    iter0 = 0,
                    verbose = false
                )
                if epoutT0.status != ChEP.CONVERGED_STATUS
                    lock(write_lock) do
                        @info(string("At threadid() = ", threadid()),
                            epoutT0.status
                        ); println()
                        D["error"] = "status = $(epoutT0.status)"
                    end
                    break
                end
            catch err
                lock(write_lock) do
                    @info(string("At threadid() = ", threadid()))
                    msg = err_str(err; max_len = 500)
                    @error(msg)
                    D["error"] = msg
                    println()
                end
                break
            end

            D["epouts"][beta] = epoutT0
        end # for (betai, beta)

        # store and save
        D["beta_range"] = collect(keys(D["epouts"])) |> unique |> sort
        D["model"] = lmodel
        D["fbaout"] = fbaout
        D["minvar"], D["maxvar"] = minvar, maxvar
        D["epsconv"] = epsconv
        D["alpha"] = alpha

        # store
        serialize(fpath, D)

        lock(write_lock) do
            println()
            @info("Finished",
                threadid(), epsconvi,
                Nepsconvs, epsconv,
                length(D["epouts"]),
                fname
            ); println()
        end

    end # @threads for var_order in var_orders
end


text:
export BoxConstraint, NormConstraint

abstract type Constraint end

repair!(c::Constraint, x) = x

repair(c::Constraint, x) = repair!(c, copy(x))

getfitness(f, c::Constraint, x) = f(x)

# NoConstraint
mutable struct NoConstraint <: Constraint
end

# BoxConstraint lo <= x <= hi
mutable struct BoxConstraint{T} <: Constraint
    lo::Vector{T} # lower bound
    hi::Vector{T} # upper bound
    α::T          # penalty scaling factor
    BoxConstraint(lo, hi, α = one(eltype(lo))) = new{eltype(lo)}(lo, hi, α)
end

repair!(c::BoxConstraint, x) = map!(clamp, x, x, c.lo, c.hi)

function getfitness(f, c::BoxConstraint, x)
    x_repair = repair(c, x)
    violation = norm(x .- x_repair)
    f(x_repair) + c.α * violation
end

# NormConstraint norm(x, p) <= θ
mutable struct NormConstraint{T} <: Constraint
    θ::T        # max norm
    p::Int      # p-norm
    α::T        # penalty scaling factor
    NormConstraint(θ, p = 2, α = one(θ)) = new{typeof(θ)}(θ, p, α)
end

function repair!(c::NormConstraint, x)
    nrm = norm(x, c.p)
    if nrm > c.θ
        rmul!(x, c.θ / nrm)
    end
    return x
end

function getfitness(f, c::NormConstraint, x)
    x_repair = repair(c, x)
    violation = norm(x .- x_repair, c.p)
    f(x_repair) + c.α * violation
end


text:
module CameraCalibrations

using Statistics, StaticArrays, CoordinateTransformations, ImageCore, ImageTransformations, Interpolations, MATLAB

export buildcalibration, calibrate

const V2 = SVector{2, Float64}

include("buildcalibrations.jl")

include("calibrate.jl")

end


text:
@testset "5.5 Inverse secant" begin
    include("5.5.1 u (a+b arcsec(c x))^n.jl")
    include("5.5.2 Inverse secant functions.jl")
end

max_stars_repo_path: test/5 Inverse trig functions/5.5 Inverse secant/.jl


text:
"""
    MOSVGP(args...; kwargs...)

Multi-Output Sparse Variational Gaussian Process

## Arguments
- `X::AbstractArray` : : Input features, if `X` is a matrix the choice of colwise/rowwise is given by the `obsdim` keyword
- `y::AbstractVector{<:AbstractVector}` : Output labels, each vector corresponds to one output dimension
- `kernel::Union{Kernel,AbstractVector{<:Kernel}` : covariance function or vector of covariance functions, can be either a single kernel or a collection of kernels for multiclass and multi-outputs models
- `likelihood::Union{AbstractLikelihood,Vector{<:Likelihood}` : Likelihood or vector of likelihoods of the model. For compatibilities, see [`Likelihood Types`](@ref likelihood_user)
- `inference` : Inference for the model, for compatibilities see the [`Compatibility Table`](@ref compat_table))
- `nLatent::Int` : Number of latent GPs
- `nInducingPoints` : number of inducing points, or collection of inducing points locations

## Keyword arguments
- `verbose::Int` : How much does the model print (0:nothing, 1:very basic, 2:medium, 3:everything)
- `optimiser` : Optimiser used for the kernel parameters. Should be an Optimiser object from the [Flux.jl](https://github.com/FluxML/Flux.jl) library, see list here [Optimisers](https://fluxml.ai/Flux.jl/stable/training/optimisers/) and on [this list](https://github.com/theogf/AugmentedGaussianProcesses.jl/tree/master/src/inference/optimisers.jl). Default is `ADAM(0.001)`
- `Zoptimiser` : Optimiser used for the inducing points locations
- `Aoptimiser` : Optimiser used for the mixing parameters.
- `atfrequency::Int=1` : Choose how many variational parameters iterations are between hyperparameters optimization
- `mean=ZeroMean()` : PriorMean object, check the documentation on it [`MeanPrior`](@ref meanprior)
- `obsdim::Int=1` : Dimension of the data. 1 : X ∈ DxN, 2: X ∈ NxD
"""
mutable struct MOSVGP{
    T<:Real,
    TLikelihood<:AbstractLikelihood,
    TInference<:AbstractInference,
    TData<:AbstractDataContainer,
    N,
    Q,
} <: AbstractGP{T,TLikelihood,TInference,N}
    data::TData
    nFeatures::Vector{Int64} # Number of features of the GP (equal to number of points)
    nf_per_task::Vector{Int64}
    f::NTuple{Q,SparseVarLatent}
    likelihood::Vector{TLikelihood}
    inference::TInference
    A::Vector{Vector{Vector{T}}}
    A_opt::Any
    verbose::Int64
    atfrequency::Int64
    trained::Bool
end

function MOSVGP(
    X::AbstractArray,
    y::AbstractVector{<:AbstractVector},
    kernel::Union{Kernel,AbstractVector{<:Kernel}},
    likelihood::Union{AbstractLikelihood,AbstractVector{<:AbstractLikelihood}},
    inference::AbstractInference,
    nLatent::Int,
    nInducingPoints::Union{Int,AbstractVector};
    verbose::Int=0,
    atfrequency::Int=1,
    mean::Union{<:Real,AbstractVector{<:Real},PriorMean}=ZeroMean(),
    variance::Real=1.0,
    optimiser=ADAM(0.01),
    Aoptimiser=ADAM(0.01),
    Zoptimiser=false,
    obsdim::Int=1,
)
    @assert length(y) > 0 "y should not be an empty vector"
    nTask = length(y)

    X, T = wrap_X(X, obsdim)

    likelihoods = if likelihood isa AbstractLikelihood
        likelihoods = [deepcopy(likelihood) for _ in 1:nTask]
    else
        likelihood
    end

    nf_per_task = zeros(Int64, nTask)
    corrected_y = Vector(undef, nTask)
    for i in 1:nTask
        corrected_y[i], nf_per_task[i], likelihoods[i] = check_data!(y[i], likelihoods[i])
    end

    inference isa AnalyticVI || error("The inference object should be of type `AnalyticVI`")
    all(implemented.(likelihood, Ref(inference))) ||
        error("The $likelihood is not compatible or implemented with the $inference")

    data = wrap_data(X, corrected_y)

    if mean isa Real
        mean = ConstantMean(mean)
    elseif mean isa AbstractVector{<:Real}
        mean = EmpiricalMean(mean)
    end

    if isa(optimiser, Bool)
        optimiser = optimiser ? ADAM(0.01) : nothing
    end

    if isa(Aoptimiser, Bool)
        Aoptimiser = Aoptimiser ? ADAM(0.01) : nothing
    end

    kernel = if kernel isa Kernel
        [kernel]
    else
        length(kernel) == nLatent ||
            error("Number of kernels should be equal to the number of tasks")
        kernel
    end
    nKernel = length(kernel)

    nInducingPoints = if nInducingPoints isa AbstractVector{<:AbstractVector{<:AbstractVector}}
        nInducingPoints
    elseif nInducingPoints isa AbstractVector{<:AbstractVector}
        [deepcopy(nInducingPoints) for _ in 1:nLatent]
    elseif nInducingPoints isa Int
        Zref = InducingPoints(KMeansAlg(nInducingPoints), X)
        [deepcopy(Zref) for _ in 1:nLatent]
    end

    nFeatures = size.(Z, 1)

    _nMinibatch = nSamples(data)
    if isStochastic(inference)
        0 < nMinibatch(inference) < nSamples || error(
            "The size of mini-batch $(nMinibatch(inference)) is incorrect (negative or bigger than number of samples), please set nMinibatch correctly in the inference object",
        )
        _nMinibatch = nMinibatch(inference)
    end

    latent_f = ntuple(
        i -> _SVGP{T}(
            nFeatures[i],
            _nMinibatch,
            Z[mod(i, nLatent) + 1],
            kernel[mod(i, nKernel) + 1],
            mean,
            optimiser,
            Zoptimiser,
        ),
        nLatent,
    )

    A = [
        [x -> x / sqrt(sum(abs2, x))(randn(T, nLatent)) for i in 1:nf_per_task[j]] for
        j in 1:nTask
    ]

    likelihoods .=
        init_likelihood.(likelihoods, inference, nf_per_task, _nMinibatch, nFeatures)
    xview = view_x(data, collect(range(1, _nMinibatch; step=1)))
    yview = view_y(likelihood, data, 1:nSamples(data))
    inference = tuple_inference(
        inference, nLatent, nFeatures, nSamples(data), _nMinibatch, xview, yview
    )

    return MOSVGP{T,eltype(likelihoods),typeof(inference),nTask,nLatent}(
        X,
        corrected_y,
        nSamples,
        nDim,
        nFeatures,
        nLatent,
        nX,
        nTask,
        nf_per_task,
        latent_f,
        likelihoods,
        inference,
        A,
        Aoptimiser,
        verbose,
        atfrequency,
        false,
    )
end

function Base.show(io::IO, model::MOSVGP)
    return print(
        io,
        "Multioutput Sparse Variational Gaussian Process with the likelihoods $(likelihood(model)) infered by $(inference(model)) ",
    )
end

@traitimpl IsMultiOutput{MOSVGP}

nOutput(::MOSVGP{<:Real,<:AbstractLikelihood,<:AbstractInference,N,Q}) where {N,Q} = Q
Zviews(m::MOSVGP) = Zview.(m.f)
objective(m::MOSVGP) = ELBO(m)


text:
using ModelingToolkit
using LinearAlgebra
using SparseArrays: sparse
using Test

@variables a,b,c,d

# test hashing
aa = a; # old a

@variables a

@test isequal(a, aa)
@test hash(a) == hash(aa)

@test isequal(get_variables(a+aa+1), [a])

@test hash(a+b ~ c+d) == hash(a+b ~ c+d)

# test some matrix operations don't throw errors
X = [a b;c d]
det(X)
lu(X)
inv(X)
qr(X)

# test operations with sparse arrays and Operations
# note `isequal` instead of `==` because `==` would give another Operation

# test that we can create a sparse array of Operation
Oarray = zeros(Num, 2,2)
Oarray[2,2] = a
@test isequal(sparse(Oarray), sparse([2], [2], [a]))

# test Operation * sparse
@test isequal(a * sparse([2], [2], [1]), sparse([2], [2], [a * 1]))

# test sparse{Operation} + sparse{Operation}
A = sparse([2], [2], [a])
B = sparse([2], [2], [b])
@test isequal(A + B, sparse([2], [2], [a+b]))

# test sparse{Operation} * sparse{Operation}
C = sparse([1, 2], [2, 1], [c, c])
D = sparse([1, 2], [2, 1], [d, d])

@test isequal(C * D, sparse([1,2], [1,2], [c * d, c * d]))

@parameters t σ ρ β
@variables x(t) y(t) z(t)
@derivatives D'~t Dx'~x Dy'~y Dz'~z
eqs = [D(x) ~ σ*(y-x),
       D(y) ~ x*(ρ-z)-y,
       D(z) ~ x*y - β*z]
J = Num[Dx(eqs[1].rhs) Dy(eqs[1].rhs) Dz(eqs[1].rhs)
 Dx(eqs[2].rhs) Dy(eqs[2].rhs) Dz(eqs[2].rhs)
 Dx(eqs[3].rhs) Dy(eqs[3].rhs) Dz(eqs[3].rhs)]

J = expand_derivatives.(J)
using LinearAlgebra
luJ = lu(J,Val(false))

using ModelingToolkit
@variables M[1:2,1:2]
inv(M)

@variables b[1:2]
M = [1 0; 0 2]
M \ b
M \ reshape(b,2,1)
M = [1 1; 0 2]
M \ reshape(b,2,1)

M = [1 a; 0 2]
M \ b
M \ [1, 2]


text:
using TidalFluxQuantities, Base.Dates
using Base.Test

ts = collect(now():now()+Day(1))

@testset "Stage" begin
    ht = Stage(ts,randn(2))
    hh = Stage(ts,randn(2))

    hs = quantity(ht)

    @test (ts,hs) == unzip(ht)

end

@testset "CrossSectionalArea" begin
    At = CrossSectionalArea(ts,randn(2))
    AA = CrossSectionalArea(ts,randn(2))

end

@testset "Velocity" begin
    vt = Velocity(ts,[(randn(3)...);(randn(3)...)])
    vv = Velocity(ts,[(randn(3)...);(randn(3)...)])

end

@testset "AlongChannelVelocity" begin
    cvt = AlongChannelVelocity(ts,randn(2))
    cvv = AlongChannelVelocity(ts,randn(2))

end

@testset "Discharge" begin
    dd = Discharge(ts,randn(2))
    dt = Discharge(ts,randn(2))

end

include("quantity_ops.jl")
include("measurements.jl")
include("fluxes.jl")
include("array_interface.jl")


text:
module Cutoffs

using JuLIP.Potentials: cutsw,
                        cutsw_d,
                        coscut,
                        coscut_d

import JuLIP
using Roots, StaticArrays

const cutsp = JuLIP.Potentials.fcut
const cutsp_d = JuLIP.Potentials.fcut_d

import JuLIP: decode_dict
import JuLIP.Potentials: cutoff
import NBodyIPs: fcut, fcut_d, NBCutoff
import Base: Dict, convert

export CosCut, CosCut2s, PolyCut, PolyCutSym, PolyCut2sA

"""
`struct CosCut{T} <: NBCutoff`

### Constructor
```
CosCut(rc1, rc2)
```

### Cut-off function
```
fcut(r) = (r < rc1) + (rc1 <= r < rc2) * (1 + cos(...))/2
```
where `cos(...)` is a rescaled cos to blend from one to zero.
"""
struct CosCut{T} <: NBCutoff
   rc1::T
   rc2::T
end
cutoff(C::CosCut) = C.rc2
fcut(C::CosCut, r::Number) = coscut(r, C.rc1, C.rc2)
fcut_d(C::CosCut, r::Number) = coscut_d(r, C.rc1, C.rc2)
Dict(C::CosCut) = Dict("__id__" => "NBodyIPs_CosCut",
                       "rc1" => C.rc1, "rc2" => C.rc2)
CosCut(D::Dict) = CosCut(D["rc1"], D["rc2"])
convert(::Val{:NBodyIPs_CosCut}, D::Dict) = CosCut(D)

"""
`struct CosCut2s{T} <: NBCutoff`

### Constructor
```
CosCut2s(ri1, ri2, ro1, ro2)
```

### Cut-off function
```
fcut(r) = (1-coscut(r, ri1, ri2)) * coscut(r, ro1, ro2)
```
"""
struct CosCut2s{T} <: NBCutoff
   ri1::T
   ri2::T
   ro1::T
   ro2::T
end
cutoff(C::CosCut2s) = C.ro2
fcut(C::CosCut2s, r::Number) = (1-coscut(r, C.ri1, C.ri2)) * coscut(r, C.ro1, C.ro2)
fcut_d(C::CosCut2s, r::Number) = (- coscut_d(r, C.ri1, C.ri2) * coscut(r, C.ro1, C.ro2)
   + (1-coscut(r, C.ri1, C.ri2)) * coscut_d(r, C.ro1, C.ro2))
Dict(C::CosCut2s) = Dict( "__id__" => "NBodyIPs_CosCut2s",
                        "ri1" => C.ri1, "ri2" => C.ri2,
                        "ro1" => C.ro1, "ro2" => C.ro2 )
CosCut2s(D::Dict) = CosCut2s(D["ri1"], D["ri2"], D["ro1"], D["ro2"])
convert(::Val{:NBodyIPs_CosCut2s}, D::Dict) = CosCut2s(D)

"""
`struct PolyCut{TI, T} <: NBCutoff`

### Constructor
```
PolyCut(p, rc)
```

### Cut-off function
```
fcut(r) = (r/rc-1)^p
```
"""
struct PolyCut{TI, T} <: NBCutoff
   p::TI
   rc::T
end
cutoff(C::PolyCut) = C.rc
fcut(C::PolyCut, r::Number) = @fastmath((r/C.rc - 1)^C.p * (r < C.rc))
fcut_d(C::PolyCut, r::Number) = @fastmath(C.p/C.rc * (r/C.rc - 1)^(C.p-1) * (r < C.rc))
Dict(C::PolyCut) = Dict( "__id__" => "NBodyIPs_PolyCut",
                        "p" => C.p, "rc" => C.rc)
PolyCut(D::Dict) = PolyCut(D["p"], D["rc"])
convert(::Val{:NBodyIPs_PolyCut}, D::Dict) = PolyCut(D)

"""
`struct PolyCutSym{TI, T} <: NBCutoff`

### Constructor
```
PolyCutSym(p, rc)
```

### Cut-off function
```
fcut(r) = (r/rc-1)^p (r/rc+1)^p
```
"""
struct PolyCutSym{TI, T} <: NBCutoff
   p::TI
   rc::T
end
cutoff(C::PolyCutSym) = C.rc
fcut(C::PolyCutSym, r::Number) = @fastmath(((r/C.rc)^2 - 1)^C.p * (r<C.rc))
fcut_d(C::PolyCutSym, r::Number) = @fastmath(2*C.p*r/C.rc^2 * ((r/C.rc)^2 - 1)^(C.p-1) * (r<C.rc))
Dict(C::PolyCutSym) = Dict( "__id__" => "NBodyIPs_PolyCutSym",
                        "p" => C.p, "rc" => C.rc)
PolyCutSym(D::Dict) = PolyCutSym(D["p"], D["rc"])
convert(::Val{:NBodyIPs_PolyCutSym}, D::Dict) = PolyCutSym(D)

"""
`struct PolyCut2sA{T} <: NBCutoff`

experimental...
"""
struct PolyCut2sA{T} <: NBCutoff
   r0::T
   rnn::T
   rc::T
   C::T
   λ::T
end

cutoff(C::PolyCut2sA) = C.rc

function PolyCut2sA(r0, rnn, rc)
   gg = λ -> exp( λ*(rc/rnn - 1) ) + exp( λ*(r0/rnn - 1) ) - 2
   λopt = find_zero(gg, -2.5)
   if abs(λopt) < 1e-5
      @error("`Roots` found a *bad* λ parameter.")
   end
   C = 1 / (exp( λopt*(rc/rnn - 1) ) - 1)
   return PolyCut2sA(r0, rnn, rc, C, λopt)
end

function fcut(C::PolyCut2sA, r::Number)
   x = @fastmath(C.C * (exp( C.λ*(r/C.rnn - 1) ) - 1))
   return (x^2-1)^2 * (C.r0 < r < C.rc)
end

function fcut_d(C::PolyCut2sA, r::Number)
   e =  @fastmath(exp( C.λ*(r/C.rnn - 1) ))
   x = C.C * (e - 1)
   dx = C.C * e * C.λ / C.rnn
   return 4*x*(x^2-1) * dx * (C.r0 < r < C.rc)
end

Dict(C::PolyCut2sA) = Dict( "__id__" => "NBodyIPs_PolyCut2sA",
                            "r0" => C.r0, "rnn" => C.rnn, "rc" => C.rc,
                            "C" => C.C, "lam" => C.λ )

PolyCut2sA(D::Dict) = PolyCut2sA(D["r0"], D["rnn"], D["rc"], D["C"], D["lam"])

convert(::Val{:NBodyIPs_PolyCut2sA}, D::Dict) = PolyCut2sA(D)

# -------------- Some General cutoff mechanics

fcut(C::NBCutoff, r::AbstractVector) = prod(fcut.(Ref(C), r))

function fcut_d(C::NBCutoff, r::AbstractVector)
   fc = fcut.(Ref(C), r)
   fctot = prod(fc)
   if fctot == 0
      return fctot, zero(typeof(r))
   end
   fc_d = fcut_d.(Ref(C), r)
   return fctot, (fctot * fc_d) ./ fc
end

end

# struct SplineCut{T} <: NBCutoff
# end
#
# struct SplineCut2s{T} <: NBCutoff
# end
#
# struct CosEnvCut{T} <: NBCutoff
# end
#
# struct CosEnvCut2s{T} <: NBCutoff
# end

max_stars_repo_name: cortner/ManyBodyIPs

max_issues_repo_name: cortner/ManyBodyIPs

max_forks_repo_name: cortner/ManyBodyIPs


text:
using ArgParse, Distributed, DrWatson

const aps = ArgParseSettings(version="1.0", add_version=true)

add_arg_group(aps, "Input and Output")
@add_arg_table aps begin
    "input"
        help = "path to the video or directory to process"
        arg_type = String
        default = datadir("videos")
end

add_arg_group(aps, "Evaluation Parameters")
@add_arg_table aps begin
    "--nperms"
        help = "number of permutations to use for significance testing"
        arg_type = Int
        default = 100000
        range_tester = p -> p > 0
end

add_arg_group(aps, "Analyses")
@add_arg_table aps begin
    "--skip-mi"
        help = "skip the mutual information analysis"
        action = :store_true
    "--skip-te"
        help = "skip the transfer entropy analysis"
        action = :store_true
end

add_arg_group(aps, "Worker Process Control")
@add_arg_table aps begin
    "--procs"
        help = "number of worker processes; must be non-negative"
        arg_type = Int
        default = 0
        range_tester = p -> p ≥ 0
    "--slurm"
        help = "use the SLURM cluster manager"
        action = :store_true
end

args = parse_args(ARGS, aps)

if args["procs"] != 0
    if args["slurm"]
        using ClusterManagers
        addprocs(SlurmManager(args["procs"]))
    else
        addprocs(args["procs"])
    end

    @everywhere begin
        using Pkg
        Pkg.activate(".")
    end
end

@everywhere begin
    using Base.Iterators, DelimitedFiles, DrWatson

    include("src/load.jl")
    include("src/info.jl")
    include("src/binning.jl")
    include("src/analysis.jl")
end

@everywhere filename(path::AbstractString) = basename(first(splitext(path)))

@everywhere function str2sym(d::Dict{String})
    e = Dict{Symbol,Any}()
    for (k, v) in d
        e[Symbol(k)] = v
    end
    e
end

@everywhere function process(filepath::AbstractString, grid::NTuple{2,Int}, nperms::Int;
                             crops=nothing, skipmi=false, skipte=false)
    if !(skipmi && skipte)
        videopath = relpath(filepath, datadir())
        parameters = str2sym(parse_savename(filepath)[2])
        merge!(parameters, Dict(:gh => first(grid), :gw => last(grid), :nperms => nperms))

        fs = frames(filepath)
        if !isnothing(crops)
            fs = crop(fs, crops...)
        end
        greenchannel = green(fs)
        greengrid = coarse(greenchannel, grid...)
        greenbinned = bin(greengrid)

        if !skipmi
            analyze(MIAnalysis(nperms=nperms), greenbinned, parameters, :lag, :mi, videopath)
        end
        if !skipte
            analyze(TEAnalysis(nperms=nperms), greenbinned, parameters, :k, :te, videopath)
        end
    end
end

function process(filepath, nperms; skipmi=false, skipte=false)
    for grid in [(1,5), (5,1), (1,10), (10,1), (30,1), (1,30), (5,5), (10,10), (30,30)]
        process(filepath, grid, nperms; crops=(400, 400, :), skipmi=skipmi, skipte=skipte)
    end
end

ismov(f) = last(splitext(f)) == ".mov"

function main(input, nperms; skipmi=false, skipte=false)
    input = abspath(input)
    if isfile(input) && ismov(input)
        process(input, nperms; skipmi=skipmi, skipte=skipte)
    elseif isdir(input)
        video_found = false
        for (root, _, files) in walkdir(input)
            for file in filter(ismov, files)
                video_found = true
                process(joinpath(root, file), nperms; skipmi=skipmi, skipte=skipte)
            end
        end
        if !video_found
            @warn "no valid movie file found in $input"
        end
    else
        @error "\"$input\" is not a path to a video file or directory"
    end
end

@time main(args["input"], args["nperms"]; skipmi=args["skip-mi"], skipte=args["skip-te"])

max_stars_repo_name: ELIFE-ASU/dglmoore

max_issues_repo_name: ELIFE-ASU/dglmoore

max_forks_repo_name: ELIFE-ASU/dglmoore


text:
struct Sets
    REGION          ::Vector{Symbol}
    FUEL            ::Vector{Symbol}
    TECH            ::Vector{Symbol}
    CLASS           ::Dict{Symbol, Vector{Symbol}}
    STORAGECLASS    ::Dict{Symbol, Vector{Symbol}}
    HOUR            ::UnitRange{Int64}
    techtype        ::Dict{Symbol,Symbol}
    techfuel        ::Dict{Symbol,Symbol}
    reservoirclass  ::Dict{Symbol, Vector{Symbol}}
    dataregions     ::Vector{Symbol}
end

# need to specify AxisArray types in more detail
struct Params
    cf                          ::AxisArray{Float64,4}
    transmissionlosses          ::AxisArray{Float64,2}
    demand                      ::AxisArray{Float64,2}
    hydrocapacity               ::AxisArray{Float64,2}
    cfhydroinflow               ::AxisArray{Float64,3}
    classlimits                 ::AxisArray{Float64,3}
    transmissionislands         ::AxisArray{Bool,2}
    efficiency                  ::AxisArray{Float64,1}
    rampingrate                 ::AxisArray{Float64,1}
    dischargetime               ::AxisArray{Float64,3}
    initialstoragelevel         ::Float64
    minflow_existinghydro       ::Float64
    emissionsCO2                ::AxisArray{Float64,1}
    fuelcost                    ::AxisArray{Float64,1}
    variablecost                ::AxisArray{Float64,1}
    smalltransmissionpenalty    ::Float64
    investcost                  ::AxisArray{Float64,2}
    crf                         ::AxisArray{Float64,1}
    fixedcost                   ::AxisArray{Float64,1}
    transmissioninvestcost      ::AxisArray{Float64,2}
    transmissionfixedcost       ::AxisArray{Float64,2}
    hydroeleccost               ::AxisArray{Float64,2}
    solarcombinedarea           ::AxisArray{Float64,3}
    pv_density                  ::Float64
    csp_density                 ::Float64
    cspsolarmultiple            ::Float64
end

struct Vars
    Systemcost                  ::JuMP.JuMPArray{JuMP.Variable,1}
    CO2emissions                ::JuMP.JuMPArray{JuMP.Variable,1}
    FuelUse                     ::JuMP.JuMPArray{JuMP.Variable,2}
    Electricity                 ::JuMP.JuMPDict{JuMP.Variable,4}
    AnnualGeneration            ::JuMP.JuMPArray{JuMP.Variable,2}
    Charging                    ::JuMP.JuMPArray{JuMP.Variable,3}
    StorageLevel                ::JuMP.JuMPDict{JuMP.Variable,4}
    Transmission                ::JuMP.JuMPArray{JuMP.Variable,3}
    TransmissionCapacity        ::JuMP.JuMPArray{JuMP.Variable,2}
    Capacity                    ::JuMP.JuMPDict{JuMP.Variable,3}
    SolarCapacity               ::JuMP.JuMPArray{JuMP.Variable,4}
end

# add type info here too
struct Constraints
    ElecCapacity
    ElecDemand
    RampingDown
    RampingUp
    StorageBalance
    MaxStorageCapacity
    InitialStorageLevel
    MaxTransmissionCapacity
    TwoWayStreet
    NoTransmission
    NoHydroCharging
    ChargingNeedsBattery
    Calculate_AnnualGeneration
    Calculate_FuelUse
    TotalCO2
    Totalcosts
end

struct HourSampling
    hours           ::Int           # unit: hours/period
    selectdays      ::Int
    skipdays        ::Int
    hoursperperiod  ::Float64       # not same as hours because cycles/year maybe not integer
    hourindexes     ::Vector{Int}   # select indexes before using hours to sample hours
end

HourSampling(options::Dict) = HourSampling(options[:hours], options[:selectdays], options[:skipdays])

function HourSampling(hours, selectdays, skipdays)
    @assert hours in (1,2,3,4,6) "hours must be one of (1,2,3,4,6)"
    hourspercycle = (selectdays + skipdays) * 24
    cyclesperyear = floor(8760/hourspercycle)
    periodsperyear = cyclesperyear * hourspercycle / hours
    hoursperperiod = 8760 / periodsperyear * (1 + skipdays/selectdays)
    hourindexes = [(c-1)*hourspercycle + h for c=1:cyclesperyear for h=1:hourspercycle if h <= selectdays*24]
    HourSampling(hours, selectdays, skipdays, hoursperperiod, hourindexes)
end

struct ModelInfo
    modelname   ::JuMP.Model
    sets        ::Sets
    params      ::Params
    vars        ::Vars
    constraints ::Constraints
    hourinfo    ::HourSampling
    options     ::Dict{Symbol,Any}
end

struct Results
    Status                      ::Symbol
    options                     ::Dict{Symbol,Any}
    hourinfo                    ::HourSampling
    sets                        ::Sets
    params                      ::Dict{Symbol,Any}
    Systemcost                  ::AxisArray{Float64,1}
    CO2emissions                ::AxisArray{Float64,1}
    FuelUse                     ::AxisArray{Float64,2}
    Electricity                 ::Dict{Tuple{Symbol,Symbol}, Array{Float64,2}}
    Charging                        ::AxisArray{Float64,3}
    StorageLevel                    ::Dict{Tuple{Symbol,Symbol}, Array{Float64,2}}
    Transmission                ::AxisArray{Float64,3}
    TransmissionCapacity        ::AxisArray{Float64,2}
    Capacity                    ::Dict{Tuple{Symbol,Symbol,Symbol}, Float64}
end

max_stars_repo_name: niclasmattsson/Supergrid

max_issues_repo_name: niclasmattsson/Supergrid

max_forks_repo_name: niclasmattsson/Supergrid


text:
module BernoulliFactory

using Random

## TODO: move this to NonUniformRandomVariateGeneration
@inline function _sampleGeometric(p::Float64, rng::RNG) where RNG<:AbstractRNG
  return floor(Int64, randexp(rng)/(-log(1.0-p)))+1
end

include("huber2016.jl")
include("huber2017.jl")
include("wastlundSqrt.jl")
include("mendoPower.jl")
include("algorithms.jl")
include("signedEstimate.jl")

end # module