text:
@testset "Indices" begin
    @test Indices() isa Indices{Any}
    @test_throws IndexError Indices([1, 1])
    @test @inferred(Indices(2*x for x in 1:10))::Indices == Indices(2:2:20)
    @test @inferred(distinct(mod(x, 5) for x in 0:9))::Indices == Indices(0:4)

    h = Indices{Int64}()

    @test isinsertable(h)
    @test length(h) == 0
    @test keys(h) === h
    @test h == h
    @test h == copy(h)
    @test isempty(h)
    @test isequal(copy(h), h)
    @test_throws IndexError h[10]
    @test length(unset!(h, 10)) == 0
    io = IOBuffer(); print(io, h); @test String(take!(io)) == "{}"
    io = IOBuffer(); show(io, MIME"text/plain"(), h); @test String(take!(io)) == "0-element Indices{Int64}"
    @test_throws IndexError delete!(h, 10)

    insert!(h, 10.0)

    @test length(h) == 1
    @test keys(h) === h
    @test unique(h) === h
    @test h == h
    @test h == copy(h)
    @test !isempty(h)
    @test isequal(copy(h), h)
    @test h[10.0] == 10
    @test_throws IndexError insert!(h, 10)
    @test length(set!(h, 10)) == 1
    @test_throws IndexError insert!(h, 10)
    io = IOBuffer(); print(io, h); @test String(take!(io)) == "{10}"
    io = IOBuffer(); show(io, MIME"text/plain"(), h); @test String(take!(io)) == "1-element Indices{Int64}\n 10"
    @test !isequal(h, empty(h))
    @test isequal(h, copy(h))
    @test isempty(empty(h))

    delete!(h, 10.0)

    @test isequal(h, Indices{Int64}())

    for i = 2:2:1000
        insert!(h, i)
    end
    @test issetequal(h, Indices(2:2:1000))
    @test all(in(i, h) == iseven(i) for i in 2:1000)
    @test isempty(empty!(h))

    # set
    @test length(set!(h, 1)) == 1
    @test length(set!(h, 2, 2)) == 2
    @test length(set!(h, 3.0, 3.0)) == 3
    @test_throws ErrorException set!(h, 4, 5)

    @testset "Comparison" begin
        i1 = Indices([1,2,3])
        i2 = Indices([1,2])
        i3 = Indices([1,2,3,4])
        i4 = Indices([3,2,1])

        @test isequal(i1, i1)
        @test hash(i1) == hash(copy(i1))

        @test isless(i2, i1)
        @test !isless(i1, i2)
        @test !isequal(i1, i2)
        @test hash(i1) != hash(i2)

        @test isless(i2, i1)
        @test !isless(i1, i2)
        @test !isequal(i1, i2)

        @test isless(i1, i3)
        @test !isless(i3, i1)
        @test !isequal(i1, i3)

        @test isless(i1, i4)
        @test !isless(i4, i1)
        @test !isequal(i1, i4)

        @test cmp(i1, i1) == 0
        @test cmp(i2, i1) == -1
        @test cmp(i1, i2) == 1
        @test cmp(i3, i1) == 1
        @test cmp(i1, i3) == -1
        @test cmp(i4, i1) == 1
        @test cmp(i1, i4) == -1

        i5 = Indices([1,2,missing])
        @test isequal(i5, i5)
        @test !isless(i5, i5)
        @test (i5 == i5) === missing
        @test (i5 == Indices([1,2,missing])) === missing
    end

    @testset "Adapated from Dict tests from Base" begin
        h = Indices{Int}()
        N = 10000

        for i in 1:N
            insert!(h, i)
        end
        for i in 1:N
            @test i in h
        end
        for i in 1:2:N
            delete!(h, i)
        end
        for i in 1:N
            @test (i in h) == iseven(i)
        end
        for i in 1:2:N
            insert!(h, i)
        end
        for i in 1:N
            @test i in h
        end
        for i in 1:N
            delete!(h, i)
        end
        @test isempty(h)
        insert!(h, 77)
        @test 77 in h
        for i in 1:N
            set!(h, i)
        end
        for i in 1:N
            @test i in h
        end
        for i in 1:2:N
            delete!(h, i)
        end
        for i in 1:N
            @test (i in h) == iseven(i)
        end
        for i in N+1:2N
            insert!(h, i)
        end
        for i in 1:2N
            @test (i in h) == (i > N || iseven(i))
        end
    end

    @testset "distinct" begin
        res = Indices([1,2,3])
        @test distinct(res) === res
        @test isequal(distinct([1,2,3]), res)
        @test isequal(distinct([1,2,3,1]), res)
        @test isequal(distinct([1,2,3]), res)
    end

    @testset "set logic" begin
        i1 = Indices([1,2])
        i2 = Indices([2,3])
        i3 = Indices([3,4])
        i4 = Indices([2])

        @test !issetequal(i1, i2)
        @test !(i1 ⊆ i2)
        @test i4 ⊆ i1

        if VERSION < v"1.5-"
            @test !disjoint(i1, i2)
            @test disjoint(i1, i3)
        else
            @test !isdisjoint(i1, i2)
            @test isdisjoint(i1, i3)
        end

        @test isequal(union(i1, i2), Indices([1,2,3]))
        @test isequal(union(i2, i1), Indices([2,3,1]))
        @test isequal(union(i1, i3), Indices([1,2,3,4]))

        @test isequal(intersect(i1, i2), Indices([2]))
        @test isequal(intersect(i1, i3), Indices([]))

        @test isequal(setdiff(i1, i2), Indices([1]))
        @test isequal(setdiff(i1, i3), Indices([1, 2]))

        @test isequal(symdiff(i1, i2), Indices([1, 3]))
        @test isequal(symdiff(i1, i3), Indices([1, 2, 3, 4]))
    end

    @testset "covert" begin
        i = Indices{Int32}([1,2,3])
        ai = ArrayIndices{Int32}([1,2,3])

        @test convert(AbstractIndices{Int32}, i) === i
        @test convert(AbstractIndices{Int64}, i)::Indices{Int64} == i

        @test convert(Indices{Int32}, i) === i
        @test convert(Indices{Int64}, i)::Indices{Int64} == i

        @test convert(Indices{Int32}, ai)::Indices{Int32} == i
        @test convert(Indices{Int64}, ai)::Indices{Int64} == i
    end

    @testset "rand" begin
        inds = Indices(["a", "b", "c", "d", "e"])
        for i = 1:100
            @test rand(inds) in ["a", "b", "c", "d", "e"]
        end

        delete!(inds, "b") # Test handling of deleted slots
        for i = 1:100
            @test rand(inds) in ["a", "c", "d", "e"]
        end
    end

    @testset "sort" begin
        inds = Indices([1, 3, 2])
        @test sort(inds)::Indices == Indices([1, 2, 3])
        @test sort(inds; rev=true)::Indices == Indices([3, 2, 1])

        indscopy = deepcopy(inds)
        sort!(indscopy)
        @test indscopy == Indices([1, 2, 3])

        indscopy = deepcopy(inds)
        sort!(indscopy; rev=true)
        @test indscopy == Indices([3, 2, 1])
    end

    # TODO: token interface
end


text:
using JuliaWeBWorK
using Plots

## --------------------------------------------------

meta=(Project="Plotting  questions",
     Question="3"
      )

qs = JuliaWeBWorK.QUESTIONS()

intro = jmt"""

[Details for students](https://webwork.maa.org/wiki/Student_Information)

----

# Including plots

This shows  how  fixed and randomized plots can be included; it also shows how to use [JSXGraph](https://jsxgraph.uni-bayreuth.de/) to create dynamic graphics (using JavaScript commands).

For randomized plots, as each  *possible*  plot is precomputed and serialized,
this should not be  used  over large event spaces.

"""

## --------------------------------------------------

## fixed plot, using string interpolation
## This is the simplest and suggested way to include a plot
p = plot(sinpi, -1/2,  5/2)
plot!(p, zero, -1/2, 5/2)
numericq(jmt"""

# Include using `(Plot(p))`

![Figure]($(Plot(p)))

(1) In the figure above, what are the zeros of the function?
""",
         () -> List((0,1,2)),
         ()
         ) |> qs

## reading a plt in from a file, `pngfile`
## here we create `pngfile` and fill it with an image; this isn't generally needed
pngfile = tempname() * ".png"
savefig(p, pngfile)

numericq(jmt"""

----
# Include file through `(File(pngfile))`

![Figure]($(File(pngfile)))

""",()->1,()) |> qs

## Reading a plot from a url
label(jmt"""

----
# Using Markdown to include from a url

![some image](https://upload.wikimedia.org/wikipedia/commons/thumb/0/0f/Tangent_to_a_curve.svg/200px-Tangent_to_a_curve.svg.png)
""") |> qs

## Randomized plot for students using the `plot` "question" type (in quotes
## as `plotq` only prings a randomized plot, the question comes after, hence
## the use of a `randomizer` to share the randomization.

r = randomizer(2:3) |> qs

label(jmt"""
----
# Randomized plots

Plot of a periodic function
""") |> qs

plotq("", (a) -> Plot(plot(sin, 0,2a*pi)), r) |> qs

numericq("(2) How many periods of the periodic function are shown in the above figure?", (a)  -> a, r) |> qs

## This shows use of [jsxgraph](https://jsxgraph.uni-bayreuth.de/) graphics for
## interactive graphics. This is not as integrated into WeBWorK as GeoGebra, but
## are easy enough to author and quite illustrative.
## To load from a file, the `JuliaWeBWorK.INCLUDE` is available to manage relative directories

label(jmt"""
----
# Using jsxgraph for interactive plots
""") |> qs

jsxgraph("""
var brd = JXG.JSXGraph.initBoard('jxgbox', {boundingbox: [-10, 10, 10, -10]});
var a = brd.create('point', [-2, 1]);
var b = brd.create('point', [-4, -5]);
var c = brd.create('point', [3, -6]);
var d = brd.create('point', [2, 3]);
var p = brd.create('polygon', [a, b, c, d], {hasInnerPoints: true});
"""; domid="jxgbox") |> qs

p = Page(intro, qs; meta...)

ext: jl


text:
using Random
using BenchmarkTools

function comprassion()
    array = []
    for i in 1:999999
        push!(array, rand())
    end
    array2 = array
    @btime sort!(array)
    @btime mergesort!(array2)
end


text:
################################################################################
# Helpers
################################################################################

function add2sub(v::SubArray, e)
	v .+= e
	return nothing
end

function add2sub(v::SubArray, e::Diagonal{T,SVector{n,T}}) where {n,T}
	for i = 1:n
		v[i,i] += e[i,i]
	end
	return nothing
end

function addI2sub(v::SubArray, e)
	n = size(v)[1]
	for i = 1:n
		v[i,i] += e
	end
	return nothing
end

function sparse_zero!(spm::SparseMatrixCSC)
	n = length(spm.nzval)
	for i = 1:n
		spm.nzval[i] = 0.0
	end
	return nothing
end

################################################################################
# Printers
################################################################################

function display_solver_header()
	@printf(
		"%-3s %-2s %-2s %-6s %-6s %-6s \n",
		"out",
		"in",
		"α",
		"Δ",
		"res",
		"reg",
		)
	return nothing
end

function display_solver_data(k, l, j, Δ, res_norm, reg)#, condi, loss, val_scale, jac_scale)
	@printf(
		"%-3s %-2s %-2s %-6s %-6s %-6s \n",
		k,
		l,
		j,
		@sprintf("%.0e", Δ),
		@sprintf("%.0e", res_norm),
		@sprintf("%.0e", reg.x),
		)
	return nothing
end

function scn(a::Number; digits::Int=1)
	@assert digits >= 0
    # a = m x 10^e
    if a == 0
        e = 0
        m = 0.0
    else
        e = Int(floor(log(abs(a))/log(10)))
        m = a*exp(-e*log(10))
    end
    m = round(m, digits=digits)
    if digits == 0
        m = Int(floor(m))
		strm = string(m)
	else
		strm = string(m)
		is_neg = m < 0.
		strm = strm*"0"^abs(2+digits+is_neg-length(strm))
    end
    sgn = a >= 0 ? " " : ""
    sgne = e >= 0 ? "+" : ""
    return "$sgn$(strm)e$sgne$e"
end

################################################################################
# Video
################################################################################

function convert_video_to_gif(video_file_path::AbstractString, output_path::AbstractString="output.gif";
    framerate::Int=30, start_time=0., duration=1e3, overwrite=false, width::Int=1080, height::Int=-2, hq_colors::Bool=false)
    output_path = abspath(output_path)

    if !isfile(video_file_path)
        error("Could not find the input file $video_file_path")
    end
    if isfile(output_path) && !overwrite
        error("The output path $output_path already exists. To overwrite that file, you can pass `overwrite=true` to this function")
    end

    mktempdir() do tmpdir
        # run(MeshCat.unpack_cmd(video_file_path, tmpdir, ".mp4", nothing)) # unpack the .tar file
        # cmd = ["-r", string(framerate), "-i", "%07d.png", "-vcodec", "libx264", "-preset", "slow", "-crf", "18"]
        color_map = hq_colors ?
            "[0:v] fps=$framerate, scale=$width:$height,split [a][b];[a] palettegen=stats_mode=single [p];[b][p] paletteuse=new=1" :
            "[0:v] fps=$framerate, scale=$width:$height,split [a][b];[a] palettegen [p];[b][p] paletteuse"
        cmd = ["-ss", string(start_time), "-t", string(duration), "-i", video_file_path, "-filter_complex", color_map]
        if overwrite
            push!(cmd, "-y")
        end
        push!(cmd, output_path)

        cd(tmpdir) do
            FFMPEG.exe(cmd...)
        end
    end
    @info("Saved output as $output_path")
    return output_path
end


text:
module Unity
import JSON
import ColorTypes
import ColorTypes: RGB, RGBA
import Colors: weighted_color_mean
import Base.TCPServer
import GeometryTypes
import GeometryTypes: SimpleRectangle, Point3f0, Face

include(joinpath("palette.jl"))
include(joinpath("colormaps.jl"))

type UnityVector3
    x::Float32
    y::Float32
    z::Float32
end

type UnityColor
    r::Float32
    g::Float32
    b::Float32
    a::Float32
end

type UnityText
    text::String
    pos::UnityVector3
    scale::UnityVector3
    rot::UnityVector3
    color::UnityColor
end

function Base.convert(::Type{UnityColor},c::ColorTypes.RGBA{Float32})
    return UnityColor(c.r,c.g,c.b,c.alpha)
end
function Base.convert(::Type{UnityVector3},x::GeometryTypes.Point3f0)
    return UnityVector3(x[1],x[2],x[3])
end
#function Base.convert(::Vector{UInt32},x::Array{GeometryTypes.Face{3,Int32},1})
#    vcat(map(xx->convert(Vector{UInt32},xx-1),x)...)
#end

#Unity mesh with c-like indexing
type UnityMesh
    id::String
    vertices::Vector{UnityVector3}
    points::Vector{UInt32}
    lines::Vector{UInt32}
    triangles::Vector{UInt32}
    colors::Vector{UnityColor}
    options::Vector{String}
    text::Vector{UnityText}
end

function UnityMesh(id::String, vertices::Vector, points::Vector,  lines::Vector, triangles::Vector, colors::Vector, options::Vector)
    return UnityMesh(id,vertices,points,lines,triangles,colors,options,UnityText[])
end

function Base.write(socket::TCPSocket, um::UnityMesh)
    jum = JSON.json(um)
    retval = write(socket, jum*"UNITY_MESH_JSON_FORMATTED")
    sleep(.1)
    return retval
end

type UnityCameraSettings
    id::String
    main_camera_position::Vector{UnityVector3} #should be all vectors of length 0 or 1
    main_camera_rotation::Vector{UnityVector3}
    main_camera_scale::Vector{UnityVector3}
    background_color::Vector{UnityColor}
    perspective::Vector{Bool}
end

function Base.write(socket::TCPSocket, ucs::UnityCameraSettings)
    jum = JSON.json(ucs)
    retval = write(socket, jum*"UNITY_CAMERA_SETTINGS")
    sleep(.1)
    return retval
end

#Unity Pyramid mesh with c-like indices
type PyramidMesh
    id::String
    vertices::Vector{UnityVector3}
    pyramids::Vector{UInt32}
    colors::Vector{ColorTypes.RGBA{Float32}}
end

begin
local const pattern = [ 0,2,1,0,3,2,2,3,1,0,1,3 ]
function Base.convert(::Type{UnityMesh},msh::PyramidMesh,dublic_vert::Bool=false)
    if dublic_vert
        return convert_and_duplicate(UnityMesh,msh,pattern)
    else
        return convert(UnityMesh,msh,pattern)
    end
end
end

function Base.convert(::Type{UnityMesh},msh::PyramidMesh,pattern::Vector{Int})
    triangles = Vector{UInt32}()
    inds = msh.pyramids
    @assert mod(length(inds),4) == 0
    for i = 1:4:length(inds)
        append!(triangles,
            [
            inds[i+pattern[1]],inds[i+pattern[2]],inds[i+pattern[3]],
            inds[i+pattern[4]],inds[i+pattern[5]],inds[i+pattern[6]],
            inds[i+pattern[7]],inds[i+pattern[8]],inds[i+pattern[9]],
            inds[i+pattern[10]],inds[i+pattern[11]],inds[i+pattern[12]]
            ]
            )
    end
    return UnityMesh(msh.id,msh.vertices,UInt32[],UInt32[],triangles,msh.colors)
end

function convert_and_duplicate(::Type{UnityMesh},msh::PyramidMesh,pattern::Vector{Int})
    id = msh.id
    verts = msh.vertices
    lines = UInt32[]
    points = UInt32[]
    inds = msh.pyramids
    clrs = msh.colors
    @assert mod(length(inds),4) == 0

    #new mesh
    #triangles = Vector{Int32}()
    vertices = Vector{UnityVector3}()
    colors = Vector{ColorTypes.RGBA{Float32}}()

    for i = 1:4:length(inds)

        tmpinds = inds[i:i+3]+1
        tmpverts = verts[tmpinds]
        tmpclrs = clrs[tmpinds]

        append!(vertices, tmpverts[pattern+1])
        append!(colors, tmpclrs[pattern+1])

    end

    triangles = UInt32[i-1 for i = 1:length(vertices)]

    return UnityMesh(id,vertices,points,lines,triangles,colors)

end

function Base.write(tcpstream::TCPSocket, um::PyramidMesh)
    jum = JSON.json(convert(UnityMesh,um))
    return write(tcpstream, jum)
end

import Combinatorics

function boundary(upm::PyramidMesh)
    verts = Vector{GeometryTypes.Point3f0}()
    pyramids = Vector{UInt32}()
    clrs = Vector{ColorTypes.RGBA{Float32}}()
    for i = 1:4:length(upm.pyramids)
        # todo
    end
end

function ColorBar(clrmap::BoundedColorMap)

    add_verts_faces_color!(vertices,faces,lines,colors,a,b,c,d,ca,cb,flip) = begin
        push!(vertices,a) #0
        push!(vertices,b) #1
        push!(vertices,c) #2
        push!(vertices,d) #3
        push!(colors,ca)
        push!(colors,cb)
        push!(colors,ca)
        push!(colors,cb)
        flip ? push!(faces,Face{3,UInt32}(length(vertices)-4,length(vertices)-1,length(vertices)-2)) : push!(faces,Face{3,UInt32}(length(vertices)-1,length(vertices)-4,length(vertices)-2))
        flip ? push!(faces,Face{3,UInt32}(length(vertices)-4,length(vertices)-3,length(vertices)-1)) : push!(faces,Face{3,UInt32}(length(vertices)-3,length(vertices)-4,length(vertices)-1))
        push!(lines,length(vertices)-4)
        push!(lines,length(vertices)-2)
        push!(lines,length(vertices)-3)
        push!(lines,length(vertices)-1)
    end

    offset = 0.
    n = 10

    rel = 10.
    sr = SimpleRectangle(0.,0.,rel,.5)

    a = map(Point3f0,convert(Vector{Float64},linspace(0,sr.w,n)),repmat([offset+sr.h],n),repmat([0],n))
    b = map(Point3f0,convert(Vector{Float64},linspace(0,sr.w,n)),repmat([offset],n),repmat([0],n))

    vals = linspace(clrmap.mindta.value, clrmap.maxdta.value, n )
    clrmap = apply(clrmap.clrmap, convert( Vector{Float32}, vals ) )

    faces = Face{3,UInt32}[]
    lines = UInt32[]
    vertices = Point3f0[]
    colors = Vector{RGBA{Float32}}()

    for i = 1:n-1
      add_verts_faces_color!(vertices,faces,lines,colors,a[i],a[i+1],b[i],b[i+1],clrmap[i],clrmap[i+1],false)
    end

    return vertices,faces,lines,colors,vals,b
end

function ColorBar(Id::String, clrmap::BoundedColorMap, alpha=0.4)

    vertices,faces,lines,colors,vals,b = ColorBar(clrmap)
    unity_vertices = map(x->UnityVector3(x[1]-8000,x[2]-8000,x[3]-105), vertices )
    unity_faces = vcat([[UInt32(x[1]),UInt32(x[2]),UInt32(x[3])] for x in faces]...)
    unity_colors = map(x->RGBA{Float32}(x.r,x.g,x.b,alpha),colors)
    text = map((x,y)->UnityText(@sprintf("%.3e",x),y-GeometryTypes.Point3f0(7999.6,8000.1,105),UnityVector3(.1,.1,.05),UnityVector3(1,180,1)), vals,b)
    #plot(ub,Id,unity_vertices,unity_faces,lines,unity_colors,ub.config[:wireframe_options],text)
    return UnityMesh(Id,unity_vertices,Int32[],Int32[],unity_faces,unity_colors,["surface_shader = transparent"],text)

end

end #module Unity


text:
# Early stopping
"""
    EarlyStopping(criteria...; kwargs...)
    EarlyStopping(n)

Stop training early when `criteria` are met. See [EarlyStopping.jl]
(https://github.com/ablaom/EarlyStopping.jl) for available stopping
criteria.

Passing an integer `n` uses the simple patience criterion: stop if the
validation loss hasn't increased for `n` epochs.

See [EarlyStopping.jl](https://github.com/JuliaAI/EarlyStopping.jl) for
documentation of all available criteria.

You can control which phases are taken to measure the out-of-sample loss
and the training loss with keyword arguments `trainphase` (default
[`AbstractTrainingPhase`](#)) and `testphase` (default [`AbstractValidationPhase`](#)).

## Examples

```julia
Learner(model, data, optimizer, lossfn, EarlyStopping(3))
```

```julia
import FluxTraining.ES: Disjunction, InvalidValue, TimeLimit

callback = EarlyStopping(Disjunction(InvalidValue(), TimeLimit(0.5)))
Learner(model, data, optimizer, lossfn, callback)
```
"""
mutable struct EarlyStopping <: Callback
    criterion::ES.StoppingCriterion
    state
    testphase::Type{<:Phase}
    trainphase::Type{<:Phase}
end

function EarlyStopping(
        criterion;
        testphase = AbstractValidationPhase,
        trainphase = AbstractTrainingPhase)
    if (testphase isa trainphase) || (trainphase isa testphase)
        error("`trainphase` and `testphase` must not be subtypes of one another")
    end
    return EarlyStopping(criterion, nothing, testphase, trainphase)
end

function EarlyStopping(n::Int; kwargs...)
    return EarlyStopping(ES.Patience(n); kwargs...)
end

Base.show(io::IO, cb::EarlyStopping) = print(io, "EarlyStopping(", cb.criterion, ")")

function on(::EpochEnd, phase::Phase, cb::EarlyStopping, learner)
    loss = last(learner.cbstate.metricsepoch[phase], :Loss)[2]

    if phase isa cb.testphase
        if isnothing(cb.state)
            cb.state = ES.EarlyStopping.update(cb.criterion, loss)
        else
            cb.state = ES.EarlyStopping.update(cb.criterion, loss, cb.state)
        end
    elseif phase isa cb.trainphase
        if isnothing(cb.state)
            cb.state = ES.EarlyStopping.update_training(cb.criterion, loss)
        else
            cb.state = ES.EarlyStopping.update_training(cb.criterion, loss, cb.state)
        end
    end
    if !isnothing(cb.state) && ES.EarlyStopping.done(cb.criterion, cb.state)
        throw(CancelFittingException(ES.EarlyStopping.message(cb.criterion, cb.state)))
    end
end

stateaccess(::EarlyStopping) = (cbstate = (metricsepoch = Read(),),)
runafter(::EarlyStopping) = (Metrics,)


text:
@testset "visualize_solution tests" begin

    file_path = joinpath(@__DIR__,"..", "examples/solver.jl")
    include(file_path)

    params = Dict{String, Any}(
        "num_nodes" => 5,
        "instance" => 1,
        "tol_zero" => 1E-4,
        "tol_psd" => 1E-3,
        "eigen_cuts_full" => true,
        "topology_flow_cuts" => true,
        "lazycuts_logging" => true
    )

    data = LOpt.get_data(params)
    model_lopt  = LOpt.build_LOModel(data)
    result_lopt = LOpt.optimize_LOModel!(model_lopt, optimizer = glpk_optimizer)

    LOpt.visualize_solution(result_lopt, data, visualizing_tool = "tikz", plot_file_format = "tex")
    LOpt.visualize_solution(result_lopt, data, visualizing_tool = "tikz", plot_file_format = "tex", display_edge_weights=true)
    LOpt.visualize_solution(result_lopt, data, visualizing_tool = "graphviz")
    LOpt.visualize_solution(result_lopt, data, visualizing_tool = "graphviz", display_edge_weights=false)

    @test result_lopt["termination_status"] == MOI.OPTIMAL
    @test result_lopt["primal_status"] == MOI.FEASIBLE_POINT

    # Nothing more to test for coverage since the visualizaiton depends on exterior packages. So, I believe the dot files generated are accurate.
end


text:
#!/usr/bin/env julia

using Luxor

using Test

using Random
Random.seed!(42)

function test_rulers(x, y, rot, w)
    gsave()
    translate(x, y)
    rotate(rot)
    box(0, 0, w, w, :clip)
    rulers()
    clipreset()
    grestore()
end

function drawing_png_buffer()
    width, height = 2000, 2000
    Drawing(width, height, :png)
    origin(1200, 1200)
    background("ivory")
    pagetiles = Tiler(width, height, 5, 5, margin=10)
    for (pos, n) in pagetiles
      test_rulers(pos.x, pos.y, rand() * 2pi, pagetiles.tilewidth)
    end

    println("...finished rulers in png buffer test")
    @test finish() == true
end

function drawing_svg_buffer()
    width, height = 2000, 2000
    Drawing(width, height, :svg)
    origin(1200, 1200)
    background("ivory")
    pagetiles = Tiler(width, height, 5, 5, margin=10)
    for (pos, n) in pagetiles
      test_rulers(pos.x, pos.y, rand() * 2pi, pagetiles.tilewidth)
    end

    println("...finished rulers in svg buffer test")
    @test finish() == true
end

function drawing(fname)
    width, height = 2000, 2000
    Drawing(width, height, fname)
    origin(1200, 1200)
    background("ivory")
    pagetiles = Tiler(width, height, 5, 5, margin=10)
    for (pos, n) in pagetiles
      test_rulers(pos.x, pos.y, rand() * 2pi, pagetiles.tilewidth)
    end

    @test finish() == true
    println("...finished rulers test: output in $(fname)")
end

drawing_png_buffer()
drawing_svg_buffer()
drawing("rulers-test.pdf")


text:
"""
Container for mocks and bookkeeping data.
Instances of this type are aware of the call depth ([`current_depth`](@ref)) and the current function/module ([`current_function`](@ref)/[`current_module`](@ref)).

All filter functions take a single argument of this type.
"""
struct Metadata{B}
    mocks::Dict{<:Tuple, <:Any}
    methods::Set{<:Tuple}
    filters::Vector{<:Function}
    mods::Vector{Module}
    funcs::Vector{Any}

    Metadata(mocks::Dict{<:Tuple}, filters::Vector{<:Function}) =
        new{!isempty(filters)}(mocks, Set(keys(mocks)), filters, [Main], [nothing])
end

"""
    current_depth(::Metadata) -> Int

Return the current call depth (the size of the call stack).
The depth is always positive, so the first function entered has a depth of 1.
"""
current_depth(m::Metadata) = length(m.funcs) - 1

"""
    current_function(::Metadata) -> Any

Return the current function (or other callable thing).
In this case, "current" refers, somewhat counterintuitively, not to the function about to be called, but to the function that is about to call it.

To illustrate:

```julia
f(x) = x
g(x) = f(x)
x = g(1)
```

In this case, when the call to `f` is reached, the function that is calling `f` is `g`.
Therefore, the "current" function is `g`.
"""
current_function(m::Metadata) = m.funcs[end-1]

"""
    current_module(m::Metadata) -> Module

Return the current module, where "current" has the same definition as in [`current_function`](@ref).
"""
current_module(m::Metadata) = m.mods[end-1]

# Ensure that the current state satisfies all filters.n
should_mock(m::Metadata{false}, method::Tuple) = method in m.methods
should_mock(m::Metadata, method::Tuple) = method in m.methods && all(f -> f(m), m.filters)

# Update the function/module stacks.
function update!(m::Metadata, ::typeof(prehook), @nospecialize(f), @nospecialize(args...))
    f_uw = unwrap_fun(f)
    Ts = Tuple{map(typeof, args)...}

    mod = if f_uw isa Builtin || !hasmethod(f_uw, Ts)
        parentmodule(f_uw)
    else
        parentmodule(f_uw, Ts)
    end

    push!(m.funcs, f_uw)
    push!(m.mods, mod)
end

function update!(m::Metadata, ::typeof(posthook), @nospecialize(args...))
    pop!(m.funcs)
    pop!(m.mods)
end

# If a function is a keyword wrapper, try to get the wrapped function.
# This garbage is the result of random experimentation and is very sketchy.
# It fails in the case of closures.
unwrap_fun(f) = f
unwrap_fun(f::Builtin) = f
function unwrap_fun(f::F) where F <: Function
    startswith(string(F.name.name), "#kw##") || return f

    name = Symbol(string(F.name.name)[6:end])
    name_hash = Symbol("#", name)
    mod = F.name.module

    return if isdefined(mod, name)
        getfield(mod, name)
    elseif isdefined(mod, name_hash)
        gf = getfield(mod, name_hash)
        isdefined(gf, :instance) ? gf.instance : f
    else
        f
    end
end


text:
using Documenter
using TensorOperations

makedocs(modules=[TensorOperations],
            sitename = "TensorOperations.jl",
            authors = "Jutho Haegeman",
            format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
            pages = [
                "Home" => ["index.md",
                            "indexnotation.md",
                            "functions.md",
                            "cache.md",
                            "implementation.md"]
            ])

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(repo = "github.com/Jutho/TensorOperations.jl.git")


text:
function default_in_partition(sites::Tuple{Int}, p::Integer, nparts::Integer)
  return p == mod1(sites[1], nparts)
end

# i ≥ j
_f(i, j) = (i - 1) * i ÷ 2 + j
function default_in_partition(sites::NTuple{2,Int}, p, nparts)
  i, j = sites
  if i ≤ j
    return p == mod1(_f(j, i), nparts)
  end
  return p == mod1(_f(i, j), nparts)
end

function default_in_partition_sorted(sites::NTuple{4,Int}, p, nparts)
  i, j, k, l = sites
  if j == k
    return default_in_partition((j,), p, nparts)
  end
  return default_in_partition((i, j), p, nparts)
end

function default_in_partition(sites::NTuple{4,Int}, p, nparts)
  return default_in_partition_sorted(sort(sites), p, nparts)
end

function _sites(t)
  return Tuple((only(ITensors.site(o)) for o in ITensors.terms(t)))
end

function partition(os::OpSum, nparts::Integer; in_partition=default_in_partition)
  oss = [OpSum() for _ in 1:nparts]
  for p in 1:nparts
    for n in 1:length(os)
      sites = _sites(os[n])
      if in_partition(sites, p, nparts)
        add!(oss[p], os[n])
      end
    end
  end
  return oss
end


text:
"""
    BiweightMidcovariance(; c=9.0, modify_sample_size=false)

The biweight midcovariance is a covariance estimator that is resilient to outliers.

The technique derives originally from astrophysics [1], and is implemented in the Python
module [Astropy](https://www.astropy.org/) [2], as well as in NIST's Dataplot [3].

Consider two random variables ``x`` and ``y``, for which we have ``n`` observations
``\\{(x_i, y_i)\\}``. The biweight midcovariance is then defined to be:
```math
n_s\\cdot\\frac{
    \\sum_{|u_i|<1,|v_i|<1}(x_i - M_x)(1 - u_i^2)^2(y_i - M_y)(1 - v_i^2)^2
}{
    \\left(\\sum_{|u_i|<1}(1 - u_i^2)(1-5u_i^2)\\right)
    \\left(\\sum_{|v_i|<1}(1 - v_i^2)(1-5v_i^2)\\right)
}
```
where ``n_s`` is the sample size, ``M_x`` and ``M_y`` are the medians of ``\\{x_i\\}`` and
``\\{y_i\\}`` respectively, and
```math
\\begin{aligned}
u_i &= \\frac{x_i - M_x}{c \\cdot \\mathrm{MAD}_x} \\\\
v_i &= \\frac{y_i - M_y}{c \\cdot \\mathrm{MAD}_y},
\\end{aligned}
```
where ``\\mathrm{MAD}`` represents the median absolute deviation,
```math
\\begin{aligned}
\\mathrm{MAD}_x &= \\mathrm{median}(\\left\\{\\left|x_i - M_x\\right|\\right\\}) \\\\
\\mathrm{MAD}_y &= \\mathrm{median}(\\left\\{\\left|y_i - M_y\\right|\\right\\}).
\\end{aligned}
```

If either ``\\mathrm{MAD}_x = 0`` or ``\\mathrm{MAD}_y = 0``, the pairwise covariance is
defined to be zero.

The parameter ``c`` is a tuning constant, for which the default is ``9.0``.
Larger values will reduce the number of outliers that are removed — i.e. reducing
robustness, but increasing sample efficiency.

# Fields
- `c::Float64`: The tuning constant corresponding to ``c`` above.
- `modify_sample_size::Bool`: If `false`, then we use a sample size ``n_s`` equal to the
    total number of observations ``n``. This is consistent with the standard definition of
    biweight midcovariance in the literature. Otherwise, we count only those elements which
    are not rejected as outliers in the numerator, i.e. those for which ``|u_i|<1``
    and ``|v_i|<1``.
    This follows the implementation in astropy [2].

# Complexity
- Space: ``O(p^2)``
- Time: ``O(np^2)``

# References
[1] Beers, Flynn, and Gebhardt (1990; AJ 100, 32)
"Measures of Location and Scale for Velocities in Clusters of Galaxies -- A Robust Approach"

[2] [Astropy biweight_midcovariance](https://docs.astropy.org/en/stable/api/astropy.stats.biweight_midcovariance.html)

[3] [NIST Dataplot biweight midcovariance](https://www.itl.nist.gov/div898/software/dataplot/refman2/auxillar/biwmidc.htm)
"""
struct BiweightMidcovariance <: CovarianceEstimator
    c::Float64
    modify_sample_size::Bool
end
function BiweightMidcovariance(; c::Real=9.0, modify_sample_size::Bool=false)
    c > 0 || throw(ArgumentError("c must be positive, got $c"))
    return BiweightMidcovariance(Float64(c), modify_sample_size)
end

function covzm(ce::BiweightMidcovariance, x::AbstractVector{<:Real})
    MADx = median!(abs.(x))
    numerator = zero(eltype(x))

    # If MADx is zero, return zero.
    iszero(MADx) && return zero(one(eltype(x)) / 1)

    denominator = zero(eltype(x))
    count = 0
    for xi in x
        ui² = (xi / (ce.c * MADx))^2
        # Benchmarking suggests that having this branch is similar or faster to a form
        # using `ifelse`, where we always compute the updates but avoid the branch.
        ui² >= 1 && continue
        count += 1
        numerator += xi^2 * (1 - ui²)^4
        denominator += (1 - ui²) * (1 - 5 * ui²)
    end
    n = ifelse(ce.modify_sample_size, count, length(x))
    return n * numerator / (denominator^2)
end

function cov(
    ce::BiweightMidcovariance,
    x::AbstractVector{<:Real};
    mean::Union{Nothing,<:Real}=nothing,
)
    Mx = mean === nothing ? median(x) : mean
    return covzm(ce, x .- Mx)
end

function covzm(
    ce::BiweightMidcovariance, x::AbstractVector{<:Real}, y::AbstractVector{<:Real}
)
    MADx = median!(abs.(x))
    MADy = median!(abs.(y))

    # Promote types between x & y for numerator
    numerator = zero(promote_type(eltype(x), eltype(y)))

    # If either of the MADs are zero, return zero.
    #   `numerator` is always zero here, and of the correct type.
    iszero(MADx) && return numerator
    iszero(MADy) && return numerator

    denominator_x = zero(eltype(x))
    denominator_y = zero(eltype(y))
    count = 0
    for (xi, yi) in zip(x, y)
        ui² = (xi / (ce.c * MADx))^2
        vi² = (yi / (ce.c * MADy))^2

        # Compute the updates, and increment if desired.
        # This way should avoid branches, and benchmarks suggest it's 1-2% faster than
        # using conditionals.
        ax = (1 - ui²) * (1 - 5 * ui²)
        ay = (1 - vi²) * (1 - 5 * vi²)
        anum = xi * (1 - ui²)^2 * yi * (1 - vi²)^2

        denominator_x += ifelse(ui² < 1, ax, zero(denominator_x))
        denominator_y += ifelse(vi² < 1, ay, zero(denominator_y))
        i_count, i_numerator = ifelse(ui² < 1 && vi² < 1, (1, anum), (0, zero(numerator)))
        count += i_count
        numerator += i_numerator
    end

    n = ifelse(ce.modify_sample_size, count, length(x))
    return n * numerator / (denominator_x * denominator_y)
end

function cov(
    ce::BiweightMidcovariance, x::AbstractVector{<:Real}, y::AbstractVector{<:Real}
)
    if size(x) != size(y)
        throw(ArgumentError("x & y have different sizes, $(size(x)) & $(size(y))"))
    end

    return covzm(ce, x .- median(x), y .- median(y))
end

function cov(
    ce::BiweightMidcovariance,
    X::AbstractMatrix{<:Real};
    mean::Union{Nothing,AbstractVector{<:Real},AbstractMatrix{<:Real}}=nothing,
    dims::Int=1,
)
    dims ∈ (1, 2) || throw(ArgumentError("Argument dims can only be 1 or 2 (given: $dims)"))

    # Ascertain the number of observations `n` and number of variables `p`.
    n, p = dims == 1 ? size(X) : (size(X, 2), size(X, 1))

    # Temporary storage to use when computing medians.
    temp = Vector{eltype(X)}(undef, n)

    # If the `mean` argument isn't provided, then we centralise with the median.
    # After this statement, Xc is standardised to be of shape (n, p)
    Xc = if mean === nothing
        # The following is equivalent to
        #   X .- median(X; dims=dims)
        # However, this implementation significantly reduces the runtime of this function.
        # This is primarily achieved by reducing the number of allocations -- we allocate a
        # single temporary buffer which we allow `median!` to mutate, and then copy the next
        # section of data into it.
        X = dims == 2 ? transpose(X) : X
        Xc = copy(X)
        @inbounds for i in 1:p
            # Avoiding copyto!(, ... @view(X[:, i])), because this can allocate a view if
            # the call to copyto! isn't inlined.
            for j in 1:n
                temp[j] = X[j, i]
            end
            Xc[:, i] .-= median!(temp)
        end
        Xc
    else
        Xc = X .- mean
        # Standardise the orientation of Xc.
        dims == 2 ? transpose(Xc) : Xc
    end

    # Compute all the median absolute deviations.
    # To avoid extra allocations, use a temporary buffer `temp` which will be populated with
    # the absolute values of the observations of each variable, and which we then allow
    # median! to mutate.
    T = typeof(one(eltype(X)) / 1)  # This is the output type of median.
    MAD = Matrix{T}(undef, 1, p)
    @inbounds for i in 1:p
        # Avoiding map!(abs, temp, @view(Xc[:, i])), because this can allocate a view if
        # the call to map! isn't inlined.
        for j in 1:n
            temp[j] = abs(Xc[j, i])
        end
        MAD[1, i] = median!(temp)
    end

    u² = @. (Xc / (ce.c * MAD))^2

    onesubu² = max.((1 .- u²), 0)
    a = @. Xc * onesubu²^2
    numerator = a' * a

    a[:, :] .= @. onesubu² * (1 - 5 * u²)  # Re-use `a` memory to avoid an allocation.
    b = sum(a; dims=1)
    denominator = b' * b

    n = if ce.modify_sample_size
        mask = onesubu² .> 0
        mask' * mask
    else
        n
    end

    # Re-use the numerator memory for the result.
    result = numerator
    result .*= n ./ denominator

    # Whereever the MAD is zero, we should set the result to zero.
    mad_zero = dropdims(iszero.(MAD); dims=1)
    result[mad_zero, :] .= zero(eltype(result))
    result[:, mad_zero] .= zero(eltype(result))

    return result
end


text:
@testset "KPLS Pediction Tests (in sample)" begin

    @testset "Test KPLS Single Non Linear Target" begin

        Random.seed!(1)

        z(x)     = 4.26 * (exp.(-x) - 4 * exp.(-2.0*x) + 3 * exp.(-3.0*x))
        x_values = Array(range(0.0,step=3.5,length=100))
        z_pure   = z(x_values)
        noise    = Random.randn(100)
        z_noisy  = z_pure + noise
        X        = MLJ.table(collect(x_values)[:,:])
        Y        = z_noisy #[:,:] #z_pure

        pls_pipe       = MLJ.@pipeline prediction_type=:deterministic MLJ.Standardizer() PLSRegressor.KPLS(n_factors=1,kernel="rbf",width=0.01)  target = MLJ.Standardizer()

		pls_machine    = MLJ.machine(pls_pipe, X, Y)

        train = range(1,stop=length(X))
        MLJ.fit!(pls_machine, rows=train,force=true)
        yhat = MLJ.predict(pls_machine, rows=train);
        @test MLJ.mae(yhat, Y[train]) |> mean < 1e-2

    end

    @testset "Test KPLS Single Target (Linear Target)" begin

        X        = MLJ.table([1 2; 2 4; 4.0 6])
        Y        = [-2; -4; -6.0] #[:,:]

        pls_pipe       = MLJ.@pipeline prediction_type=:deterministic MLJ.Standardizer() PLSRegressor.KPLS(n_factors=1,kernel="rbf",width=0.01)  target = MLJ.Standardizer()

		pls_machine    = MLJ.machine(pls_pipe, X, Y)

        train = range(1,stop=length(X))
        MLJ.fit!(pls_machine, rows=train,force=true)
        yhat = MLJ.predict(pls_machine, rows=train);
        @test MLJ.mae(yhat, Y[train]) |> mean < 1e-2

        X        = MLJ.table([1 2; 2 4; 4.0 6])
        Y        = [2; 4; 6.0]#[:,:]

        pls_pipe       = MLJ.@pipeline prediction_type=:deterministic MLJ.Standardizer() PLSRegressor.KPLS(n_factors=1,kernel="rbf",width=0.01)  target = MLJ.Standardizer()

		pls_machine    = MLJ.machine(pls_pipe, X, Y)

        train = range(1,stop=length(X))
        MLJ.fit!(pls_machine, rows=train,force=true)
        yhat = MLJ.predict(pls_machine, rows=train);
        @test MLJ.mae(yhat, Y[train]) |> mean < 1e-2

    end

    @testset "Test KPLS Multiple Target (Linear Target)" begin

        X        = MLJ.table([1; 2; 3.0][:,:])
        Y        = MLJ.table([1 1; 2 2; 3 3.0]) #[:,:]

        pls_pipe       = MLJ.@pipeline prediction_type=:deterministic MLJ.Standardizer() PLSRegressor.KPLS(n_factors=1,kernel="rbf",width=0.01)  target = MLJ.Standardizer()

		pls_machine    = MLJ.machine(pls_pipe, X, Y)

        train = range(1,stop=length(X))
        MLJ.fit!(pls_machine, rows=train,force=true)
        yhat = MLJ.predict(pls_machine, rows=train);
        @test abs.(MLJ.matrix(yhat) .- MLJ.matrix(Y)[train,:]) |> mean < 1e-6

        X        = MLJ.table([1; 2; 3.0][:,:])
        Y        = MLJ.table([1 -1; 2 -2; 3 -3.0])#[:,:]

        pls_pipe       = MLJ.@pipeline prediction_type=:deterministic MLJ.Standardizer() PLSRegressor.KPLS(n_factors=1,kernel="rbf",width=0.01)  target = MLJ.Standardizer()

		pls_machine    = MLJ.machine(pls_pipe, X, Y)

        train = range(1,stop=length(X))
        MLJ.fit!(pls_machine, rows=train,force=true)
        yhat = MLJ.predict(pls_machine, rows=train);
        @test abs.( MLJ.matrix(yhat) .- MLJ.matrix(Y)[train,:]) |> mean < 1e-6

        @testset "Linear Prediction Tests " begin

        X        = MLJ.table([1 2; 2 4; 4 6.0])
        Y        = MLJ.table([4 2;6 4;8 6.0])#[:,:]

        pls_pipe       = MLJ.@pipeline prediction_type=:deterministic MLJ.Standardizer() PLSRegressor.KPLS(n_factors=1,kernel="rbf",width=0.01)  target = MLJ.Standardizer()

		pls_machine    = MLJ.machine(pls_pipe, X, Y)

        train = range(1,stop=length(X))
        MLJ.fit!(pls_machine, rows=train,force=true)
        yhat = MLJ.predict(pls_machine, rows=train);
        @test abs.( MLJ.matrix(yhat) .- MLJ.matrix(Y)[train,:]) |> mean < 1e-6

        X           = MLJ.table([1 -2; 2 -4; 4 -6.0])
        Y           = MLJ.table([-4 -2;-6 -4;-8 -6.0])#[:,:]

        pls_pipe       = MLJ.@pipeline prediction_type=:deterministic MLJ.Standardizer() PLSRegressor.KPLS(n_factors=1,kernel="rbf",width=0.01)  target = MLJ.Standardizer()

		pls_machine    = MLJ.machine(pls_pipe, X, Y)

        train = range(1,stop=length(X))
        MLJ.fit!(pls_machine, rows=train,force=true)
        yhat = MLJ.predict(pls_machine, rows=train);
        @test abs.(MLJ.matrix(yhat) .- MLJ.matrix(Y)[train,:]) |> mean < 1e-6

        end

    end

end;


text:
# Reference: Temporal Distance Metrics for Social Network Analysis,
#  John Tang, Mirco Musolesi, Cecilia Mascolo and Vito Latora,
# 2009.

# find the shortest temporal path via depth first search
function _DFS_shortest_temporal_path(g::AbstractEvolvingGraph,
                                     v1::Tuple,  # v1 and v2 are active nodes (v,t)
                                     v2::Tuple,
                                     path = TemporalPath(),
                                     shortest =  Union{};
                                     verbose::Bool = false)
    path = deepcopy(path)
    path.walks = [path.walks; v1]

    if verbose
        println("current:", path)
    end

    if v1 == v2
        return path
    end
    for node in forward_neighbors(g, v1)
        if !(has_node(path, node)) # avoid cycles
            if node[2] <= v2[2] # avoid searching nodes at timestamps > v2[2]
                if (shortest ==  Union{}) || (spatial_length(path) < spatial_length(shortest))
                    newPath = _DFS_shortest_temporal_path(g, node, v2, path, shortest, verbose = verbose)
                    if newPath !=  Union{}
                        shortest = newPath
                    end
                end
            end
        end
    end
    return shortest
end

"""
        shortest_temporal_path(g, v1, t1, v2, t2 [, verbose = false])

Find the shortest temporal path from node `v1` at timestamp `t1`
to node `v2` at timestamp `t2` on the evolving graph `g`.
If `verbose = true`, prints the current path at each search step.
"""
shortest_temporal_path(g::AbstractEvolvingGraph, v1, t1, v2, t2; verbose::Bool = false) =    _DFS_shortest_temporal_path(g, (v1, t1), (v2, t2), verbose = verbose)
shortest_temporal_path{V,  T}(g::EvolvingGraph{V, T},
                                                        v1::V, t1::T, v2::V, t2::T;
                                                        verbose::Bool = false) =
   _DFS_shortest_temporal_path(g, (v1, t1), (v2, t2), verbose = verbose)

function shortest_temporal_path{V, T}(g::EvolvingGraph{V, T}, v1, t1, v2, t2;
                                                                        verbose::Bool = false)
    v1 = find_node(g, v1)
    v2 = find_node(g, v2)
    shortest_temporal_path(g, v1, T(t1), v2, T(t2), verbose = verbose)
end

"""
        shortest_temporal_distance(g, v1, t1, v2, t2)

Find the shortest temporal distance from node `v1` at timestamp `t1` to
node `v2` at timestamp `t2` on the evolving graph `g`.
"""
function shortest_temporal_distance(g::AbstractEvolvingGraph, v1, t1, v2, t2)
    if shortest_temporal_path(g, v1, t1, v2, t2) == Union{}
        return Inf
    else
        return length(shortest_temporal_path(g, v1, t1, v2, t2)) - 1
    end
end

"""
        temporal_efficiency(g, v1, t1, v2, t2)

Find the temporal efficiency from node `v1` at timestamp `t1` to node `v2` at timestamp
`t2` on the evolving graph `g`.
"""
temporal_efficiency(g::AbstractEvolvingGraph, v1, t1, v2, t2) =
         one(Float64)/shortest_temporal_distance(g, v1, t1, v2, t2)

"""
        global_temporal_efficiency(g, t1, t2)

Return the global temporal efficiency of the evolving graph `g` between
timestamp `t1` and `t2`. The global temporal efficiency is a measure of
how well information flow between two given timestamps.
"""
function global_temporal_efficiency(g::AbstractEvolvingGraph, t1, t2)
    result = zero(Float64)
    for j in nodes(g)
        for i in nodes(g)
            result += temporal_efficiency(g, i,t1, j,t2)
        end
    end
    result
end


text:
module SpringCollab2020NonBadelineMovingBlock

using ..Ahorn, Maple

@pardef NonBadelineMovingBlock(x1::Integer, y1::Integer, x2::Integer=x1+16, y2::Integer=y1, width::Integer=Maple.defaultBlockWidth, height::Integer=Maple.defaultBlockHeight, tiletype::String="g", highlightTiletype::String="G") = Entity("SpringCollab2020/nonBadelineMovingBlock", x=x1, y=y1, nodes=Tuple{Int, Int}[(x2, y2)], width=width, height=height, tiletype=tiletype, highlightTiletype=highlightTiletype)

const placements = Ahorn.PlacementDict(
    "Non-Badeline Boss Moving Block (Spring Collab 2020)" => Ahorn.EntityPlacement(
        NonBadelineMovingBlock,
        "rectangle",
        Dict{String, Any}(),
        function(entity)
            entity.data["nodes"] = [(Int(entity.data["x"]) + Int(entity.data["width"]) + 8, Int(entity.data["y"]))]
        end
    )
)

Ahorn.editingOptions(entity::NonBadelineMovingBlock) = Dict{String, Any}(
    "tiletype" => Ahorn.tiletypeEditingOptions(),
    "highlightTiletype" => Ahorn.tiletypeEditingOptions()
)

Ahorn.nodeLimits(entity::NonBadelineMovingBlock) = 1, 1
Ahorn.minimumSize(entity::NonBadelineMovingBlock) = 8, 8
Ahorn.resizable(entity::NonBadelineMovingBlock) = true, true

function Ahorn.selection(entity::NonBadelineMovingBlock)
    if entity.name == "SpringCollab2020/nonBadelineMovingBlock"
        x, y = Ahorn.position(entity)
        nx, ny = Int.(entity.data["nodes"][1])

        width = Int(get(entity.data, "width", 8))
        height = Int(get(entity.data, "height", 8))

        return [Ahorn.Rectangle(x, y, width, height), Ahorn.Rectangle(nx, ny, width, height)]
    end
end

function Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::NonBadelineMovingBlock, room::Maple.Room)
    Ahorn.drawTileEntity(ctx, room, entity, material=get(entity.data, "tiletype", "g")[1], blendIn=false)
end

function Ahorn.renderSelectedAbs(ctx::Ahorn.Cairo.CairoContext, entity::NonBadelineMovingBlock, room::Maple.Room)
    x, y = Ahorn.position(entity)
    nodes = get(entity.data, "nodes", ())

    width = Int(get(entity.data, "width", 8))
    height = Int(get(entity.data, "height", 8))

    if !isempty(nodes)
        nx, ny = Int.(nodes[1])
        cox, coy = floor(Int, width / 2), floor(Int, height / 2)

        # Use 'G' instead of 'g', as that is the highlight color of the block (the active color)
        fakeTiles = Ahorn.createFakeTiles(room, nx, ny, width, height, get(entity.data, "highlightTiletype", "G")[1], blendIn=false)
        Ahorn.drawFakeTiles(ctx, room, fakeTiles, room.objTiles, true, nx, ny, clipEdges=true)
        Ahorn.drawArrow(ctx, x + cox, y + coy, nx + cox, ny + coy, Ahorn.colors.selection_selected_fc, headLength=6)
    end
end

end


text:
abstract type AbstractGroupRegressor <: MMI.Deterministic end
abstract type AbstractGroupRidgeRegressor <: AbstractGroupRegressor end

function StatsBase.fit(grp_ridge::AbstractGroupRidgeRegressor, X, Y, grp::GroupedFeatures)
    decomposition = grp_ridge.decomposition
    tuning = _main_hyperparameter_value(grp_ridge)
    nobs = length(Y)
    if decomposition === :default
        decomposition = (nfeatures(grp) <= 4*nobs) ? :cholesky : :woodbury
    end

    if decomposition === :cholesky
        pred = CholeskyRidgePredictor(X)
    elseif decomposition === :woodbury
        pred = WoodburyRidgePredictor(X)
    else
        "Only :default, :cholesky and :woodbury currently supported"
    end
    workspace = BasicGroupRidgeWorkspace(X = X, Y = Y, groups = grp, XtXpΛ_chol = pred)
    StatsBase.fit!(workspace, tuning)
    workspace
end


text:
abstract type RpTAlgorithm end

function get_rpt_alg()
    implemented = [RpTa()]
    N = Options.get("pt_alg")
    try
        return implemented[N]
    catch BoundsError
        throw(FermiException("implementation number $N not available for RCCSD(T)."))
    end
end
"""
    Fermi.CoupledCluster.RCCSDpT

Fermi struct that holds information about RCCSD(T) wavefunctions

_struct tree:_

**RCCSD(T)** <: AbstractCCWavefunction <: AbstractCorrelatedWavefunction <: AbstractWavefunction
"""
struct RCCSDpT{T} <: AbstractCCWavefunction
    CCSD::RCCSD{T}
    energy::T
    correction::T
end

function RCCSDpT(x...)
    if !any(i-> i isa RpTAlgorithm, x)
        RCCSDpT(x..., get_rpt_alg())
    else
        # Print the type of arguments given for a better feedback
        args = "("
        for a in x[1:end-1]
            args *= "$(typeof(a)), "
        end
        args = args[1:end-2]*")"
        throw(FermiException("invalid arguments for RCCSD(T) method: $args"))
    end
end

#implementations
struct RpTa <: RpTAlgorithm end
include("RpTa.jl")


text:
```
PowerSpectra()
```

max_stars_repo_name: yukki-jpn/JuliaTSA

max_issues_repo_name: yukki-jpn/JuliaTSA

max_forks_repo_name: yukki-jpn/JuliaTSA


text:
"""
    cartesian2jacobi!(r, r_j, m, m_sum)

Transforms cartesian coordinates into jacobi coordinates.

# Args

* `r`: NxD matrix where N is the number of bodies and D is the number of dimensions, each row contains the positions of the bodies in cartesian coordinates.
* `r_j`: NxD matrix in which the results will be saved.
* `m`: array of N elements, where element number i contains the mass of i-th body.
* `m_sum`: array of N elements, where element with index i contains the mass of the bodies from 1 to i
```math
m_{sum_{i}} = \\sum_{j = 1}^{i} m_{j} \\\\
```
"""
function cartesian2jacobi!(r, r_j, m, m_sum)
    R = m[1]*r[1,:]
    for i in 2:size(r)[1]
        r_j[i,:] = r[i,:] - (R/m_sum[i-1])
        R = R*(1.0+m[i]/m_sum[i-1]) + m[i]*r_j[i,:]
    end
    r_j[1,:] = R/m_sum[size(r)[1]]
end

"""
    jacobi2cartesian!(r, r_j, m, m_sum)

Transforms jacobi coordinates into cartesian coordinates.

# Args

* `r`: NxD matrix in which the results will be saved.
* `r_j`: NxD matrix where N is the number of bodies and D is the number of dimensions, each row contains the positions of the bodies in jacobi coordinates.
* `m`: array of N elements, where element number i contains the mass of i-th body.
* `m_sum`: array of N elements, where element with index i contains the mass of the bodies from 1 to i
```math
m_{sum_{i}} = \\sum_{j = 1}^{i} m_{j} \\\\
```
"""
function jacobi2cartesian!(r, r_j, m, m_sum)
    R = r_j[1,:]*m_sum[size(r)[1]]
    for i in size(r)[1]:-1:2
        R = (R - m[i]*r_j[i,:])/m_sum[i]
        r[i,:] = r_j[i,:] + R
        R *= m_sum[i-1]
    end
    r[1,:] = R/m[1]
end

max_stars_repo_name: salanueva/KeplerFlow

max_issues_repo_name: salanueva/KeplerFlow

max_forks_repo_name: salanueva/KeplerFlow


text:
using ITensors
using Printf

function ITensors.op(::OpName"expτSS", ::SiteType"S=1/2", s1::Index, s2::Index; τ)
  h =
    1 / 2 * op("S+", s1) * op("S-", s2) +
    1 / 2 * op("S-", s1) * op("S+", s2) +
    op("Sz", s1) * op("Sz", s2)
  return exp(τ * h)
end

function main(; N=10, cutoff=1E-8, δτ=0.1, beta_max=2.0)

  # Make an array of 'site' indices
  s = siteinds("S=1/2", N; conserve_qns=true)

  # Make gates (1,2),(2,3),(3,4),...
  gates = ops([("expτSS", (n, n + 1), (τ=-δτ / 2,)) for n in 1:(N - 1)], s)
  # Include gates in reverse order to complete Trotter formula
  append!(gates, reverse(gates))

  # Initial state is infinite-temperature mixed state
  rho = MPO(s, "Id") ./ √2

  # Make H for measuring the energy
  terms = OpSum()
  for j in 1:(N - 1)
    terms += 1 / 2, "S+", j, "S-", j + 1
    terms += 1 / 2, "S-", j, "S+", j + 1
    terms += "Sz", j, "Sz", j + 1
  end
  H = MPO(terms, s)

  # Do the time evolution by applying the gates
  # for Nsteps steps
  for β in 0:δτ:beta_max
    energy = inner(rho, H)
    @printf("β = %.2f energy = %.8f\n", β, energy)
    rho = apply(gates, rho; cutoff)
    rho = rho / tr(rho)
  end

  return nothing
end


text:
function accept_rparen(ps)
    if kindof(ps.nt) === Tokens.RPAREN
        return EXPR(next(ps))
    else
        return mErrorToken(ps, EXPR(:RPAREN, 0, 0), UnexpectedToken)
    end
end
accept_rparen(ps::ParseState, args) = push!(args, accept_rparen(ps))

function accept_rsquare(ps)
    if kindof(ps.nt) === Tokens.RSQUARE
        return EXPR(next(ps))
    else
        return mErrorToken(ps, EXPR(:RSQUARE, 0, 0), UnexpectedToken)
    end
end
accept_rsquare(ps::ParseState, args) = push!(args, accept_rsquare(ps))

function accept_rbrace(ps)
    if kindof(ps.nt) === Tokens.RBRACE
        return EXPR(next(ps))
    else
        return mErrorToken(ps, EXPR(:RBRACE, 0, 0), UnexpectedToken)
    end
end
accept_rbrace(ps::ParseState, args) = push!(args, accept_rbrace(ps))

function accept_end(ps::ParseState)
    if kindof(ps.nt) === Tokens.END
        return EXPR(next(ps))
    else
        return mErrorToken(ps, EXPR(:END, 0, 0), UnexpectedToken)
    end
end
accept_end(ps::ParseState, args) = push!(args, accept_end(ps))

function accept_comma(ps)
    if iscomma(ps.nt)
        return EXPR(next(ps))
    else
        return EXPR(:RPAREN, 0, 0)
    end
end
accept_comma(ps::ParseState, args) = push!(args, accept_comma(ps))

function requires_ws(x, ps)
    if x.span == x.fullspan
        return mErrorToken(ps, x, Unknown)
    else
        return x
    end
end

function requires_no_ws(x, ps)
    if !(kindof(ps.nt) === Tokens.RPAREN ||
        kindof(ps.nt) === Tokens.RBRACE ||
        kindof(ps.nt) === Tokens.RSQUARE || kindof(ps.nt) === Tokens.END) && x.span != x.fullspan
        return mErrorToken(ps, x, UnexpectedWhiteSpace)
    else
        return x
    end
end


text:
using Statistics

# include("sindyc.jl")
# include("../utils/regression_utils.jl")

"""
Initialize the control signal using the residual of a linear fit model, dmd.
By default, removes small entries from the residual that may be noise

Normal use will be to run this to get the initial guess, then use 'sra_dmdc()'
to further refine the guess.
"""
function sra_dmd(dat, dat_grad;
                 quantile_threshold=0.95, add_const=false)
    n, m = size(dat)
    #   Note: add a constant field
    if add_const
        dat_and_ctr = vcat(dat, ctr);
    else
        dat_and_ctr = dat;
    end
    AB_guess = dat_grad / dat_and_ctr
    A_guess = AB_guess[:, 1:n]
    B_guess = AB_guess[:, n+1:end]
    # Stay in gradient space
    U_guess = calculate_residual(dat_and_ctr, dat_grad, AB_guess)
    # grad_guess = AB_guess*dat_and_ctr;
    # U_guess = dat_grad .- grad_guess
    # if using_previous_ctr
    #     U_guess = vcat(U_guess, ctr)
    # end
    # plot(U_func.(ts), label="Control Signal")
    # plot!(U_guess[2,:], label="Control signal guess")

    # Remove the smaller elements of U_guess
    sparsify_signal!(U_guess, quantile_threshold=quantile_threshold)

    return (U=U_guess, A=A_guess, B=B_guess)
end

"""
Initialize the control signal using the residual of a NONLINEAR fit model, SINDy.
By default, removes small entries from the residual that may be noise

Normal use will be to run this to get the initial guess, then use 'sra_sindy()'
to further refine the guess.
"""
function sra_sindy(dat, dat_grad;
                 quantile_threshold=0.95, library=nothing)
    n, m = size(dat)
    model = sindyc(dat, dat_grad,
                   library=library,
                   use_lasso=true)
    # Stay in gradient space
    U_guess = calculate_residual(model, dat, dat_grad)

    # Remove the smaller elements of U_guess
    sparsify_signal!(U_guess, quantile_threshold=quantile_threshold)

    return (U_guess=U_guess, model=model)
end

"""
Calculates the residual of the learned DMD dynamics
"""
function calculate_residual(dat, dat_grad, A)
    return dat_grad .- A*dat
end

"""
Calculates the residual of the learned SINDy dynamics
"""
function calculate_residual(m::sindycModel, dat, dat_grad)
    # return dat_grad .- m(dat)
    return dat_grad .- intrinsic_dynamics(m, dat)
end

export sra_dmd, sra_sindy, calculate_residual, sparsify_signal!


text:
module NOAu111

using JuLIP: JuLIP
using JuLIP.Potentials: z2i, ZList
using StaticArrays: SMatrix, SVector
using LinearAlgebra: norm, Hermitian, dot
using DataStructures: DefaultDict
using NQCBase: au_to_ang, eV_per_ang_to_au, eV_to_au
using NQCModels: NQCModels, DiabaticModels

export NOAu

include("parameters.jl")
include("au_au.jl")
include("diabatic_elements.jl")
include("h11.jl")

struct NOAu{V1,V2,V3,V4,V5,A} <: DiabaticModels.DiabaticModel
    n_states::Int
    H00::V1
    H11::V2
    H01::V3
    AuAu::V4
    image_potential::V5
    atoms::A
    Nindex::Int
    Oindex::Int
end

NQCModels.nstates(::NOAu) = 2

NOAu(jatoms, Nindex, Oindex) = NOAu(2, H00(), H11(jatoms), H01(), AuAu(), image(D, C, zimage), jatoms, Nindex, Oindex)

function NOAu(symbols::AbstractVector{Symbol}, cell, R)
    jatoms = JuLIP.Atoms(
        X=zeros(3,length(symbols)),
        M=JuLIP.atomic_mass.(symbols),
        Z=JuLIP.atomic_number.(symbols),
        cell=au_to_ang.(cell.vectors'),
        pbc=cell.periodicity
        )

    JuLIP.set_positions!(jatoms, au_to_ang.(R))
    Nindex = findfirst(isequal(7), jatoms.Z)
    Oindex = findfirst(isequal(8), jatoms.Z)
    NOAu(jatoms, Nindex, Oindex)
end

function NQCModels.potential(model::NOAu, R::AbstractMatrix)

    JuLIP.set_positions!(model.atoms, au_to_ang.(R))
    Au = eV_to_au(JuLIP.energy(model.AuAu, model.atoms))
    V11 = eV_to_au(JuLIP.energy(model.H00, model.atoms)) + Au
    V22 = eV_to_au(JuLIP.energy(model.H11, model.atoms) + ϕ - Eₐ) + Au
    V12 = eV_to_au(JuLIP.energy(model.H01, model.atoms))

    return Hermitian(SMatrix{2,2}(V11, V12, V12, V22))
end

function NQCModels.derivative!(model::NOAu, D::AbstractMatrix{<:Hermitian}, R::AbstractMatrix)

    JuLIP.set_positions!(model.atoms, au_to_ang.(R))
    Au = -JuLIP.forces(model.AuAu, model.atoms)
    D11 = -JuLIP.forces(model.H00, model.atoms) + Au
    D22 = -JuLIP.forces(model.H11, model.atoms) + Au
    D12 = -JuLIP.forces(model.H01, model.atoms)

    for i=1:length(model.atoms)
        for j=1:3
            d11 = eV_per_ang_to_au(D11[i][j])
            d12 = eV_per_ang_to_au(D12[i][j])
            d22 = eV_per_ang_to_au(D22[i][j])
            D[j,i] = Hermitian(SMatrix{2,2}(d11, d12, d12, d22))
        end
    end

    return D
end

end