text:
module MCOPlots

    # ----------------------------------------------------------------------------------
    #   Imports
    # ----------------------------------------------------------------------------------
    using MonteCarloObservable; const MCO = MonteCarloObservable;
    using StatsBase
    using PyPlot

    import PyPlot: plot

    # ----------------------------------------------------------------------------------
    #   Plots
    # ----------------------------------------------------------------------------------
    """
        plot(obs::Observable{T}[; errors=true, digits=3])

    Plot the observable's time series.
    """
    function plot(obs::Observable{T}; errors=true, digits=3) where T
        ts = timeseries(obs)
        Xmean = mean(obs)
        err = error(obs)
        Xstd = std(ts)

        fig, ax = subplots(1,1)
        ax[:plot](ts, ".-")
        ax[:set_ylabel](name(obs))
        ax[:set_xlabel]("Monte Carlo time \$ t \$")
        # ax[:set_yticks]([])
        ax[:axhline](Xmean, color="black", label="\$ $(round.(Xmean, digits=digits)) \$ (mean)", linewidth=2.0)

        if errors
            ax[:axhline](Xmean+err, color="r", label="\$ \\pm $(round.(err, digits=digits)) \$ (σ error)", linewidth=2.0)
            ax[:axhline](Xmean-err, color="r", linewidth=2.0)

            ax[:axhline](Xmean+2*err, color="r", alpha=.3, label="\$ \\pm $(round.(2*err, digits=digits)) \$ (2σ error)", linewidth=2.0)
            ax[:axhline](Xmean-2*err, color="r", alpha=.3, linewidth=2.0)

            ax[:axhline](Xmean+Xstd, color="g", label="\$ \\pm $(round.(Xstd, digits=digits)) \$ (std)", linewidth=2.0)
            ax[:axhline](Xmean-Xstd, color="g", linewidth=2.0)
        end

        ax[:legend](frameon=true, loc="best")

        tight_layout()
        nothing
    end

    ################################################
    """
        hist(obs::Observable{T}[; errors=true, digits=3])

    Plot a histogram of the observable's time series.
    """
    function hist(obs::Observable{T}; errors=true, digits=3, customtitle="") where T
        ts = timeseries(obs)
        Xmean = mean(obs)
        err = error(obs)
        Xstd = std(ts)

        fig, ax = subplots(1,1)
        ax[:hist](ts, 50, color="gray", alpha=.5, density=1)
        ax[:set_ylabel]("Frequency")
        ax[:set_xlabel](name(obs))
        ax[:set_yticks]([])
        ax[:axvline](Xmean, color="black", label="\$ $(round.(Xmean, digits=digits)) \$ (mean)", linewidth=2.0)

        if errors
            ax[:axvline](Xmean+err, color="r", label="\$ \\pm $(round.(err, digits=digits)) \$ (σ error)", linewidth=2.0)
            ax[:axvline](Xmean-err, color="r", linewidth=2.0)

            ax[:axvline](Xmean+2*err, color="r", alpha=.3, label="\$ \\pm $(round.(2*err, digits=digits)) \$ (2σ error)", linewidth=2.0)
            ax[:axvline](Xmean-2*err, color="r", alpha=.3, linewidth=2.0)

            ax[:axvline](Xmean+Xstd, color="g", label="\$ \\pm $(round.(Xstd, digits=digits)) \$ (std)", linewidth=2.0)
            ax[:axvline](Xmean-Xstd, color="g", linewidth=2.0)
        end

        ax[:legend](frameon=true, loc="best")

        if customtitle != ""
            title(customtitle)
        end

        tight_layout()
        nothing
    end

    ################################################
    """
        corrplot(obs::Observable{T})

    Plot the autocorrelation function of the observable.
    """
    function corrplot(obs::Observable{T}; showtau=false) where T
        ts = timeseries(obs)

        fig, ax = subplots(1,1)
        ax[:plot](autocor(ts), "-", color="k", linewidth=2.0)
        ax[:set_xlabel]("Monte Carlo time \$ t \$")
        ax[:set_ylabel]("Autocorrelation of $(name(obs))")

        if showtau
            t = tau(obs)
            tfinder = tau(obs, MCO.Rplateaufinder(obs)[2])
            tsum = sum(StatsBase.autocor(ts))
            ax[:axvline](x=t, color="red", label="tau ≈ $(round(t,2))")
            ax[:axvline](x=tfinder, color="orange", label="tau (finder) ≈ $(round(tfinder,2))")
            ax[:axvline](x=tsum, color="green", label="tau (sum) ≈ $(round(tsum,2))")
            # @show sum(StatsBase.autocor(ts))
            # @show Rplateaufinder(obs)[2]
            # @show tau(obs)
            ax[:axhline](y=exp(-1), color="gray", alpha=.2, label="1/e")
            ax[:legend]()
        end

        tight_layout()
        nothing
    end

    ################################################
    """
        binningplot(obs::Observable{T}[; min_nbins=32])

    Creates a plot of the binning error coefficient `R` as a function of bin size.

    The coefficient `R` should (up to statistical fluctuations) show a plateau for larger bin sizes,
    indicating that the bin averages have become independent.
    For correlated data one has `R>≈1` and `sqrt(R)` quantifies how much one would have underestimated
    the one-sigma errorbar.

    See [`binning_error`](@ref).
    """
    function binningplot(obs::Observable{T}; min_nbins=50) where T
        ts = timeseries(obs)

        bss, R, means = MCO.R_function(ts, min_nbins=min_nbins)

        figure()
        plot(bss, R, "m.-", label="R")
        plot(bss, means, ".-", label="R means")
        ylabel("error coefficient")
        legend()
        xlabel("bin size")
        tight_layout();
        nothing
    end

    ################################################
    """
        errorplot(obs::Observable{T})

    Creates a plot of the statistical standard error as a function of bin size.

    The standard error should (up to statistical fluctuations) show a plateau for larger bin sizes,
    indicating that the bin averages have become independent.

    See [`binning_error`](@ref), [`error`](@ref).
    """
    function errorplot(obs::Observable; customtitle::String="")
        ts = timeseries(obs)

        bss, R, means = MCO.R_function(ts, min_nbins=50)
        rawerrors = [MCO.binning_error_from_R(ts, r) for r in R]
        meanerrors = [MCO.binning_error_from_R(ts, r) for r in means]

        figure()
        plot(bss, rawerrors, "m.-", label="raw")
        plot(bss, meanerrors, ".-", label="means")
        ylabel("error")
        legend()
        xlabel("bin size")
        if customtitle != ""
            title(customtitle)
        end
        tight_layout();
        nothing
    end

    # ----------------------------------------------------------------------------------
    #   Exports
    # ----------------------------------------------------------------------------------
    export plot, hist, binningplot, errorplot, corrplot

end

    # # --------------------------------------
    # #           Playground
    # # --------------------------------------
    # # using Plots instead of PyPlot

    # function plot_histogram_Plots(obs::Observable{T}; errors=true, digits=3) where T
    #     ts = timeseries(obs)
    #     Xmean = mean(obs)
    #     err = std(obs)
    #     Xstd = std(ts)

    #     histogram(ts, framestyle=:box, grid=false, normed=true, label="", color="lightgrey")
    #     # ax[:hist](ts, 50, color="gray", alpha=.5, normed=1)
    #     ylabel!("Frequency")
    #     xlabel!(name(obs))
    #     yticks!(Float64[])
    #     vline!([Xmean], color="black", linewidth=2.0, label="\$ $(round.(Xmean, digits=digits)) \$ (mean)")

    #     if errors
    #         vline!([Xmean+err, Xmean-err], color="red", linewidth=2.0, label="\$ \\pm $(round.(err, digits=digits)) \$ (σ error)")
    #         vline!([Xmean+2*err, Xmean-2*err], color="red", linewidth=2.0, label="\$ \\pm $(round.(2*err, digits=digits)) \$ (2σ error)", alpha=.3)
    #         vline!([Xmean+Xstd, Xmean-Xstd], color="green", linewidth=2.0, label="\$ \\pm $(round.(Xstd, digits=digits)) \$ (std)", alpha=.3)
    #     end
    #     nothing
    # end

max_stars_repo_name: crstnbr/MCOPlotting

max_issues_repo_name: crstnbr/MCOPlotting

max_forks_repo_name: crstnbr/MCOPlotting


text:
filename() = "day13.txt"

earliest() = parse(Int, readlines(filename())[1])

buses() = parse.(Int, filter(x -> x != "x", split(readlines(filename())[2], ",")))

# Find the earliest time that I can take a bus, at or after time `target`.
# The buses departs at different but consistent times independently. The
# bus schedule is stored in an array of integers `bs`.
#
# Algorithm:
#
# The next available departure time for each bus can be determined by
# a simple formula: (target ÷ b + 1) * b
#
# Just loop through the buses and find the best one. Keep the best
# solution so far in two variables:
# 1. `ride`: which bus? Just an index into the `bs`.
# 2. `best`: best time to board that bus.
#
# Thoughts:
# This could be improved by some kind of `argmin` algorithm.
function part1()
    target = earliest()
    bs = buses()
    best = typemax(Int)
    ride = -1
    for (i, b) in enumerate(bs)
        next = ((target ÷ b) + 1) * b
        if next < best
            best = next
            ride = i
        end
    end
    return bs[ride] * (best - target)
end

# ------------------------------------------
#=
The meaning of the bus schedule input has changed.
- First line is no longer used
- The positions of second line's numbers represents a constraint about how much
  time they have to be aligned against the first bus.

Brute force algorithm (with some math intuition):

Start with the first and second bus. The intuition is that if the first bus departs
at a multiple of `x`, then we must solve a `y` such that the second bus departs at
that time and yet it's exactly `z` minutes after. For example:

    3x + z = 5y  (x_id=3, y_id=5, z=1)

So we can write a function that is given a value of `x`, we can find the *next*
`x` & `y` value such that the equation still holds.
=#
function search_xy(apart, x_id, y_id, x)
    while true
        x += 1
        y, rem = divrem(apart + x_id * x, y_id)
        rem == 0 && return x, y
    end
end

#=
Now, it turns out that the next `x`'s are predictable. In the same example above,
the values are (x,y)'s are as follows:

julia> search_xy(1, 3, 5, 0)
(3, 2)

julia> search_xy(1, 3, 5, 3)
(8, 5)

julia> search_xy(1, 3, 5, 8)
(13, 8)

So it's clear that the value of `x` increments by 5 every time. Hence, once we
have bootstrapped the first candidate value of `x`, the following values are
readily available.

Given these information, the algorithm works as follows:
1. Find the first

=#
read_buses2() = [(pos = i - 1, val = parse(Int, v))
    for (i,v) in enumerate(split(readlines(filename())[2], ",")) if v != "x"]

function part2(atleast=1)
    schedule = read_buses2()
    buses = [s.val for s in schedule]
    apart = diff([s.pos for s in schedule])

    # Find a sensible starting point for first bus such that there's a solution
    # for second bus.
    main_start, _ = search_xy(apart[1], buses[1], buses[2], atleast)

    # The increment for first bus is just the cycle time for second bus according
    # to our math intution above.
    incr = buses[2]

    iter = 0
    while true
        start = main_start
        n = length(buses) - 1
        for i in 1:n
            x = buses[i] * start + apart[i]
            d, r = divrem(x, buses[i+1])
            if r == 0 && i == n  # no remainder all along?
                return main_start * buses[1]
            elseif r != 0        # has remainder, not a solution
                break
            end
            start = d            # starting point for next bus
        end
        iter += 1
        main_start += incr
        # iter > 30 && break     # circut breaking for debugging purpose
    end
end

# https://shainer.github.io/crypto/math/2017/10/22/chinese-remainder-theorem.html

function chinese_remainder_gauss(n, a)
    N = prod(n)
    result = 0
    for i in 1:length(n)
        ai = a[i]
        ni = n[i]
        bi = N ÷ ni
        result += ai * bi * invmod(bi, ni)
    end
    return result % N
end

function part2_chinese_remainder()
    schedule = read_buses2()
    a = [i == 0 ? 0 : v.val - v.pos for (i,v) in enumerate(schedule)]
    n = [v.val for v in schedule]
    return chinese_remainder_gauss(n, a)
end

ext: jl


text:
using libcint_jll

function intPtrFunc!(::Val{:cint1e_ovlp_cart}, buf::Ptr{Cdouble}, shls::Ptr{Cint},
                     atm::Ptr{Cint}, natm::Cint, bas::Ptr{Cint}, nbas::Cint,
                     env::Ptr{Cdouble}, opt::Ptr=C_NULL)
    ccall((:cint1e_ovlp_cart, libcint), Cvoid,
          (Ptr{Cdouble}, Ptr{Cint}, Ptr{Cint}, Cint, Ptr{Cint}, Cint, Ptr{Cdouble}),
          buf, shls, atm, natm, bas, nbas, env)
end

function intPtrFunc!(::Val{:cint1e_nuc_cart}, buf::Ptr{Cdouble}, shls::Ptr{Cint},
                     atm::Ptr{Cint}, natm::Cint, bas::Ptr{Cint}, nbas::Cint,
                     env::Ptr{Cdouble}, opt::Ptr=C_NULL)
    ccall((:cint1e_nuc_cart, libcint), Cvoid,
          (Ptr{Cdouble}, Ptr{Cint}, Ptr{Cint}, Cint, Ptr{Cint}, Cint, Ptr{Cdouble}),
          buf, shls, atm, natm, bas, nbas, env)
end

function intPtrFunc!(::Val{:cint1e_kin_cart}, buf::Ptr{Cdouble}, shls::Ptr{Cint},
                     atm::Ptr{Cint}, natm::Cint, bas::Ptr{Cint}, nbas::Cint,
                     env::Ptr{Cdouble}, opt::Ptr=C_NULL)
    ccall((:cint1e_kin_cart, libcint), Cvoid,
          (Ptr{Cdouble}, Ptr{Cint}, Ptr{Cint}, Cint, Ptr{Cint}, Cint, Ptr{Cdouble}),
          buf, shls, atm, natm, bas, nbas, env)
end

function intPtrFunc!(::Val{:cint2e_cart}, buf::Ptr{Cdouble}, shls::Ptr{Cint},
                     atm::Ptr{Cint}, natm::Cint, bas::Ptr{Cint}, nbas::Cint,
                     env::Ptr{Cdouble}, opt::Ptr=C_NULL)
    ccall((:cint2e_cart, libcint), Cvoid,
          (Ptr{Cdouble}, Ptr{Cint}, Ptr{Cint}, Cint, Ptr{Cint}, Cint, Ptr{Cdouble}, Ptr{Cvoid}),
          buf, shls, atm, natm, bas, nbas, env, opt)
end


text:
struct TrackProperties{T<:AbstractFloat}
    "Width of each rail track"
    track_gauge::T
    "Distance between tracks measured between the inner faces of the rails"
    track_width::T
    "A vector of the orientation, in radians, of the trackbed"
    trackbed_orientation::Vector{T}
    "Width of each sleeper or railroad tie"
    sleeper_width::T
    "Length of each sleeper or railroad tie"
    sleeper_length::T
    "Distance between sleepers or railroad ties"
    sleeper_distance::T
end

function TrackProperties(; track_gauge::T = 1.435, track_width::T = 10e-3, trackbed_orientation::Vector{T} = [0.0], sleeper_width::T = 0.25, sleeper_length::T = 2.6, sleeper_distance::T = 0.70) where T
    return TrackProperties{T}(track_gauge,track_width,trackbed_orientation,sleeper_width,sleeper_length,sleeper_distance)
end

"""
    TracksAhead

Represents everything about the track ahead of the train car. The origin of the coordinate system is on the track ballast and between the two tracks immediately infront of the train car.

The ``(x,y,z)`` of the coordinate system has ``x`` along the rail axis, ``z`` points up towards the sky, and ``y`` then has to point towards the left rail.
"""
struct TracksAhead{T<:AbstractFloat}
    "The positions of the left track when facing forward"
    left_track::Vector{SVector{3,T}}
    "The positions of the right track when facing forward"
    right_track::Vector{SVector{3,T}}
    "the physical properties of the tracks"
    trackproperties::TrackProperties{T}
    function TracksAhead{T}(left_track::Vector{V}, right_track::Vector{V}, trackproperties::TrackProperties{T} = TrackProperties(track_gauge = mean(right_track - left_track)[2])) where {T, V <: AbstractVector{T}}

        if trackproperties.track_gauge < 0
            # throw(DimensionMismatch(""))
            @error "The left_track should be to the left of the right track when ahead of the train. Note the coordinate system used."
        # elseif std(right_track - left_track)[2] > 0.01
        #     @error "The distance between the left_track and right track vary too much."
        else
            new{T}(SVector{3}.(left_track), SVector{3}.(right_track), trackproperties)
        end
    end
end

function TracksAhead(left_track::Vector{V}, right_track::Vector{V}; trackproperties::TrackProperties{T} = TrackProperties(track_gauge = mean(right_track - left_track)[2])) where {T, V <: AbstractVector{T}}

    return TracksAhead{T}(left_track, right_track, trackproperties)

end

"""
    OpticalProperties

Represents everything about the optical properties of the camera.

For the possible fields we have

    - `focal_length`
    - `sensor_width`
    - `sensor_height`
    - `pixelspermeter`
"""
struct OpticalProperties{T<:AbstractFloat}
    "The focal length in meters"
    focal_length::T
    "Sensor width in pixels"
    sensor_width::Int
    "Sensor height in pixels"
    sensor_height::Int
    "Meters in one pixel length"
    pixelspermeter::T
end

function OpticalProperties(focal_length::T; sensor_width::Int = 6500, sensor_height::Int = 5200, pixelspermeter::T = 1 / 5.5e-6) where T
    return OpticalProperties{T}(focal_length,sensor_width,sensor_height,pixelspermeter)
end

"""
    TrainCar

Represents everything about the train car.

For the possible fields we have

    - `xyz` is the centre of mass of the train car in the form ``[x,y,z]`` .
    - `ψθφ` is the orientation of the train car relative to the tracks and is in the form [ψ,θ,φ], where ψ is the roll, θ the pitch, and φ the yaw of the camera.

"""
struct TrainCar{T<:AbstractFloat}
    ψθφ::SVector{3,T}
    xyz::SVector{3,T}
end

TrainCar(ψθφ::AbstractVector{T}, xyz::AbstractVector{T}) where T = TrainCar{T}( SVector{3,T}(ψθφ),SVector{3,T}(xyz))

TrainCar(ψθφ::AbstractVector{T} = zeros(T,3); xyz::AbstractVector{T} = zeros(T,3)) where T = TrainCar(ψθφ,xyz)

"""
    VideoCamera

Represents everything about the camera position and the camera's [`OpticalProperties`](@ref).

For the possible fields we have

    - `xyz` is a vector of the form ``[x,y,z]`` which gives the position of the camera if the track had no defects and the train was still.
    - `ψθφ` is a vector of the form [ψ,θ,φ] where ψ is the roll, θ the pitch, and φ the yaw of the camera.

"""
struct VideoCamera{T<:AbstractFloat}
    xyz::SVector{3,T}
    ψθφ::SVector{3,T}
    opticalproperties::OpticalProperties{T}
end

VideoCamera(xyz::AbstractVector{T}, ψθφ::AbstractVector{T}, opticalproperties::OpticalProperties{T}) where T = VideoCamera{T}(SVector{3,T}(xyz), SVector{3,T}(ψθφ), opticalproperties)

function VideoCamera(cameraposition::AbstractVector{T};
    focalpoint::AbstractVector{T} = [T(30),T(0),T(0)],
    ψθφ::AbstractVector{T} = [T(0); focalθφ(focalpoint,cameraposition)],
    focal_length::T = 6e-3,
    kws...) where T

    return VideoCamera(cameraposition, ψθφ, OpticalProperties(focal_length; kws...))
end

"""
    VideoCamera(camera::VideoCamera{T}, distortion::Dict)

Returns a `VideoCamera` with the same properties as camera, except with the added distortion.
"""
function VideoCamera(camera::VideoCamera{T}, distortions::Dict) where T
    δxyz = [get(distortions, s, zero(T)) for s in [:X,:Y,:Z]]
    δψθφ = [get(distortions, s, zero(T)) for s in [:ψ,:θ,:φ]]
    δf = get(distortions, :f, zero(T))

    dict = Dict(
        n => getfield(camera.opticalproperties,n)
    for n in fieldnames(OpticalProperties))

    f = dict[:focal_length] + δf
    delete!(dict,:focal_length)

    opt = OpticalProperties(f; dict...)

    return VideoCamera(camera.xyz - δxyz, camera.ψθφ + δψθφ, opt)
end

focal_length(camera::VideoCamera{T}) where T = camera.opticalproperties.focal_length

"""
    focalθφ(focalpoint::Vector, cameraposition::Vector)

Calculate the pitch and yaw of the camera (θ, φ) that would make the focalpoint in 3D be centered in the camera's focus / sensor.
"""
function focalθφ(focalpoint::Vector{T}, cameraposition::Vector{T}) where T
# focalpoint = [xyz[1] + xyz[3]*cot(θ)*cos(φ), xyz[2] + xyz[3]*cot(θ)*sin(φ), 0.0]

    cx, cy, cz = cameraposition

    x = focalpoint[1] - cx
    y = focalpoint[2] - cy

    φ = atan(y,x)
    θ = atan(cz * cos(φ) / x)

    if θ < 0 || abs(φ) > pi
        @warn "unexpected angles θ = $θ or φ = $φ for a camera mounted on a train"
    end

    return [θ, φ]
end

Rψ(ψ::T) where T = SMatrix{3,3}(
    [T(1) T(0)    T(0);
     T(0) cos(ψ)  sin(ψ);
     T(0) -sin(ψ) cos(ψ)]
)

Rθ(θ::T) where T = SMatrix{3,3}(
    [cos(θ) T(0) -sin(θ);
     T(0)   T(1) T(0);
     sin(θ) T(0) cos(θ)]
)

Rφ(φ::T) where T = SMatrix{3,3}(
    [cos(φ)  sin(φ) T(0);
     -sin(φ) cos(φ) T(0);
     T(0)    T(0)   T(1)]
)

"""
    track_image_u(camera::VideoCamera, Y::T, vs::AbstractVector)

returns the us, corresponding to vs, in the camera's image of a track which is aligned with the camera's X axis. This track is assumed to be completely straight, and the camera properties are all known.
"""
function track_image_u(camera::VideoCamera{T}, Y::T, vs::AbstractVector{T}) where T
    Z = - camera.xyz[3]

    f = camera.opticalproperties.focal_length * camera.opticalproperties.pixelspermeter

    ψ, θ, φ = camera.ψθφ

    a1 = Y * cos(θ) * cos(ψ) - Z * sin(θ) * cos(ψ) * sin(φ) + Z * cos(φ) * sin(ψ)
    a2 = Z * cos(φ) * cos(ψ) + (Z * sin(θ) * sin(φ) - Y * cos(θ)) * sin(ψ)
    a = a1 / a2

    b1 = - f * (Y * sin(θ) + Z * cos(θ) * sin(φ))
    b = b1 / a2

    us = b .+ a .* vs
    return us
end

function left_track_image_u(camera::VideoCamera{T}, trackproperties::TrackProperties{T}, vs::AbstractVector{T}) where T
    # traincar::TrainCar{T} = TrainCar()) where T

    return track_image_u(camera, -trackproperties.track_gauge / T(2) - camera.xyz[2], vs)
end

function right_track_image_u(camera::VideoCamera{T}, trackproperties::TrackProperties{T}, vs::AbstractVector{T}) where T

    return track_image_u(camera, trackproperties.track_gauge / T(2) - camera.xyz[2], vs)
end

function camera_image(camera::VideoCamera{T}, spatial_points::Vector{V}, traincar::TrainCar{T} = TrainCar(zeros(T,3))) where {T, V<:AbstractVector{T}}

    ψ, θ, φ = camera.ψθφ
    RC = Rψ(ψ) * Rθ(θ) * Rφ(φ)

    ψT, θT, φT = traincar.ψθφ
    RT = Rψ(ψT) * Rθ(θT) * Rφ(φT)

    R = RC * RT

    ξηζs = [ R * (xyz - camera.xyz) for xyz in spatial_points]

    f = camera.opticalproperties.focal_length * camera.opticalproperties.pixelspermeter

    uvs = map(ξηζs) do ξηζ
        scaletopixel = f / ξηζ[1]
        SVector{2}(ξηζ[2:3] .* scaletopixel)
    end

   inside_sensor(uv::SVector{2,T}) = abs(uv[1]) <= camera.opticalproperties.sensor_width / T(2) && abs(uv[2]) <= camera.opticalproperties.sensor_height / T(2)

   return filter(inside_sensor, uvs)

end


text:
# SEP-RANK
## Greeting and welcome to my code for computing lowerboounds on the separable rank.

Pkg.activate(".")
    srcDir  = dirname(@__FILE__)*"\\";
    include(srcDir*"Examples\\Examples.jl")
    include(srcDir*"Model\\R_sep_Model.jl")
    include(srcDir*"Model\\C_sep_Model.jl")
    include(srcDir*"sep_Compute.jl")

    using .Examples
    using .R_sep_Model
    using .C_sep_Model
    using .sep_Compute
### Quick start
#### Manually load data and run the model:
ρ  =  [1. 0. 0. 0.
       0. 0. 0. 0.
       0. 0. 0. 0.
       0. 0. 0. 1.]
d₁,d₂  = (2,2) ; cl = "S1sG" ; t  = (2,2)
sep_Cmod = C_sep_Model.Modelξₜˢᵉᵖ(ρ,(d₁,d₂),t;con_list = cl)
sep_Cmod_sol = sep_Compute.Computeξₜˢᵉᵖ(sep_Cmod)

Lx = sep_Cmod[:Lx]

#### Run the **ℝeal** model with the same data as above:
sep_Rmod = R_sep_Model.Modelξₜˢᵉᵖ(ρ,(d₁,d₂),t;con_list = cl)
sep_Rmod_sol = sep_Compute.Computeξₜˢᵉᵖ(sep_Rmod)

#### Load example data and run complex model
ex = "E7_i"
d = ex.get_example_meta(ex).Size[1] ; ρ = ex.get_example(ex) ; cl = "S1sG" ; t  = (2,2)
sep_Cmod = C_sep_Model.Modelξₜˢᵉᵖ(ρ,(d₁,d₂),t;con_list = cl)
sep_Cmod_sol = sep_Compute.Computeξₜˢᵉᵖ(sep_Cmod)

### Details ***(optional)***
#### Preloaded Examples
df = Examples.get_example_meta()
ρ_dict = Examples.get_example()

df["E7_i",:]
ρ_dict["E7_i"]
#Look up example in appendix.

#### Running a batch of examples
include(srcDir*"batch.jl")
using .batch

batch.batch_model(t,ρ_dic,df;sdir = save_directory,con_list= ["S1sG","S2sG","S3sG"])

#### Loading the constraint matrices
include(srcDir*"Constraints\\R_constraints.jl")
include(srcDir*"Constraints\\C_constraints.jl")
using .C_constraints
using .R_constraints
using JuMP

model = JuMP.Model() ; @variable(model, Lx[ccon.make_mon_expo_keys(d,t[1])])

#### The **ℂomplex** moment matrix, 4ᵗʰ-order constraints, localizing maps and matrix localizing map.
C_PSD_con = C_constraints.make_PSD_con(d,t,Lx)
C_ord4_con = C_constraints.make_ord4_con(d,Lx)
C_loc_cons_S1 = C_constraints.make_loc_cons_S1(ρ,d,t,Lx)
C_loc_cons_S2 = C_constraints.make_loc_cons_S2(ρ,d,t,Lx)
C_loc_cons_S3 = C_constraints.make_loc_cons_S3(ρ,d,t,Lx)
C_Gᴿ_con = C_constraints.make_Gᴿ_con(ρ,d,t,Lx)

#### The **ℝeal** moment matrix, 4ᵗʰ-order constraints, localizing maps and matrix localizing map.
R_PSD_con = R_constraints.make_PSD_con(d,t,Lx)
R_ord4_con = R_constraints.make_ord4_con(d,Lx)
R_loc_cons_S1 = R_constraints.make_loc_cons_S1(ρ,d,t,Lx)
R_loc_cons_S2 = R_constraints.make_loc_cons_S2(ρ,d,t,Lx)
R_loc_cons_S3 = R_constraints.make_loc_cons_S3(ρ,d,t,Lx)
R_G_con = R_constraints.make_G_con(ρ,d,t,Lx)

max_stars_repo_name: JAndriesJ/sep-rank

max_issues_repo_name: JAndriesJ/sep-rank

max_forks_repo_name: JAndriesJ/sep-rank


text:
using Test
using ClimateMachine
using ClimateMachine.Mesh.Grids
using ClimateMachine.Mesh.Topologies
using ClimateMachine.MPIStateArrays
using Impero, Printf, MPI, LinearAlgebra, Statistics, GaussQuadrature
using ImperoHooks
import ImperoHooks: launch_volume_gradient!, launch_interface_gradient!
import ImperoHooks: launch_volume_divergence!, launch_interface_divergence!

include(pwd() * "/test/gradient_test_utils.jl")

ClimateMachine.init()
const ArrayType = ClimateMachine.array_type()
const mpicomm = MPI.COMM_WORLD
const FT = Float64
Ω = Circle(-1,1) × Circle(-1,1) × Circle(-1,1)
dims = ndims(Ω)

ClimateMachine.gpu_allowscalar(true)

polynomialorder = (3,3,3)
elements = (2,1,1)
nrealelem = prod(elements)
grid = DiscontinuousSpectralElementGrid(Ω, elements = elements, polynomialorder = polynomialorder, array = ArrayType)
println(nrealelem - length(grid.interiorelems))
x, y, z = coordinates(grid)
device = array_device(x)
##

# Define Arrays
# Divergence Arrays
ijksize = prod(polynomialorder .+ 1)
Φ = MPIStateArray{FT}(mpicomm, ArrayType, ijksize, nrealelem , 3)
v_flux_divergence = MPIStateArray{FT}(mpicomm, ArrayType, ijksize, nrealelem , 1)
s_flux_divergence = MPIStateArray{FT}(mpicomm, ArrayType, ijksize, nrealelem , 1)
# Gradient Arrays
Q    = MPIStateArray{FT}(mpicomm, ArrayType, ijksize, nrealelem, 1)
v_∇Q = MPIStateArray{FT}(mpicomm, ArrayType, ijksize, nrealelem, 3)
s_∇Q = MPIStateArray{FT}(mpicomm, ArrayType, ijksize, nrealelem, 3)
##
# Test 1: Continuous function
a = 1
b = 1
c = 1
ClimateMachine.gpu_allowscalar(true)
@. Φ.realdata[:,:,1] = a * sin(π*x)
@. Φ.realdata[:,:,2] = b * sin(π*y)
@. Φ.realdata[:,:,3] = c * cos(π*z)
@. Q = a * sin(π*x) + b * sin(π*y) + c * cos(π*z)

# Divergence
event = launch_volume_divergence!(v_flux_divergence, Φ, grid)
wait(event)
event = launch_interface_divergence!(s_flux_divergence, Φ, grid)
wait(event)

# Gradient
event = launch_volume_gradient!(v_∇Q, Q, grid)
wait(event)
event = launch_interface_gradient!(s_∇Q, Q, grid)
wait(event)

tol = eps(1e5)
L∞(x) = maximum(abs.(x))
@testset "Continuous: Divergence = Trace(Gradient)" begin
    total_divergence = -v_flux_divergence + s_flux_divergence
    total_trace_gradient = sum(v_∇Q - s_∇Q, dims=3)
    @test L∞(s_∇Q  ) < tol
    @test L∞(total_divergence - ArrayType(total_trace_gradient)) < tol
end

##
# Test 2: Discontinuous function
a = 1
b = 0
c = 0
ClimateMachine.gpu_allowscalar(true)
# Heaviside Step Function Check
@. Φ[:,1,1] = 0.0
@. Φ[:,2,1] = 1.0
@. Φ[:,:,2] = 0.0
@. Φ[:,:,3] = 0.0
@. Q = Φ[:,:,1]

# Divergence
event = launch_volume_divergence!(v_flux_divergence, Φ, grid)
wait(event)
event = launch_interface_divergence!(s_flux_divergence, Φ, grid)
wait(event)

# Gradient
event = launch_volume_gradient!(v_∇Q, Q, grid)
wait(event)
event = launch_interface_gradient!(s_∇Q, Q, grid)
wait(event)

tol = eps(1e5)
L∞(x) = maximum(abs.(x))
@testset "Discontinuous: Divergence = Trace(Gradient)" begin
    total_divergence = -v_flux_divergence + s_flux_divergence
    total_gradient =  v_∇Q + s_∇Q
    total_trace_gradient = sum(total_gradient, dims=3)
    @test L∞(v_∇Q) < tol # since piecewise constant
    @test L∞(total_divergence - ArrayType(total_trace_gradient)) < tol
end

max_stars_repo_name: CliMA/ImperoHooks

max_issues_repo_name: CliMA/ImperoHooks

max_forks_repo_name: CliMA/ImperoHooks


text:
using MPI

function runmpi(file; ntasks = 1)
    MPI.Initialized() &&
    !MPI.Finalized() &&
    error(
        "runmpi does not work if MPI has been " *
        "Initialized but not Finalized",
    )

    if MPI.MPI_LIBRARY == MPI.OpenMPI
        oversubscribe = `--oversubscribe`
    else
        oversubscribe = ``
    end

    @info "Running MPI test..." file ntasks
    # Running this way prevents:
    #   Balance Law Solver | No tests
    # since external tests are not returned as passed/fail
    @time @test mpiexec() do cmd
        run(`$cmd $oversubscribe -n $ntasks $(Base.julia_cmd()) --startup-file=no --project=$(Base.active_project()) $file`)
        true
    end
end


text:
abstract type StageCost end

struct StageCosts <: Objective
    cost::Vector{StageCost}
    T::Int
end

g(obj::StageCosts, x, u, t) = 0.0

function objective(obj::StageCosts, x, u)
    T = obj.T
    J = 0.0
    for t = 1:T-1
        J += g(obj, x[t], u[t], t)
    end
    J += g(obj, x[T], nothing, T)
    return J
end

"""
    nonlinear stage cost
"""
struct NonlinearCost <: StageCost
end

NonlinearCosts = Vector{NonlinearCost}

"""
    quadratic stage cost
"""
struct QuadraticCost <: StageCost
    Q
    q
    R
    r
end

QuadtraticCosts = Vector{QuadraticCost}

max_issues_repo_name: thowell/DirectMotionPlanning


text:
"""
    get_pedfile_trios(pedfile.txt)

Read a standard format pedigree file and extract
complete trios and return a dictionary of pedigrees
with a list of trios (child,father,mother) in
each pedigree.
"""
function get_pedfile_trios(pedfile)

ptrios = Dict{String, Array}()

open(pedfile) do file

    for i in eachline(file)

        if occursin(r"^[#|P|p]", i)
            continue
        end

        line = split(i, r"\t|\,|\s+")

        if (line[3] != "0") && (line[4] != "0")

            t = join([ line[2], line[3], line[4] ], ",")

            if !haskey(ptrios, line[1])
                ptrios[ line[1] ] = [t]
            else
                push!(ptrios[ line[1] ], t)
            end

        end
    end
end
    peds = collect(keys(ptrios))
    return(ptrios, peds)
end

max_stars_repo_name: ScottWatkins/STRdiff

max_issues_repo_name: ScottWatkins/STRdiff

max_forks_repo_name: ScottWatkins/STRdiff


text: function f() :: Int64 return 1 end

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
@testset "Local Linear Trend Model" begin
    finland_fatalities = CSV.File(StateSpaceModels.VEHICLE_FATALITIES) |> DataFrame
    log_finland_fatalities = log.(finland_fatalities.ff)

    @test has_fit_methods(LocalLinearTrend)
    # Based on
    # https://www.statsmodels.org/stable/examples/notebooks/generated/statespace_local_linear_trend.html?highlight=linear%20trend
    # Fitting the model with mod = sm.tsa.UnobservedComponents(df['lff'], 'local linear trend') gives 26.740
    model = LocalLinearTrend(log_finland_fatalities)
    fit!(model)
    @test loglike(model) ≈ 26.740 atol = 1e-5 rtol = 1e-5

    # forecasting
    forec = forecast(model, 10)
    @test monotone_forecast_variance(forec)

    predicted_mean = [
        5.933257,
        5.897660,
        5.862062,
        5.826465,
        5.790867,
        5.755270,
        5.719672,
        5.684074,
        5.648477,
        5.612879
    ]
    @test maximum(abs.(predicted_mean .- vcat(forec.expected_value...))) < 1.5e-3
end


text:
@testset "Wilcoxon rank sum exact test" begin
    x1 = [0.80, 0.83, 1.89, 1.04, 1.45, 1.38, 1.91, 1.64, 0.73, 1.46]
    x2 = [1.15, 0.88, 0.90, 0.74, 1.21]
    res_wtest = wilcox_test(x1, x2)
    test_show(res_wtest)
    @test res_wtest.statistic ≈ 35
    @test res_wtest.pval ≈ 0.254412254412254
end

@testset "Wilcoxon rank sum test with continuity correction" begin
    df_sleep = dataset("datasets", "sleep")
    res_wtest = wilcox_test(@formula(Extra ~ Group), df_sleep)
    test_show(res_wtest)
    @test res_wtest.statistic ≈ 25.5
    @test res_wtest.pval ≈ 0.0693275754336266
end

@testset "Wilcoxon signed rank sum exact test" begin
    x1 = [1.83,  0.50,  1.62,  2.48, 1.68, 1.88, 1.55, 3.06, 1.30]
    x2 = [0.878, 0.647, 0.598, 2.05, 1.06, 1.29, 1.06, 3.14, 1.29]
    res_wtest = wilcox_test(x1, x2, paired=true)
    test_show(res_wtest)
    @test res_wtest.statistic ≈ 40
    @test res_wtest.pval ≈ 0.0390625
end

@testset "Wilcoxon signed rank sum test with continuity correction" begin
    x1 = [-2.0, -3.2, 1.2, 5.8, 19.6, 19.9, 28.9, 31.1, 27.9, 21.9, 8.8, 7.5]
    x2 = [-1.8, -2.9, 0.9, 1.5, 11.5, 18.8, 23.9, 24.9, 23.8, 19.4, 8.8, 4.8]
    res_wtest = wilcox_test(x1, x2, paired=true)
    test_show(res_wtest)
    @test res_wtest.statistic ≈ 62
    @test res_wtest.pval ≈ 0.011278190116074
end


text:
<<<<<<< HEAD
=======
include("globVar.jl")

using .errVar
using .glfwVar

>>>>>>> 5f9517fcfc507bc9c83bc0ca5755ab4a4d01f653
function errExit()
  global fstream

  flush(fstream)
  close(fstream)
end

function mainExit()
  errExit()

  GLFW.DestroyWindow(win)
end


text:
@testset "Test function update_coeffs I" begin

    Nx = 2
    Ne = 8
    m = 5
    X = zeros(Nx, Ne)

    X .=  [0.267333   1.43021;
              0.364979   0.607224;
             -1.23693    0.249277;
             -2.0526     0.915629;
             -0.182465   0.415874;
              0.412907   1.01672;
              1.41332   -0.918205;
              0.766647  -1.00445]';
    # Initialize map with zero coefficients
    C_old = HermiteMapComponent(m, Nx; α = 1e-6, b= "CstProHermiteBasis");

    setcoeff!(C_old, [1.5])
    reduced_margin = getreducedmargin(getidx(C_old))
    @show reduced_margin
    idx_new = vcat(getidx(C_old), reduced_margin)

    C_new = HermiteMapComponent(m, Nx, idx_new, vcat(getcoeff(C_old), zeros(2)));

    coeff_new, coeff_idx_added, idx_added = update_coeffs(C_old, C_new)

    @test norm(coeff_new - [1.5; 0.0; 0.0])<1e-8
    @test coeff_idx_added == [2; 3]
    @test idx_added == [0 1; 1 0]
end

@testset "Test function update_coeffs II" begin

    Nx = 2
    Ne = 8
    m = 5
    X = zeros(Nx, Ne)

    X    .=  [0.267333   1.43021;
              0.364979   0.607224;
             -1.23693    0.249277;
             -2.0526     0.915629;
             -0.182465   0.415874;
              0.412907   1.01672;
              1.41332   -0.918205;
              0.766647  -1.00445]';
    # Initialize map with zero coefficients
    C_old = HermiteMapComponent(m, Nx, [0 0; 1 0; 2 0; 0 1; 0 2; 1 1], randn(6); α = 1e-6);

    reduced_margin = getreducedmargin(getidx(C_old))
    idx_new = vcat(getidx(C_old), reduced_margin)
    C_new = HermiteMapComponent(m, Nx, idx_new, vcat(getcoeff(C_old), zeros(4)));
    coeff_new, coeff_idx_added, idx_added = update_coeffs(C_old, C_new)

    @test norm(coeff_new - vcat(getcoeff(C_old), zeros(4)))<1e-8
    @test coeff_idx_added == [7; 8; 9; 10]
    @test idx_added == [0 3; 1 2; 2 1; 3 0]
end

@testset "Test function update_coeffs III" begin
    Nx = 3
    m = 5
    # Initialize map with zero coefficients
    C_old = HermiteMapComponent(m, Nx, [0 0 0; 0 0 1; 0 0 2; 0 1 0; 0 2 0], randn(5); α = 1e-6);

    reduced_margin = getreducedmargin(getidx(C_old))
    idx_new = vcat(getidx(C_old), reduced_margin)
    C_new = HermiteMapComponent(m, Nx, idx_new, vcat(getcoeff(C_old), zeros(4)));
    coeff_new, coeff_idx_added, idx_added = update_coeffs(C_old, C_new)

    @test norm(coeff_new - vcat(getcoeff(C_old), zeros(4)))<1e-8
    @test coeff_idx_added == [6; 7; 8; 9]
    @test idx_added == reduced_margin
end

@testset "Test update_component I" begin
    Nx = 2
    Ne = 8
    m = 5
    X = zeros(Nx, Ne)

    X .=     [0.267333   1.43021;
              0.364979   0.607224;
             -1.23693    0.249277;
             -2.0526     0.915629;
             -0.182465   0.415874;
              0.412907   1.01672;
              1.41332   -0.918205;
              0.766647  -1.00445]';
    X = X
    # Initialize map with zero coefficients
    C_old = HermiteMapComponent(m, Nx; α = 1e-6);

    S = Storage(C_old.I.f, X)

    dJ_old = zeros(1)
    J_old = 0.0
    J_old = negative_log_likelihood!(J_old, dJ_old, getcoeff(C_old), S, C_old, X)

    @test abs(J_old - 1.317277788545110)<1e-10

    @test norm(dJ_old - [0.339034875000000])<1e-10

    # setcoeff!(C_old, [1.5])
    reduced_margin0 = getreducedmargin(getidx(C_old))
    idx_old0 = getidx(C_old)
    idx_new0 = vcat(idx_old0, reduced_margin0)

    # Define updated map
    f_new = ExpandedFunction(C_old.I.f.B, idx_new0, vcat(getcoeff(C_old), zeros(size(reduced_margin0,1))))
    C_new = HermiteMapComponent(f_new; α = 1e-6)
    idx_new, reduced_margin = update_component!(C_old, X, reduced_margin0, S)

    dJ_new = zeros(3)
    J_new = 0.0
    S = update_storage(S, X, reduced_margin0)
    J_new = negative_log_likelihood!(J_new, dJ_new, getcoeff(C_new), S, C_new, X)

    @test abs(J_new -    1.317277788545110)<1e-8

    @test norm(dJ_new - [0.339034875000000;   0.228966410748600;   0.192950409869679])<1e-8

    @test findmax(abs.(dJ_new[2:3]))[2] == 2-1

    @test idx_new0[3,:] == reduced_margin0[2,:]

    @test reduced_margin == updatereducedmargin(getidx(C_old), reduced_margin0, 1)[2]
    @test idx_new == updatereducedmargin(getidx(C_old), reduced_margin0, 1)[1]
end

@testset "Test update_component II" begin
    Nx = 2
    Ne = 8
    m = 5
    X = zeros(Nx, Ne)

    X .=     [0.267333   1.43021;
              0.364979   0.607224;
             -1.23693    0.249277;
             -2.0526     0.915629;
             -0.182465   0.415874;
              0.412907   1.01672;
              1.41332   -0.918205;
              0.766647  -1.00445]';
    X = X

    idx = [0 0; 0 1; 1 0; 0 2; 2 0; 1 1]

    coeff = [ -0.9905841755746164;
          0.6771992097558741;
         -2.243695806805015;
         -0.34591297359447354;
         -1.420159186008486;
         -0.5361337327704369]

    C_old = HermiteMapComponent(m, Nx, idx, coeff; α = 1e-6);

    S = Storage(C_old.I.f, X)

    dJ_old = zeros(6)
    J_old = 0.0
    J_old = negative_log_likelihood!(J_old, dJ_old, getcoeff(C_old), S, C_old, X)

    @test abs(J_old - 3.066735373495125)<1e-5

    dJt_old  = [-1.665181693060798;
                -0.919442416275569;
                -0.900922248322526;
                -0.119522662066802;
                -0.472794468187974;
                -0.496862354324658]

    @test norm(dJ_old - dJt_old)<1e-5

    reduced_margin0 = getreducedmargin(getidx(C_old))
    idx_old0 = getidx(C_old)
    idx_new0 = vcat(idx_old0, reduced_margin0)

    # Define updated map
    f_new = ExpandedFunction(C_old.I.f.B, idx_new0, vcat(getcoeff(C_old), zeros(size(reduced_margin0,1))))
    C_new = HermiteMapComponent(f_new; α = 1e-6)
    idx_new, reduced_margin = update_component!(C_old, X, reduced_margin0, S)

    dJ_new = zeros(10)
    J_new = 0.0
    S = update_storage(S, X, reduced_margin0)
    J_new = negative_log_likelihood!(J_new, dJ_new, getcoeff(C_new), S, C_new, X)

     Jt_new  = 3.066735373495125
    dJt_new  = [-1.665181693060798;
                -0.919442416275569;
                -0.900922248322526;
                -0.119522662066802;
                -0.472794468187974;
                -0.496862354324658;
                 0.338831623096620;
                -0.083502724667314;
                -0.281130533740555;
                 0.271587190562667]

    @test abs(J_new -    Jt_new)<1e-5
    @test norm(dJ_new - dJt_new)<1e-5

    @test findmax(abs.(dJ_new[7:10]))[2] == 7-6

    @test idx_new0[8,:] == reduced_margin0[2,:]

    @test reduced_margin == updatereducedmargin(getidx(C_old), reduced_margin0, 1)[2]
    @test idx_new == updatereducedmargin(getidx(C_old), reduced_margin0, 1)[1]
end

@testset "Test greedy optimization on training set only without QR " begin

Nx = 2
m = 15

X = Matrix([0.267333   1.43021;
          0.364979   0.607224;
         -1.23693    0.249277;
         -2.0526     0.915629;
         -0.182465   0.415874;
          0.412907   1.01672;
          1.41332   -0.918205;
          0.766647  -1.00445]')
L = LinearTransform(X)
# transform!(L, X)
C_new, error = greedyfit(m, Nx, X, 6; withconstant = false, withqr = false, verbose = false)

@test norm(error -  [1.317277788545110;
                     1.240470006245235;
                     1.239755181933741;
                     1.191998917636778;
                     0.989821300578785;
                     0.956587763340890;
                     0.848999846549310])<1e-4

@test norm(getcoeff(C_new) - [ 2.758739330966050;
                               4.992420909046952;
                               0.691952693078807;
                               5.179459837781270;
                              -3.837569097976292;
                               1.660285776681347])<1e-4

@test norm(getidx(C_new) -   [0  1;
                              0  2;
                              0  3;
                              0  4;
                              1  0;
                              2  0])<1e-10
end

@testset "Test greedy optimization on training set only with QR " begin

Nx = 2
m = 15

X = Matrix([0.267333   1.43021;
          0.364979   0.607224;
         -1.23693    0.249277;
         -2.0526     0.915629;
         -0.182465   0.415874;
          0.412907   1.01672;
          1.41332   -0.918205;
          0.766647  -1.00445]')
L = LinearTransform(X)
# transform!(L, X)
C_new, error = greedyfit(m, Nx, X, 6; withconstant = false, withqr = true, verbose = false)

@test norm(error -  [1.317277788545110;
                     1.240470006245235;
                     1.239755181933741;
                     1.191998917636778;
                     0.989821300578785;
                     0.956587763340890;
                     0.848999846549310])<1e-4

@test norm(getcoeff(C_new) - [ 2.758739330966050;
                               4.992420909046952;
                               0.691952693078807;
                               5.179459837781270;
                              -3.837569097976292;
                               1.660285776681347])<1e-4

@test norm(getidx(C_new) -   [0  1;
                              0  2;
                              0  3;
                              0  4;
                              1  0;
                              2  0])<1e-10
end

@testset "Test greedy optimization on train dataset without/with QR" begin

Nx = 2
m = 10
Ne = 10^3
X = randn(2, Ne).^2 + 0.1*randn(2, Ne)
L = LinearTransform(X)
transform!(L, X)

C_noqr, error_noqr = greedyfit(m, Nx, X, 10; withconstant = true, withqr = false, verbose = false);
C_qr, error_qr = greedyfit(m, Nx, X, 10; withconstant = true, withqr = true, verbose = false);

train_error_noqr = error_noqr
train_error_qr   = error_qr

C_test = deepcopy(C_noqr)
setcoeff!(C_test, zero(getcoeff(C_noqr)));

S_test = Storage(C_test.I.f, X)
coeff_test = getcoeff(C_test)

res = Optim.optimize(Optim.only_fg!(negative_log_likelihood(S_test, C_test, X)), coeff_test, Optim.BFGS())

@test norm(getcoeff(C_noqr)-Optim.minimizer(res))<5e-5
@test norm(getcoeff(C_qr)-Optim.minimizer(res))<5e-5

@test norm(train_error_noqr[end] - res.minimum)<5e-5
@test norm(train_error_qr[end]   - res.minimum)<5e-5

end

@testset "Test greedy optimization on train/valid dataset without/with QR" begin

Nx = 2
m = 10
Ne = 10^3
X = randn(2, Ne).^2 + 0.1*randn(2, Ne)
L = LinearTransform(X)
transform!(L, X)

X_train = X[:,1:800]
X_valid = X[:,801:1000]

C_noqr, error_noqr = greedyfit(m, Nx, X_train, X_valid, 10; withconstant = true, withqr = false, verbose = false);
C_qr, error_qr = greedyfit(m, Nx, X_train, X_valid, 10; withconstant = true, withqr = true, verbose = false);

train_error_noqr, valid_error_noqr = error_noqr
train_error_qr, valid_error_qr = error_qr

C_test = deepcopy(C_noqr)
setcoeff!(C_test, zero(getcoeff(C_noqr)));

S_test = Storage(C_test.I.f, X_train)
coeff_test = getcoeff(C_test)

res = Optim.optimize(Optim.only_fg!(negative_log_likelihood(S_test, C_test, X_train)), coeff_test, Optim.BFGS())

@test norm(getcoeff(C_noqr) - Optim.minimizer(res))<5e-5
@test norm(getcoeff(C_qr)  - Optim.minimizer(res))<5e-5

@test norm(train_error_noqr[end] - res.minimum)<5e-5
@test norm(train_error_qr[end] - res.minimum)<5e-5
@test norm(valid_error_noqr[end] - valid_error_qr[end])<5e-5

end

@testset "Test greedyfit routine with/without QR, with/without preconditioner" begin

    Nx = 2
    Ne = 100

    X = Matrix([   0.096932076386797  -0.062113573952594
       1.504907056818036  -1.270761688825393
       0.122588515870727   0.138604054973047
       0.095040305032907  -0.097989669097228
       0.056360334325948  -0.098106117790437
       0.034757567319988  -0.152864122871058
       0.260831995265119  -0.564465060524812
       0.191690179163063   0.097123146831467
       0.131835643270215  -0.303532750891873
       0.121404994711150  -0.414223062306829
      -0.046027217149733   0.267918560101165
       0.222891154723706   0.054734106232808
       0.264238681579693  -0.561044970267127
       0.021345749974379   0.023020769850203
       0.395423632217738  -0.619284581169702
       0.231607199656522   0.319089487692097
       0.389505867820041   0.572610562725449
      -0.146845693861564  -0.007959126279965
       0.138489687840940   0.330481872227252
       1.617286792209003  -1.272765331253319
       0.013078993247482   0.006243503157725
       0.304266481947197  -0.508127858160603
      -0.134224572142272  -0.057784550290700
       0.058245459428341  -0.388163133932979
       0.143773871017238  -0.569911233902741
      -0.020852128497605   0.159888469685458
       0.042924979621615  -0.285726823344180
       0.566514904142701  -0.815010017214964
       0.022782788493659  -0.468191560991914
      -0.072005082647238  -0.139317736738881
      -0.088861617714439  -0.042323531263036
       0.146019714866274  -0.441835803838725
       0.062827433058116   0.255887684932700
       0.024887528661752  -0.036316627235643
       0.079288422911736   0.009625876863572
       0.124929172288805  -0.196711126223702
       0.886085483369468   1.012837297288688
      -0.062868708331948   0.065147939699519
       0.117047537973873  -0.301180703268764
       0.147117619773736   0.466610279392604
       0.690523534483231   0.722717906410571
      -0.110919157441807   0.000008657527041
       2.055034349742512   1.423365847218695
       0.091130503033388  -0.192933472923564
       0.075510323340039   0.075087088602955
      -0.181121769738736  -0.158594842114751
      -0.080866378396288  -0.081393537160189
       0.184561993178819   0.409693133491638
      -0.032018126656783   0.040851821749292
       0.702808256935611   0.978185567019995
      -0.029041940950512   0.097268795862442
       1.406479923500753  -1.166340130567743
      -0.109758856235453  -0.297180107753266
       0.236494055349260   0.548223372779986
      -0.108550945852467  -0.097740533128956
       0.908726493797006  -0.957752761353643
       0.178092997560103  -0.359297653971968
       1.103683368372646  -1.064482193617671
       1.907065203740822   1.349542162065178
       0.523722913693736   0.725695151794695
       0.261131020580618   0.576363120484564
       0.118044539009197   0.196304662470752
       0.289261334786348   0.399639383890177
       0.902906400981006  -0.957301599045371
      -0.054657884786803  -0.292760625531884
      -0.021735291600425   0.029650166664389
       0.065200888050752  -0.295894159582647
       1.486186253485923  -1.217814186300608
       0.889545420155124   0.939789761164950
      -0.174386606410644  -0.092037014472893
      -0.065037226616579   0.009771974040525
       0.074486430996400  -0.287910597788305
       0.174336742307535  -0.400464726305446
       0.096781997529499  -0.153178844921250
       0.796408810115516   0.881930856163525
       0.005874471676882   0.067750993468743
       0.156654113730103  -0.239182272065197
       0.333688106104332  -0.629954291766549
       0.086388606764696   0.305488995071947
       0.211268899950691  -0.299878322704640
       0.104223240397571   0.199354790284364
       0.336858958710283  -0.620166113933000
       0.145071750152222   0.250136305618056
       0.032242317206686  -0.233223578816564
       0.064558616046395  -0.007577632839606
      -0.055022872335109  -0.190212128003969
      -0.169436992810515  -0.206948615170099
       0.150088939478557   0.090560438773547
       0.256256842094403   0.598874371523965
       0.340882741111244  -0.516594535669757
       0.278186325120137  -0.547968005821968
       0.645979568807173  -0.827959899268083
       0.436535833804569   0.689957746461832
       0.268437499571141   0.341281325011944
       0.120485843238972  -0.301999984465851
       0.160365386980321  -0.202012022307596
       0.154560496835611   0.244912011403144
       0.117966782622546   0.342990099492354
       0.280465408470057  -0.526206432878627
       1.174002195932550   1.136790584581798]')

    L = LinearTransform(X)
    X = transform(L, X)

    nterms = 15

	# For b = "ProHermiteBasis"

    Cnoqrnoprecond, error_noqrnoprecond = greedyfit(10, 2, X, nterms; α = 0.0, withconstant = false,
	                   withqr = false,
                       hessprecond = false, b = "ProHermiteBasis",
                       ATMcriterion = "gradient")

	Cqrnoprecond, error_qrnoprecond = greedyfit(10, 2, X, nterms; α = 0.0, withconstant = false,
	                  withqr = true,
	                  hessprecond = false, b = "ProHermiteBasis",
	                  ATMcriterion = "gradient")

	Cnoqrprecond, error_noqrprecond = greedyfit(10, 2, X, nterms; α = 0.0, withconstant = false,
  	                     withqr = false,
                         hessprecond = true, b = "ProHermiteBasis",
                         ATMcriterion = "gradient")

  	Cqrprecond, error_qrprecond = greedyfit(10, 2, X, nterms; α = 0.0, withconstant = false,
  	                  withqr = true,
  	                  hessprecond = true, b = "ProHermiteBasis",
  	                  ATMcriterion = "gradient")

	error_true =   [1.413939;
					1.410124;
					1.403666;
					1.398928;
					1.398540;
					1.374076;
					1.374064;
					1.368187;
					1.308003;
					1.233256;
					1.230376;
					1.177314;
					1.155683;
					1.138491;
					1.117992]

	idx_true =      [0     0;
				     0     1;
				     0     2;
				     1     0;
				     0     3;
				     0     4;
				     2     0;
				     3     0;
				     4     0;
				     5     0;
				     0     5;
				     0     6;
				     0     7;
				     0     8;
				     6     0]
	coeff_true = [0.105975175137040; -64.829184645818373;   3.644736570204142;  34.333722919244586; -42.641654439858939;  10.491202988540742;   6.170682598059267;
				  15.243539772142837; 13.033494780747679;  3.833643021385606; -25.133002059654828;  12.504145228058636;  -7.625230653938331;   4.946600220335136;   7.531771217294497]

	# For b = "CstProHermiteBasis"
end


text:
"""
    SimpleTestProgram

    Executable test program for ImarisWriter.jl
    Using or importing this module will execute generate three IMS files:
    0. out_0.ims generated by `testConvert(0)` with progress callback
    1. out_1.ims generated by `testConvert(1)` with progress callback
    2. out_2.ims generated by `testConvertActual(2)` without progress callback

    Based loosely upon https://github.com/imaris/ImarisWriterTest/blob/master/testC/bpImarisWriter96TestProgram.c
"""
module SimpleTestProgram
    using ImarisWriter
    # Manually import shorthand types
    import ImarisWriter: Options, Size5D, Parameter, ParameterSection, DimensionSequence5D
    import ImarisWriter: Color, ColorInfo, ImageExtent, Index5D, TimeInfo

    # For progress callback
    struct bpCallbackData
        handle::Ptr{Nothing} # Base.AsyncCondition.handle
        aProgress::bpConverterTypesC_Float # Passed in from callback
        aTotalBytesWritten::bpConverterTypesC_UInt64 # Passed in from callback
        mImageIndex::UInt32 # User parameter
        mProgress::Int32 # User parameter
    end

    """
        progressCallbackActual(vCallbackData::bpCallbackData)

        Actually handles the callback from ImarisWriter and updates status
    """
    function progressCallbackActual(vCallbackData::bpCallbackData)
        # Extract data passed in from callback
        aProgress = vCallbackData.aProgress
        aTotalBytesWritten = vCallbackData.aTotalBytesWritten

        vProgress = round(Int,aProgress * 100)

        if(vProgress - vCallbackData.mProgress < 5)
            return nothing
        end

        vImageIndex = vCallbackData.mImageIndex
        if aTotalBytesWritten < 10 * 1024 * 1024
            println("Progress image $vImageIndex: $vProgress% [$(aTotalBytesWritten ÷ 1024) KB]")
        else
            println("Progress image $vImageIndex: $vProgress% [$(aTotalBytesWritten ÷ (1024*1024)) MB]")
        end
        # vCallbackData.mProgress[] = vProgress
    end

    """
        progressCallback(aProgress, aTotalBytesWritten, aUserData)

        ImarisWriter calls back from an outside thread, and this is not thread safe.
        aUserData is a Ptr{Nothing} that points to a struct that must be isbits && isimmutable
    """
    function progressCallback(aProgress::Float32, aTotalBytesWritten::UInt64, aUserData::Ptr{Nothing})
        # Based on https://github.com/JuliaGPU/OpenCL.jl/blob/master/src/event.jl
        ptr = convert(Ptr{bpCallbackData}, aUserData)
        val = unsafe_load(ptr)
        handle = val.handle

        val = bpCallbackData(val.handle, aProgress, aTotalBytesWritten, val.mImageIndex, val.mProgress)
        unsafe_store!(ptr, val)
        ccall(:uv_async_send, Nothing, (Ptr{Nothing},), handle)
        # ccall(:printf, Cint, (Cstring,Ptr{Nothing}), "Hello World %d \n", handle)
        return nothing
    end

    # C function pointer for callback
    const aProgressCallback = @cfunction(progressCallback, Nothing, (bpConverterTypesC_Float, bpConverterTypesC_UInt64, Ptr{Nothing}) )

    # Callback Task
    function callbackTask(cond, aCallbackUserData)
        run = true
        while(run)
            try
                #println("Waiting... $(cond.handle)")
                wait(cond)
                data = aCallbackUserData[]
                progressCallbackActual(data)
            catch
                run = false
                # rethrow()
            end
        end
        close(cond)
        println("Exiting callback loop for test $aTestIndex")
    end

    """
        testConvert(aTestIndex = 0)

        Test with progress callback. Complex.
    """
    function testConvert(aTestIndex = 0; minimalMetadata = false)
        println("testConvert($aTestIndex; minimalMetadata = $minimalMetadata)")

        # Setup callback
        cond = Base.AsyncCondition()
        # println("This is the beginning $aTestIndex: $(cond.handle)")
        aCallbackUserData = Ref( bpCallbackData(cond.handle, 0.0, 0, aTestIndex, -5) )
        ptr = Base.unsafe_convert(Ptr{Nothing}, aCallbackUserData)

        # Callback task to receive callback from foreign thread
        t = @async callbackTask(cond, aCallbackUserData)

        #Start actual program, preserving cond and aCallbackUserData from garbage collection
        GC.@preserve cond aCallbackUserData testConvertActual(aTestIndex, aCallbackUserData; minimalMetadata)

        # Clean up callback
        progressCallbackActual(aCallbackUserData[])
        close(cond)
        t
    end

    """
        testConvertActual(aTestIndex = 0, aCallbackuserData = C_NULL)

        Test program without progress callback. Simpler.
    """
    function testConvertActual(aTestIndex = 0, aCallbackUserData = C_NULL; minimalMetadata = false)
        ### Setup

        aDataType = UInt8
        aImageSize = Size5D(X = 512, Y = 512, Z = 32, C = 4, T = 2)
        vVoxelSizeXY = 1.4
        vVoxelSizeZ = 5.4
        aImageExtent = ImageExtent(
            X = (0, aImageSize.mValueX * vVoxelSizeXY),
            Y = (0, aImageSize.mValueY * vVoxelSizeXY),
            Z = (0, aImageSize.mValueZ)
        )
        aSample = Size5D()
        aDimensionSequence = DimensionSequence5D(:X, :Y, :Z, :C, :T)
        aBlockSize = Size5D(X = 256, Y = 256, Z = 8)
        aOutputFile = "./out_$aTestIndex.ims"

        # ImarisWriter Options
        aOptions = Options(
            mEnableLogProgress = aCallbackUserData != C_NULL, # Disable callback logging if second argument is a null pointer
            mNumberOfThreads = length(Sys.cpu_info()) ÷ 2, # Uses half of the logical processors, which should be # of cores
            mDisablePyramid = false # Create multiresolution pyramid
        )

        aApplicationName = "simpleTestProgram.jl"
        aApplicationVersion = "1.0.0"

        ### Create File and Handle

        vConverter = ImarisWriter.create(
            aDataType, aImageSize, aSample,
            aDimensionSequence, aBlockSize,
            aOutputFile, aOptions,
            aApplicationName, aApplicationVersion,
            aProgressCallback, aCallbackUserData
        )
        ImarisWriter.checkErrors(vConverter)

        ### Write Blocks

        vBlockSize = aBlockSize.mValueX * aBlockSize.mValueY * aBlockSize.mValueZ * aBlockSize.mValueC * aBlockSize.mValueT
        vData = Vector{UInt8}(undef, vBlockSize)
        for vIndex in eachindex(vData)
            vData[vIndex] = mod(vIndex-1, 256)
        end

        vNBlocksX = ImarisWriter.numBlocks(aImageSize.mValueX, aBlockSize.mValueX)
        vNBlocksY = ImarisWriter.numBlocks(aImageSize.mValueY, aBlockSize.mValueY)
        vNBlocksZ = ImarisWriter.numBlocks(aImageSize.mValueZ, aBlockSize.mValueZ)
        vNBlocksC = ImarisWriter.numBlocks(aImageSize.mValueC, aBlockSize.mValueC)
        vNBlocksT = ImarisWriter.numBlocks(aImageSize.mValueT, aBlockSize.mValueT)

        aBlockIndex = Index5D()

        for aBlockIndex.mValueC = 0:vNBlocksC-1,
            aBlockIndex.mValueT = 0:vNBlocksT-1,
            aBlockIndex.mValueZ = 0:vNBlocksZ-1,
            aBlockIndex.mValueY = 0:vNBlocksY-1,
            aBlockIndex.mValueX = 0:vNBlocksX-1

            ImarisWriter.copyBlock(vConverter, vData, aBlockIndex)
            ImarisWriter.checkErrors(vConverter)
            yield()
        end

        ### Write Metadata

        if !minimalMetadata

            vNumberOfOtherSections = 1
            vNumberOfSections = vNumberOfOtherSections + aImageSize.mValueC
            vParameterSections = Vector{ParameterSection}(undef,vNumberOfSections)
            vUnitParameter = Parameter("Unit", "um")
            vImageSection = ParameterSection("Image", vUnitParameter)
            vParameterSections[1] = vImageSection
            vNumberOfParametersPerChannel = 3
            for vC = 1:aImageSize.mValueC
                vChannelParameters = Vector{Parameter}(undef, vNumberOfParametersPerChannel)
                vChannelParameters[1] = Parameter("Name", vC == 0 ? "First channel" : vC == 1 ? "Second channel" : vC == 2 ? "Third channel" : "Other channel")
                vChannelParameters[2] = Parameter("LSMEmissionWavelength", "700")
                vChannelParameters[3] = Parameter("OtherChannelParameter", "OtherChannelValue")
                vParameterSections[vNumberOfOtherSections + vC] = ParameterSection("Channel $vC", vChannelParameters)
            end

            #aParameters = bpConverterTypesC_Parameters(vParameterSections, vNumberOfSections)

            vTimeInfos = Vector{TimeInfo}(undef,aImageSize.mValueT)
            for vT = 1:aImageSize.mValueT
                vSeconds = vT + 4 + 60 * (27 + 60 *15)
                vTimeInfos[vT] = TimeInfo( 2458885, vSeconds * 1_000_000_000 )
            end
            # Convert vTimesInfos to a pointer
            #aTimeInfoPerTimePoint = bpConverterTypesC_TimeInfos(vTimeInfos,  aImageSize.mValueT)

            vColorInfos = Vector{ColorInfo}(undef, aImageSize.mValueC)
            for vC = 1:aImageSize.mValueC
                vColorInfos[vC] = ColorInfo(
                    Color(
                        (vC - 1 % 3) == 0 ? 1 : 0, # mRed
                        (vC - 1 % 3) == 1 ? 1 : 0, # mGreen
                        (vC - 1 % 3) == 2 ? 1 : 0, # mBlue
                        1
                    ) # mBaseColor
                )
            end
        else
            # Minimal Metadata requires the next four lines
            vParameterSections = ParameterSection[]
            vTimeInfos = TimeInfo[]
            vColorInfos = fill(ColorInfo(), aImageSize.mValueC)
        end

        aAutoAdjustColorRange = true

        ImarisWriter.finish(vConverter,
            aImageExtent, vParameterSections, vTimeInfos,
            vColorInfos, aAutoAdjustColorRange
        )
        ImarisWriter.checkErrors(vConverter)

        # This is really important, file will not be written until this is called
        ImarisWriter.destroy(vConverter)

    end

    function __init__()
        testConvert(0)
        testConvert(1; minimalMetadata = true)
        println("No callback test")
        testConvertActual(2)
    end

end


text:
function _is_cell_done(cell)
    if cell.metadata["disabled"]
        return true
    else
        return !cell.queued && !cell.running
    end
end

function nothingstring(x::Union{Nothing,AbstractString})::Union{Nothing,String}
    return x isa Nothing ? x : string(x)::String
end

@enum OutputFormat begin
    documenter_output
    franklin_output
    html_output
end

const WRITE_FILES_DEFAULT = true
const PREVIOUS_DIR_DEFAULT = nothing
const OUTPUT_FORMAT_DEFAULT = html_output
const ADD_DOCUMENTER_CSS_DEFAULT = true
const USE_DISTRIBUTED_DEFAULT = true
const MAX_CONCURRENT_RUNS_DEFAULT = 4

"""
    BuildOptions(
        dir::AbstractString;
        write_files::Bool=$WRITE_FILES_DEFAULT,
        previous_dir::Union{Nothing,AbstractString}=$PREVIOUS_DIR_DEFAULT,
        output_format::OutputFormat=$OUTPUT_FORMAT_DEFAULT,
        add_documenter_css::Bool=$ADD_DOCUMENTER_CSS_DEFAULT,
        use_distributed::Bool=$USE_DISTRIBUTED_DEFAULT,
        max_concurrent_runs::Int=$MAX_CONCURRENT_RUNS_DEFAULT
    )

Arguments:

- `dir`:
    Directory in which the Pluto notebooks are stored.
- `write_files`:
    Write files to `joinpath(dir, "\$file.html")`.
- `previous_dir::Union{Nothing,AbstractString}=Nothing`:
    Use the output from the previous run as a cache to speed up running time.
    To use the cache, specify a directory `previous_dir::AbstractString` which contains HTML or Markdown files from a previous run.
    Specifically, files are expected to be at `joinpath(previous_dir, "\$file.html")`.
    The output from the previous run may be embedded in a larger HTML or Markdown file.
    This package will extract the original output from the full file contents.
    By default, caching is disabled.
-  `output_format`:
    What file to write the output to.
    By default this is `html_output::OutputFormat` meaning that the output of the HTML method is pure HTML.
    To generate Franklin or Documenter files, use respectively `franklin_output` or `documenter_output`.
    When `BuildOptions.write_files == true` and `output_format != html_output`, the output file has a ".md" extension instead of ".html".
- `add_documenter_css` whether to add a CSS style to the HTML when `documenter_output=true`.
- `use_distributed`:
    Whether to build the notebooks in different processes.
    By default, this is enabled just like in Pluto and the notebooks are build in parallel.
    The benefit of different processes is that things are more independent of each other.
    Unfortunately, the drawback is that compilation has to happen for each process.
    By setting this option to `false`, all notebooks are built sequentially in the same process which avoids recompilation.
    This is likely quicker in situations where there are few threads available such as GitHub Runners depending on the notebook contents.
    Beware that `use_distributed=false` will not work with Pluto's built-in package manager.
- `max_concurrent_runs`:
    Maximum number of notebooks to evaluate concurrently when `use_distributed=true`.
    Note that each notebook starts in a different thread and can start multiple threads, so don't set this number too high.
"""
struct BuildOptions
    dir::String
    write_files::Bool
    previous_dir::Union{Nothing,String}
    output_format::OutputFormat
    add_documenter_css::Bool
    use_distributed::Bool
    max_concurrent_runs::Int

    function BuildOptions(
            dir::AbstractString;
            write_files::Bool=WRITE_FILES_DEFAULT,
            previous_dir::Union{Nothing,AbstractString}=PREVIOUS_DIR_DEFAULT,
            output_format::OutputFormat=OUTPUT_FORMAT_DEFAULT,
            add_documenter_css::Bool=ADD_DOCUMENTER_CSS_DEFAULT,
            use_distributed::Bool=USE_DISTRIBUTED_DEFAULT,
            max_concurrent_runs::Int=MAX_CONCURRENT_RUNS_DEFAULT
        )
        return new(
            string(dir)::String,
            write_files,
            nothingstring(previous_dir),
            output_format,
            add_documenter_css,
            use_distributed,
            max_concurrent_runs
        )
    end
end

"""
    _is_notebook_done(notebook::Notebook)

Return whether all cells in the `notebook` have executed.
This method is more reliable than using `notebook.executetoken` because Pluto.jl drops that lock also after installing packages.
"""
function _notebook_done(notebook::Notebook)
    cells = [last(elem) for elem in notebook.cells_dict]
    return all(_is_cell_done, cells)
end

function _wait_for_notebook_done(nb::Notebook)
    while !_notebook_done(nb)
        sleep(1)
    end
end

function extract_previous_output(html::AbstractString)::String
    start_range = findfirst(BEGIN_IDENTIFIER, html)
    @assert !isnothing(start_range)
    start = first(start_range)

    stop_range = findfirst(END_IDENTIFIER, html)
    @assert !isnothing(stop_range)
    stop = last(stop_range)

    return html[start:stop]
end

"""
    Previous(state::State, html::String)

- `state`: Previous state.
- `text`: Either
    - HTML starting with "$BEGIN_IDENTIFIER" and ending with "$END_IDENTIFIER" or
    - Contents of a Franklin/Documenter Markdown file.
"""
struct Previous
    state::Union{State,Nothing}
    text::String
end

function Previous(text::String)
    state = contains(text, STATE_IDENTIFIER) ?
        extract_state(text) :
        nothing
    text = contains(text, BEGIN_IDENTIFIER) ?
        extract_previous_output(text) :
        ""
    return Previous(state, text)
end

function Previous(bopts::BuildOptions, in_file)
    prev_dir = bopts.previous_dir
    if isnothing(prev_dir)
        return Previous(nothing, "")
    end
    name, _ = splitext(in_file)
    ext = bopts.output_format == html_output ? ".html" : ".md"
    prev_path = joinpath(prev_dir, "$name$ext")
    if !isfile(prev_path)
        return Previous(nothing, "")
    end
    text = read(prev_path, String)
    return Previous(text)
end

function reuse_previous(previous::Previous, dir, in_file)::Bool
    in_path = joinpath(dir, in_file)
    curr = path2state(in_path)

    prev = previous.state
    isnothing(prev) && return false

    sha_match = prev.input_sha == curr.input_sha
    julia_match = prev.julia_version == curr.julia_version
    reuse = sha_match && julia_match
    return reuse
end

"""
Write to a ".html" or ".md" file depending on `HTMLOptions.output_format`.
The output file is always a sibling to the file at `in_path`.
"""
function _write_main_output(in_path, text, bopts::BuildOptions, hopts::HTMLOptions)
    ext = bopts.output_format == html_output ? ".html" : ".md"
    if bopts.write_files
        dir = dirname(in_path)
        in_file = basename(in_path)
        without_extension, _ = splitext(in_file)
        out_file = "$(without_extension)$(ext)"
        out_path = joinpath(dir, out_file)
        write(out_path, text)
    end
    return nothing
end

"Used when creating the page for the first time and to restore the cache."
function _wrap_franklin_output(html)
    return "~~~\n$(html)\n~~~"
end

"Used when creating the page for the first time and when restoring the cache."
function _wrap_documenter_output(html::String, bopts::BuildOptions, in_path::String)
    if bopts.add_documenter_css
        html = _add_documenter_css(html)
    end
    editurl = _editurl_text(bopts, in_path)
    return """
        ```@raw html
        $html
        ```
        $editurl
        """
end

function _outcome2text(session, prev::Previous, in_path::String, bopts, hopts)::String
    text = prev.text
    if bopts.output_format == franklin_output
        text = _wrap_franklin_output(text)
    end
    if bopts.output_format == documenter_output
        text = _wrap_documenter_output(text, bopts, in_path)
    end
    _write_main_output(in_path, text, bopts, hopts)
    return text
end

function _inject_script(html, script)
    l = length(END_IDENTIFIER)
    without_end = html[1:end-l]
    return string(without_end, '\n', script, '\n', END_IDENTIFIER)
end

function _outcome2text(session, nb::Notebook, in_path::String, bopts, hopts)::String
    _throw_if_error(session, nb)

    # Grab output before changing binds via `_run_dynamic!`.
    # Otherwise, the outputs look wrong when opening a page for the first time.
    html = notebook2html(nb, in_path, hopts)

    if bopts.output_format == franklin_output
        html = _wrap_franklin_output(html)
    end
    if bopts.output_format == documenter_output
        html = _wrap_documenter_output(html, bopts, in_path)
    end

    _write_main_output(in_path, html, bopts, hopts)

    # The sleep avoids `AssertionError: will_run_code(notebook)`
    @async begin
        sleep(5)
        SessionActions.shutdown(session, nb; verbose=false)
    end

    return string(html)::String
end

const TimeState = Dict{String,DateTime}

_time_init!(time_state::TimeState, in_file::String) = setindex!(time_state, now(), in_file)
_time_elapsed(time_state::TimeState, in_file::String) = now() - time_state[in_file]

"Return `... minutes and ... seconds`."
function _pretty_elapsed(t::Millisecond)
    value = t.value
    seconds = Float64(value) / 1000
    min, sec = divrem(seconds, 60)
    rmin = round(Int, min)
    min_text = rmin == 0 ? "" :
        rmin == 1 ? "$rmin minute and " :
        "$rmin minutes and "
    rsec = round(Int, sec)
    sec_text = rsec == 1 ? "$rsec second" : "$rsec seconds"
    return string(min_text, sec_text)
end

function _evaluate_file(bopts::BuildOptions, hopts::HTMLOptions, session, in_file, time_state)
    dir = bopts.dir
    in_path = joinpath(dir, in_file)::String
    @assert isfile(in_path) "Expected .jl file at $in_path"

    previous = Previous(bopts, in_file)
    if reuse_previous(previous, dir, in_file)
        @info "Using cache for Pluto notebook at $in_file"
        return previous
    else
        @info "Starting evaluation of Pluto notebook $in_file"
        if bopts.use_distributed
            run_async = true
        else
            run_async = false
            # `use_distributed` means mostly "whether to run in a new process".
            session.options.evaluation.workspace_use_distributed = false
        end
        nb = run_notebook!(in_path, session; hopts, run_async)
        if bopts.use_distributed
            # The notebook is running in a distributed process, but we still need to wait to
            # avoid spawning too many processes in `_evaluate_parallel`.
            _wait_for_notebook_done(nb)
        end
        elapsed = _time_elapsed(time_state, in_file)
        pretty_elapsed = _pretty_elapsed(elapsed)
        @info "Finished evaluation of Pluto notebook $in_file in $pretty_elapsed"
        return nb
    end
end

"""
Evaluate `files` in parallel.

Using asynchronous tasks instead of multi-threading since Downloads is not thread-safe on Julia 1.6/1.7.
https://github.com/JuliaLang/Downloads.jl/issues/110.
"""
function _evaluate_parallel(bopts, hopts, session, files, time_state)
    ntasks = bopts.max_concurrent_runs
    asyncmap(files; ntasks) do in_file
        _time_init!(time_state, in_file)
        _evaluate_file(bopts, hopts, session, in_file, time_state)
    end
end

function _evaluate_sequential(bopts, hopts, session, files, time_state)
    return map(files) do in_file
        _time_init!(time_state, in_file)
        _evaluate_file(bopts, hopts, session, in_file, time_state)
    end
end

"""
    build_notebooks(
        bopts::BuildOptions,
        files,
        hopts::HTMLOptions=HTMLOptions();
        session=ServerSession()
    ) -> Vector{String}

Build all `files` in `dir` in parallel.
"""
function build_notebooks(
        bopts::BuildOptions,
        files,
        hopts::HTMLOptions=HTMLOptions();
        session=ServerSession()
    )::Vector{String}

    time_state = TimeState()
    func = bopts.use_distributed ? _evaluate_parallel : _evaluate_sequential
    X = func(bopts, hopts, session, files, time_state)

    H = map(zip(files, X)) do (in_file, x)
        in_path = joinpath(bopts.dir, in_file)
        html = _outcome2text(session, x, in_path, bopts, hopts)
        return html
    end

    return H
end

function _is_pluto_file(path::AbstractString)::Bool
    first(eachline(string(path))) == "### A Pluto.jl notebook ###"
end

"""
    build_notebooks(
        bopts::BuildOptions,
        hopts::HTMLOptions=HTMLOptions()
    ) -> Vector{String}

Build all ".jl" files in `dir`.
"""
function build_notebooks(
        bopts::BuildOptions,
        hopts::HTMLOptions=HTMLOptions()
    )::Vector{String}
    dir = bopts.dir
    files = filter(readdir(dir)) do file
        path = joinpath(dir, file)
        endswith(file, ".jl") && _is_pluto_file(path) && !startswith(file, TMP_COPY_PREFIX)
    end
    return build_notebooks(bopts, files, hopts)
end


text:
function test_utils(::Type{T}) where T <: CasadiSymbolicObject
    @testset "$( string("Function 'substitute' for ", T, "                     ") )" begin
        x = T("x")
        y = T("y")
        z = T("z")
        f = x + 2y + z

        @test casadi.is_equal(substitute(f, x, y), y + 2y + z, 3)
        @test casadi.is_equal(substitute(f, z, 1), x + 2y + 1, 3)
        @test to_julia( substitute(f, [y z x], [-1 0 1]) ) == -1
        @test to_julia( substitute(f, [y; z; x], [2; 0; -4]) ) == 0
    end
end


text:
using Pkg

Pkg.activate(joinpath(@__DIR__, "..")); Pkg.instantiate()
Pkg.activate(); Pkg.instantiate()

pushfirst!(LOAD_PATH, joinpath(@__DIR__, ".."))

using Documenter, LSHFunctions

makedocs(
    sitename = "LSHFunctions.jl",
    format   = Documenter.HTML(),
    modules  = [LSHFunctions],
    pages    = ["Home" => "index.md",
                "The LSHFunction API" => "lshfunction_api.md",
                "Similarity statistics" => [
                    "Cosine similarity" => joinpath("similarities", "cosine.md"),
                    "``\\ell^p`` distance" => joinpath("similarities", "lp_distance.md"),
                    "Jaccard similarity" => joinpath("similarities", "jaccard.md"),
                    "Inner product similarity" => joinpath("similarities", "inner_prod.md")],
                "Function-space hashing" => "function_hashing.md",
                "Performance tips" => "performance.md",
                "FAQ" => "faq.md",
                "Glossary" => "glossary.md",
                "API reference" => "full_api.md",
               ]
)

deploydocs(
    repo = "github.com/kernelmethod/LSHFunctions.jl.git"
)


text:
export data,
       TensorStorage,
       randn!,
       scale!,
       norm

abstract type TensorStorage{ElT} <: AbstractVector{ElT} end

data(S::TensorStorage) = S.data

Base.eltype(::TensorStorage{ElT}) where {ElT} = ElT

Base.eltype(::Type{<:TensorStorage{ElT}}) where {ElT} = ElT

Base.iterate(S::TensorStorage,args...) = iterate(data(S),args...)

dense(S::TensorStorage) = S

# This is necessary since for some reason inference doesn't work
# with the more general definition (eltype(Nothing) === Any)
Base.eltype(::TensorStorage{Nothing}) = Nothing

Base.length(S::TensorStorage) = length(data(S))

Base.size(S::TensorStorage) = size(data(S))

Base.@propagate_inbounds Base.getindex(S::TensorStorage,
                                       i::Integer) = data(S)[i]
Base.@propagate_inbounds Base.setindex!(S::TensorStorage,v,
                                        i::Integer) = (setindex!(data(S),v,i); S)

# Needed for passing Tensor{T,2} to BLAS/LAPACK
function Base.unsafe_convert(::Type{Ptr{ElT}},
                             T::TensorStorage{ElT}) where {ElT}
  return Base.unsafe_convert(Ptr{ElT},data(T))
end

# This may need to be overloaded, since storage types
# often have other fields besides data

Base.conj!(S::TensorStorage) = (conj!(data(S)); return S)

function Base.conj(S::T;always_copy = false) where {T<:TensorStorage}
  if always_copy
    return conj!(copy(S))
  end
  return T(conj(data(S)))
end

Base.complex(S::T) where {T<:TensorStorage} = complex(T)(complex(data(S)))

Base.copyto!(S1::TensorStorage,
             S2::TensorStorage) = (copyto!(data(S1),data(S2)); S1)

Random.randn!(S::TensorStorage) = (randn!(data(S)); S)

Base.fill!(S::TensorStorage,v) = (fill!(data(S),v); S)

LinearAlgebra.rmul!(S::TensorStorage,v::Number) = (rmul!(data(S),v); S)
scale!(S::TensorStorage,v::Number) = rmul!(S,v)

LinearAlgebra.norm(S::TensorStorage) = norm(data(S))

Base.convert(::Type{T},S::T) where {T<:TensorStorage} = S

blockoffsets(S::TensorStorage) = S.blockoffsets

"""
nzblocks(T::TensorStorage)

Return a vector of the non-zero blocks of the BlockSparse storage.
"""
nzblocks(T::TensorStorage) = nzblocks(blockoffsets(T))

nnzblocks(S::TensorStorage) = length(blockoffsets(S))
nnz(S::TensorStorage) = length(S)

offset(S::TensorStorage,
       block) = offset(blockoffsets(S), block)

nzblock(S::TensorStorage,
        n::Int) = nzblock(blockoffsets(S), n)


text:
"""
    find_global_maximum_dfs(graph::Dict{Int64,Array{Int64,1}})

Find the total global maximum based on a comparsion between the intial vertex
(`global_maximum = 0`), the current dictionary key, and the accesible vertexes from
the array-list (value).

# Arguments
- `graph::Dict{Int64,Array{Int64,1}}`: Graph of the connected nodes
"""
function find_global_maximum_dfs(graph::Dict{Int64,Array{Int64,1}})
    global_maximum = 0
    # Going through the dictionary
    for (key, value) in graph
        # Check the keys as new global maximum
        if key > global_maximum
            global_maximum = key
        end

        # Check the values as new local maximum
        local_maximum = maximum(value)

        # Check if new local maximum becomes global maximum
        if local_maximum > global_maximum
            global_maximum = local_maximum
        end
    end
    return global_maximum
end

"""
    dfs_recrusive(
        graph::Dict{Int64,Array{Int64,1}},
        start::Int64,
        visited::Array{Bool,1},
        solution::Array{Int64,1},
    )

The recursively call of the dfs_recrusive is essential for exploring each single branch
of the graph.

# Arguments
- `graph::Dict{Int64,Array{Int64,1}}`: Graph of the connected nodes
- `start::Int64`: Startpoint (first selected vertex) of the graph-traveling process
- `visited::Array{Bool,1}`: Visited nodes of the graph
- `solution::Array{Int64,1}`: Solution of the raph-traveling proces
"""
function dfs_recrusive(
    graph::Dict{Int64,Array{Int64,1}},
    start::Int64,
    visited::Array{Bool,1},
    solution::Array{Int64,1},
)
    # Mark the current node as visited and enqueue it
    visited[start] = true
    append!(solution, start)

    # Recur for all the vertices adjacent to this vertex
    for i in graph[start]
        if visited[i] == false
            dfs_recrusive(graph, i, visited, solution)
        end
    end
    return solution

end

"""
    depth_first_search(graph::Dict{Int64,Array{Int64,1}}, start::Int64)

The depth-first search (DFS) is an algorithm dedicated to traversing or searching for
tree or graph data structures. It starts at a specified tree root (start) for exploring
as far as possible at each branch. Afterthat the BFS starts automatically backtracking.
DFS is stack-based.
For more information see: [https://en.wikipedia.org/wiki/Depth-first_search](https://en.wikipedia.org/wiki/Depth-first_search)

# Arguments
- `graph::Dict{Int64,Array{Int64,1}}`: Graph of the connected nodes
- `start::Int64`: Startpoint (first selected vertex) of the graph-traveling process

# Examples
```julia-repl
julia> import ClassicAlgorithmsCollections
julia> graph = Dict(1=> [2, 3], 2=> [3], 3=> [1, 4], 4=> [4])
julia> ClassicAlgorithmsCollections.depth_first_search(graph, 3)
[3, 1, 2, 4]
```
"""
function depth_first_search(graph::Dict{Int64,Array{Int64,1}}, start::Int64)

    # Mark all the vertices as not visited
    n = find_global_maximum_dfs(graph)
    visited = zeros(Bool, n)

    # Create array-list for the searching solution
    solution = zeros(Int64, false)

    result = dfs_recrusive(graph, start, visited, solution)
    return result
end

max_stars_repo_name: Anselmoo/ClassicAlgorithmsCollections

max_issues_repo_name: Anselmoo/ClassicAlgorithmsCollections

max_forks_repo_name: Anselmoo/ClassicAlgorithmsCollections