text:
using AbstractPlotting.Keyboard
using AbstractPlotting.Mouse
using AbstractPlotting: ispressed, is_mouseinside

function axislines!(scene, rect)
    points = lift(rect) do r
        p1 = Point2(r.origin[1], r.origin[2] + r.widths[2])
        p2 = Point2(r.origin[1], r.origin[2])
        p3 = Point2(r.origin[1] + r.widths[1], r.origin[2])
        [p1, p2, p3]
    end
    lines!(scene, points, linewidth = 2, show_axis = false)
end

"""
A cheaper function that tries to come up with usable tick locations for a given value range
"""
function locateticks(xmin, xmax)
    # whats the distance?
    d = xmax - xmin
    # which order of magnitude is the distance in?
    ex = log10(d)
    # round the exponent to the closest one
    exrounded = round(ex)
    # this factor is used so we can find integer steps that then map back to
    # nice numbers in the given value range
    factor = 1 / 10 ^ (exrounded - 1)

    # minimum needs to be at least the lower value times the scaling factor
    xminf = ceil(xmin * factor)
    # maximum needs to be at most the higher value times the scaling factor
    xmaxf = floor(xmax * factor)

    # xminf and xmaxf are now both integers that are in an order of magnitude around ten steps apart
    df = xmaxf - xminf

    # step sizes we like
    steps = [5, 4, 2, 1]

    # from the highest to the lowest step size, choose the first that fits at least
    # two times between the end values (gives three ticks including the end value)
    for s in steps
        n, remainder = divrem(df, s)
        if n >= 2
            rang = 1:n
            # calculate the ticks by dividing with the factor from above
            ticks = [xminf; [xminf + x * s for x in rang]] ./ factor
            return ticks
        end
    end
end

# struct LimitCamera <: AbstractCamera end

function add_pan!(scene::SceneLike, limits)
    startpos = Base.RefValue((0.0, 0.0))
    pan = Mouse.right
    xzoom = Keyboard.x
    yzoom = Keyboard.y
    e = events(scene)
    on(
        camera(scene),
        # Node.((scene, cam, startpos))...,
        Node.((scene, startpos))...,
        e.mousedrag
    ) do scene, startpos, dragging
        # pan = cam.panbutton[]
        mp = e.mouseposition[]
        if ispressed(scene, pan) && is_mouseinside(scene)
            window_area = pixelarea(scene)[]
            if dragging == Mouse.down
                startpos[] = mp
            elseif dragging == Mouse.pressed && ispressed(scene, pan)
                diff = startpos[] .- mp
                startpos[] = mp
                pxa = scene.px_area[]
                diff_fraction = Vec2f0(diff) ./ Vec2f0(widths(pxa))

                diff_limits = diff_fraction .* widths(limits[])

                if ispressed(scene, xzoom)
                    limits[] = FRect(
                        Vec2f0(limits[].origin) .+ Vec2f0(diff_limits[1], 0), widths(limits[])
                    )
                elseif ispressed(scene, yzoom)
                    limits[] = FRect(
                        Vec2f0(limits[].origin) .+ Vec2f0(0, diff_limits[2]), widths(limits[])
                    )
                else
                    limits[] = FRect(Vec2f0(limits[].origin) .+ Vec2f0(diff_limits), widths(limits[]))
                end
            end
        end
        return
    end
end

function add_zoom!(scene::SceneLike, limits)

    e = events(scene)
    cam = camera(scene)
    on(cam, e.scroll) do x
        # @extractvalue cam (zoomspeed, zoombutton, area)
        zoomspeed = 0.10f0
        zoombutton = nothing
        zoom = Float32(x[2])
        if zoom != 0 && ispressed(scene, zoombutton) &&        AbstractPlotting.is_mouseinside(scene)
            pa = pixelarea(scene)[]

            # don't let z go negative
            z = max(0.1f0, 1f0 + (zoom * zoomspeed))

            # limits[] = FRect(limits[].origin..., (limits[].widths .* 0.99)...)
            mp_fraction = (Vec2f0(e.mouseposition[]) - minimum(pa)) ./ widths(pa)

            mp_data = limits[].origin .+ mp_fraction .* limits[].widths

            xorigin = limits[].origin[1]
            yorigin = limits[].origin[2]

            xwidth = limits[].widths[1]
            ywidth = limits[].widths[2]
            newxwidth = xwidth * z
            newywidth = ywidth * z

            newxorigin = xorigin + mp_fraction[1] * (xwidth - newxwidth)
            newyorigin = yorigin + mp_fraction[2] * (ywidth - newywidth)

            if AbstractPlotting.ispressed(scene, AbstractPlotting.Keyboard.x)
                limits[] = FRect(newxorigin, yorigin, newxwidth, ywidth)
            elseif AbstractPlotting.ispressed(scene, AbstractPlotting.Keyboard.y)
                limits[] = FRect(xorigin, newyorigin, xwidth, newywidth)
            else
                limits[] = FRect(newxorigin, newyorigin, newxwidth, newywidth)
            end
        end
        return
    end
end

function LayoutedAxis(parent::Scene)
    scene = Scene(parent, Node(IRect(0, 0, 100, 100)), raw = true)
    limits = Node(FRect(0, 0, 100, 100))
    xlabel = Node("x label")
    ylabel = Node("y label")

    add_pan!(scene, limits)
    add_zoom!(scene, limits)

    campixel!(scene)

    ticksnode = Node(Point2f0[])
    ticks = linesegments!(
        parent, ticksnode, linewidth = 2, show_axis = false
    )[end]

    # the algorithm from above seems to not give more than 7 ticks with the step sizes I chose
    nmaxticks = 7

    xticklabelnodes = [Node("0") for i in 1:nmaxticks]
    xticklabelposnodes = [Node(Point(0.0, 0.0)) for i in 1:nmaxticks]
    xticklabels = map(1:nmaxticks) do i
        text!(
            parent,
            xticklabelnodes[i],
            position = xticklabelposnodes[i],
            align = (:center, :top),
            textsize = 20,
            show_axis = false
        )[end]
    end

    yticklabelnodes = [Node("0") for i in 1:nmaxticks]
    yticklabelposnodes = [Node(Point(0.0, 0.0)) for i in 1:nmaxticks]
    yticklabels = map(1:nmaxticks) do i
        text!(
            parent,
            yticklabelnodes[i],
            position = yticklabelposnodes[i],
            align = (:center, :bottom),
            rotation = pi/2,
            textsize = 20,
            show_axis = false
        )[end]
    end

    on(camera(scene), pixelarea(scene), limits) do pxa, lims

        nearclip = -10_000f0
        farclip = 10_000f0

        limox, limoy = Float32.(lims.origin)
        limw, limh = Float32.(widths(lims))
        l, b = Float32.(pxa.origin)
        w, h = Float32.(widths(pxa))
        # projection = AbstractPlotting.orthographicprojection(0f0, w * 2f0, 0f0, h, nearclip, farclip)
        projection = AbstractPlotting.orthographicprojection(limox, limox + limw, limoy, limoy + limh, nearclip, farclip)
        camera(scene).projection[] = projection
        camera(scene).projectionview[] = projection

        # pxa = scene.px_area[]
        px_aspect = pxa.widths[1] / pxa.widths[2]

        # @printf("cam %.1f, %.1f, %.1f, %.1f\n", a.origin..., a.widths...)
        # @printf("pix %.1f, %.1f, %.1f, %.1f\n", pxa.origin..., pxa.widths...)

        width = lims.widths[1]
        # width = px_aspect > 1 ? a.widths[1] / px_aspect : a.widths[1]
        xrange = (lims.origin[1], lims.origin[1] + width)

        if width == 0 || !isfinite(xrange[1]) || !isfinite(xrange[2])
            return
        end

        xtickvals = locateticks(xrange...)
        # xtickvals, vminbest, vmaxbest = optimize_ticks(xrange...)

        # this code here tries to transform between values given in pixels of the
        # scene and the camera area, but this is incorrect right now and everything
        # should just be determined by the now unused limits that are saved in the
        # LayoutedAxis object
        xfractions = (xtickvals .- xrange[1]) ./ width
        xrange_scene = (pxa.origin[1], pxa.origin[1] + pxa.widths[1])
        width_scene = xrange_scene[2] - xrange_scene[1]
        xticks_scene = xrange_scene[1] .+ width_scene .* xfractions
        ticksize = 10 # px
        y = pxa.origin[2]
        xtickstarts = [Point(x, y) for x in xticks_scene]
        xtickends = [t + Point(0.0, -ticksize) for t in xtickstarts]

        # height = px_aspect < 1 ? a.widths[2] * px_aspect : a.widths[2]
        height = lims.widths[2]
        yrange = (lims.origin[2], lims.origin[2] + height)

        ytickvals = locateticks(yrange...)
        # ytickvals, vminbest, vmaxbest = optimize_ticks(yrange...)
        yfractions = (ytickvals .- yrange[1]) ./ height
        yrange_scene = (pxa.origin[2], pxa.origin[2] + pxa.widths[2])
        height_scene = yrange_scene[2] - yrange_scene[1]
        yticks_scene = yrange_scene[1] .+ height_scene .* yfractions
        ticksize = 10 # px
        x = pxa.origin[1]
        ytickstarts = [Point(x, y) for y in yticks_scene]
        ytickends = [t + Point(-ticksize, 0.0) for t in ytickstarts]

        # set and position tick labels
        xtickstrings = Showoff.showoff(xtickvals, :plain)
        nxticks = length(xtickvals)
        for i in 1:nmaxticks
            if i <= nxticks
                xticklabelnodes[i][] = xtickstrings[i]
                xticklabelposnodes[i][] = xtickends[i] + Point(0.0, -10.0)
                xticklabels[i].visible = true
            else
                xticklabels[i].visible = false
            end
        end

        ytickstrings = Showoff.showoff(ytickvals, :plain)
        nyticks = length(ytickvals)
        for i in 1:nmaxticks
            if i <= nyticks
                yticklabelnodes[i][] = ytickstrings[i]
                yticklabelposnodes[i][] = ytickends[i] + Point(-10.0, 0.0)
                yticklabels[i].visible = true
            else
                yticklabels[i].visible = false
            end
        end

        # set tick mark positions
        ticksnode[] = collect(Iterators.flatten(zip(
            [xtickstarts; ytickstarts],
            [xtickends; ytickends]
        )))
    end

    labelgap = 50

    xlabelpos = lift(scene.px_area) do a
        Point2(a.origin[1] + a.widths[1] / 2, a.origin[2] - labelgap)
    end

    ylabelpos = lift(scene.px_area) do a
        Point2(a.origin[1] - labelgap, a.origin[2] + a.widths[2] / 2)
    end

    tx = text!(
        parent, xlabel, textsize = 20, position = xlabelpos, show_axis = false
    )[end]
    tx.align = [0.5, 1]
    ty = text!(
        parent, ylabel, textsize = 20,
        position = ylabelpos, rotation = pi/2, show_axis = false
    )[end]
    ty.align = [0.5, 0]

    axislines!(parent, scene.px_area)

    LayoutedAxis(parent, scene, xlabel, ylabel, limits)
end

function applylayout(sg::SolvedGridLayout)
    for c in sg.content
        applylayout(c.al)
    end
end

function applylayout(sa::SolvedAxisLayout)
    sa.axis.scene.px_area[] = IRect2D(sa.inner)
end

function shrinkbymargin(rect, margin)
    IRect((rect.origin .+ margin), (rect.widths .- 2 .* margin))
end

function linkxaxes!(a::LayoutedAxis, b::LayoutedAxis)
    on(a.limits) do alim
        blim = b.limits[]

        ao = alim.origin[1]
        bo = blim.origin[1]
        aw = alim.widths[1]
        bw = blim.widths[1]

        if ao != bo || aw != bw
            b.limits[] = FRect(ao, blim.origin[2], aw, blim.widths[2])
        end
    end

    on(b.limits) do blim
        alim = a.limits[]

        ao = alim.origin[1]
        bo = blim.origin[1]
        aw = alim.widths[1]
        bw = blim.widths[1]

        if ao != bo || aw != bw
            a.limits[] = FRect(bo, alim.origin[2], bw, alim.widths[2])
        end
    end
end

function linkyaxes!(a::LayoutedAxis, b::LayoutedAxis)
    on(a.limits) do alim
        blim = b.limits[]

        ao = alim.origin[2]
        bo = blim.origin[2]
        aw = alim.widths[2]
        bw = blim.widths[2]

        if ao != bo || aw != bw
            b.limits[] = FRect(blim.origin[1], ao, blim.widths[1], aw)
        end
    end

    on(b.limits) do blim
        alim = a.limits[]

        ao = alim.origin[2]
        bo = blim.origin[2]
        aw = alim.widths[2]
        bw = blim.widths[2]

        if ao != bo || aw != bw
            a.limits[] = FRect(alim.origin[1], bo, alim.widths[1], bw)
        end
    end
end


text:
# Approaches v/norm(v) for large v, but approaches zero smoothly as v -> 0
function smooth_normalization(v, α=1.e3)
    h = norm(v) + 1/α * log(1 + exp(-2*α*norm(v)))
    return v/h
end

# Smooth approximation of Euclidean norm
function smooth_norm(x, α=1.e-3)
    sqrt(sum(x.^2) + α^2) - α
end

# 1D Barrier functions
# Finite exponential barrier with value specification
# Inputs:
#   x : Eval location
#   x_lim : Barrier location
#   x_ctl : Control location
#   f_lim : Value at x_lim
#   ctl_ratio : Control ratio (value at x_ctl = f_lim*ctl_ratio)
function exp_barrier_value(x, x_lim, x_ctl, f_lim, ctl_rat)
    σ = (x_ctl - x_lim)/log(1/ctl_rat)
    f_lim*exp(-(x - x_lim)/σ)
end

# Finite exponential barrier with slope specification
# Works for both min and max barrier
# Inputs:
#   x : Eval location
#   x_lim : Barrier location
#   x_ctl : Control location
#   fdot_lim : Slope at x_lim
#   ctl_ratio : Control ratio (slope at x_γ = fdot_lim*ctl_ratio)
function exp_barrier_slope(x, x_lim, x_ctl, fdot_lim, ctl_ratio)
    σ = (x_lim - x_ctl)/log(1/ctl_ratio)
    fdot_lim*σ*exp((x - x_lim)/σ)
end

# Softplus with slope specification
# Inputs:
# Works for both min and max barrier
#   x : Eval location
#   x_lim : Barrier location
#   x_ctl : Control location
#   fdot_lim : Slope at x_lim
#   ctl_ratio : Control ratio (slope at x_ctl = fdot_lim*ctl_ratio)
function softplus_slope(x, x_lim, x_ctl, fdot_lim, γ)
    σ = (x_ctl - x_lim)*log(1 - 2*ctl_ratio)
    2*γ*fdot_lim*σ*log(1 + exp((x - x_lim)/σ))
end

# Smooth approximation of max and min functions
smooth_max(x::AbstractVector, α=1.e2) = sum(@. x*exp(α*x))/sum(@. exp(α*x))
smooth_max(x, y, α=1.e2) = (x*exp(α*x) + y*exp(α*y))/(exp(α*x) + exp(α*y))

smooth_min(x::AbstractVector, α=1.e2) = -smooth_max(-x, α)
smooth_min(x, y, α=1.e2) = -smooth_max(-x, -y, α)

# Smooth approximation of clamp function
# Inputs:
#  x : Eval location
#  x_trans : Clamp transition start point
#  x_lim : Clamp transition ending point
#  slope_lim : Slope at x_lim
function smooth_clamp(x, x_trans, x_lim, slope_lim=1.e-3)
    if x < -x_trans
        c = 2. - 4. / slope_lim
        d = (-c + sqrt(c^2 - 4.))/2.
        k = -log(d)/(-x_lim + x_trans)
        L = 4. / k
        h = -x_trans - 0.5*L

        return L/(1+exp(-k*(x + x_trans))) + h

    elseif -x_trans <= x < x_trans
        return x

    else
        c = 2. - 4. / slope_lim
        d = (-c + sqrt(c^2 - 4.))/2.
        k = -log(d)/(x_lim - x_trans)
        L = 4. / k
        h = x_trans - 0.5*L

        return L/(1+exp(-k*(x - x_trans))) + h
    end
end

# Smooth absolute value function
function smooth_abs(x, ε=1.e-3)
    sqrt(x^2 + ε)
end

# Logistic function with value specification
# Inputs:
#  x : Eval location
#  xmin : Location of specified near-minimum value
#  xmax : Location of specified near-maximum value
#  fmin : Minimum value at limit
#  fmax : Maximum value at limit
#  ε : Gap from fmin/fmax at xmin/xmax
function logistic_val(x, xmin, xmax, fmin, fmax, ε)
    h = fmin
    L = fmax - fmin
    k = log(((fmin-fmax+ε)/ε)^2)/(xmax - xmin)
    x0 = 1. / k * log((fmax-fmin)/ε - 1.) + xmin

    L/(1. + exp(-k*(x-x0))) + h
end

# Angle wrapping functions
wrapTo2Pi(x) = mod(2*π + mod(x, 2*π), 2*π)
wrapToPi(x) = -π + wrapTo2Pi(x + π)

# Distance to line segment
# Inputs:
#  x : Eval location
#  p1 : Endpoint 1 of line segment
#  p2 : Endpoint 2 of line segment
function line_segment_distance(x, p1, p2)
    p2p1 = p2 - p1
    xp1 = x - p1
    xp2 = x - p2
    if dot(xp1, p2p1) > 0.
        if dot(xp2, p2p1) < 0.
            normal_vec = xp1 - p2p1*dot(xp1, p2p1)/norm(p2p1)^2
            return norm(normal_vec)
        else
            return norm(xp2)
        end
    else
        return smooth_norm(xp1)
    end
end

# Distance from point with positive x to box bordering y axis
# Inputs:
#  x : Eval location
#  p_lo : Low corner location
#  p_hi : High corner location
#  α : smoothing parameter
function aligned_box_distance(x, p_lo, p_hi, α=1.e2)
    if x[1] <= p_lo[1]
        if p_lo[2] <= x[2] <= p_hi[2]
            return 0.
        else
            return smooth_max(smooth_max(0., p_lo[2] - x[2], α), smooth_max(0., x[2] - p_hi[2], α), α)
        end
    else
        return line_segment_distance(x, p_lo, p_hi)
    end
end

## Plotting utiliies
function html_video(filename)
    base64_video = Base64.base64encode(open(read, filename))
    """<video controls src="data:video/x-m4v;base64,$base64_video">"""
end

ext: jl


text:
using DiffEqFlux, Flux, Zygote

Zygote.@nograd diagm
include("eig.jl")

""" Regularization terms as function of the parameters """

function no_regularization(ode, link, mod, feedthrough, p_ode, p_link, p_mod, p_feed)
    return 0.0
end

""" Regularize magnitude of Lur'e-type Lyapunov function parameters """
function lds_condition_V(ode, link, mod, feedthrough, p_ode, p_link, p_mod, p_feed)
    return mean(mod.get_Ω(p_mod))
end

""" Regularize eigenvalues of stability matrix """
function lds_condition_ctrnn0(ode, link, mod, feedthrough, p_ode, p_link, p_mod, p_feed)
    Ω = mod.get_Ω(p_mod)
    τ = exp.(view(p_ode, ode.index_dict["logτ"]))
    A = reshape(view(p_ode, ode.index_dict["A"]), (ode.state_dim, ode.state_dim))

    Ωmat = Ω .* diagm(ones(length(Ω)))
    mat = 0.5 * τ .* (Ω .* A)
    return softplus_eigreg(
        0.5 * τ .* Ωmat * A + 0.5 * τ .* A' * Ωmat - Ω .* diagm(ones(length(Ω))),
    )
end

""" Regularize max eigenvalue of stability matrix """
function lds_condition_ctrnn1(ode, link, mod, feedthrough, p_ode, p_link, p_mod, p_feed)
    Ω = mod.get_Ω(p_mod)
    τ = exp.(view(p_ode, ode.index_dict["logτ"]))
    W = reshape(view(p_ode, ode.index_dict["W"]), ode.state_dim, ode.hidden_dim)
    A = reshape(view(p_ode, ode.index_dict["A"]), ode.hidden_dim, ode.state_dim)

    Ωmat = Ω .* diagm(ones(length(Ω)))
    mat = 0.5 * τ .* (Ω .* A) * W
    return softplus_eigreg(mat + mat' - Ωmat)
end

max_stars_repo_name: syanga/ISSNeuralODE

max_issues_repo_name: syanga/ISSNeuralODE

max_forks_repo_name: syanga/ISSNeuralODE


text:
for N in [Float64, Rational{Int}, Float32]
    P = BallInf(zeros(N, 2), N(1//10))
    d = N[1//10, 1//10]
    answer = is_interior_point(d, P)
    @test answer == (N <: Rational)
    @test !is_interior_point(d, P; ε=N(1))
    @test is_interior_point(d, P; ε=N(0))

    d = N[1//10, 1//10] .- LazySets._rtol(N)
    @test is_interior_point(d, P)
end

# tests that do not work with Rational{Int}
for N in [Float64, Float32]
    P = BallInf(zeros(N, 2), N(1//10))
    d = N[1//10, 1//10]
    @test !is_interior_point(d, P; p=N(2))

    d = N[1//10, 1//10] .- LazySets._rtol(N)
    @test is_interior_point(d, P; p=N(2))
end


text:
# ModifiedBinaryLayer and ModifiedBinaryMERA types, and methods thereof.
# To be `included` in MERAKit.jl.

"""
    ModifiedBinaryLayer{ST, ET, Tan} <: SimpleLayer

The type for layers of a modified binary MERA.

Each layer consists of three tensors, a 2-to-2 disentangler, often called `u`, and two
2-to-1 isometries, often called `wl` and `wr`, for left and right. Their relative locations
are
```
│     │
wl   wr
│ ╲ ╱ │
│  u  │
│ ╱ ╲ │
```

The type parameters are `ST` for space type, e.g. `ComplexSpace` or `SU2Space`; `ET` for
element type, e.g. `Complex{Float64}`; and `Tan` for whether this layer is a tangent layer
or not.  If `Tan = false`, the layer in question is an actual MERA layer. If `Tan = true` it
consists, instead of the actual tensors, of Stiefel/Grassmann tangent vectors of these
tensors.

`ModifiedBinaryLayer` implements the iteration interface, returning first the disentangler,
then the left isometry, and finally the right isometry.

Index numbering convention is as follows, where the physical indices are at the bottom:
Disentangler:
```
 3│    4│
 ┌┴─────┴┐
 │   u   │
 └┬─────┬┘
 1│    2│
```

Isometries:
```
    3│
 ┌───┴───┐
 │ wl/wr │
 └┬─────┬┘
 1│    2│
```
"""
struct ModifiedBinaryLayer{ST, ET, Tan} <: SimpleLayer
    # Even though the types are here marked as any, they are restricted to be specific,
    # concrete types dependent on ST, ET and Tan, both in the constructor and in
    # getproperty.
    disentangler::Any
    isometry_left::Any
    isometry_right::Any

    function ModifiedBinaryLayer{ST, ET, Tan}(
        disentangler, isometry_left, isometry_right
    ) where {ST, ET, Tan}
        DisType = disentangler_type(ST, ET, Tan)
        IsoType = binaryisometry_type(ST, ET, Tan)
        disconv = convert(DisType, disentangler)::DisType
        isoleftconv = convert(IsoType, isometry_left)::IsoType
        isorightconv = convert(IsoType, isometry_right)::IsoType
        return new{ST, ET, Tan}(disconv, isoleftconv, isorightconv)
    end
end

function ModifiedBinaryLayer(
    disentangler::DisType, isometry_left::IsoType, isometry_right::IsoType
) where {ST, DisType <: AbstractTensorMap{ST, 2, 2}, IsoType <: AbstractTensorMap{ST, 2, 1}}
    ET = eltype(DisType)
    @assert eltype(IsoType) === ET
    return ModifiedBinaryLayer{ST, ET, false}(disentangler, isometry_left, isometry_right)
end

function ModifiedBinaryLayer(
    disentangler::DisTanType, isometry_left::IsoTanType, isometry_right::IsoTanType
) where {
    ST,
    DisType <: AbstractTensorMap{ST, 2, 2},
    IsoType <: AbstractTensorMap{ST, 2, 1},
    DisTanType <: Stiefel.StiefelTangent{DisType},
    IsoTanType <: Grassmann.GrassmannTangent{IsoType}
}
    ET = eltype(DisType)
    @assert eltype(IsoType) === ET
    return ModifiedBinaryLayer{ST, ET, true}(disentangler, isometry_left, isometry_right)
end

function Base.getproperty(
    l::ModifiedBinaryLayer{ST, ET, Tan}, sym::Symbol
) where {ST, ET, Tan}
    if sym === :disentangler
        T = disentangler_type(ST, ET, Tan)
    elseif sym === :isometry_left || sym === :isometry_right
        T = binaryisometry_type(ST, ET, Tan)
    else
        T = Any
    end
    return getfield(l, sym)::T
end

"""
    ModifiedBinaryMERA{N}

A modified binary MERA is a MERA consisting of `ModifiedBinaryLayer`s.
"""
ModifiedBinaryMERA{N} = GenericMERA{N, T, O} where {T <: ModifiedBinaryLayer, O}
layertype(::Type{ModifiedBinaryMERA}) = ModifiedBinaryLayer
#Base.show(io::IO, ::Type{ModifiedBinaryMERA}) = print(io, "ModifiedBinaryMERA")
#function Base.show(io::IO, ::Type{ModifiedBinaryMERA{N}}) where {N}
#    return print(io, "ModifiedBinaryMERA{($N)}")
#end

# Implement the iteration and indexing interfaces. Allows things like `u, wl, wr = layer`.
# See simplelayer.jl for details.
function _tuple(layer::ModifiedBinaryLayer)
    return (layer.disentangler, layer.isometry_left, layer.isometry_right)
end

function operatortype(::Type{ModifiedBinaryLayer{ST, ET, false}}) where {ST, ET}
    return ModifiedBinaryOp{tensormaptype(ST, 2, 2, ET)}
end
operatortype(::Type{ModifiedBinaryLayer{ST, ET, true}}) where {ST, ET} = Nothing

TensorKit.spacetype(::Type{ModifiedBinaryLayer{ST, ET, Tan}}) where {ST, ET, Tan} = ST

scalefactor(::Type{<:ModifiedBinaryLayer}) = 2

causal_cone_width(::Type{<:ModifiedBinaryLayer}) = 2

Base.eltype(::Type{ModifiedBinaryLayer{ST, ET, Tan}}) where {ST, ET, Tan} = ET

outputspace(layer::ModifiedBinaryLayer) = space(layer.disentangler, 1)
inputspace(layer::ModifiedBinaryLayer) = space(layer.isometry_left, 3)'
internalspace(layer::ModifiedBinaryLayer) = space(layer.isometry_right, 1)
internalspace(m::ModifiedBinaryMERA, depth) = internalspace(getlayer(m, depth))

function expand_inputspace(layer::ModifiedBinaryLayer, V_new)
    u, wl, wr = layer
    wl = pad_with_zeros_to(wl, 3 => V_new')
    wr = pad_with_zeros_to(wr, 3 => V_new')
    return ModifiedBinaryLayer(u, wl, wr)
end

function expand_outputspace(layer::ModifiedBinaryLayer, V_new)
    u, wl, wr = layer
    u = pad_with_zeros_to(u, 1 => V_new, 2 => V_new)
    wl = pad_with_zeros_to(wl, 1 => V_new)
    wr = pad_with_zeros_to(wr, 2 => V_new)
    return ModifiedBinaryLayer(u, wl, wr)
end

function expand_internalspace(layer::ModifiedBinaryLayer, V_new)
    u, wl, wr = layer
    u = pad_with_zeros_to(u, 3 => V_new', 4 => V_new')
    wl = pad_with_zeros_to(wl, 2 => V_new)
    wr = pad_with_zeros_to(wr, 1 => V_new)
    return ModifiedBinaryLayer(u, wl, wr)
end

function randomlayer(
    ::Type{ModifiedBinaryLayer}, ::Type{T}, Vin, Vout, Vint = Vout;
    random_disentangler = false
) where {T}
    wl = randomisometry(T, Vout ⊗ Vint, Vin)
    # We make the initial guess be reflection symmetric, since that's often true of the
    # desired MERA too (at least if random_disentangler is false, but we do it every time
    # any way).
    # TODO For anyons, should there be a twist on the top leg too, and does the braid make
    # sense?
    wr = braid(wl, (1, 2, 3), (2, 1), (3,); copy = true)
    u = initialize_disentangler(T, Vout, Vint, random_disentangler)
    return ModifiedBinaryLayer(u, wl, wr)
end

function ascending_fixedpoint(layer::ModifiedBinaryLayer)
    width = causal_cone_width(layer)
    Vtotal = ⊗(ntuple(n->inputspace(layer), Val(width))...)
    eye = id(storagetype(operatortype(layer)), Vtotal)
    return ModifiedBinaryOp(sqrt(8.0/5.0) * eye, sqrt(2.0/5.0) * eye)
end

function scalingoperator_initialguess(l::ModifiedBinaryLayer, irreps...)
    width = causal_cone_width(l)
    V = inputspace(l)
    interlayer_space = ⊗(ntuple(n->V, Val(width))...)
    outspace = interlayer_space
    inspace = interlayer_space
    for irrep in irreps
        inspace = inspace ⊗ spacetype(V)(irrep => 1)
    end
    typ = eltype(l)
    t = TensorMap(randn, typ, outspace ← inspace)
    return ModifiedBinaryOp(t)
end

function gradient(
    layer::ModifiedBinaryLayer, env::ModifiedBinaryLayer; metric = :euclidean
)
    u, wl, wr = layer
    uenv, wlenv, wrenv = env
    # The environment is the partial derivative. We need to turn that into a tangent vector
    # of the Stiefel manifold point u or w.
    # The factor of two is from the partial_x + i partial_y derivative of the cost function,
    # and how it depends on both v and v^dagger.
    ugrad = Stiefel.project!(2*uenv, u; metric = metric)
    wlgrad = Grassmann.project!(2*wlenv, wl)
    wrgrad = Grassmann.project!(2*wrenv, wr)
    return ModifiedBinaryLayer(ugrad, wlgrad, wrgrad)
end

function precondition_tangent(layer::ModifiedBinaryLayer, tan::ModifiedBinaryLayer, rho)
    u, wl, wr = layer
    utan, wltan, wrtan = tan
    @planar rho_wl_mid[-1; -11] := rho.mid[-1 1; -11 1]
    @planar rho_wl_gap[-1; -11] := rho.gap[1 -1; 1 -11]
    rho_wl = (rho_wl_mid + rho_wl_gap) / 2.0
    @planar rho_wr_mid[-1; -11] := rho.mid[1 -1; 1 -11]
    @planar rho_wr_gap[-1; -11] := rho.gap[-1 1; -11 1]
    rho_wr = (rho_wr_mid + rho_wr_gap) / 2.0
    @planar(
        rho_u[-1 -2; -11 -12] :=
        wl'[12; 1 -11] * wr'[22; -12 2] *
        rho.mid[11 21; 12 22] *
        wl[1 -1; 11] * wr[-2 2; 21]
    )
    utan_prec = precondition_tangent(utan, rho_u)
    wltan_prec = precondition_tangent(wltan, rho_wl)
    wrtan_prec = precondition_tangent(wrtan, rho_wr)
    return ModifiedBinaryLayer(utan_prec, wltan_prec, wrtan_prec)
end

# # # Invariants

function space_invar_intralayer(layer::ModifiedBinaryLayer)
    u, wl, wr = layer
    matching_bonds = (
        (space(u, 3)', space(wl, 2)),
        (space(u, 4)', space(wr, 1)),
    )
    allmatch = all(pair -> ==(pair...), matching_bonds)
    # Check that the dimensions are such that isometricity can hold.
    allmatch &= all((u, wl, wr)) do v
        codom, dom = fuse(codomain(v)), fuse(domain(v))
        infimum(dom, codom) == dom
    end
    return allmatch
end

function space_invar_interlayer(layer::ModifiedBinaryLayer, next_layer::ModifiedBinaryLayer)
    u, wl, wr = layer
    unext, wlnext, wrnext = next_layer
    matching_bonds = (
        (space(wl, 3)', space(unext, 1)),
        (space(wl, 3)', space(unext, 2)),
        (space(wr, 3)', space(unext, 1)),
        (space(wr, 3)', space(unext, 2)),
    )
    allmatch = all(pair -> ==(pair...), matching_bonds)
    return allmatch
end

# # # Ascending and descending superoperators

"""
Ascend a two-site `op` from the bottom of the given layer to the top.
"""
function ascend(op::ModifiedBinaryOp{T}, layer::ModifiedBinaryLayer) where {
    S1,
    T <: Union{SquareTensorMap{2},
    AbstractTensorMap{S1, 2, 3},
    AbstractTensorMap{S1, 2, 4}}
}
    l = ascend_left(op, layer)
    r = ascend_right(op, layer)
    m = ascend_mid(op, layer)
    b = ascend_between(op, layer)
    scaled_op_mid = (l+r+m) / 2
    scaled_op_gap = b / 2
    scaled_op = ModifiedBinaryOp(scaled_op_mid, scaled_op_gap)
    return scaled_op
end

function ascend(op::ModifiedBinaryOp{T}, layer::ModifiedBinaryLayer) where {
    T <: SquareTensorMap{1}
}
    op = expand_support(op, causal_cone_width(ModifiedBinaryLayer))
    return ascend(op, layer)
end

function ascend(op::AbstractTensorMap, layer::ModifiedBinaryLayer)
    op = ModifiedBinaryOp(op)
    return ascend(op, layer)
end

# TODO Think about how to best remove the code duplication of having the separate methods
# for ordinary, charged, and double charged operators.
# Operators without charge legs.
function ascend_left(op::ModifiedBinaryOp{T}, layer::ModifiedBinaryLayer) where {
    T <: SquareTensorMap{2}
}
    u, wl, wr = layer
    op_mid, op_gap = op
    # Cost: 2X^7 + 3X^6 + 1X^5
    @planar(
        scaled_op[-100 -200; -300 -400] :=
        wl[7 8; -300] * wr[10 1; -400] *
        u[6 5; 8 10] *
        op_gap[3 4; 7 6] *
        u'[2 9; 4 5] *
        wl'[-100; 3 2] * wr'[-200; 9 1]
    )
    return scaled_op
end

function ascend_right(op::ModifiedBinaryOp{T}, layer::ModifiedBinaryLayer) where {
    T <: SquareTensorMap{2}
}
    u, wl, wr = layer
    op_mid, op_gap = op
    # Cost: 2X^7 + 3X^6 + 1X^5
    @planar(
        scaled_op[-100 -200; -300 -400] :=
        wl[1 10; -300] * wr[8 7; -400] *
        u[5 6; 10 8] *
        op_gap[4 3; 6 7] *
        u'[9 2; 5 4] *
        wl'[-100; 1 9] * wr'[-200; 2 3]
    )
    return scaled_op
end

function ascend_mid(op::ModifiedBinaryOp{T}, layer::ModifiedBinaryLayer) where {
    T <: SquareTensorMap{2}
}
    u, wl, wr = layer
    op_mid, op_gap = op
    # Cost: 4X^6 + 2X^5
    @planar(
        scaled_op[-100 -200; -300 -400] :=
        wl[12 23; -300] * wr[21 11; -400] *
        u[3 4; 23 21] *
        op_mid[1 2; 3 4] *
        u'[24 22; 1 2] *
        wl'[-100; 12 24] * wr'[-200; 22 11]
    )
    return scaled_op
end

function ascend_between(op::ModifiedBinaryOp{T}, layer::ModifiedBinaryLayer) where {
    T <: SquareTensorMap{2}
}
    u, wl, wr = layer
    op_mid, op_gap = op
    # Cost: 2X^6 + 2X^5
    @planar(
        scaled_op[-100 -200; -300 -400] :=
        wr[12 23; -300] * wl[21 11; -400] *
        op_mid[24 22; 23 21] *
        wr'[-100; 12 24] * wl'[-200; 22 11]
    )
    return scaled_op
end

# Operators with one charge leg.
function ascend_left(op::ModifiedBinaryOp{T}, layer::ModifiedBinaryLayer) where {
    S1, T <: AbstractTensorMap{S1, 2, 3}
}
    u, wl, wr = layer
    op_mid, op_gap = op
    # Cost: 2X^7 + 3X^6 + 1X^5
    @tensor(
        scaled_op[-100 -200; -300 -400 -1000] :=
        wl[7 8; -300] * wr[10 1; -400] *
        u[6 5; 8 10] *
        op_gap[3 4; 7 6 -1000] *
        u'[2 9; 4 5] *
        wl'[-100; 3 2] * wr'[-200; 9 1]
    )
    return scaled_op
end

function ascend_right(op::ModifiedBinaryOp{T}, layer::ModifiedBinaryLayer) where {
    S1, T <: AbstractTensorMap{S1, 2, 3}
}
    u, wl, wr = layer
    op_mid, op_gap = op
    # Cost: 2X^7 + 3X^6 + 1X^5
    @tensor(
        scaled_op[-100 -200; -300 -400 -1000] :=
        wl[1 10; -300] * wr[8 7; -400] *
        u[5 6; 10 8] *
        op_gap[4 3; 6 7 -1000] *
        u'[9 2; 5 4] *
        wl'[-100; 1 9] * wr'[-200; 2 3]
    )
    return scaled_op
end

function ascend_mid(op::ModifiedBinaryOp{T}, layer::ModifiedBinaryLayer) where {
    S1, T <: AbstractTensorMap{S1, 2, 3}
}
    u, wl, wr = layer
    op_mid, op_gap = op
    # Cost: 4X^6 + 2X^5
    @tensor(
        scaled_op[-100 -200; -300 -400 -1000] :=
        wl[12 23; -300] * wr[21 11; -400] *
        u[3 4; 23 21] *
        op_mid[1 2; 3 4 -1000] *
        u'[24 22; 1 2] *
        wl'[-100; 12 24] * wr'[-200; 22 11]
    )
    return scaled_op
end

function ascend_between(op::ModifiedBinaryOp{T}, layer::ModifiedBinaryLayer) where {
    S1, T <: AbstractTensorMap{S1, 2, 3}
}
    u, wl, wr = layer
    op_mid, op_gap = op
    # Cost: 2X^6 + 2X^5
    @tensor(
        scaled_op[-100 -200; -300 -400 -1000] :=
        wr[12 23; -300] * wl[21 11; -400] *
        op_mid[24 22; 23 21 -1000] *
        wr'[-100; 12 24] * wl'[-200; 22 11]
    )
    return scaled_op
end

# Operators with two charge legs.
function ascend_left(op::ModifiedBinaryOp{T}, layer::ModifiedBinaryLayer) where {
    S1, T <: AbstractTensorMap{S1, 2, 4}
}
    u, wl, wr = layer
    op_mid, op_gap = op
    # Cost: 2X^7 + 3X^6 + 1X^5
    @planar(
        temp1[-1 -2 -3; -11 -12 -13 -14] :=
        op_gap[3 4; -1 -11 -12 -13] *
        u'[2 -3; 4 -14] *
        wl'[-2; 3 2]
    )
    temp2 = braid(temp1, (11, 12, 13, 14, 1, 100, 15), (1, 2, 5, 6, 3), (4, 7))
    @planar(
        temp3[-100 -1000 -2000 -200; -300 -400] :=
        wl[7 8; -300] * wr[10 1; -400] *
        u[6 5; 8 10] *
        temp2[7 -100 -1000 -2000 9; 6 5] *
        wr'[-200; 9 1]
    )
    scaled_op = braid(temp3, (11, 1, 100, 12, 13, 14), (1, 4), (5, 6, 2, 3))
    return scaled_op
end

function ascend_right(op::ModifiedBinaryOp{T}, layer::ModifiedBinaryLayer) where {
    S1, T <: AbstractTensorMap{S1, 2, 4}
}
    u, wl, wr = layer
    op_mid, op_gap = op
    # Cost: 2X^7 + 3X^6 + 1X^5
    # This happens to not need any braiding.
    @planar(
        scaled_op[-100 -200; -300 -400 -1000 -2000] :=
        wl[1 10; -300] * wr[8 7; -400] *
        u[5 6; 10 8] *
        op_gap[4 3; 6 7 -1000 -2000] *
        u'[9 2; 5 4] *
        wl'[-100; 1 9] * wr'[-200; 2 3]
    )
    return scaled_op
end

function ascend_mid(op::ModifiedBinaryOp{T}, layer::ModifiedBinaryLayer) where {
    S1, T <: AbstractTensorMap{S1, 2, 4}
}
    u, wl, wr = layer
    op_mid, op_gap = op
    # Cost: 4X^6 + 2X^5
    @planar(
        temp1[-1 -2; -3 -4 -5 -6] :=
        u[3 4; -3 -4] *
        op_mid[1 2; 3 4 -5 -6] *
        u'[-1 -2; 1 2]
    )
    temp2 = braid(temp1, (11, 12, 13, 14, 1, 100), (1, 2), (3, 5, 6, 4))
    @planar(
        temp3[-100 -200; -300 -1000 -2000 -400] :=
        wl[12 23; -300] * wr[21 11; -400] *
        temp2[24 22; 23 -1000 -2000 21] *
        wl'[-100; 12 24] * wr'[-200; 22 11]
    )
    scaled_op = braid(temp3, (11, 12, 13, 1, 100, 14), (1, 2), (3, 6, 4, 5))
    return scaled_op
end

function ascend_between(op::ModifiedBinaryOp{T}, layer::ModifiedBinaryLayer) where {
    S1, T <: AbstractTensorMap{S1, 2, 4}
}
    u, wl, wr = layer
    op_mid, op_gap = op
    # Cost: 2X^6 + 2X^5
    temp1 = braid(op_mid, (11, 12, 13, 14, 1, 100), (1, 2), (3, 5, 6, 4))
    @planar(
        temp2[-100 -200; -300 -1000 -2000 -400] :=
        wr[12 23; -300] * wl[21 11; -400] *
        temp1[24 22; 23 -1000 -2000 21] *
        wr'[-100; 12 24] * wl'[-200; 22 11]
    )
    scaled_op = braid(temp2, (11, 12, 13, 1, 100, 14), (1, 2), (3, 6, 4, 5))
    return scaled_op
end

function descend_left(rho::ModifiedBinaryOp, layer::ModifiedBinaryLayer)
    u, wl, wr = layer
    rho_mid, rho_gap = rho
    # Cost: 2X^7 + 3X^6 + 1X^5
    @planar(
        scaled_rho[-100 -200; -300 -400] :=
        u'[5 6; -400 7] *
        wl'[4; -300 5] * wr'[2; 6 1] *
        rho_mid[9 3; 4 2] *
        wl[-100 10; 9] * wr[8 1; 3] *
        u[-200 7; 10 8]
    )
    return scaled_rho
end

function descend_right(rho::ModifiedBinaryOp, layer::ModifiedBinaryLayer)
    u, wl, wr = layer
    rho_mid, rho_gap = rho
    # Cost: 2X^7 + 3X^6 + 1X^5
    @planar(
        scaled_rho[-100 -200; -300 -400] :=
        u'[6 5; 7 -300] *
        wl'[2; 1 6] * wr'[4; 5 -400] *
        rho_mid[3 9; 2 4] *
        wl[1 8; 3] * wr[10 -200; 9] *
        u[7 -100; 8 10]
    )
    return scaled_rho
end

function descend_mid(rho::ModifiedBinaryOp, layer::ModifiedBinaryLayer)
    u, wl, wr = layer
    rho_mid, rho_gap = rho
    # Cost: 4X^6 + 2X^5
    @planar(
        scaled_rho[-100 -200; -300 -400] :=
        u'[9 10; -300 -400] *
        wl'[5; 4 9] * wr'[2; 10 1] *
        rho_mid[6 3; 5 2] *
        wl[4 7; 6] * wr[8 1; 3] *
        u[-100 -200; 7 8]
    )
    return scaled_rho
end

function descend_between(rho::ModifiedBinaryOp, layer::ModifiedBinaryLayer)
    u, wl, wr = layer
    rho_mid, rho_gap = rho
    # Cost: 2X^6 + 2X^5
    @planar(
        scaled_rho[-100 -200; -300 -400] :=
        wr'[5; 4 -300] * wl'[2; -400 1] *
        rho_gap[6 3; 5 2] *
        wr[4 -100; 6] * wl[-200 1; 3]
    )
    return scaled_rho
end

function descend(rho::ModifiedBinaryOp, layer::ModifiedBinaryLayer)
    u, wl, wr = layer
    rho_mid, rho_gap = rho
    l = descend_left(rho, layer)
    r = descend_right(rho, layer)
    m = descend_mid(rho, layer)
    b = descend_between(rho, layer)
    scaled_rho_mid = (m + b) / 2.0
    scaled_rho_gap = (l + r) / 2.0
    scaled_rho = ModifiedBinaryOp(scaled_rho_mid, scaled_rho_gap)
    return scaled_rho
end

function descend(op::AbstractTensorMap, layer::ModifiedBinaryLayer)
    return descend(ModifiedBinaryOp(op), layer)
end

# # # Optimization

function environment(op, layer::ModifiedBinaryLayer, rho; vary_disentanglers = true)
    if vary_disentanglers
        env_u = environment_disentangler(op, layer, rho)
    else
        # The adjoint is just for type stability.
        env_u = zero(layer.disentangler')'
    end
    env_wl = environment_isometry_left(op, layer, rho)
    env_wr = environment_isometry_right(op, layer, rho)
    return ModifiedBinaryLayer(env_u, env_wl, env_wr)
end

function minimize_expectation_ev(
    layer::ModifiedBinaryLayer, env::ModifiedBinaryLayer; vary_disentanglers = true
)
    u = if vary_disentanglers
        projectisometric(env.disentangler; alg = Polar())
    else
        layer.disentangler
    end
    wl = projectisometric(env.isometry_left; alg = Polar())
    wr = projectisometric(env.isometry_right; alg = Polar())
    return ModifiedBinaryLayer(u, wl, wr)
end

function environment_disentangler(
    h::ModifiedBinaryOp, layer::ModifiedBinaryLayer, rho::ModifiedBinaryOp
)
    u, wl, wr = layer
    h_mid, h_gap = h
    rho_mid, rho_gap = rho
    # Cost: 2X^7 + 3X^6 + 1X^5
    @planar(
        envl[-100 -200; -300 -400] :=
        rho_mid[4 2; 10 3] *
        wl[7 5; 4] * wr[6 1; 2] *
        u[8 -200; 5 6] *
        h_gap[9 -100; 7 8] *
        wl'[10; 9 -300] * wr'[3; -400 1]
    )

    # Cost: 2X^7 + 3X^6 + 1X^5
    @planar(
        envr[-100 -200; -300 -400] :=
        rho_mid[2 4; 3 10] *
        wl[1 6; 2] * wr[5 7; 4] *
        u[-100 8; 6 5] *
        h_gap[-200 9; 8 7] *
        wl'[3; 1 -300] * wr'[10; -400 9]
    )

    # Cost: 4X^6 + 2X^5
    @planar(
        envm[-100 -200; -300 -400] :=
        rho_mid[6 4; 5 3] *
        wl[1 7; 6] * wr[8 2; 4] *
        u[9 10; 7 8] *
        h_mid[-100 -200; 9 10] *
        wl'[5; 1 -300] * wr'[3; -400 2]
    )

    env = (envl + envr + envm) / 4.0
    return env
end

function environment_disentangler(h::SquareTensorMap{2}, layer::ModifiedBinaryLayer, rho)
    h = ModifiedBinaryOp(expand_support(h, causal_cone_width(ModifiedBinaryLayer)))
    return environment_disentangler(h, layer, rho)
end

function environment_disentangler(h::SquareTensorMap{1}, layer::ModifiedBinaryLayer, rho)
    h = ModifiedBinaryOp(expand_support(h, causal_cone_width(ModifiedBinaryLayer)))
    return environment_disentangler(h, layer, rho)
end

function environment_isometry_left(h::ModifiedBinaryOp, layer, rho::ModifiedBinaryOp)
    u, wl, wr = layer
    h_mid, h_gap = h
    rho_mid, rho_gap = rho
    # Cost: 2X^7 + 3X^6 + 1X^5
    @planar(
        envl[-100 -200; -300] :=
        rho_mid[4 3; -300 2] *
        wl[7 5; 4] * wr[6 1; 3] *
        u[8 10; 5 6] *
        h_gap[-100 9; 7 8] *
        u'[-200 11; 9 10] *
        wr'[2; 11 1]
    )

    # Cost: 2X^7 + 3X^6 + 1X^5
    @planar(
        envr[-100 -200; -300] :=
        rho_mid[10 8; -300 9] *
        wl[-100 11; 10] * wr[1 2; 8] *
        u[5 3; 11 1] *
        h_gap[4 6; 3 2] *
        u'[-200 7; 5 4] *
        wr'[9; 7 6]
    )

    # Cost: 4X^6 + 2X^5
    @planar(
        envm[-100 -200; -300] :=
        rho_mid[10 9; -300 8] *
        wl[-100 11; 10] * wr[7 5; 9] *
        u[3 4; 11 7] *
        h_mid[1 2; 3 4] *
        u'[-200 6; 1 2] *
        wr'[8; 6 5]
    )

    # Cost: 2X^6 + 2X^5
    @planar(
        envb[-100 -200; -300] :=
        rho_gap[5 7; 4 -300] *
        wr[1 3; 5] * wl[6 -200; 7] *
        h_mid[2 -100; 3 6] *
        wr'[4; 1 2]
    )

    env = (envl + envr + envm + envb) / 4.0
    return env
end

function environment_isometry_left(
    h::SquareTensorMap{2}, layer::ModifiedBinaryLayer, rho::ModifiedBinaryOp
)
    h = ModifiedBinaryOp(expand_support(h, causal_cone_width(ModifiedBinaryLayer)))
    return environment_isometry_left(h, layer, rho)
end

function environment_isometry_left(
    h::SquareTensorMap{1}, layer::ModifiedBinaryLayer, rho::ModifiedBinaryOp
)
    h = ModifiedBinaryOp(expand_support(h, causal_cone_width(ModifiedBinaryLayer)))
    return environment_isometry_left(h, layer, rho)
end

function environment_isometry_right(h::ModifiedBinaryOp, layer, rho::ModifiedBinaryOp)
    u, wl, wr = layer
    h_mid, h_gap = h
    rho_mid, rho_gap = rho
    # Cost: 2X^7 + 3X^6 + 1X^5
    @planar(
        envl[-100 -200; -300] :=
        rho_mid[8 10; 9 -300] *
        wl[2 1; 8] * wr[11 -200; 10] *
        u[3 5; 1 11] *
        h_gap[6 4; 2 3] *
        u'[7 -100; 4 5] *
        wl'[9; 6 7]
    )

    # Cost: 2X^7 + 3X^6 + 1X^5
    @planar(
        envr[-100 -200; -300] :=
        rho_mid[3 4; 2 -300] *
        wl[1 6; 3] * wr[5 7; 4] *
        u[10 8; 6 5] *
        h_gap[9 -200; 8 7] *
        u'[11 -100; 10 9] *
        wl'[2; 1 11]
    )

    # Cost: 4X^6 + 2X^5
    @planar(
        envm[-100 -200; -300] :=
        rho_mid[9 10; 8 -300] *
        wl[5 7; 9] * wr[11 -200; 10] *
        u[4 3; 7 11] *
        h_mid[2 1; 4 3] *
        u'[6 -100; 2 1] *
        wl'[8; 5 6]
    )

    # Cost: 2X^6 + 2X^5
    @planar(
        envb[-100 -200; -300] :=
        rho_gap[7 5; -300 4] *
        wr[-100 6; 7] * wl[3 1; 5] *
        h_mid[-200 2; 6 3] *
        wl'[4; 2 1]
    )

    env = (envl + envr + envm + envb) / 4.0
    return env
end

function environment_isometry_right(
    h::SquareTensorMap{2}, layer::ModifiedBinaryLayer, rho::ModifiedBinaryOp
)
    h = ModifiedBinaryOp(expand_support(h, causal_cone_width(ModifiedBinaryLayer)))
    return environment_isometry_right(h, layer, rho)
end

function environment_isometry_right(
    h::SquareTensorMap{1}, layer::ModifiedBinaryLayer, rho::ModifiedBinaryOp
)
    h = ModifiedBinaryOp(expand_support(h, causal_cone_width(ModifiedBinaryLayer)))
    return environment_isometry_right(h, layer, rho)
end


text:
# mirror of jpeglib.h
module LibJpeg

import JpegTurbo_jll: libjpeg
using CEnum

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

# Generated from jpeglib.h using Clang.jl with:
#=
    using Clang, JpegTurbo_jll

    jconfig_h = joinpath(JpegTurbo_jll.artifact_dir, "include", "jconfig.h")
    jmorecfg_h = joinpath(JpegTurbo_jll.artifact_dir, "include", "jmorecfg.h")
    jpeglib_h = joinpath(JpegTurbo_jll.artifact_dir, "include", "jpeglib.h")
    wc = init(; headers = [jconfig_h, jmorecfg_h, jpeglib_h],
            output_file = joinpath(@__DIR__, "libjpeg_out.jl"),
            common_file = joinpath(@__DIR__, "libjpeg_common.jl"),
            header_wrapped = (root, current)->root == current,
            header_library = x->"libjpeg"
            )
    run(wc)
=#
# ... plus some manual cleanup.

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

const JSAMPLE = UInt8

# incomplete internal types, defined in jpegint.h, which
# we only need pointers to and hence can define as Cvoid
const jpeg_comp_master = Cvoid
const jpeg_c_main_controller = Cvoid
const jpeg_c_prep_controller = Cvoid
const jpeg_c_coef_controller = Cvoid
const jpeg_marker_writer = Cvoid
const jpeg_color_converter = Cvoid
const jpeg_downsampler = Cvoid
const jpeg_forward_dct = Cvoid
const jpeg_entropy_encoder = Cvoid
const jpeg_decomp_master = Cvoid
const jpeg_d_main_controller = Cvoid
const jpeg_d_coef_controller = Cvoid
const jpeg_d_post_controller = Cvoid
const jpeg_input_controller = Cvoid
const jpeg_marker_reader = Cvoid
const jpeg_entropy_decoder = Cvoid
const jpeg_inverse_dct = Cvoid
const jpeg_upsampler = Cvoid
const jpeg_color_deconverter = Cvoid
const jpeg_color_quantizer = Cvoid

######################################################################################################
# Automatically generated using Clang.jl

const JPEG_LIB_VERSION = 62
const LIBJPEG_TURBO_VERSION = ".1."
const LIBJPEG_TURBO_VERSION_NUMBER = 2001000
const BITS_IN_JSAMPLE = 8
const MAX_COMPONENTS = 10

# Skipping MacroDefinition: GETJOCTET ( value ) ( value )

const JPEG_MAX_DIMENSION = Int32(65500)

# Skipping MacroDefinition: METHODDEF ( type ) static type
# Skipping MacroDefinition: LOCAL ( type ) static type
# Skipping MacroDefinition: GLOBAL ( type ) type
# Skipping MacroDefinition: EXTERN ( type ) extern type
# Skipping MacroDefinition: JMETHOD ( type , methodname , arglist ) type ( * methodname ) arglist

const JCOEF = Int16
const JOCTET = Cuchar
const UINT8 = Cuchar
const UINT16 = UInt32
const INT16 = Int16
const INT32 = Int32
const JDIMENSION = UInt32
const boolean = Sys.iswindows() ? Cuchar : Cint
const DCTSIZE = 8
const DCTSIZE2 = 64
const NUM_QUANT_TBLS = 4
const NUM_HUFF_TBLS = 4
const NUM_ARITH_TBLS = 16
const MAX_COMPS_IN_SCAN = 4
const MAX_SAMP_FACTOR = 4
const C_MAX_BLOCKS_IN_MCU = 10
const D_MAX_BLOCKS_IN_MCU = 10
const JCS_EXTENSIONS = 1
const JCS_ALPHA_EXTENSIONS = 1

@cenum J_DCT_METHOD::UInt32 begin
    JDCT_ISLOW = 0
    JDCT_IFAST = 1
    JDCT_FLOAT = 2
end

const JDCT_DEFAULT = JDCT_ISLOW
const JDCT_FASTEST = JDCT_IFAST

# Skipping MacroDefinition: jpeg_common_fields struct jpeg_error_mgr * err ; /* Error handler module */ struct jpeg_memory_mgr * mem ; /* Memory manager module */ struct jpeg_progress_mgr * progress ; /* Progress monitor, or NULL if none */ void * client_data ; /* Available for use by application */ boolean is_decompressor ; /* So common code can tell which is which */ int global_state

const JMSG_LENGTH_MAX = 200
const JMSG_STR_PARM_MAX = 80
const JPOOL_PERMANENT = 0
const JPOOL_IMAGE = 1
const JPOOL_NUMPOOLS = 2

# Skipping MacroDefinition: JPP ( arglist ) arglist
# Skipping MacroDefinition: jpeg_create_compress ( cinfo ) jpeg_CreateCompress ( ( cinfo ) , JPEG_LIB_VERSION , ( size_t ) sizeof ( struct jpeg_compress_struct ) )
# Skipping MacroDefinition: jpeg_create_decompress ( cinfo ) jpeg_CreateDecompress ( ( cinfo ) , JPEG_LIB_VERSION , ( size_t ) sizeof ( struct jpeg_decompress_struct ) )

const JPEG_SUSPENDED = 0
const JPEG_HEADER_OK = 1
const JPEG_HEADER_TABLES_O<NLY = 2
const JPEG_REACHED_SOS = 1
const JPEG_REACHED_EOI = 2
const JPEG_ROW_COMPLETED = 3
const JPEG_SCAN_COMPLETED = 4
const JPEG_RST0 = 0xd0
const JPEG_EOI = 0xd9
const JPEG_APP0 = 0xe0
const JPEG_COM = 0xfe
const JSAMPROW = Ptr{JSAMPLE}
const JSAMPARRAY = Ptr{JSAMPROW}
const JSAMPIMAGE = Ptr{JSAMPARRAY}
const JBLOCK = NTuple{64, JCOEF}
const JBLOCKROW = Ptr{JBLOCK}
const JBLOCKARRAY = Ptr{JBLOCKROW}
const JBLOCKIMAGE = Ptr{JBLOCKARRAY}
const JCOEFPTR = Ptr{JCOEF}

struct JQUANT_TBL
    quantval::NTuple{64, UINT16}
    sent_table::boolean
end

struct JHUFF_TBL
    bits::NTuple{17, UINT8}
    huffval::NTuple{256, UINT8}
    sent_table::boolean
end

struct jpeg_component_info
    component_id::Cint
    component_index::Cint
    h_samp_factor::Cint
    v_samp_factor::Cint
    quant_tbl_no::Cint
    dc_tbl_no::Cint
    ac_tbl_no::Cint
    width_in_blocks::JDIMENSION
    height_in_blocks::JDIMENSION
    DCT_scaled_size::Cint
    downsampled_width::JDIMENSION
    downsampled_height::JDIMENSION
    component_needed::boolean
    MCU_width::Cint
    MCU_height::Cint
    MCU_blocks::Cint
    MCU_sample_width::Cint
    last_col_width::Cint
    last_row_height::Cint
    quant_table::Ptr{JQUANT_TBL}
    dct_table::Ptr{Cvoid}
end

struct jpeg_scan_info
    comps_in_scan::Cint
    component_index::NTuple{4, Cint}
    Ss::Cint
    Se::Cint
    Ah::Cint
    Al::Cint
end

struct jpeg_marker_struct
    next::Ptr{jpeg_marker_struct}
    marker::UINT8
    original_length::UInt32
    data_length::UInt32
    data::Ptr{JOCTET}
end

const jpeg_saved_marker_ptr = Ptr{jpeg_marker_struct}

@cenum J_COLOR_SPACE::UInt32 begin
    JCS_UNKNOWN = 0
    JCS_GRAYSCALE = 1
    JCS_RGB = 2
    JCS_YCbCr = 3
    JCS_CMYK = 4
    JCS_YCCK = 5
    JCS_EXT_RGB = 6
    JCS_EXT_RGBX = 7
    JCS_EXT_BGR = 8
    JCS_EXT_BGRX = 9
    JCS_EXT_XBGR = 10
    JCS_EXT_XRGB = 11
    JCS_EXT_RGBA = 12
    JCS_EXT_BGRA = 13
    JCS_EXT_ABGR = 14
    JCS_EXT_ARGB = 15
    JCS_RGB565 = 16
end

@cenum J_DITHER_MODE::UInt32 begin
    JDITHER_NONE = 0
    JDITHER_ORDERED = 1
    JDITHER_FS = 2
end

struct ANONYMOUS1_msg_parm
    s::NTuple{80, UInt8}
end

struct jpeg_error_mgr
    error_exit::Ptr{Cvoid}
    emit_message::Ptr{Cvoid}
    output_message::Ptr{Cvoid}
    format_message::Ptr{Cvoid}
    reset_error_mgr::Ptr{Cvoid}
    msg_code::Cint
    msg_parm::ANONYMOUS1_msg_parm
    trace_level::Cint
    num_warnings::Clong
    jpeg_message_table::Ptr{Cstring}
    last_jpeg_message::Cint
    addon_message_table::Ptr{Cstring}
    first_addon_message::Cint
    last_addon_message::Cint
end

struct jpeg_memory_mgr
    alloc_small::Ptr{Cvoid}
    alloc_large::Ptr{Cvoid}
    alloc_sarray::Ptr{Cvoid}
    alloc_barray::Ptr{Cvoid}
    request_virt_sarray::Ptr{Cvoid}
    request_virt_barray::Ptr{Cvoid}
    realize_virt_arrays::Ptr{Cvoid}
    access_virt_sarray::Ptr{Cvoid}
    access_virt_barray::Ptr{Cvoid}
    free_pool::Ptr{Cvoid}
    self_destruct::Ptr{Cvoid}
    max_memory_to_use::Clong
    max_alloc_chunk::Clong
end

struct jpeg_progress_mgr
    progress_monitor::Ptr{Cvoid}
    pass_counter::Clong
    pass_limit::Clong
    completed_passes::Cint
    total_passes::Cint
end

struct jpeg_common_struct
    err::Ptr{jpeg_error_mgr}
    mem::Ptr{jpeg_memory_mgr}
    progress::Ptr{jpeg_progress_mgr}
    client_data::Ptr{Cvoid}
    is_decompressor::boolean
    global_state::Cint
end

const j_common_ptr = Ptr{jpeg_common_struct}

struct jpeg_destination_mgr
    next_output_byte::Ptr{JOCTET}
    free_in_buffer::Cint
    init_destination::Ptr{Cvoid}
    empty_output_buffer::Ptr{Cvoid}
    term_destination::Ptr{Cvoid}
end

Base.@kwdef mutable struct jpeg_compress_struct
    err::Ptr{jpeg_error_mgr} = C_NULL
    mem::Ptr{jpeg_memory_mgr} = C_NULL
    progress::Ptr{jpeg_progress_mgr} = C_NULL
    client_data::Ptr{Cvoid} = C_NULL
    is_decompressor::boolean = false
    global_state::Cint = 0
    dest::Ptr{jpeg_destination_mgr} = C_NULL
    image_width::JDIMENSION = 0
    image_height::JDIMENSION = 0
    input_components::Cint = 0
    in_color_space::J_COLOR_SPACE = JCS_UNKNOWN
    input_gamma::Cdouble = 0
    data_precision::Cint = 0
    num_components::Cint = 0
    jpeg_color_space::J_COLOR_SPACE = JCS_UNKNOWN
    comp_info::Ptr{jpeg_component_info} = C_NULL
    quant_tbl_ptrs::NTuple{4, Ptr{JQUANT_TBL}} = (C_NULL,C_NULL,C_NULL,C_NULL)
    dc_huff_tbl_ptrs::NTuple{4, Ptr{JHUFF_TBL}} = (C_NULL,C_NULL,C_NULL,C_NULL)
    ac_huff_tbl_ptrs::NTuple{4, Ptr{JHUFF_TBL}} = (C_NULL,C_NULL,C_NULL,C_NULL)
    arith_dc_L::NTuple{16, UINT8} = ntuple(i -> 0x00, Val(16))
    arith_dc_U::NTuple{16, UINT8} = ntuple(i -> 0x00, Val(16))
    arith_ac_K::NTuple{16, UINT8} = ntuple(i -> 0x00, Val(16))
    num_scans::Cint = 0
    scan_info::Ptr{jpeg_scan_info} = C_NULL
    raw_data_in::boolean = false
    arith_code::boolean = false
    optimize_coding::boolean = false
    CCIR601_sampling::boolean = false
    smoothing_factor::Cint = 0
    dct_method::J_DCT_METHOD = JDCT_DEFAULT
    restart_interval::UInt32 = 0
    restart_in_rows::Cint = 0
    write_JFIF_header::boolean = false
    JFIF_major_version::UINT8 = 0
    JFIF_minor_version::UINT8 = 0
    density_unit::UINT8 = 0
    X_density::UINT16 = 0
    Y_density::UINT16 = 0
    write_Adobe_marker::boolean = false
    next_scanline::JDIMENSION = 0
    progressive_mode::boolean = false
    max_h_samp_factor::Cint = 0
    max_v_samp_factor::Cint = 0
    total_iMCU_rows::JDIMENSION = 0
    comps_in_scan::Cint = 0
    cur_comp_info::NTuple{4, Ptr{jpeg_component_info}} = (C_NULL,C_NULL,C_NULL,C_NULL)
    MCUs_per_row::JDIMENSION = 0
    MCU_rows_in_scan::JDIMENSION = 0
    blocks_in_MCU::Cint = 0
    MCU_membership::NTuple{10, Cint} = ntuple(i -> Cint(0), Val(10))
    Ss::Cint = 0
    Se::Cint = 0
    Ah::Cint = 0
    Al::Cint = 0
    master::Ptr{jpeg_comp_master} = C_NULL
    main::Ptr{jpeg_c_main_controller} = C_NULL
    prep::Ptr{jpeg_c_prep_controller} = C_NULL
    coef::Ptr{jpeg_c_coef_controller} = C_NULL
    marker::Ptr{jpeg_marker_writer} = C_NULL
    cconvert::Ptr{jpeg_color_converter} = C_NULL
    downsample::Ptr{jpeg_downsampler} = C_NULL
    fdct::Ptr{jpeg_forward_dct} = C_NULL
    entropy::Ptr{jpeg_entropy_encoder} = C_NULL
    script_space::Ptr{jpeg_scan_info} = C_NULL
    script_space_size::Cint = 0
end

const j_compress_ptr = Ref{jpeg_compress_struct}

struct jpeg_source_mgr
    next_input_byte::Ptr{JOCTET}
    bytes_in_buffer::Cint
    init_source::Ptr{Cvoid}
    fill_input_buffer::Ptr{Cvoid}
    skip_input_data::Ptr{Cvoid}
    resync_to_restart::Ptr{Cvoid}
    term_source::Ptr{Cvoid}
end

Base.@kwdef mutable struct jpeg_decompress_struct
    err::Ptr{jpeg_error_mgr} = C_NULL
    mem::Ptr{jpeg_memory_mgr} = C_NULL
    progress::Ptr{jpeg_progress_mgr} = C_NULL
    client_data::Ptr{Cvoid} = C_NULL
    is_decompressor::boolean = false
    global_state::Cint = 0
    src::Ptr{jpeg_source_mgr} = C_NULL
    image_width::JDIMENSION = 0
    image_height::JDIMENSION = 0
    num_components::Cint = 0
    jpeg_color_space::J_COLOR_SPACE = JCS_UNKNOWN
    out_color_space::J_COLOR_SPACE = JCS_UNKNOWN
    scale_num::UInt32 = 0
    scale_denom::UInt32 = 0
    output_gamma::Cdouble = 0
    buffered_image::boolean = false
    raw_data_out::boolean = false
    dct_method::J_DCT_METHOD = JDCT_DEFAULT
    do_fancy_upsampling::boolean = false
    do_block_smoothing::boolean = false
    quantize_colors::boolean = false
    dither_mode::J_DITHER_MODE = JDITHER_NONE
    two_pass_quantize::boolean = false
    desired_number_of_colors::Cint = 0
    enable_1pass_quant::boolean = false
    enable_external_quant::boolean = false
    enable_2pass_quant::boolean = false
    output_width::JDIMENSION = 0
    output_height::JDIMENSION = 0
    out_color_components::Cint = 0
    output_components::Cint = 0
    rec_outbuf_height::Cint = 0
    actual_number_of_colors::Cint = 0
    colormap::JSAMPARRAY = C_NULL
    output_scanline::JDIMENSION = 0
    input_scan_number::Cint = 0
    input_iMCU_row::JDIMENSION = 0
    output_scan_number::Cint = 0
    output_iMCU_row::JDIMENSION = 0
    coef_bits::Ptr{NTuple{64, Cint}} = C_NULL
    quant_tbl_ptrs::NTuple{4, Ptr{JQUANT_TBL}} = (C_NULL,C_NULL,C_NULL,C_NULL)
    dc_huff_tbl_ptrs::NTuple{4, Ptr{JHUFF_TBL}} = (C_NULL,C_NULL,C_NULL,C_NULL)
    ac_huff_tbl_ptrs::NTuple{4, Ptr{JHUFF_TBL}} = (C_NULL,C_NULL,C_NULL,C_NULL)
    data_precision::Cint = 0
    comp_info::Ptr{jpeg_component_info} = C_NULL
    progressive_mode::boolean = false
    arith_code::boolean = false
    arith_dc_L::NTuple{16, UINT8} = ntuple(i -> 0x00, Val(16))
    arith_dc_U::NTuple{16, UINT8} = ntuple(i -> 0x00, Val(16))
    arith_ac_K::NTuple{16, UINT8} = ntuple(i -> 0x00, Val(16))
    restart_interval::UInt32 = 0
    saw_JFIF_marker::boolean = false
    JFIF_major_version::UINT8 = 0
    JFIF_minor_version::UINT8 = 0
    density_unit::UINT8 = 0
    X_density::UINT16 = 0
    Y_density::UINT16 = 0
    saw_Adobe_marker::boolean = false
    Adobe_transform::UINT8 = 0
    CCIR601_sampling::boolean = false
    marker_list::jpeg_saved_marker_ptr = C_NULL
    max_h_samp_factor::Cint = 0
    max_v_samp_factor::Cint = 0
    min_DCT_scaled_size::Cint = 0
    total_iMCU_rows::JDIMENSION = 0
    sample_range_limit::Ptr{JSAMPLE} = C_NULL
    comps_in_scan::Cint = 0
    cur_comp_info::NTuple{4, Ptr{jpeg_component_info}} = (C_NULL,C_NULL,C_NULL,C_NULL)
    MCUs_per_row::JDIMENSION = 0
    MCU_rows_in_scan::JDIMENSION = 0
    blocks_in_MCU::Cint = 0
    MCU_membership::NTuple{10, Cint} = ntuple(i -> Cint(0), Val(10))
    Ss::Cint = 0
    Se::Cint = 0
    Ah::Cint = 0
    Al::Cint = 0
    unread_marker::Cint = 0
    master::Ptr{jpeg_decomp_master} = C_NULL
    main::Ptr{jpeg_d_main_controller} = C_NULL
    coef::Ptr{jpeg_d_coef_controller} = C_NULL
    post::Ptr{jpeg_d_post_controller} = C_NULL
    inputctl::Ptr{jpeg_input_controller} = C_NULL
    marker::Ptr{jpeg_marker_reader} = C_NULL
    entropy::Ptr{jpeg_entropy_decoder} = C_NULL
    idct::Ptr{jpeg_inverse_dct} = C_NULL
    upsample::Ptr{jpeg_upsampler} = C_NULL
    cconvert::Ptr{jpeg_color_deconverter} = C_NULL
    cquantize::Ptr{jpeg_color_quantizer} = C_NULL
end

const j_decompress_ptr = Ref{jpeg_decompress_struct}
const jvirt_sarray_control = Cvoid
const jvirt_sarray_ptr = Ptr{jvirt_sarray_control}
const jvirt_barray_control = Cvoid
const jvirt_barray_ptr = Ptr{jvirt_barray_control}
const jpeg_marker_parser_method = Ptr{Cvoid}
# Julia wrapper for header: jconfig.h
# Automatically generated using Clang.jl

# Julia wrapper for header: jmorecfg.h
# Automatically generated using Clang.jl

# Julia wrapper for header: jpeglib.h
# Automatically generated using Clang.jl

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

function jpeg_std_error(err)
    ccall((:jpeg_std_error, libjpeg), Ptr{jpeg_error_mgr}, (Ptr{jpeg_error_mgr},), err)
end

function jpeg_CreateCompress(cinfo, version, structsize)
    ccall((:jpeg_CreateCompress, libjpeg), Cvoid, (j_compress_ptr, Cint, Cint), cinfo, version, structsize)
end

function jpeg_CreateDecompress(cinfo, version, structsize)
    ccall((:jpeg_CreateDecompress, libjpeg), Cvoid, (j_decompress_ptr, Cint, Cint), cinfo, version, structsize)
end

function jpeg_destroy_compress(cinfo)
    ccall((:jpeg_destroy_compress, libjpeg), Cvoid, (j_compress_ptr,), cinfo)
end

function jpeg_destroy_decompress(cinfo)
    ccall((:jpeg_destroy_decompress, libjpeg), Cvoid, (j_decompress_ptr,), cinfo)
end

function jpeg_stdio_dest(cinfo, outfile)
    ccall((:jpeg_stdio_dest, libjpeg), Cvoid, (j_compress_ptr, Ptr{Libc.FILE}), cinfo, outfile)
end

function jpeg_stdio_src(cinfo, infile)
    ccall((:jpeg_stdio_src, libjpeg), Cvoid, (j_decompress_ptr, Ptr{Libc.FILE}), cinfo, infile)
end

function jpeg_mem_dest(cinfo, outbuffer, outsize)
    ccall((:jpeg_mem_dest, libjpeg), Cvoid, (j_compress_ptr, Ptr{Ptr{Cuchar}}, Ptr{Culong}), cinfo, outbuffer, outsize)
end

function jpeg_mem_src(cinfo, inbuffer, insize)
    ccall((:jpeg_mem_src, libjpeg), Cvoid, (j_decompress_ptr, Ptr{Cuchar}, Culong), cinfo, inbuffer, insize)
end

function jpeg_set_defaults(cinfo)
    ccall((:jpeg_set_defaults, libjpeg), Cvoid, (j_compress_ptr,), cinfo)
end

function jpeg_set_colorspace(cinfo, colorspace)
    ccall((:jpeg_set_colorspace, libjpeg), Cvoid, (j_compress_ptr, J_COLOR_SPACE), cinfo, colorspace)
end

function jpeg_default_colorspace(cinfo)
    ccall((:jpeg_default_colorspace, libjpeg), Cvoid, (j_compress_ptr,), cinfo)
end

function jpeg_set_quality(cinfo, quality, force_baseline)
    ccall((:jpeg_set_quality, libjpeg), Cvoid, (j_compress_ptr, Cint, boolean), cinfo, quality, force_baseline)
end

function jpeg_set_linear_quality(cinfo, scale_factor, force_baseline)
    ccall((:jpeg_set_linear_quality, libjpeg), Cvoid, (j_compress_ptr, Cint, boolean), cinfo, scale_factor, force_baseline)
end

function jpeg_add_quant_table(cinfo, which_tbl, basic_table, scale_factor, force_baseline)
    ccall((:jpeg_add_quant_table, libjpeg), Cvoid, (j_compress_ptr, Cint, Ptr{UInt32}, Cint, boolean), cinfo, which_tbl, basic_table, scale_factor, force_baseline)
end

function jpeg_quality_scaling(quality)
    ccall((:jpeg_quality_scaling, libjpeg), Cint, (Cint,), quality)
end

function jpeg_simple_progression(cinfo)
    ccall((:jpeg_simple_progression, libjpeg), Cvoid, (j_compress_ptr,), cinfo)
end

function jpeg_suppress_tables(cinfo, suppress)
    ccall((:jpeg_suppress_tables, libjpeg), Cvoid, (j_compress_ptr, boolean), cinfo, suppress)
end

function jpeg_alloc_quant_table(cinfo)
    ccall((:jpeg_alloc_quant_table, libjpeg), Ptr{JQUANT_TBL}, (j_common_ptr,), cinfo)
end

function jpeg_alloc_huff_table(cinfo)
    ccall((:jpeg_alloc_huff_table, libjpeg), Ptr{JHUFF_TBL}, (j_common_ptr,), cinfo)
end

function jpeg_start_compress(cinfo, write_all_tables)
    ccall((:jpeg_start_compress, libjpeg), Cvoid, (j_compress_ptr, boolean), cinfo, write_all_tables)
end

function jpeg_write_scanlines(cinfo, scanlines, num_lines)
    ccall((:jpeg_write_scanlines, libjpeg), JDIMENSION, (j_compress_ptr, JSAMPARRAY, JDIMENSION), cinfo, scanlines, num_lines)
end

function jpeg_finish_compress(cinfo)
    ccall((:jpeg_finish_compress, libjpeg), Cvoid, (j_compress_ptr,), cinfo)
end

function jpeg_write_raw_data(cinfo, data, num_lines)
    ccall((:jpeg_write_raw_data, libjpeg), JDIMENSION, (j_compress_ptr, JSAMPIMAGE, JDIMENSION), cinfo, data, num_lines)
end

function jpeg_write_marker(cinfo, marker, dataptr, datalen)
    ccall((:jpeg_write_marker, libjpeg), Cvoid, (j_compress_ptr, Cint, Ptr{JOCTET}, UInt32), cinfo, marker, dataptr, datalen)
end

function jpeg_write_m_header(cinfo, marker, datalen)
    ccall((:jpeg_write_m_header, libjpeg), Cvoid, (j_compress_ptr, Cint, UInt32), cinfo, marker, datalen)
end

function jpeg_write_m_byte(cinfo, val)
    ccall((:jpeg_write_m_byte, libjpeg), Cvoid, (j_compress_ptr, Cint), cinfo, val)
end

function jpeg_write_tables(cinfo)
    ccall((:jpeg_write_tables, libjpeg), Cvoid, (j_compress_ptr,), cinfo)
end

function jpeg_write_icc_profile(cinfo, icc_data_ptr, icc_data_len)
    ccall((:jpeg_write_icc_profile, libjpeg), Cvoid, (j_compress_ptr, Ptr{JOCTET}, UInt32), cinfo, icc_data_ptr, icc_data_len)
end

function jpeg_read_header(cinfo, require_image)
    ccall((:jpeg_read_header, libjpeg), Cint, (j_decompress_ptr, boolean), cinfo, require_image)
end

function jpeg_start_decompress(cinfo)
    ccall((:jpeg_start_decompress, libjpeg), boolean, (j_decompress_ptr,), cinfo)
end

function jpeg_read_scanlines(cinfo, scanlines, max_lines)
    ccall((:jpeg_read_scanlines, libjpeg), JDIMENSION, (j_decompress_ptr, JSAMPARRAY, JDIMENSION), cinfo, scanlines, max_lines)
end

function jpeg_skip_scanlines(cinfo, num_lines)
    ccall((:jpeg_skip_scanlines, libjpeg), JDIMENSION, (j_decompress_ptr, JDIMENSION), cinfo, num_lines)
end

function jpeg_crop_scanline(cinfo, xoffset, width)
    ccall((:jpeg_crop_scanline, libjpeg), Cvoid, (j_decompress_ptr, Ptr{JDIMENSION}, Ptr{JDIMENSION}), cinfo, xoffset, width)
end

function jpeg_finish_decompress(cinfo)
    ccall((:jpeg_finish_decompress, libjpeg), boolean, (j_decompress_ptr,), cinfo)
end

function jpeg_read_raw_data(cinfo, data, max_lines)
    ccall((:jpeg_read_raw_data, libjpeg), JDIMENSION, (j_decompress_ptr, JSAMPIMAGE, JDIMENSION), cinfo, data, max_lines)
end

function jpeg_has_multiple_scans(cinfo)
    ccall((:jpeg_has_multiple_scans, libjpeg), boolean, (j_decompress_ptr,), cinfo)
end

function jpeg_start_output(cinfo, scan_number)
    ccall((:jpeg_start_output, libjpeg), boolean, (j_decompress_ptr, Cint), cinfo, scan_number)
end

function jpeg_finish_output(cinfo)
    ccall((:jpeg_finish_output, libjpeg), boolean, (j_decompress_ptr,), cinfo)
end

function jpeg_input_complete(cinfo)
    ccall((:jpeg_input_complete, libjpeg), boolean, (j_decompress_ptr,), cinfo)
end

function jpeg_new_colormap(cinfo)
    ccall((:jpeg_new_colormap, libjpeg), Cvoid, (j_decompress_ptr,), cinfo)
end

function jpeg_consume_input(cinfo)
    ccall((:jpeg_consume_input, libjpeg), Cint, (j_decompress_ptr,), cinfo)
end

function jpeg_calc_output_dimensions(cinfo)
    ccall((:jpeg_calc_output_dimensions, libjpeg), Cvoid, (j_decompress_ptr,), cinfo)
end

function jpeg_save_markers(cinfo, marker_code, length_limit)
    ccall((:jpeg_save_markers, libjpeg), Cvoid, (j_decompress_ptr, Cint, UInt32), cinfo, marker_code, length_limit)
end

function jpeg_set_marker_processor(cinfo, marker_code, routine)
    ccall((:jpeg_set_marker_processor, libjpeg), Cvoid, (j_decompress_ptr, Cint, jpeg_marker_parser_method), cinfo, marker_code, routine)
end

function jpeg_read_coefficients(cinfo)
    ccall((:jpeg_read_coefficients, libjpeg), Ptr{jvirt_barray_ptr}, (j_decompress_ptr,), cinfo)
end

function jpeg_write_coefficients(cinfo, coef_arrays)
    ccall((:jpeg_write_coefficients, libjpeg), Cvoid, (j_compress_ptr, Ptr{jvirt_barray_ptr}), cinfo, coef_arrays)
end

function jpeg_copy_critical_parameters(srcinfo, dstinfo)
    ccall((:jpeg_copy_critical_parameters, libjpeg), Cvoid, (j_decompress_ptr, j_compress_ptr), srcinfo, dstinfo)
end

function jpeg_abort_compress(cinfo)
    ccall((:jpeg_abort_compress, libjpeg), Cvoid, (j_compress_ptr,), cinfo)
end

function jpeg_abort_decompress(cinfo)
    ccall((:jpeg_abort_decompress, libjpeg), Cvoid, (j_decompress_ptr,), cinfo)
end

function jpeg_abort(cinfo)
    ccall((:jpeg_abort, libjpeg), Cvoid, (j_common_ptr,), cinfo)
end

function jpeg_destroy(cinfo)
    ccall((:jpeg_destroy, libjpeg), Cvoid, (j_common_ptr,), cinfo)
end

function jpeg_resync_to_restart(cinfo, desired)
    ccall((:jpeg_resync_to_restart, libjpeg), boolean, (j_decompress_ptr, Cint), cinfo, desired)
end

function jpeg_read_icc_profile(cinfo, icc_data_ptr, icc_data_len)
    ccall((:jpeg_read_icc_profile, libjpeg), boolean, (j_decompress_ptr, Ptr{Ptr{JOCTET}}, Ptr{UInt32}), cinfo, icc_data_ptr, icc_data_len)
end

######################################################################################################
# macros from jpeglib.h

jpeg_create_compress(cinfo) =
    jpeg_CreateCompress(cinfo, JPEG_LIB_VERSION, sizeof(jpeg_compress_struct))
jpeg_create_decompress(cinfo) =
    jpeg_CreateDecompress(cinfo, JPEG_LIB_VERSION, sizeof(jpeg_decompress_struct))

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

end # module


text:
# morlet wavelet
function morlet(F::T, fs::T; dtype=Float32) where {T}
    wtime = range(-1, stop=1, length=Int(2 * fs))
    s = 6 / ( 2 * pi * F)
    wavelet = exp.(2im * pi * F .* wtime) .* exp.(-wtime.^2 ./ (2*s^2));
    return wavelet |> Vector{Complex{dtype}}
end

"""
function dwt(data, fs, frange; reflection=false, wavelet=morlet)

Arguments:
- data
- fs
- frange

Keyword Arguments:
- reflection
- wavelet

Return result:

Note:
discrete wavelet transformation by convolution between the fft results of
original values and the wavelet function.
"""
function dwt(data::Vector{D}, fs::T, frange::Vector{T};
             reflection::Bool=false, wavelet::Function=morlet) where {D, T}
    data_fft = data
    if reflection
        data_flip = reverse(data_fft, dims=1)
        data_fft = [data_flip;data_fft;data_flip]
    end

    data_fft = data_fft |> Vector{Complex{Float32}}
    result = zeros(Complex{Float32}, size(data,1), length(frange))

    for (fidx, ftarget) in enumerate(frange)
        w = wavelet(ftarget, fs)
        if reflection
            result[:, fidx] = DSP.conv(data_fft, w)[Int(fs)+size(data,1):end-Int(fs)-size(data, 1)]
        else
            result[:, fidx] = DSP.conv(data_fft, w)[Int(fs):end-Int(fs)]
        end
    end

    result
end


text:
# Regex
regx_bands 		= r"(.*)?what bands were born in (\d{2}s)(.*)?\?"
regx_recommend  = r"(.*)?can you recommend me (.*)?\?"
regx_tell_me    = r"(.*)?tell me about (.*)?"
regx_do_you 	= r"(.*)?what (.*)? do you (.*)?\?"

# Dictionaries
bands 	= Dict("60s"  => ["the beatles","the rolling stones","the beach boys"],
				"70s" => ["ramones","black sabbath","eagles"],
				"80s" => ["depeche mode","the smiths", "nine inch nails", "blur"],
				"90s" => ["the strokes", "weezer","oasis","the smiths"],
				"2000s" => ["artic monkeys","arcade fire","franz ferdinand","the smiths"])
genres 	= Dict("rock" 	   => ["the rolling stones","arctic monkeys"],
			   "pop"  	   => ["dua lipa","sam smith","lorde"],
			   "metal"	   => ["iron maiden","megadeth","black sabbath"],
			   "hip hop"   => ["kendrick lamar","mac miller","wiz khalifa"],
			   "reggaeton" => ["j balvin","bad bunny","nicky jam"])

# Reponses
resp_spot 	 = "Spotyman: "
resp_hello 	 = "What do you wanna know of music?"
resp_bye	 = "Good times!, see you later!"
resp_user_hello = s"\2, What you wanna know of music?"
resp_you 	 = s"I \3 \2 like ... Is this question relevant?"
resp_what 	 = "Some things I don't know"
resp_default = "Go on, ask me!"

# search argument element in values of dictionary argument dic
# return key of the value array where element was found otherwise return ""
function find(dic, element)
	for (key, value) in dic
			index = findall(b -> cmp(b, element) == 0, value)
			if !isempty(index)
				return key
			end
		end
	return ""
end

println("Hello I am Spotyman an expert chatbot of music. What is your name?")
while true
	s = readline()
	s = lowercase(string(s))
	if occursin(r".*hello.*", s)
		println(resp_spot * response_hello)
	elseif occursin(r"(.*)?my name is (.*)?", s)
		println(resp_spot * replace(s,r"(.*)?my name is(.*)?" => resp_user_hello))
	elseif occursin(regx_bands, s)
		println( resp_spot * join(get(bands, replace(s, regx_bands => s"\2"), 0),", ", " and ") * " were born.")
	elseif occursin(regx_recommend, s)
		r = join(get(genres, replace(s, regx_recommend => s"\2"),""),", ", " and ")
		println( resp_spot * r * " recommend you to hear.")
	elseif occursin(regx_tell_me, s)
		band = replace(s, regx_tell_me=> s"\2")
		time = find(bands, band)
		if cmp(time, "") != 0
			println(resp_spot * band * " were born in " * time)
		else
			println(resp_spot * " I don't know about " * band)
		end
	elseif occursin(regx_do_you , s)
		println(resp_spot * replace(s, regx_do_you => resp_you))
	elseif occursin(r"(.*)?what(.*)?", s)
		println(resp_what)
	elseif occursin(r"(.*)?bye(.*)?", s)
		println(resp_bye)
		break
	else
		println(resp_default)
	end
end

max_stars_repo_name: ablarry/Spotyman

max_issues_repo_name: ablarry/Spotyman

max_forks_repo_name: ablarry/Spotyman


text:
function doubling(a::Array{T, 2}, b::Array{T, 2}, g::Array{T, 2}, tol::T, maxiters::S) where {T <: AbstractFloat, S <: Int}

  retcode = 0
  iters = 1
  len = Inf

  gn = similar(g)

  while len > tol && iters <= maxiters
    gn = g+a*g*b
    len = maximum(abs,gn-g)
    a = a*a
    b = b*b
    g = copy(gn)
    iters += 1
  end

  if iters > maxiters
    retcode = 1
  end

  return gn,retcode

end


text:
using Zygote

f(x, y) = (x+exp(y), y)
invf(x, y) = (x-exp(y), y)

# ∂L/∂x2 = ∂L/∂x1*∂x1/∂x2 + ∂L/∂y1*∂y1/∂y2 = ∂L/∂x1*invf'(x2) + ∂L/∂y1*invf'(y2)
x1, y1 = 1.4, 4.4
x2, y2 = f(x,y)
function gf(x, y, gx, gy)
    x2, y2 = f(x, y)
    invJ1 = gradient((x2, y2)->invf(x2, y2)[1], x2, y2)
    invJ2 = gradient((x2, y2)->invf(x2, y2)[2], x2, y2)
    return (x2, y2, gx, gy)
end

gradient((x, y)->invf(x, y)[1], x, y)

mutable struct A{T}
    x::T
end

Base.:*(x1::A, x2::A) = A(x1.x*x2.x)
Base.:+(x1::A, x2::A) = A(x1.x+x2.x)
Base.zero(::A{T}) where T = A(T(0))

struct A2{T}
    x::T
end

Base.:*(x1::A2, x2::A2) = A2(x1.x*x2.x)
Base.:+(x1::A2, x2::A2) = A2(x1.x+x2.x)
Base.zero(::A2{T}) where T = A2(T(0))

struct BG{T}
    x::T
    g::B{T}
    BG(x::T) where T = new{T}(x)
end

struct BG{T}
    x::T
    g::BG{T}
    BG(x::T) where T = new{T}(x)
end

mutable struct AG{T}
    x::T
    g::AG{T}
    AG(x::T) where T = new{T}(x)
    AG(x::T, g::TG) where {T,TG} = new{T}(x, T(g))
end
Base.:*(x1::AG, x2::AG) = AG(x1.x*x2.x)
Base.:+(x1::AG, x2::AG) = AG(x1.x+x2.x)
Base.zero(::AG{T}) where T = AG(T(0))
init(ag::AG{T}) where T = (ag.g = AG(T(0)))

using BenchmarkTools
ma = fill(A(1.0), 100,100)
ma2 = fill(A2(1.0), 100,100)
function f(ma, mb)
    M, N, K = size(ma, 1), size(mb, 2), size(ma, 2)
    res = fill(zero(ma[1]), M, N)
    for i=1:M
        for j=1:N
            for k=1:K
                @inbounds res[i,j] += ma[i,k]*mb[k,j]
            end
        end
    end
    return res
end

@benchmark f(ma, ma)
@benchmark f(ma2, ma2)
ma = fill(AG(1.0), 100,100)
@benchmark ma*ma

a = A(0.4)
ag = AG(0.4)
using NiLangCore
@benchmark isdefined($ag, :g)
@benchmark $ag + $ag
ag.g = AG(0.0)
@benchmark $a + $a

struct SG{T}
    x::T
    g::Ref{T}
    SG(x::T) where T = new{T}(x)
end
Base.:*(x1::SG, x2::SG) = SG(x1.x*x2.x)
Base.:+(x1::SG, x2::SG) = SG(x1.x+x2.x)
Base.zero(::SG{T}) where T = SG(T(0))
init(ag::AG{T}) where T = (ag.g = AG(T(0)))

using BenchmarkTools
ma = fill(SG(1.0), 100,100)
@benchmark ma*ma

a = A(0.4)
ag = AG(0.4)
using NiLangCore
@benchmark isdefined($ag, :g)
@benchmark $ag + $ag
ag.g = AG(0.0)
@benchmark $a + $a

using NiLang, NiLang.AD
@i function test(x, one, N::Int)
    for i = 1:N
        x += one
    end
end

invcheckon(true)
@benchmark test'(Loss(0.0), 1.0, 1000000)


text:
export
    bom_to_hamiltonian

I(::Type{T}) where T<:Real = ones(T, 2)
Z(::Type{T}) where T<:Real = [one(T), -one(T)]

function _string_to_vec(T, s::String)
    result = [one(T)]
    for b = s
        if b == 'I'
            result = kron(result, I(T))
        elseif b == 'Z'
            result = kron(result, Z(T))
        else
            error("Bad string: chars have to be \'I\' or \'Z\'")
        end
    end
    result
end

function bom_to_hamiltonian(ising::Ising{T}) where T<:Real
    variables = used_variables(ising)
    degoffreedom = length(variables)

    result = zeros(T, 2^degoffreedom)
    for (key, val) = values(ising)
        pauli = ones(T, 1)
        for i=1:degoffreedom
            pauli = kron(pauli, variables[i] ∈ key ? Z(T) : I(T))
        end
        result += val * pauli
    end
    result
end

bom_to_hamiltonian(hobo::HOBO) = bom_to_hamiltonian(Ising(hobo))


text:
# j: counts dimension of Y
# k: counts component

# function expert_ll_pos_list(Y, model)
#     # result = [[hcat([expert_ll_pos.(model[k, j], Y[:, 4*(k-1)+1], Y[:, 4*(k-1)+2], Y[:, 4*(k-1)+3], Y[:, 4*(k-1)+4]) for j in 1:size(model)[2]]...) for k in 1:size(model)[1]]...]
#     return [[hcat([expert_ll_pos.(model[j, k], Y[:, 4*(j-1)+1], Y[:, 4*(j-1)+2], Y[:, 4*(j-1)+3], Y[:, 4*(j-1)+4]) for k in 1:size(model)[2]]...) for j in 1:size(model)[1]]...]
# end

function expert_ll_ind_mat(Y, model)
    result = fill(NaN, size(model)[1], size(model)[2])
    for j in 1:size(model)[1]
        for k in 1:size(model)[2]
            result[j, k] = LRMoE.expert_ll.(model[j, k], Y[4*(j-1)+1], Y[4*(j-1)+2], Y[4*(j-1)+3], Y[4*(j-1)+4])
        end
    end
    return result
end

function expert_ll_list(Y, model)
    result = fill(NaN, size(model)[1], size(model)[2], size(Y)[1])
    for i in 1:size(Y)[1]
        result[:,:,i] = expert_ll_ind_mat(Y[i,:], model[:,:,i])
    end
    return result
    # return [[hcat([expert_ll.(model[j, k], Y[:, 4*(j-1)+1], Y[:, 4*(j-1)+2], Y[:, 4*(j-1)+3], Y[:, 4*(j-1)+4]) for k in 1:size(model)[2]]...) for j in 1:size(model)[1]]...]
end

# function expert_tn_pos_list(Y, model)
#     return [[hcat([expert_tn_pos.(model[j, k], Y[:, 4*(j-1)+1], Y[:, 4*(j-1)+2], Y[:, 4*(j-1)+3], Y[:, 4*(j-1)+4]) for k in 1:size(model)[2]]...) for j in 1:size(model)[1]]...]
# end

function expert_tn_ind_mat(Y, model)
    result = fill(NaN, size(model)[1], size(model)[2])
    for j in 1:size(model)[1]
        for k in 1:size(model)[2]
            result[j, k] = LRMoE.expert_tn.(model[j, k], Y[4*(j-1)+1], Y[4*(j-1)+2], Y[4*(j-1)+3], Y[4*(j-1)+4])
        end
    end
    return result
end

function expert_tn_list(Y, model)
    result = fill(NaN, size(model)[1], size(model)[2], size(Y)[1])
    for i in 1:size(Y)[1]
        result[:,:,i] = expert_tn_ind_mat(Y[i,:], model[:,:,i])
    end
    return result
    # return [[hcat([expert_tn.(model[j, k], Y[:, 4*(j-1)+1], Y[:, 4*(j-1)+2], Y[:, 4*(j-1)+3], Y[:, 4*(j-1)+4]) for k in 1:size(model)[2]]...) for j in 1:size(model)[1]]...]
end

# function expert_tn_bar_pos_list(Y, model)
#     return [[hcat([expert_tn_bar_pos.(model[j, k], Y[:, 4*(j-1)+1], Y[:, 4*(j-1)+2], Y[:, 4*(j-1)+3], Y[:, 4*(j-1)+4]) for k in 1:size(model)[2]]...) for j in 1:size(model)[1]]...]
# end

function expert_tn_bar_ind_mat(Y, model)
    result = fill(NaN, size(model)[1], size(model)[2])
    for j in 1:size(model)[1]
        for k in 1:size(model)[2]
            result[j, k] = LRMoE.expert_tn_bar.(model[j, k], Y[4*(j-1)+1], Y[4*(j-1)+2], Y[4*(j-1)+3], Y[4*(j-1)+4])
        end
    end
    return result
end

function expert_tn_bar_list(Y, model)
    result = fill(NaN, size(model)[1], size(model)[2], size(Y)[1])
    for i in 1:size(Y)[1]
        result[:,:,i] = expert_tn_bar_ind_mat(Y[i,:], model[:,:,i])
    end
    return result
    # return [[hcat([expert_tn_bar.(model[j, k], Y[:, 4*(j-1)+1], Y[:, 4*(j-1)+2], Y[:, 4*(j-1)+3], Y[:, 4*(j-1)+4]) for k in 1:size(model)[2]]...) for j in 1:size(model)[1]]...]
end

# function gate_dim_expert_pos(gate, Y, model)
#     return gate + sum(expert_ll_pos_list(Y, model), dims = 1)[1]
# end

# function loglik_aggre_all(ll_mat)
#     return sum(ll_mat)
# end

function loglik_np(Y, gate, model)

    # By dimension, then by component
    # expert_ll_pos_dim_comp = expert_ll_pos_list(Y, model)
    # expert_tn_pos_dim_comp = expert_tn_pos_list(Y, model)
    # expert_tn_bar_pos_dim_comp = expert_tn_bar_pos_list(Y, model)

    expert_ll_dim_comp = expert_ll_list(Y, model)
    expert_tn_dim_comp = expert_tn_list(Y, model)
    expert_tn_bar_dim_comp = expert_tn_bar_list(Y, model)

    # Aggregate by dimension
    # expert_ll_pos_comp = loglik_aggre_dim(expert_ll_pos_dim_comp)
    # expert_tn_pos_comp = loglik_aggre_dim(expert_tn_pos_dim_comp)
    # expert_tn_bar_pos_comp = loglik_aggre_dim(expert_tn_bar_pos_dim_comp)

    expert_ll_comp = loglik_aggre_dim(expert_ll_dim_comp)
    expert_tn_comp = loglik_aggre_dim(expert_tn_dim_comp)
    expert_tn_bar_comp = loglik_aggre_dim(expert_tn_bar_dim_comp)

    # Adding the gating function
    # gate_expert_ll_pos_comp = loglik_aggre_gate_dim(gate, expert_ll_pos_comp)
    # gate_expert_tn_pos_comp = loglik_aggre_gate_dim(gate, expert_tn_pos_comp)
    # gate_expert_tn_bar_pos_comp = gate + log1mexp.(expert_tn_pos_comp)
    # gate_expert_tn_bar_pos_comp_k = gate + expert_tn_bar_pos_comp

    gate_expert_ll_comp = loglik_aggre_gate_dim(gate, expert_ll_comp)
    gate_expert_tn_comp = loglik_aggre_gate_dim(gate, expert_tn_comp)
    gate_expert_tn_bar_comp = gate + log1mexp.(expert_tn_comp)
    gate_expert_tn_bar_comp_k = gate + log1mexp.(expert_tn_bar_comp)
    gate_expert_tn_bar_comp_z_lat = loglik_aggre_gate_dim(gate, expert_tn_bar_comp)

    # Aggregate by component
    # gate_expert_ll_pos = loglik_aggre_gate_dim_comp(gate_expert_ll_pos_comp)
    # gate_expert_tn_pos = loglik_aggre_gate_dim_comp(gate_expert_tn_pos_comp)
    # gate_expert_tn_bar_pos = loglik_aggre_gate_dim_comp(gate_expert_tn_bar_pos_comp)
    # gate_expert_tn_bar_pos_k = loglik_aggre_gate_dim_comp(gate_expert_tn_bar_pos_comp_k)

    gate_expert_ll = loglik_aggre_gate_dim_comp(gate_expert_ll_comp)
    gate_expert_tn = loglik_aggre_gate_dim_comp(gate_expert_tn_comp)
    gate_expert_tn_bar = loglik_aggre_gate_dim_comp(gate_expert_tn_bar_comp)
    gate_expert_tn_bar_k = loglik_aggre_gate_dim_comp(gate_expert_tn_bar_comp_k)
    gate_expert_tn_bar_z_lat = loglik_aggre_gate_dim_comp(gate_expert_tn_bar_comp_z_lat)

    # Normalize by tn & tn_bar
    # norm_gate_expert_ll_pos = gate_expert_ll_pos - gate_expert_tn_pos
    # norm_gate_expert_ll = gate_expert_ll - gate_expert_tn
    norm_gate_expert_ll = gate_expert_ll - gate_expert_tn_bar_k

    # Sum over all observations
    ll = sum(norm_gate_expert_ll)

    return (# expert_ll_pos_dim_comp = expert_ll_pos_dim_comp,
            # expert_tn_pos_dim_comp = expert_tn_pos_dim_comp,
            # expert_tn_bar_pos_dim_comp = expert_tn_bar_pos_dim_comp,

            expert_ll_dim_comp = expert_ll_dim_comp,
            expert_tn_dim_comp = expert_tn_dim_comp,
            expert_tn_bar_dim_comp = expert_tn_bar_dim_comp,

            # expert_ll_pos_comp = expert_ll_pos_comp,
            # expert_tn_pos_comp = expert_tn_pos_comp,
            # expert_tn_bar_pos_comp = expert_tn_bar_pos_comp,

            expert_ll_comp = expert_ll_comp,
            expert_tn_comp = expert_tn_comp,
            expert_tn_bar_comp = expert_tn_bar_comp,

            # gate_expert_ll_pos_comp = gate_expert_ll_pos_comp,
            # gate_expert_tn_pos_comp = gate_expert_tn_pos_comp,
            # gate_expert_tn_bar_pos_comp = gate_expert_tn_bar_pos_comp,
            # gate_expert_tn_bar_pos_comp_k = gate_expert_tn_bar_pos_comp_k,

            gate_expert_ll_comp = gate_expert_ll_comp,
            gate_expert_tn_comp = gate_expert_tn_comp,
            gate_expert_tn_bar_comp = gate_expert_tn_bar_comp,
            gate_expert_tn_bar_comp_k = gate_expert_tn_bar_comp_k,
            gate_expert_tn_bar_comp_z_lat = gate_expert_tn_bar_comp_z_lat,

            # gate_expert_ll_pos = gate_expert_ll_pos,
            # gate_expert_tn_pos = gate_expert_tn_pos,
            # gate_expert_tn_bar_pos = gate_expert_tn_bar_pos,
            # gate_expert_tn_bar_pos_k = gate_expert_tn_bar_pos_k,

            gate_expert_ll = gate_expert_ll,
            gate_expert_tn = gate_expert_tn,
            gate_expert_tn_bar = gate_expert_tn_bar,
            gate_expert_tn_bar_k = gate_expert_tn_bar_k,
            gate_expert_tn_bar_z_lat = gate_expert_tn_bar_z_lat,

            # norm_gate_expert_ll_pos = norm_gate_expert_ll_pos,
            norm_gate_expert_ll = norm_gate_expert_ll,

            ll = ll
        )

end

max_stars_repo_name: sparktseung/JuliaLRMoE

max_issues_repo_name: sparktseung/JuliaLRMoE

max_forks_repo_name: sparktseung/JuliaLRMoE


text:
@testset "$TEST $G" begin

g = CompleteDiGraph(5, DG)
@test nv(g) == 5 && ne(g) == 20
@test typeof(g) == DG

g = CompleteGraph(5, G)
@test nv(g) == 5 && ne(g) == 10
@test typeof(g) == G

g = CompleteBipartiteGraph(5, 8, G)
@test nv(g) == 13 && ne(g) == 40
@test typeof(g) == G
g = StarDiGraph(5, DG)
@test nv(g) == 5 && ne(g) == 4
g = StarGraph(5, G)
@test typeof(g) == G
@test nv(g) == 5 && ne(g) == 4
g = StarGraph(1, G)
@test typeof(g) == G
@test nv(g) == 1 && ne(g) == 0

g = PathDiGraph(5, DG)
@test typeof(g) == DG
@test nv(g) == 5 && ne(g) == 4
g = PathGraph(5, G)
@test typeof(g) == G
@test nv(g) == 5 && ne(g) == 4

g = CycleDiGraph(5, DG)
@test typeof(g) == DG
@test nv(g) == 5 && ne(g) == 5
g = CycleGraph(5, G)
@test typeof(g) == G
@test nv(g) == 5 && ne(g) == 5

g = WheelDiGraph(5, DG)
@test typeof(g) == DG
@test nv(g) == 5 && ne(g) == 8
g = WheelGraph(5, G)
@test typeof(g) == G
@test nv(g) == 5 && ne(g) == 8

g = Grid([3,3,4], G)
@test typeof(g) == G
@test nv(g) == 3*3*4
@test ne(g) == 75
@test maximum(degree(g)) == 6
@test minimum(degree(g)) == 3

g = CliqueGraph(3,5, G)
@test typeof(g) == G
@test nv(g) == 15 && ne(g) == 20
@test g[1:3] == CompleteGraph(3, G)

g = crosspath(BinaryTree(2, G), 3)
@test typeof(g) == G
@test nv(g) == 9
@test ne(g) == 12

g = DoubleBinaryTree(3, G)
@test typeof(g) == G
@test nv(g) == 14
@test ne(g) == 13

rg3 = RoachGraph(3, G)
@test typeof(g) == G
@test nv(g) == 14
@test ne(g) == 13

end # testset


text:
function gauleg(a,b,n)
#==
    1D Gauss-Legendre quadrature rule.

    quadrule = gauleg(a,b,n) computes the n-point Gauss-Legendre
    quadrature rule on the interval [a,b] up to mashine precision.

    Note that all quadrature rules obtained from gauleg are of order 2*n-1.

   The output named tuple quadrule contains the following fields:
    w n-by-1 vector specifying the weights of the quadrature rule.
    X n-by-1 vector specifying the abscissae of the quadrature rule.

   Example:
    quad = gauleg(0,1,10);
    quad.x
    quad.w
=#
    if n < 0
        throw(DomainError(n, "Input n must be a non-negative integer"))
    elseif n==1
        x = 0; w = 2;
    else
        d = zeros(n)
        c = zeros(n-1)
        for i in 1:(n-1)
            c[i]=i/sqrt(4*i*i-1)
        end
        J=LinearAlgebra.SymTridiagonal(d,c)
        x=LinearAlgebra.eigvals(J)
        vec=LinearAlgebra.eigvecs(J)
        w=2*vec[1,:].*vec[1,:]

    end

    # The above rule computed a quadrature rule for the reference interval.
    # The reference interval is always [-1,1], and if [a,b] != [-1,1] then we
    # re-map the abscissas (x) and rescale the weights.

    xm = (b+a)/2; # midpoint
    xl = (b-a)/2; # area of the requested interval / area of the reference interval

    x = xm .+ xl*x;
    w = w * xl;

    return (x=x,w=w) #named tuple
end


text:
# Autogenerated wrapper script for Hwloc_jll for i686-w64-mingw32
export libhwloc

JLLWrappers.@generate_wrapper_header("Hwloc")
JLLWrappers.@declare_library_product(libhwloc, "libhwloc-15.dll")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libhwloc,
        "bin\\libhwloc-15.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module MagicUnderscores

export @_, @lazy_

import Base: tail

#-------------------------------------------------------------------------------
# Machinery for lowered form
struct Placeholder ; end
Base.show(io::IO, p::Placeholder) = print(io, "_")

(::Placeholder)(x) = x

struct PartialApply{Func,Args} <: Function
    f::Func
    args::Args
end
Base.show(io::IO, p::PartialApply) = print(io, "$(p.f)($(join(sprint.(show,p.args), ",")))")

struct TightPartialApply{Func,Arg1,Args} <: Function
    f::Func
    arg1::Arg1
    args::Args
end
Base.show(io::IO, p::TightPartialApply) = print(io, "$(p.f)($(p.arg1), $(join(sprint.(show,p.args), ",")))")

struct TightPartialApply2{Func,Arg1,Arg2,Args} <: Function
    f::Func
    arg1::Arg1
    arg2::Arg2
    args::Args
end
Base.show(io::IO, p::TightPartialApply2) = print(io, "$(p.f)($(p.arg1), $(p.arg2), $(join(sprint.(show,p.args), ",")))")

for T in [PartialApply,TightPartialApply,TightPartialApply2]
    @eval (p::$T)(args...) = fillplaceholders((p,), args)[1][1]
end

#
#   fillplaceholders(tofill, args)
#
# Replaces any `Placeholder`s in the tuple `to_fill` with elements of `args`,
# from left to right. Returns the tuple `(filled,extras)`, where `extras` are
# any remaining elements of `args` which weren't consumed by a placeholder.
#
# When an element of `tofill` is a partially applied function, execute it
# eagerly, consuming elements from `args` and inlining the result in the
# `filled` tuple.
#
fillplaceholders(tofill::Tuple{}, extra) = (),extra
function fillplaceholders(tofill::Tuple{Placeholder,Vararg}, extra)
    inner, extra2 = fillplaceholders(tail(tofill), tail(extra))
    (extra[1], inner...), extra2
end
function fillplaceholders(tofill, extra)
    inner, extra2 = fillplaceholders(tail(tofill), extra)
    (tofill[1], inner...), extra2
end

# Versions for evaluating partial applications
function fillplaceholders(tofill::Tuple{PartialApply,Vararg}, extra)
    pa = tofill[1]
    inner, extra2 = fillplaceholders(pa.args, extra)
    func_result = pa.f(inner...)
    inner2, extra3 = fillplaceholders(tail(tofill), extra2)
    (func_result, inner2...), extra3
end

function fillplaceholders(tofill::Tuple{TightPartialApply,Vararg}, extra)
    pa = tofill[1]
    inner, extra2 = fillplaceholders(pa.args, extra)
    func_result = pa.f(pa.arg1, inner...)
    inner2, extra3 = fillplaceholders(tail(tofill), extra2)
    (func_result, inner2...), extra3
end

function fillplaceholders(tofill::Tuple{TightPartialApply2,Vararg}, extra)
    pa = tofill[1]
    inner, extra2 = fillplaceholders((pa.arg1, pa.args...), extra)
    func_result = pa.f(inner[1], pa.arg2, tail(inner)...)
    inner2, extra3 = fillplaceholders(tail(tofill), extra2)
    (func_result, inner2...), extra3
end

# Bind underscore; loose binding by default
function ubind(f, args...)
    PartialApply(f, args)
end

# Execute partial application, assuming it's fully applied
materialize(a) = a
materialize(p::Union{PartialApply,TightPartialApply,TightPartialApply2}) = p()

#-------------------------------------------------------------------------------
# Lowering
function has_underscores(ex)
    if ex === :_
        return true
    elseif !(ex isa Expr)
        return false
    else
        return any(has_underscores, ex.args)
    end
end

function lower_underscores(ex)
    # Note that this is more complicated than it should be because it has to
    # detect and handle special lowerings like getproperty
    if ex === :_
        return MagicUnderscores.Placeholder()
    elseif !(ex isa Expr)
        return ex
    elseif has_underscores(ex)
        if ex.head == :call
            return Expr(:call, MagicUnderscores.ubind, map(lower_underscores, ex.args)...)
        elseif ex.head == :. && length(ex.args) == 2 && ex.args[2] isa QuoteNode
            # getproperty
            return lower_underscores(Expr(:call, :getproperty, ex.args...))
        end
    end
    return ex
end

"""
    @_ expression

Transform uses of `_` in `expression` into arguments of a lambda, and evaluate
the resulting expression eagerly.
"""
macro _(ex)
    esc(Expr(:call, MagicUnderscores.materialize, lower_underscores(ex)))
end

"""
    @lazy_ expression

Lazy (unevaluated) version of `@_`.
"""
macro lazy_(ex)
    esc(lower_underscores(ex))
end

#-------------------------------------------------------------------------------
# Tight binding behaviour for for some Base functions
for func in [:map, :filter, :findall]
    @eval function ubind(::typeof($func), f, a, args...)
        TightPartialApply($func, f, (a, args...))
    end
end

for func in [:|>]
    @eval function ubind(::typeof($func), v, f, args...)
        TightPartialApply2($func, v, f, args)
    end
end

end