text:
module WeirdDetector

    using Interpolations
    using DataFrames
    using FITSIO
    using Optim
    using Statistics
    using Distributed
    using Random

    export Point, periodogram, aliases, optimal_periods, getFITS, loadFITS, pointsify, flatten, scrambled_periodogram, movingstd

    struct Point
        t :: Float32
        modt :: Float32
        F :: Float32
        sigmaF :: Float32
        smoothedF :: Float32
    end
    Point(t::Real, F::Real, sigmaF::Real) = Point(Float32(t), NaN, Float32(F), Float32(sigmaF), NaN)
    Point(t::Real, F::Real) = Point(t, F, 5f-5)

    """
    convert a DataFrame in the format returned by `loadFITS`
    to an array of `Points`
    """
    function pointsify(df; keep_interpolated::Bool=false) :: Vector{Point}
        if keep_interpolated
            df[df[:interpolated] .== true, :sigmaF] = mean(df[df[:interpolated] .== false, :sigmaF])
        else
           df = df[df[:interpolated] .== false, :]
        end
        Point.((df[:t]), (df[:F]), (df[:sigmaF]))
    end

    function fold!(period::Float32, data::Array{Point})
        for (i,p) in enumerate(data)
            data[i] = Point(p.t, p.t%period, p.F, p.sigmaF, p.smoothedF)
        end
        sort!(data, by=(x -> x.modt))
    end

    "helper function for smoothing"
    function access(data::Vector{A}, i::Int) :: A where A
        npoints = length(data)
        j = mod(i,npoints) == 0 ? npoints : mod(i,npoints)
        data[j]
    end

    "helper function for smoothing"
    function wrappedtime(data::Vector{Point}, i::Int, P::Float32) :: Float32
        npoints = length(data)
        modt = access(data, i).modt
        if i < 1
            return modt - P
        elseif i > npoints
            return modt + P
        else
            return modt
        end
    end

    function smooth2!(data::Vector{Point}, P::Float32, kw::Float32)
        kw = Int(round(kw * length(data)))
        s = map(collect(-kw:kw)) do i
            p = access(data, i)
            (p.F / p.sigmaF, 1/p.sigmaF)
        end
        sumFoverSigma = sum(first.(s))
        sumInverseSigma = sum(last.(s))
        lb = -kw
        ub = kw
        for i = 1:length(data)
            point = access(data, lb)
            sumFoverSigma -= point.F / point.sigmaF
            sumInverseSigma -= 1e0 / point.sigmaF
            lb += 1

            ub += 1
            point = access(data, ub)
            sumFoverSigma += point.F / point.sigmaF
            sumInverseSigma += 1e0 / point.sigmaF

            data[i] = Point(data[i].t, data[i].modt, data[i].F, data[i].sigmaF,
                            Float32(sumFoverSigma/sumInverseSigma))
        end
    end

    "smooth the data with a rolling mean"
    function smooth!(data::Vector{Point}, P::Float32, kw::Float32)
        kw = kw * P
        sumFoverSigma = 0e0
        sumInverseSigma = 0e0
        lb = 2
        ub = 1
        while wrappedtime(data, lb, P) > data[1].modt - kw
            lb -= 1
            point = access(data, lb)
            sumFoverSigma += point.F / point.sigmaF
            sumInverseSigma += 1e0 / point.sigmaF
        end
        for i = 1:length(data)
            while wrappedtime(data, lb, P) < data[i].modt - kw
                point = access(data, lb)
                sumFoverSigma -= point.F / point.sigmaF
                sumInverseSigma -= 1e0 / point.sigmaF
                lb += 1
            end
            while wrappedtime(data, ub, P) < data[i].modt + kw
                ub += 1
                point = access(data, ub)
                sumFoverSigma += point.F / point.sigmaF
                sumInverseSigma += 1e0 / point.sigmaF
            end
            data[i] = Point(data[i].t, data[i].modt, data[i].F, data[i].sigmaF,
                            Float32(sumFoverSigma/sumInverseSigma))
        end
    end

    "calculate reduced chi squared on data that has already been folded and smoothed"
    function chi2(data::Vector{Point}) :: Float32
        chi2 = 0e0
        for p in data
            chi2 += ((p.F - p.smoothedF)/p.sigmaF)^2
        end
        Float32(chi2)
    end

    "fold, smooth, and calculate reduced chi squared"
    function chi2(data::Vector{Point}, period::Float32; kw::Float32=0.002f0)
        fold!(period, data)
        smooth!(data, period, kw)
        chi2(data)
    end

    function kurtosis(data::Vector{Point}) :: Float32
        fs = (x->x.smoothedF).(data)
        mu = mean(fs)
        sigma = std(fs)
        mean((fs .- mu).^4)/sigma^4
    end

    "returns dataframe containing chi-squared and kurtosis by period"
    function periodogram(data::Vector{Point}, periods::Vector{Float32}; kw=0.002f0,
                         parallel=true, datakw=false, postprocess=true)
        df = DataFrame(chi2=Float32[], kurtosis=Float32[])
        s = div(length(periods), nworkers()*3)
        results = pmap(periods, distributed=parallel, batch_size=s) do p
            fold!(p, data)
            row = Vector{Any}()
            if datakw
                smooth2!(data, p, kw)
            else
                smooth!(data, p, kw)
            end
            push!(row, chi2(data))
            push!(row, kurtosis(data))
            push!(df, row)
        end
        df
    end

    function aliases(downto::Int=5; upperBound=nothing) :: Vector{Rational}
        lines = Set{Rational}()
        for n in 1:downto
            for m in 1:n
                if upperBound == nothing
                    push!(lines, m//n)
                else
                    i = 0
                    while i + m//n < upperBound
                        push!(lines, i + m//n)
                        i += 1
                    end
                end
            end
        end
        collect(lines)
    end

    function flatten(periods::Vector{Float32}, chi2s::Vector{Float32}; stepwidth::Float32=10.21f0,
                     preflipped=false, tess=false)
        nchi2s = similar(chi2s)
        if (!tess)
            steps = Int(ceil(periods[end]/stepwidth))
            for n in 1:steps
                lb = first(searchsorted(periods, (n-1)*stepwidth))
                ub = last(searchsorted(periods, n*stepwidth))
                line(x, p) = p[1]*x + p[2]
                c2(p) = sum((line(periods[lb:ub], p) .- chi2s[lb:ub]).^2)
                p0 = [0e0, 0e0]
                p0[1] = (chi2s[ub] - chi2s[lb])/(periods[ub] - periods[lb])
                p0[2] = chi2s[lb] - p0[1]*periods[lb]
                res = Optim.optimize(c2, p0, Newton())
                p = Optim.minimizer(res)
                nchi2s[lb:ub] .= Float32.(line(periods[lb:ub], p) - chi2s[lb:ub])
            end
        else
            medianchi2s = fill(median(chi2s), length(chi2s))
            nchi2s = medianchi2s - chi2s
        end
        if preflipped
            -nchi2s
        else
            nchi2s
        end
    end

    function movingstd(npowers::Vector{Float32}, kw=200) :: Vector{Float32}
        sigmas = similar(npowers)
        for i in 1:length(npowers)
            lb = i - kw < 1 ? 1 : i - kw
            ub = i + kw > length(npowers) ? length(npowers) : i + kw
            sigmas[i] = std(npowers[lb:ub])
        end
        sigmas
    end

    function scrambled_periodogram(df::DataFrame, periods::Vector{Float32}; tess=false, kwargs...)
        df = deepcopy(df)
        dropmissing!(df, :F, disallowmissing=false)
        df[:F] = df[randperm(size(df)[1]), :F]
        df[:sigmaF] = df[randperm(size(df)[1]), :sigmaF]
        data = pointsify(df)
        periodogram(data, periods; kwargs...)
    end

    function optimal_periods(pmin=0.25f0, pmax=5f1; n=6)
        pmin = Float32(pmin)
        pmax = Float32(pmax)
       exp.((log(pmin) : (0.001f0/n) : log(pmax)))
    end

    "download FITS data for KIC"
    getFITS(KIC::Int; fitsdir="fitsfiles/") = getFITS(lpad(KIC,9,0); fitsdir=fitsdir)
    function getFITS(KIC::String; fitsdir="fitsfiles/", force=false, tic_id="", ete6=false)
        if tic_id == ""
            #info("Dowloading data for KIC $(KIC)")
            ftpfolder = "http://archive.stsci.edu/pub/kepler/lightcurves//" * KIC[1:4] * "/" * KIC * "/"
        elseif ete6
            tic_id = lpad(tic_id, 9, 0)
            ftpfolder = "http://archive.stsci.edu/missions/tess/ete-6/tid/00/000/00" * tic_id[1:1] * "/" * tic_id[2:4] * "/tess2019128220341-" * lpad(tic_id,16,0) * "-0016-s_lc.fits"
        else
            tic_id = lpad(tic_id, 9, 0)     #good for sector 1 data only
            sector = lpad(sector, 2, 0)
            ftpfolder = "https://mast.stsci.edu/api/v0.1/Download/file/?uri=mast:TESS/product/tess2018206045859-s0001-" * lpad(tic_id,16,0) * "-0120-s_lc.fits"
        end
        command = `wget -P $fitsdir -nH --cut-dirs=6 -r -c -N -np -q -R '*.tar' -R 'index*' -erobots=off $ftpfolder`
        run(command)
    end

    #window width is kw*2+1
    "remove outliers from time series.  This is a hleper function for loadFITS"
    function prune(df::DataFrame, threshold_sigma::Int, kw::Int; usephoto=false) :: DataFrame
        nrows = size(df,1)
        med = map(1:nrows) do i
            lb = i-kw < 1 ? 1 : i-kw
            ub = i+kw < nrows ? i+kw : nrows
            m = median(df[:F][lb:ub])
            m, median(abs.(m .- df[:F][lb:ub]))
        end
        df = df[1+kw : end-kw, :]
        med = med[1+kw : end-kw]
        sigma = (m->m[2]).(med)
        med = (m->m[1]).(med)
        if usephoto
            df[abs.(df[:F] .- med) .< threshold_sigma.*df[:sigmaF], :]
        else
            df[abs.(df[:F] .- med) .< threshold_sigma.*sigma.*1.4826, :]
        end
    end

    "detrend with moving median.  This is a helper function for loadFITS"
    function detrend(df::DataFrame, P, trim=true, func=median)
        m = func(df[:F])
        df[:F] .-= m
        fs = 60*24 / 29.4 #sampling frequency [days^-1]
        kw = Int(round(P*fs/2))
        nrows = size(df, 1)
        smoothedF = Vector{Float32}(undef, nrows)
        for i in 1:nrows
            lb = i-kw < 1 ? 1 : i-kw
            ub = i+kw < nrows ? i+kw : nrows
            smoothedF[i] = func(df[:F][lb:ub])
        end
        df[:F] .= df[:F] .+ m .- smoothedF
        if trim
            df[1+kw: end-kw, :]
        else
            df
        end
    end

    "Interpolate missing points from LC.  This is a helper function for loadFITS"
    function interpolate_missing!(df::DataFrame; tess::Bool=false)
        #collect points to interpolate
        df[:interpolated] = false
        newTimes = Vector{Float64}()
        dt = 0.020416 #kepler cadence in days
        if (tess)
            dt = 0.000694444 #tess cadence in days
        end
        epsilon = 0.005
        for i in 2:size(df,1)
            t = df[i-1, :t]
            while df[i,:t] - t > dt + epsilon
                t += dt
                push!(newTimes, t)
            end
        end
        itp = Interpolations.interpolate((Vector{Float32}(df[:t]),),
                                         Vector{Float32}(df[:F]), Gridded(Linear()))
        for t in newTimes
            if (!tess)
                row = [t, itp(t), missing, true]
            else
                row = [t, itp(t), -1, true]
            end
            push!(df, row)
        end
        #sort in interpolated points
        sort!(df, [:t])
    end

    """
    returns a DataFrame containing the lightcurve across all quarters for KIC
    with outliers removed and missing and bad chi2s = chi2(data, periods) points interpolated.
    """
    loadFITS(KIC::Int; kwargs...) = loadFITS(lpad(string(KIC),9,"0"); kwargs...)
    function loadFITS(KIC::String; usetimecorr=true, prune_kw=5, prune_theshold=5,
                      prune_usephoto=false, usephasmaP=0, cadence="llc",
                      nodetrend=false, detrend_kw=2, splitwidth=0.3, trim=true, detrend_with=median,
                      quarters=1:16, fitsdir="fitsfiles/", usePDC=true, tic_id="", ete6=false) :: DataFrame

	filenames = []

        if (tic_id=="")
            filenames = [fitsdir*fn for fn in readdir(fitsdir)
                         if occursin(KIC, fn) && occursin(cadence, fn)]
        elseif ete6
            tic_id = lpad(tic_id, 16, "0")
            fitsdir = fitsdir * "/" * tic_id[6:8] * "/" * tic_id[9:11] * "/"
            filenames = [fitsdir*fn for fn in readdir(fitsdir)
                        if occursin(tic_id, fn) && occursin("-s_lc", fn)]
        else
            tic_id = lpad(tic_id, 16, "0")
	    for fn in readdir(fitsdir)
		if occursin(tic_id, fn) && occursin("-s_lc", fn)
			filenames = [fitsdir*fn]
		end
	    end
        end
        if filenames == []
            return DataFrame()
        end
        dfs = []

        for fn in filenames
            f = FITS(fn)
            quarter = 0
            if tic_id == ""
                quality = "SAP_QUALITY"
                quarter = read_key(f[1], "QUARTER")[1]
            else
                quality = "QUALITY"
            end

            if (tic_id != "") || (quarter in quarters)
                df = DataFrame()
                fluxerrcol = "PDCSAP_FLUX_ERR"
                if usePDC
                    fluxcol = "PDCSAP_FLUX"
                else
                    fluxcol = "SAP_FLUX"
                    if (tic_id != "")
                        fluxerrcol = "SAP_FLUX_ERR"
                    end
                end
                names = [("TIME", :t), ("TIMECORR", :tcorr), (fluxcol , :F),
                        (fluxerrcol, :sigmaF), (quality, :QUALITY)]
                for (oldname, newname) in names
                     df[newname] = [isnan(x) ? missing : x for x in read(f[2], oldname)]
                end

                if usetimecorr
                    df[:t] .+= df[:tcorr]
                end
                deletecols!(df, :tcorr)
                dropmissing!(df, :F, disallowmissing=false)

                #remove outliers
                df = prune(df, 5, 5, usephoto=false)
                #drop bad points
                df = df[df[:QUALITY] .== 0,:]
                deletecols!(df, :QUALITY)

                ts = df[:t]
                delts = [ts[i+1] - ts[i] for i in 1:(length(ts)-1)]
                boundaries = vcat([0], findall((dt->dt>splitwidth), delts), [length(ts)])
                if boundaries == [0,0]
                    continue
                end
                for i in 1:(length(boundaries)-1)
                    sdf = df[boundaries[i]+1:boundaries[i+1], :]
                    if sdf[end, :t] - sdf[1, :t] < splitwidth
                        continue
                    end
                    push!(dfs, sdf)
                end
            end
            close(f)
        end

        #detrend and normalize each segment
        for i in 1:length(dfs)
            dfs[i][:F] = convert(Vector{Float32}, dfs[i][:F])

            if (tic_id == "")
                interpolate_missing!(dfs[i])
            else
                interpolate_missing!(dfs[i], tess = true)
            end
            if usephasmaP != 0
                dfs[i] = phasma(dfs[i], usephasmaP)
            elseif !nodetrend
                dfs[i] = detrend(dfs[i], detrend_kw, trim, detrend_with)
            end
            #normalized each quarter individually
            m = mean(dfs[i][:F])
            dfs[i][:F] ./= m
            dfs[i][:sigmaF] ./= m
            dfs[i][:F] .-= 1
        end

        if length(dfs) == 0
            return DataFrame()
        elseif length(dfs) == 1
            return dfs[1]
        else
            #combine quarters
            return vcat(dfs...)
        end
    end
end

max_stars_repo_name: ajwheeler/weirddetector

max_issues_repo_name: ajwheeler/weirddetector

max_forks_repo_name: ajwheeler/weirddetector


text:
seats_check = zeros(Bool, 2^10)

seats = []
open("day05/day05.in") do io
    global map_str, passports
    for line in readlines(io)
        push!(seats, line)
    end
end

function seatID(s::String)::Int64
    id = 0
    for chi in 1:10
        id *= 2
        if s[chi] == 'B' || s[chi] == 'R'
            id += 1
        end
    end
    return id
end

for seat_s in seats
    global seats_check
    seats_check[seatID(seat_s)+1] = true
end

for seat in 100:800
    global seats_check
    if !seats_check[seat]
        println(seat-1)
    end
end

ext: jl

max_stars_repo_name: detrin/advent-of-code

max_issues_repo_name: detrin/advent-of-code

max_forks_repo_name: detrin/advent-of-code


text:
using SuccessiveConvexProgrammings
using LinearAlgebra

function my_func_array(x, u, t, k)
    A = [1 2;
         3 4;
         5 6]
    B = Matrix(I, 3, 3)
    return A*x + B*u
end

x = [1, 2]
u = [3, 4, 5]
t = [6, 7]
k = [8, 9]
jacob = get_jacobian(my_func_array, x, u, t, k)


text:
"""
    solve(prob::FODESystem, h, FLMMBDF())

Use [Backward Differentiation Formula](https://en.wikipedia.org/wiki/Backward_differentiation_formula) generated weights fractional linear multi-steps method to solve system of FODE.

### References

```tex
@article{Garrappa2015TrapezoidalMF,
  title={Trapezoidal methods for fractional differential equations: Theoretical and computational aspects},
  author={Roberto Garrappa},
  journal={ArXiv},
  year={2015},
  volume={abs/1912.09878}
}
```
"""
struct FLMMBDF <: FractionalDiffEqAlgorithm end

function solve(prob::FODESystem, h, ::FLMMBDF; reltol=1e-6, abstol=1e-6)
    @unpack f, α, u0, tspan = prob
    t0 = tspan[1]; tfinal = tspan[2]
    fdefun, alphas, y0 = f, α, u0
    alpha = alphas[1]
    itmax = 100

    # generate jacobian of input function
    Jfdefun(t, u) = jacobian_of_fdefun(fdefun, t, u)

    m_alpha = ceil.(Int, alpha)
    m_alpha_factorial = factorial.(collect(0:m_alpha-1))
    # Structure for storing information of the problem

    problem_size = size(y0, 1)

    # Number of points in which to evaluate the solution or the BDFWeights
    r::Int = 16
    N::Int = ceil(Int, (tfinal-t0)/h)
    Nr::Int = ceil(Int, (N+1)/r)*r
    Q::Int = ceil(Int, log2((Nr)/r))-1
    global NNr = 2^(Q+1)*r

    # Preallocation of some variables
    y = zeros(problem_size, N+1)
    fy = zeros(problem_size, N+1)
    zn = zeros(problem_size, NNr+1)

    # Evaluation of convolution and starting BDFWeights of the FLMM
    (omega, w, s) = BDFWeights(alpha, NNr+1)
    halpha = h^alpha

    # Initializing solution and proces of computation
    t = collect(0:N)*h
    y[:, 1] = y0[:, 1]
    temp = zeros(length(y0[:, 1]))
    fdefun(temp, y0[:, 1], nothing, t0)
    #fy[:, 1] = BDFf_vectorfield(t0, y0[:, 1], fdefun)
    fy[:, 1] = temp
    (y, fy) = BDFFirstApproximations(t, y, fy, abstol, itmax, s, halpha, omega, w, problem_size, fdefun, Jfdefun, y0, m_alpha, t0, m_alpha_factorial)
    (y, fy) = BDFTriangolo(s+1, r-1, 0, t, y, fy, zn, N, abstol, itmax, s, w, omega, halpha, problem_size, fdefun, Jfdefun, y0, m_alpha, t0, m_alpha_factorial)

    # Main process of computation by means of the FFT algorithm
    nx0::Int = 0; ny0::Int = 0
    ff = zeros(1, 2^(Q+2), 1)
    ff[1:2] = [0 2]
    for q = 0:Q
        L::Int = 2^q
        (y, fy) = BDFDisegnaBlocchi(L, ff, r, Nr, nx0+L*r, ny0, t, y, fy, zn, N, abstol, itmax, s, w, omega, halpha, problem_size, fdefun, Jfdefun, y0, m_alpha, t0, m_alpha_factorial) ;
        ff[1:4*L] = [ff[1:2*L]; ff[1:2*L-1]; 4*L]
    end
    # Evaluation solution in TFINAL when TFINAL is not in the mesh
    if tfinal < t[N+1]
        c = (tfinal - t[N])/h
        t[N+1] = tfinal
        y[:, N+1] = (1-c)*y[:, N] + c*y[:, N+1]
    end
    t = t[1:N+1]; y = y[:, 1:N+1]
    return t, y
end

function BDFDisegnaBlocchi(L, ff, r, Nr, nx0, ny0, t, y, fy, zn, N , abstol, itmax, s, w, omega, halpha, problem_size, fdefun, Jfdefun, y0, m_alpha, t0, m_alpha_factorial)
    nxi::Int = copy(nx0); nxf::Int = copy(nx0 + L*r - 1)
    nyi::Int = copy(ny0); nyf::Int = copy(ny0 + L*r - 1)
    is::Int = 1
    s_nxi = zeros(N)
    s_nxf = zeros(N)
    s_nyi = zeros(N)
    s_nyf = zeros(N)
    s_nxi[is] = nxi; s_nxf[is] = nxf; s_nyi[is] = nyi; s_nyf[is] = nyf
    i_triangolo = 0;  stop = false
    while ~stop
        stop = (nxi+r-1 == nx0+L*r-1) || (nxi+r-1>=Nr-1)
        zn = BDFQuadrato(nxi, nxf, nyi, nyf, fy, zn, omega, problem_size)
        (y, fy) = BDFTriangolo(nxi, nxi+r-1, nxi, t, y, fy, zn, N, abstol, itmax, s, w, omega, halpha, problem_size, fdefun, Jfdefun, y0, m_alpha, t0, m_alpha_factorial) ;#fy出问题了
        i_triangolo = i_triangolo + 1

        if ~stop
            if nxi+r-1 == nxf   # Il triangolo finisce dove finisce il quadrato -> si scende di livello
                i_Delta = ff[i_triangolo]
                Delta = i_Delta*r
                nxi = s_nxf[is]+1; nxf = s_nxf[is] + Delta
                nyi = s_nxf[is] - Delta +1; nyf = s_nxf[is]
                s_nxi[is] = nxi; s_nxf[is] = nxf; s_nyi[is] = nyi; s_nyf[is] = nyf
            else # Il triangolo finisce prima del quadrato -> si fa un quadrato accanto
                nxi = nxi + r; nxf = nxi + r - 1 ; nyi = nyf + 1 ; nyf = nyf + r
                is = is + 1
                s_nxi[is] = nxi; s_nxf[is] = nxf ; s_nyi[is] = nyi ; s_nyf[is] = nyf
            end
        end

    end
    return y, fy
end

function BDFQuadrato(nxi, nxf, nyi, nyf, fy, zn, omega, problem_size)
    coef_beg::Int = nxi-nyf; coef_end::Int = nxf-nyi+1
    funz_beg::Int = nyi+1; funz_end::Int = nyf+1
    vett_coef = omega[coef_beg+1:coef_end+1]
    vett_funz = [fy[:, funz_beg:funz_end]  zeros(problem_size, funz_end-funz_beg+1)]
    zzn = real(FastConv(vett_coef, vett_funz))
    zn[:, nxi+1:nxf+1] = zn[:, nxi+1:nxf+1] + zzn[:, nxf-nyf:end-1]
    return zn
end

function BDFTriangolo(nxi, nxf, j0, t, y, fy, zn, N, abstol, itmax, s, w, omega, halpha, problem_size, fdefun, Jfdefun, y0, m_alpha, t0, m_alpha_factorial)
    for n = nxi:min(N, nxf)
        n1::Int = n+1
        St = ABMStartingTerm(t[n1], y0, m_alpha, t0, m_alpha_factorial)

        Phi = zeros(problem_size, 1)
        for j = 0:s
            Phi = Phi + w[j+1, n1]*fy[:, j+1]
        end
        for j = j0:n-1
            Phi = Phi + omega[n-j+1]*fy[:, j+1]
        end
        Phi_n = St + halpha*(zn[:, n1] + Phi)

        yn0 = y[:, n]
        temp = zeros(length(yn0))
        fdefun(temp, yn0, nothing, t[n1])
        #fn0 = BDFf_vectorfield(t[n1], yn0, fdefun)
        fn0 = temp
        Jfn0 = Jf_vectorfield(t[n1], yn0, Jfdefun)
        Gn0 = yn0 - halpha*omega[1]*fn0 - Phi_n
        stop = false; it = 0
        while ~stop
            JGn0 = zeros(problem_size, problem_size)+I - halpha*omega[1]*Jfn0
            global yn1 = yn0 - JGn0\Gn0
            global fn1 = zeros(length(yn1)) #BDFf_vectorfield(t[n1], yn1, fdefun)
            fdefun(fn1, yn1, nothing, t[n1])
            Gn1 = yn1 - halpha*omega[1]*fn1 - Phi_n
            it = it + 1

            stop = (norm(yn1-yn0, Inf) < abstol) || (norm(Gn1, Inf)<abstol)
            if it > itmax && ~stop
                @warn "Non Convergence"
                stop = true
            end

            yn0 = yn1; Gn0 = Gn1
            if ~stop
                Jfn0 = Jf_vectorfield(t[n1], yn0, Jfdefun)
            end

        end
        y[:, n1] = yn1
        fy[:, n1] = fn1
    end
    return y, fy
end

function BDFFirstApproximations(t, y, fy, abstol, itmax, s, halpha, omega, w, problem_size, fdefun, Jfdefun, y0, m_alpha, t0, m_alpha_factorial)
    m = problem_size
    Im = zeros(m, m)+I; Ims = zeros(m*s, m*s)+I
    Y0 = zeros(s*m, 1); F0 = copy(Y0); B0 = copy(Y0)
    for j = 1 : s
        Y0[(j-1)*m+1:j*m, 1] = y[:, 1]
        temp = zeros(length(y[:, 1]))
        fdefun(temp, y[:, 1], nothing, t[j+1])
        F0[(j-1)*m+1:j*m, 1] = temp#BDFf_vectorfield(t[j+1], y[:, 1], fdefun)
        St = ABMStartingTerm(t[j+1], y0, m_alpha, t0, m_alpha_factorial)
        B0[(j-1)*m+1:j*m, 1] = St + halpha*(omega[j+1]+w[1, j+1])*fy[:, 1]
    end
    W = zeros(s, s)
    for i = 1:s
        for j = 1:s
            if i >= j
                W[i, j] = omega[i-j+1] + w[j+1, i+1]
            else
                W[i, j] = w[j+1, i+1]
            end
        end
    end
    W = halpha*kron(W, Im)
    G0 = Y0 - B0 - W*F0
    JF = zeros(s*m, s*m)
    for j = 1:s
        JF[(j-1)*m+1:j*m, (j-1)*m+1:j*m] = Jf_vectorfield(t[j+1], y[:, 1], Jfdefun)
    end
    stop = false; it = 0
    F1 = zeros(s*m, 1)
    while ~stop
        JG = Ims - W*JF
        global Y1 = Y0 - JG\G0

        for j = 1 : s
            temp = zeros(length(Y1[(j-1)*m+1:j*m, 1]))
            fdefun(temp, Y1[(j-1)*m+1:j*m, 1], nothing, t[j+1])
            F1[(j-1)*m+1:j*m, 1] = temp#BDFf_vectorfield(t[j+1], Y1[(j-1)*m+1:j*m, 1], fdefun)
        end
        G1 = Y1 - B0 - W*F1

        it = it + 1

        stop = (norm(Y1-Y0, Inf) < abstol) || (norm(G1, Inf) <  abstol)
        if it > itmax && ~stop
            @warn "Non Convergence"
            stop = 1
        end

        Y0 = Y1 ; G0 = G1
        if ~stop
            for j = 1 : s
                JF[(j-1)*m+1:j*m, (j-1)*m+1:j*m] = Jf_vectorfield(t[j+1], Y1[(j-1)*m+1:j*m, 1], Jfdefun)
            end
        end

    end
    for j = 1 : s
        y[:, j+1] = Y1[(j-1)*m+1:j*m, 1]
        fy[:, j+1] = F1[(j-1)*m+1:j*m, 1]
    end
    return y, fy
end

function BDFWeights(alpha, N)
    # BDF-2 with generating function (2/3/(1-4x/3+x^2/3))^alpha
    omega = zeros(1, N+1)
    onethird = 1/3; fourthird = 4/3
    twothird_oneminusalpha = 2/3*(1-alpha)
    fourthird_oneminusalpha = 4/3*(1-alpha)
    omega[1] = 1; omega[2] = fourthird*alpha*omega[1]
    for n = 2:N
        omega[n+1] = (fourthird - fourthird_oneminusalpha/n)*omega[n] + (twothird_oneminusalpha/n - onethird)*omega[n-1]
    end
    omega = omega*((2/3)^(alpha))

    k = floor(1/abs(alpha))
    if abs(k - 1/alpha) < 1.0e-12
        A = collect(0:k)*abs(alpha)
    else
        A = [collect(0:k)*abs(alpha); 1]
    end
    s = length(A) - 1
    # Generation of the matrix and the right hand--side vectors of the system
    nn = collect(0:N)
    V = zeros(s+1, s+1); jj_nu = zeros(s+1, N+1); nn_nu_alpha = copy(jj_nu)
    for i = 0:s
        nu = A[i+1]
        V[i+1, :] = collect(0:s).^nu
        jj_nu[i+1, :] = nn.^nu
        if alpha > 0
            nn_nu_alpha[i+1, :] = gamma(nu+1)/gamma(nu+1+alpha)*nn.^(nu+alpha)
        else
            if i == 0
                nn_nu_alpha[i+1, :] = zeros(1, N+1)
            else
                nn_nu_alpha[i+1, :] = gamma(nu+1)/gamma(nu+1+alpha)*nn.^(nu+alpha)
            end
        end
    end

    temp = FastConv([omega  zeros(size(omega))], [jj_nu zeros(size(jj_nu))])
    temp = real.(temp)
    b = nn_nu_alpha - temp[:, 1:N+1]
    # Solution of the linear system with multiple right-hand side
    w = real.(V\b)
    return omega, w, s
end

Jf_vectorfield(t, y, Jfdefun) = Jfdefun(t, y)

function ABMStartingTerm(t,y0, m_alpha, t0, m_alpha_factorial)
    ys = zeros(size(y0, 1), 1)
    for k = 1:Int64(m_alpha)
        ys = ys + (t-t0)^(k-1)/m_alpha_factorial[k]*y0[:, k]
    end
    return ys
end

function jacobian_of_fdefun(f, t, y)
    ForwardDiff.jacobian(y) do y # Thanks Bagge Carlson for helping me with this
    du = similar(y)
    f(du, y, 0, t)
    du
    end
end


text:
function write_doc(doc, rendered, out_path)
    return write_doc(doc.format, doc, rendered, out_path)
end

function write_doc(::WeaveFormat, doc, rendered, out_path)
    write(out_path, rendered)
    return out_path
end

include("pandoc.jl")
include("latex.jl")


text:
@testset "spherical harmonics" begin
  ɛ = 4*eps(Float64)
  @polyvar x y z

  #@test_throws(ErrorException(BoundsError), ylm(0,1,x,y,z))
  @test_throws DomainError ylm(0,1,x,y,z)
  # l = 0
  #Y_{0,0} = √{\frac{1}{4π}}
  @test isapprox(ylm(0,0,x,y,z),sqrt(1/(4pi))+0*x,atol=ɛ)
  # l = 1
  #Y_{1,-1} = √{\frac{3}{4π}}\frac{y}{r}
  #Y_{1,0} = √{\frac{3}{4π}}\frac{z}{r}
  #Y_{1,1} = √{\frac{3}{4π}}\frac{x}{r}
  @test isapprox(ylm(1,-1,x,y,z),0+sqrt(3/(4pi))*y,atol=ɛ)
  @test isapprox(ylm(1,0,x,y,z),0+sqrt(3/(4pi))*z,atol=ɛ)
  @test isapprox(ylm(1,1,x,y,z),0+sqrt(3/(4pi))*x,atol=ɛ)
  # l = 2
  #Y_{2,-2} = √{\frac{15}{4π}}\frac{xy}{r²}
  #Y_{2,-1} = √{\frac{15}{4π}}\frac{yz}{r²}
  #Y_{2,0} = √{\frac{5}{16π}}\frac{3z² - 1}{r²}
  #Y_{2,1} = √{\frac{15}{4π}}\frac{zx}{r²}
  #Y_{2,2} = √{\frac{15}{16π}}\frac{x²-y²}{r²}
  @test isapprox(ylm(2,-2,x,y,z),0+sqrt(15/(4pi))*x*y,atol=ɛ)
  @test isapprox(ylm(2,-1,x,y,z),0+sqrt(15/(4pi))*y*z,atol=ɛ)
  @test isapprox(ylm(2,0,x,y,z),sqrt(5/(16pi))*(3*z^2 - 1),atol=ɛ)
  @test isapprox(ylm(2,1,x,y,z),0+sqrt(15/(4pi))*z*x,atol=ɛ)
  @test isapprox(ylm(2,2,x,y,z),sqrt(15/(16pi))*(x^2-y^2),atol=ɛ)
  # l = 3
  #Y_{3,-3} = √{\frac{35}{32π}}\frac{(3x²-y²)y}{r³}
  #Y_{3,-2} = √{\frac{105}{4π}}\frac{xyz}{r³}
  #Y_{3,-1} = √{\frac{21}{32π}}\frac{y(5z²-1)}{r³}
  #Y_{3,0} = √{\frac{7}{16π}}\frac{5z³-3z}{r³}
  #Y_{3,1} = √{\frac{21}{32π}}\frac{x(5z²-1)}{r³}
  #Y_{3,2} = √{\frac{105}{16π}}\frac{(x²-y²)z}{r³}
  #Y_{3,3} = √{\frac{35}{32π}}\frac{(x²-3y²)x}{r³}
  @test isapprox(ylm(3,-3,x,y,z),sqrt(35/(32pi))*(3*x^2-y^2)*y,atol=ɛ)
  @test isapprox(ylm(3,-2,x,y,z),0+sqrt(105/(4pi))*x*y*z,atol=ɛ)
  @test isapprox(ylm(3,-1,x,y,z),sqrt(21/(32pi))*y*(5*z^2-1),atol=ɛ)
  @test isapprox(ylm(3,0,x,y,z),sqrt(7/(16pi))*(5*z^3-3*z),atol=ɛ)
  @test isapprox(ylm(3,1,x,y,z),sqrt(21/(32pi))*x*(5*z^2-1),atol=ɛ)
  @test isapprox(ylm(3,2,x,y,z),sqrt(105/(16pi))*(x^2-y^2)*z,atol=ɛ)
  @test isapprox(ylm(3,3,x,y,z),sqrt(35/(32pi))*(x^2-3*y^2)*x,atol=ɛ)
end

@testset "evaluate r^l*ylm" begin
	ε = 32*eps(Float64)
	@polyvar x y z

	@test isapprox(SphericalHarmonics.rlylm(0,0,x,y,z),sqrt(1/(4pi))+0*x,atol=ɛ)

	@test isapprox(SphericalHarmonics.rlylm(1,-1,x,y,z),0+sqrt(3/(4pi))*y,atol=ɛ)
	@test isapprox(SphericalHarmonics.rlylm(1,0,x,y,z),0+sqrt(3/(4pi))*z,atol=ɛ)
	@test isapprox(SphericalHarmonics.rlylm(1,1,x,y,z),0+sqrt(3/(4pi))*x,atol=ɛ)

	@test isapprox(SphericalHarmonics.rlylm(2,-2,x,y,z),0+sqrt(15/(4pi))*x*y,atol=ɛ)
	@test isapprox(SphericalHarmonics.rlylm(2,-1,x,y,z),0+sqrt(15/(4pi))*y*z,atol=ɛ)
	@test isapprox(SphericalHarmonics.rlylm(2,0,x,y,z),sqrt(5/(16pi))*(2*z^2 - x^2 - y^2),atol=ɛ)
	@test isapprox(SphericalHarmonics.rlylm(2,1,x,y,z),0+sqrt(15/(4pi))*z*x,atol=ɛ)
	@test isapprox(SphericalHarmonics.rlylm(2,2,x,y,z),sqrt(15/(16pi))*(x^2-y^2),atol=ɛ)

	@test isapprox(SphericalHarmonics.rlylm(3,-3,x,y,z),sqrt(35/(32pi))*(3*x^2-y^2)*y,atol=ɛ)
	@test isapprox(SphericalHarmonics.rlylm(3,-2,x,y,z),0+sqrt(105/(4pi))*x*y*z,atol=ɛ)
	@test isapprox(SphericalHarmonics.rlylm(3,-1,x,y,z),sqrt(21/(32pi))*y*(4*z^2-x^2-y^2),atol=ɛ)
	@test isapprox(SphericalHarmonics.rlylm(3,0,x,y,z),sqrt(7/(16pi))*(5*z^3-3*z*(x^2+y^2+z^2)),atol=ɛ)
	@test isapprox(SphericalHarmonics.rlylm(3,1,x,y,z),sqrt(21/(32pi))*x*(4*z^2-x^2-y^2),atol=ɛ)
	@test isapprox(SphericalHarmonics.rlylm(3,2,x,y,z),sqrt(105/(16pi))*(x^2-y^2)*z,atol=ɛ)
	@test isapprox(SphericalHarmonics.rlylm(3,3,x,y,z),sqrt(35/(32pi))*(x^2-3*y^2)*x,atol=ɛ)

end

@testset "solid harmonic" begin
	ε = 32*eps(Float64)
	@polyvar x y z

	# test for l = 0,1,2,3

	# l = 0
	@test isapprox(rlm(0,0,x,y,z),1+0*x,atol=ε)

	# l = 1
	@test isapprox(rlm(1,-1,x,y,z),y+0,atol=ε)
	@test isapprox(rlm(1,0,x,y,z),z+0,atol=ε)
	@test isapprox(rlm(1,1,x,y,z),x+0,atol=ε)

	# l = 2
	@test isapprox(rlm(2,-2,x,y,z),sqrt(3)*x*y+0,atol=ε)
	@test isapprox(rlm(2,-1,x,y,z),sqrt(3)*y*z+0,atol=ε)
	@test isapprox(rlm(2,0,x,y,z),1/2*(2*z^2-x^2-y^2),atol=ε)
	@test isapprox(rlm(2,1,x,y,z),sqrt(3)*x*z+0,atol=ε)
	@test isapprox(rlm(2,2,x,y,z),sqrt(3/4)*(x^2-y^2),atol=ε)

	# l = 3
	@test isapprox(rlm(3,-3,x,y,z),sqrt(5/8)*(3*x^2*y-y^3),atol=ε)
	@test isapprox(rlm(3,-2,x,y,z),sqrt(15)*x*y*z+0,atol=ε)
	@test isapprox(rlm(3,-1,x,y,z),sqrt(3/8)*y*(4*z^2-x^2-y^2),atol=ε)
	@test isapprox(rlm(3,0,x,y,z),1/2*(5*z^3-3*z*(x^2+y^2+z^2)),atol=ε)
	@test isapprox(rlm(3,1,x,y,z),sqrt(3/8)*x*(4*z^2-x^2-y^2),atol=ε)
	@test isapprox(rlm(3,2,x,y,z),sqrt(15/4)*(x^2-y^2)*z,atol=ε)
	@test isapprox(rlm(3,3,x,y,z),sqrt(5/8)*(x^3-3*x*y^2),atol=ε)
end


text:
using TangentMps
using Test

@testset "TangentMps.jl" begin
    # Write your own tests here.
end


text:
# ==================
# Dependencies
# ==================
using ReachabilityAnalysis, CarlemanLinearization
using Plots, LaTeXStrings, LinearAlgebra, SparseArrays
using Plots.PlotMeasures
using DifferentialEquations

using CarlemanLinearization: _error_bound_specabs_R

include("../utils.jl")

# ==================
# Model definition
# ==================

# the following are constants of the model
# the domain of x is [-L0/2, L0/2]
L0 = 1.0     # domain length
ν = 0.05     # kinematic viscosity
U0 = 1.0
Re = 20
ν = U0 * L0 / Re
Tnl = L0 / U0
Tmax = Tnl/2

# ================================
# Spatial flowpipe visualization
# ================================

# no-op
ReachabilityAnalysis.convexify(R::ReachabilityAnalysis.AbstractReachSet) = R

function flowpipe_spatial(sol, t, n)
    solt = sol(t) |> convexify

    xdom = range(-L0/2, L0/2, length=n)
    xdomred = xdom[2:end-1]

    X = [set(overapproximate(Projection(solt, i:i), Interval)) for i in 1:dim(sol)]
    R = ReachSet(CartesianProductArray(X), tspan(solt))
    arr = [Interval(xdomred[i]) × CartesianProductArray(X).array[i] for i in 1:length(xdomred)]
    arr = vcat(Interval(xdom[1]) × Interval(0), arr, Interval(xdom[end]) × Interval(0))
    U = UnionSetArray(arr)
    Uch = [ConvexHull(U.array[i], U.array[i+1]) for i in 1:length(U.array)-1] |> UnionSetArray

    return Uch
end

# model for solution with TM
@taylorize function burgers!(du, u, p, t)
    local nf = length(u)
    n = nf + 2
    local Δx = L0/(n-1)
    local c1 = ν/Δx^2
    local c2 = 1/(4Δx)

    du[1] = c1*(-2*u[1]+u[2]) - c2*(u[2]^2)
    du[nf] = c1*(u[nf-1]-2*u[nf]) - c2*(-u[nf-1]^2)

    for i in 2:nf-1
        du[i] = c1*(u[i-1]-2*u[i]+u[i+1]) - c2*(u[i+1]^2 - u[i-1]^2)
    end
end

# ====================
# Carleman method
# ====================

using DynamicPolynomials

# we remove the endpoints
function burgers_carlin(nf)
    n = nf + 2
    Δx = L0/(n-1)
    c1 = ν/Δx^2
    c2 = 1/(4Δx)

    # fix the endpoints
    n = nf # rename to simplify the rest of the code
    dv = c1 * fill(-2.0, n)
    ev = c1 * fill(1.0, n-1)

    F1 = SymTridiagonal(dv, ev)

    F2 = zeros(n, n^2)
    x, = @polyvar x[1:n]
    y = kron(x, x)

    idx = findfirst(y, x[2]^2)
    F2[1, idx] = 1.0

    idx = findfirst(y, x[n]^2)
    F2[n, idx] = -1.0

    for i in 2:n-1
        idx = findfirst(y, x[i-1]^2)
        F2[i, idx] = -1.0

        idx = findfirst(y, x[i+1]^2)
        F2[i, idx] = 1.0
    end
    F2 .*= -c2
    F2 = sparse(F2)

    return F1, F2
end

# ===================
# Solution method
# ===================

function _solve_burgers_carlin(; n, N, width0=0, δ=0.02)

    nf = n-2
    F1, F2 = burgers_carlin(nf)

    A = build_matrix(F1, F2, N)

    xdom = range(-L0/2, L0/2, length=n)
    c0 = Singleton(-U0*sin.(2*π*xdom[2:end-1]/L0))
    if width0 == 0
        X0 = overapproximate(c0, Hyperrectangle);
    else
        H0 = Hyperrectangle(zeros(nf), fill(width0, nf))
        X0 = overapproximate(c0 ⊕ H0, Hyperrectangle)
    end

    dirs = _template(n=nf, N=N)
    # amodel = Forward(inv=true, setops=:lazy, sih=:concrete) # FIXME set options
    alg = LGG09(δ=δ, template=dirs, approx_model=Forward())

    sol_carlin = _solve_CARLIN(X0, F1, F2; alg=alg, N=N, T=Tmax, bloat=false, kron_pow_algorithm="symbolic");
    sol_carlin_spatial = flowpipe_spatial(sol_carlin, Tnl/2, n);
end

# ===================
# TM solution method
# ===================

function _solve_burgers_TM(; n, width0=0)

    nf = n - 2

    xdom = range(-L0/2, L0/2, length=n)
    c0 = Singleton(-U0*sin.(2*π*xdom[2:end-1]/L0))
    if width0 == 0
        X0 = overapproximate(c0, Hyperrectangle)
    else
        H0 = Hyperrectangle(zeros(nf), fill(width0, nf))
        X0 = overapproximate(c0 ⊕ H0, Hyperrectangle)
    end

    prob = IVP(BlackBoxContinuousSystem(burgers!, nf), X0)
    solTM = solve(prob, tspan=(0.0, Tmax), alg=TMJets20(), ensemble=true, trajectories=1);
    solzTM = overapproximate(solTM, Zonotope);

    sol_spatial = flowpipe_spatial(solzTM, Tnl/2, n);
end

function _get_R(; n)
    nf = n-2
    F1, F2 = burgers_carlin(nf);
    xdom = range(-L0/2, L0/2, length=n)
    x0 = -U0*sin.(2*π*xdom[2:end-1]/L0);
    R, Re_λ1 = _error_bound_specabs_R(x0, Matrix(F1), Matrix(F2); check=true)
end

# ================================
# Results
# ================================

n = 10

# Initial point
_solve_burgers_TM(n=n); tt = @elapsed _solve_burgers_TM(n=n)
print(io, "Burgers, TMJets, initial point, -, $(tt)\n")

_solve_burgers_carlin(n=n, N=2, δ=0.01); tt = @elapsed _solve_burgers_carlin(n=n, N=2, δ=0.01)
print(io, "Burgers, Carleman, initial point, N=2, $(tt)\n")

_solve_burgers_carlin(n=n, N=3, δ=0.01); tt = @elapsed _solve_burgers_carlin(n=n, N=3, δ=0.01)
print(io, "Burgers, Carleman, initial point, N=3, $(tt)\n")

# Initial set
width0 = 0.03
_solve_burgers_TM(n=n, width0=width0); tt = @elapsed _solve_burgers_TM(n=n, width0=width0)
print(io, "Burgers, TMJets, initial set, -, $(tt)\n")

_solve_burgers_carlin(n=n, N=2, δ=0.01, width0=width0); tt = @elapsed _solve_burgers_carlin(n=n, N=2, δ=0.01, width0=width0)
print(io, "Burgers, Carleman, initial set, N=2, $(tt)\n")

_solve_burgers_carlin(n=n, N=3, δ=0.01, width0=width0); tt = @elapsed _solve_burgers_carlin(n=n, N=3, δ=0.01, width0=width0)
print(io, "Burgers, Carleman, initial set, N=3, $(tt)\n")

# ==================
# Figures
# ==================

# initial point
function figure_burger_point()

    fig = plot(legend=:topright, xlab=L"x", ylab=L"u(x, 0.5)",
               legendfontsize=25,
               tickfont=font(25, "Times"),
               guidefontsize=25,
               xguidefont=font(35, "Times"),
               yguidefont=font(35, "Times"),
               xtick = xticklatex([-0.50, -0.25, 0.0, 0.25, 0.50], 2),
               ytick = xticklatex([-1.0, -0.5, 0.0, 0.5, 1.0], 2),
               xlims=(-0.5, 0.5),
               ylims=(-1.05, 1.0),
               bottom_margin=10mm, left_margin=6mm, right_margin=6mm, top_margin=3mm, size=(900, 600))

    # model values
    n = 10

    # initial condition (smooth)
    xdom = range(-L0/2, L0/2, length=200)
    c0 = Singleton(-U0*sin.(2*π*xdom/L0))
    plot!(fig, xdom, element(c0), seriestype=:path, c=:black, marker=:none, ls=:dash, lab=L"u(x, 0)", lw=2.0, alpha=1.)

    # initial cond, samples
    xdom = range(-L0/2, L0/2, length=n)
    c0 = Singleton(-U0*sin.(2*π*xdom/L0))
    plot!(fig, xdom, element(c0), seriestype=:path, c=:grey, marker=:square, ls=:dash, lab="", lw=2.0, alpha=1.)

    solN2_carlin = _solve_burgers_carlin(n=n, N=2, δ=0.01)
    solN3_carlin = _solve_burgers_carlin(n=n, N=3, δ=0.01);

    sol_TM = _solve_burgers_TM(n=n);

    plot!(fig, solN3_carlin, c=:orange, lab=L"N = 3", alpha=1., lc=:orange)
    plot!(fig, solN2_carlin, c=:red, lab=L"N = 2", alpha=1., lc=:red)
    plot!(fig, sol_TM, c=:blue, lab=L"\textrm{TM}", alpha=1., lc=:blue)

    return fig
end

fig = figure_burger_point()
savefig(fig, joinpath(TARGET_FOLDER, "figure_burger_point.pdf"))

# ==================
# initial set
# ==================

function figure_burger_set()

    fig = plot(legend=:topright, xlab=L"x", ylab=L"u(x, 0.5)",
               legendfontsize=25,
               tickfont=font(25, "Times"),
               guidefontsize=25,
               xguidefont=font(35, "Times"),
               yguidefont=font(35, "Times"),
               xtick = xticklatex([-0.50, -0.25, 0.0, 0.25, 0.50], 2),
               ytick = xticklatex([-1.0, -0.5, 0.0, 0.5, 1.0], 2),
               xlims=(-0.5, 0.5),
               ylims=(-1.05, 1.0),
               bottom_margin=10mm, left_margin=6mm, right_margin=6mm, top_margin=3mm, size=(900, 600))

    # model values
    n = 10
    width0 = 0.03

    # initial condition (smooth)
    xdom = range(-L0/2, L0/2, length=200)
    c0 = Singleton(-U0*sin.(2*π*xdom/L0))
    plot!(fig, xdom, element(c0), seriestype=:path, c=:black, marker=:none, ls=:dash, lab=L"u(x, 0)", lw=2.0, alpha=1.)

    # initial cond, samples
    xdom = range(-L0/2, L0/2, length=n)
    c0 = Singleton(-U0*sin.(2*π*xdom/L0))
    plot!(fig, xdom, element(c0), seriestype=:path, c=:grey, marker=:square, ls=:dash, lab="", lw=2.0, alpha=1.)

    solN2_carlin = _solve_burgers_carlin(n=n, N=2, width0=width0, δ=0.01)
    solN3_carlin = _solve_burgers_carlin(n=n, N=3, width0=width0, δ=0.01);

    sol_TM = _solve_burgers_TM(n=n, width0=width0);

    plot!(fig, solN3_carlin, c=:orange, lab=L"N = 3", alpha=1., lc=:orange)
    plot!(fig, solN2_carlin, c=:red, lab=L"N = 2", alpha=1., lc=:red)
    plot!(fig, sol_TM, c=:blue, lab=L"\textrm{TM}", alpha=1., lc=:blue)

    return fig
end

fig = figure_burger_set()
savefig(fig, joinpath(TARGET_FOLDER, "figure_burger_set.pdf"))


text:
"""
    BlockIndex{N}

A `BlockIndex` is an index which stores a global index in two parts: the block
and the offset index into the block.
"""
immutable BlockIndex{N}
    I::NTuple{N, Int}
    α::NTuple{N, Int}
end

@inline function _find_block(block_sizes::BlockSizes, dim::Int, i::Int)
    accum = 0
    block = 0
    @inbounds while accum < i
        block += 1
        accum += block_sizes[dim, block]
    end
    @inbounds accum -= block_sizes[dim, block]
    return block, i - accum
end

"""
    global2blockindex{N}(block_sizes::BlockSizes{N}, inds...) -> BlockIndex{N}

Converts from global indices `inds` to a `BlockIndex`.
"""
@generated function global2blockindex{N}(block_sizes::BlockSizes{N}, i::NTuple{N, Int})
    # TODO: Try get rid of @generated
    block_index_ex = Expr(:tuple, [:(_find_block(block_sizes, $k, i[$k])) for k = 1:N]...)
    I_ex = Expr(:tuple, [:(block_index[$k][1]) for k = 1:N]...)
    α_ex = Expr(:tuple, [:(block_index[$k][2]) for k = 1:N]...)
    return quote
        @inbounds block_index = $block_index_ex
        @inbounds I = $I_ex
        @inbounds α = $α_ex
        return BlockIndex(I, α)
    end
end

"""
    blockindex2global{N}(block_sizes::BlockSizes{N}, block_index::BlockIndex{N}) -> inds

Converts from a block index to a tuple containing the global indices
"""
@generated function blockindex2global{N}(block_sizes::BlockSizes{N}, block_index::BlockIndex{N})
    # TODO: Try get rid of @generated
    ex = Expr(:tuple, [:(_sumiter(block_sizes[$k], block_index.I[$k]-1) + block_index.α[$k]) for k = 1:N]...)
    return quote
        @inbounds v = $ex
        return v
    end
end


text:
__precompile__()
module EventBasedIntegrator
using Statistics, StatsBase, SpecialFunctions, Compat, Distributed, LinearAlgebra

# Macro for common field sharing
macro def(name, definition)
    return quote
        macro $(esc(name))()
            esc($(Expr(:quote, definition)))
        end
    end
end

struct LIFParams{Ttau,Treset,Tthresh,Tdelay,Tref,Tstrace}
    tau::Ttau
    reset::Treset
    threshold::Tthresh
    delay::Tdelay
    refractory_period::Tref
    spiketrace_timeconstant::Tstrace
end
export LIFParams
"""
    LIFParams(; tau = 10., reset = 0., threshold = 20., delay = .5,
                 refractory_period = 2., spiketrace_timeconstant = 20.)

Each keyword argument accepts either a float or an array of floats with
individual values for each neuron (except `spiketrace_timeconstant` only float).
Choosing individual values for `tau` may increase simulation time a lot.
"""
LIFParams(; tau = 25., reset = 0., threshold = 20., delay = 1e-1,
            refractory_period = 0., spiketrace_timeconstant = 20.) =
LIFParams(tau, reset, threshold, delay, refractory_period,
           spiketrace_timeconstant)

@enum Event SpikeInsert RefractEnd Spike NoEvent

abstract type AbstractNet{Tplasticity} end
mutable struct Net{Tv,Tt,Tw,Tstrace,Tparams,Trecorder,Tplasticity,Tnoise} <: AbstractNet{Tplasticity}
    v::Tv
    n_of_neurons::Int64
    t::Tt
    posts::Array{Array{Tuple{Int64,Int64}, 1}, 1}
    pres::Array{Array{Tuple{Int64,Int64}, 1}, 1}
    weights::Tw
    nextevents::Array{Tuple{Event, Int64, Tt}, 1}
    spikecandidates::Array{Int64, 1}
    recorder::Trecorder
    spiketrace::Tstrace
    parameters::Tparams
    plasticityrule::Tplasticity
    noise::Tnoise
end
export Net
function Net(n_of_neurons, pres, posts, weights, recorder, parameters,
             plasticityrule, noise)
    Net(Float64.(rand(n_of_neurons) .* (parameters.threshold .- parameters.reset) .+
                 parameters.reset),
        n_of_neurons,
        0.,
        pres,
        posts,
        Float64.(weights),
        Tuple{Event, Int64, Float64}[],
        Int64[],
        recorder,
        zeros(n_of_neurons),
        parameters,
        plasticityrule,
        noise)
end

include("escapenoise.jl")
include("eulerforward.jl")
include("recorders.jl")
include("plasticity.jl")
include("ratemodel.jl")
include("helpers.jl")

"""
    function BalancedNet(n_of_neurons, J, g, parameters = LIFParams();
                         f = 3/4,
                         recorder = AllSpikesRecorder(),
                         plasticityrule = NoPlasticity())
"""
function BalancedNet(n_of_neurons, J, g; parameters = LIFParams(),
                     f = 3/4, recorder = AllSpikesRecorder(),
                     plasticityrule = NoPlasticity(),
                     noise = NoNoise())
    posts = [Tuple{Int64,Int64}[] for i in 1:n_of_neurons]
    pres = [Tuple{Int64,Int64}[] for i in 1:n_of_neurons]
    n_of_excitatory = floor(Int64,f*n_of_neurons)
    n_of_inhibitory = floor(Int64, (1-f)*n_of_neurons)
    weights = Float64[J; -g*J]
    for post in 1:n_of_neurons
        for pre in sample(1:n_of_excitatory, div(n_of_excitatory, 10),
                            replace = false)
            push!(posts[pre], (post, 1))
            push!(pres[post], (pre, 1))
        end
        for pre in sample(n_of_excitatory+1:n_of_neurons, div(n_of_inhibitory, 10),
                            replace = false)
            push!(posts[pre], (post, 2))
            push!(pres[post], (pre, 2))
        end
    end
    Net(n_of_neurons, posts, pres, weights, recorder, parameters, plasticityrule, noise)
end
export BalancedNet
"""
    FeedForwardNet(layersizes, parameters = LIFParams();
                   sampleinitw =  (lpre, i, lpost, j) -> .1 * rand(),
                   recorder = AllSpikesRecorder(),
                   plasticityrule = NoPlasticity())
"""
function FeedForwardNet(layersizes, parameters = LIFParams();
                        sampleinitw =  (lpre, i, lpost, j) -> randn() .* 20 ./ sqrt(layersizes[lpre]),
                        recorder = AllSpikesRecorder(),
                        plasticityrule = NoPlasticity(),
                        noise = NoNoise())
    ns = [0; cumsum(layersizes)]
    posts = [Tuple{Int64,Int64}[] for i in 1:ns[end]]
    pres = [Tuple{Int64,Int64}[] for i in 1:ns[end]]
    weights = Float64[]
    weightsindex = 1
    for i in 1:length(layersizes) - 1
        for pre in ns[i] + 1:ns[i+1]
            for post in ns[i+1] + 1:ns[i+2]
                if (w = sampleinitw(i, pre - ns[i], i+1, post - ns[i+1])) != 0
                    push!(posts[pre], (post, weightsindex))
                    push!(pres[post], (pre, weightsindex))
                    push!(weights, w)
                    weightsindex += 1
                end
            end
        end
    end
    Net(ns[end], posts, pres, weights, recorder, parameters, plasticityrule, noise)
end
export FeedForwardNet
"""
    deepcopynet(net;
                recorder = net.recorder,
                plasticityrule = net.plasticityrule)

Copy a net and optionally replace the `recorder` or the `plasticityrule`.
"""
function deepcopynet(net;
                     recorder = deepcopy(net.recorder),
                     plasticityrule = deepcopy(net.plasticityrule),
                     parameters = deepcopy(net.parameters))
    Net([field == :recorder ? recorder :
         field == :plasticityrule ? plasticityrule :
         field == :parameters ? parameters :
         deepcopy(getfield(net, field)) for field in fieldnames(typeof(net))]...)
end
export deepcopynet

let basetypes = (Number, AbstractArray)
    for vargt in basetypes; for inargt in basetypes; for targt in basetypes;
        v = vargt == Number ? :v : :(v[i])
        input = inargt == Number ? :input : :(input[i])
        threshold = targt == Number ? :threshold : :(threshold[i])
        eval(:(@inline arglog(v::$vargt, input::$inargt, threshold::$targt, i::Int64) =
               arglog($v, $input, $threshold)))
    end; end; end
end
@inline function arglog(v, input, threshold)
    threshold >= input && return typeof(v)(Inf)
    (v - input)/(threshold - input)
end

for func in (:gettau, :getrefractory_period, :getdelay, :getreset)
    eval(:(@inline $func(x::Number, i) = x))
    eval(:(@inline $func(x::AbstractArray, i) = x[i]))
end

@inline function nextspiketime(arglog, parameters, i)
    gettau(parameters.tau, i) * log(arglog)
end

function insert_event!(nextevents, event)
    if length(nextevents) == 0 || nextevents[end][3] <= event[3]
        push!(nextevents, event)
        return
    end
    i = length(nextevents)
    while i > 0 && nextevents[i][3] > event[3]
        i -= 1
    end
    insert!(nextevents, i+1, event)
end

@inline function addspikepostprocess!(net, t, i)
    recordspike!(net.recorder, net, (i, t))
    if net.parameters.spiketrace_timeconstant > 0
        net.spiketrace[i] += 1/net.parameters.spiketrace_timeconstant
    end
    net.v[i] = -typeof(net.v[i])(Inf)
end
function addspike!(net, t, i)
    insert_event!(net.nextevents, (RefractEnd, i, t +
                    getrefractory_period(net.parameters.refractory_period, i)))
    (typeof(net.noise) == RectLinEscapeNoise) &&
    (net.noise.endofrefrperiod[i] = t + getrefractory_period(net.parameters.refractory_period, i))
    insert_event!(net.nextevents, (SpikeInsert, i, t +
                                   getdelay(net.parameters.delay, i)))
    addspikepostprocess!(net, t, i)
end

@inline addinput!(gamma, input::Number, v) = v .+= (1 .- gamma) .* input
@inline addinput!(gamma::Number, input::AbstractArray, v) =
    axpy!((1 - gamma), input, v)
@inline function addinput!(gamma::AbstractArray, input::AbstractArray, v)
    v .+= (1 .- gamma) .* input
end

function updatevars!(net, eventtime, input)
    gamma = exp.(-(eventtime - net.t)./net.parameters.tau)
    net.v .*= gamma
    addinput!(gamma, input, net.v)
    if net.parameters.spiketrace_timeconstant > 0
        rmul!(net.spiketrace,
              exp(-(eventtime - net.t)/net.parameters.spiketrace_timeconstant))
    end
    net.t = eventtime
end

@inline function getinext(net, input::Number, threshold::Number)
    if typeof(net.noise) == NoNoise
      v, i = findmax(net.v)
      arglog(v, input, threshold), i
    elseif typeof(net.noise) == RectLinEscapeNoise
      inext = 0; tmin = typeof(net.v[1])(Inf)
      @inbounds for i in 1:net.n_of_neurons
          if (t = getspiketimesample(net.v[i], input, net.parameters.threshold, net.t,
             net.parameters.tau, net.noise.β, net.parameters.reset, (net.v[i] == -Inf), net.noise.endofrefrperiod[i])) < tmin
              tmin = t; inext = i
          end
      end
      return tmin, inext
    end
end
@inline function getvmin(v, input, threshold, range)
    inext = 0; vmin = typeof(v[1])(Inf)
    @inbounds for i in range
        if (vi = arglog(v[i], input, threshold, i)) < vmin
            vmin = vi; inext = i
        end
    end
    return vmin, inext
end
@inline function getvminparallel(v, th, inp)
    nthreads = Threads.nthreads()
    vmin = Inf*ones(nthreads); imin = zeros(Int64, nthreads)
    N = div(length(v), nthreads)
    Threads.@threads for i in 1:nthreads
        if i == nthreads
            range = (i-1)*N + 1:length(v)
        else
            range = (i-1)*N + 1:i * N
        end
        vmin[i], imin[i] = getvmin(v, th, inp, range)
    end
    i0 = findmin(vmin)[2]
    return vmin[i0], imin[i0]
end
@inline function getinext(net, input, threshold)
    if typeof(net.noise) == NoNoise
        return getvminparallel(net.v, input, threshold)
    elseif typeof(net.noise) == RectLinEscapeNoise
      inext = 0; tmin = typeof(net.v[1])(Inf)
      @inbounds for i in 1:net.n_of_neurons
          if (t = getspiketimesample(net.v[i], input[i], net.parameters.threshold[i], net.t,
            net.parameters.tau, net.noise.β, net.parameters.reset, (net.v[i] == -Inf), net.noise.endofrefrperiod[i])) < tmin
              tmin = t; inext = i
          end
      end
      return tmin, inext
    end
end

@inline abovethreshold(v, th::Number, i) = v[i] >= th
@inline abovethreshold(v, th::AbstractArray, i) = v[i] >= th[i]

function broadcastspike!(net, i)
    @inbounds @simd for post in net.posts[i]
        net.v[post[1]] += net.weights[post[2]] # * 1000./net.parameters.tau to match with LIF rate model
        if abovethreshold(net.v, net.parameters.threshold, post[1])
            push!(net.spikecandidates, post[1])
        end
    end
end

nexteventatsamemoment(net) = false
nexteventatsamemoment(net::Net) = (length(net.nextevents) > 0 && net.t == net.nextevents[1][3])
function processspikecandidates!(net)
    if nexteventatsamemoment(net) return end
    @inbounds for i in net.spikecandidates
        if abovethreshold(net.v, net.parameters.threshold, i)
            addspike!(net, net.t, i)
        end
    end
    empty!(net.spikecandidates)
end

# here input = external_input + u_reset
"""
    function integratenet!(net, input, tmax)
"""
function integratenet!(net, input, tmax)
    net.t >= tmax && return
    while updatenet!(net, input, tmax) end
end
export integratenet!

function nextevent(net, input, tmax)
    if length(net.spikecandidates) > 0
        processspikecandidates!(net)
    end
    if length(net.nextevents) > 0 && net.t == net.nextevents[1][3]
        return net.nextevents[1]
    end
    if typeof(net.noise) == NoNoise
      arglog, inext = getinext(net, input, net.parameters.threshold)
    elseif typeof(net.noise) == RectLinEscapeNoise
      tnext, inext = getinext(net, input, net.parameters.threshold)
    end
    if inext == 0
        if length(net.nextevents) == 0
            return (NoEvent, 0, tmax)
        else
            return net.nextevents[1]
        end
    else
        if typeof(net.noise) == NoNoise
          t = net.t + nextspiketime(arglog, net.parameters, inext)
        elseif typeof(net.noise) == RectLinEscapeNoise
          t = tnext
        end
        if length(net.nextevents) == 0 || t < net.nextevents[1][3]
            if t > tmax
                return (NoEvent, 0, tmax)
            else
                return (Spike, inext, t)
            end
        else
            return net.nextevents[1]
        end
    end
end

updatenet!(net, input) = updatenet!(net, input, Inf)
function updatenet!(net, input, tmax)
    updatenetevent!(net, nextevent(net, input, tmax), input, tmax)
end
function updatenetevent!(net, event, input, tmax, doshift = true)
    if event[1] == NoEvent && event[3] == tmax || event[3] > tmax
        updatevars!(net, tmax, input)
        return false
    elseif net.t < event[3]
        updatevars!(net, event[3], input)
    end
    if event[1] == NoEvent
        return true
    elseif event[1] == Spike
        addspike!(net, event[3], event[2])
        return true
    else
        if event[1] == SpikeInsert
            broadcastspike!(net, event[2])
            updateweights!(net, event[2])
        else
            net.v[event[2]] = getreset(net.parameters.reset, event[2])
        end
        if doshift; popfirst!(net.nextevents); end
        return true
    end
end

end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    pairwise(γ, domain)

Evaluate variogram `γ` between all elements in the `domain`.
"""
function pairwise(γ::Variogram, domain)
  Dim = embeddim(domain)
  T = coordtype(domain)
  n = nelements(domain)
  x = rand(Point{Dim,T})
  R = result_type(γ, x, x)
  Γ = Matrix{R}(undef, n, n)
  pairwise!(Γ, γ, domain)
end

function pairwise!(Γ, γ::Variogram, domain)
  @inbounds for j=1:nelements(domain)
    xj = centroid(domain, j)
    for i=j+1:nelements(domain)
      xi = centroid(domain, i)
      Γ[i,j] = γ(xi, xj)
    end
    Γ[j,j] = γ(xj, xj)
    for i=1:j-1
      Γ[i,j] = Γ[j,i] # leverage the symmetry
    end
  end
  Γ
end

"""
    pairwise(γ, domain₁, domain₂)

Evaluate variogram `γ` between all elements of `domain₁` and `domain₂`.
"""
function pairwise(γ::Variogram, domain₁, domain₂)
  Dim = embeddim(domain₁)
  T = coordtype(domain₁)
  m = nelements(domain₁)
  n = nelements(domain₂)
  x = rand(Point{Dim,T})
  R = result_type(γ, x, x)
  Γ = Array{R}(undef, m, n)
  pairwise!(Γ, γ, domain₁, domain₂)
end

function pairwise!(Γ, γ::Variogram, domain₁, domain₂)
  @inbounds for j=1:nelements(domain₂)
    xj = centroid(domain₂, j)
    for i=1:nelements(domain₁)
      xi = centroid(domain₁, i)
      Γ[i,j] = γ(xi, xj)
    end
  end
  Γ
end


text:
"""
Abstract type for storing Neuroimaing data.

All other neuroimaging types inherit from this type.
All neuroimaing types support the following functions:

* `samplingrate()`
* `channelnames()`
* `remove_channel!()`
* `keep_channel!()`
* `trim_channel()`
* `highpass_filter()`
* `lowpass_filter()`
* `rereference()`

# Examples
```julia
data = # load your neuroimaging data
samplingrate(data)  # Returns the sampling rate
channelnames(data)  # Returns the channel names
```
"""
abstract type NeuroimagingMeasurement end


text:
using DelimitedFiles
using CUDA, CuYao
device!(parse(Int, ARGS[1]))

using Random
using TropicalTensors

function run(::Type{T}, n::Int; nrepeat, usecuda) where T
    saveto = joinpath(@__DIR__, "potts_n$(n)_elsl.dat")
    elsl = zeros(T, 3, nrepeat)
    println("The number of repetition is $nrepeat.")
    for i=1:nrepeat
        lt = SquareLattice(n, n)
        t = @elapsed res = solve_potts(CountingTropical{T}, Val(3), lt, TropicalTensors.potts_randpm_J(lt); usecuda=usecuda)
        println("Δt = $t, maximum energy = $(res.n), degeneracy = $(res.c)")
        elsl[1,i] = res.n/n
        elsl[2,i] = log(res.c)/n
        elsl[3,i] = t
    end
    writedlm(saveto, elsl)
    return elsl
end

L = parse(Int, ARGS[2])
run(Float32, L; nrepeat=100, usecuda=true)


text:
# TODO: FIXME for general case
# wrong for transposed matrices, and certain views/SubArrays.
unitstride(op::Operation, s) = first(op.ref.ref) === s

function cost(op::Operation, unrolled::Symbol, Wshift::Int, size_T::Int = op.elementbytes)
    isconstant(op) && return 0.0, 0, 1
    # Wshift == dependson(op, unrolled) ? Wshift : 0
    # c = first(cost(instruction(op), Wshift, size_T))::Int
    instr = instruction(op)
    opisunrolled = dependson(op, unrolled)
    srt, sl, srp = opisunrolled ? vector_cost(instr, Wshift, size_T) : scalar_cost(instr)
    if accesses_memory(op)
        # either vbroadcast/reductionstore, vmov(a/u)pd, or gather/scatter
        # @show instr, unrolled, loopdependencies(op), unitstride(op, unrolled)
        if opisunrolled
            if !unitstride(op, unrolled)# || !isdense(op) # need gather/scatter
                r = (1 << Wshift)
                srt *= r
                sl *= r
            # else # vmov(a/u)pd
            end
        elseif instr === :setindex! # broadcast or reductionstore; if store we want to penalize reduction
            srt *= 2
            sl *= 2
        end
    end
    srt, sl, srp
end

    # Base._return_type()

function biggest_type_size(ls::LoopSet)
    maximum(elsize, operations(ls))
end
function VectorizationBase.pick_vector_width(ls::LoopSet, u::Symbol)
    VectorizationBase.pick_vector_width(length(ls, u), biggest_type_size(ls))
end
function VectorizationBase.pick_vector_width_shift(ls::LoopSet, u::Symbol)
    VectorizationBase.pick_vector_width_shift(length(ls, u), biggest_type_size(ls))
end
function hasintersection(a, b)
    for aᵢ ∈ a, bᵢ ∈ b
        aᵢ === bᵢ && return true
    end
    false
end

# evaluates cost of evaluating loop in given order
# heuristically, could simplify analysis by just unrolling outer loop?
function evaluate_cost_unroll(
    ls::LoopSet, order::Vector{Symbol}, max_cost = typemax(Float64), vectorized::Symbol = first(order)
)
    # included_vars = Set{UInt}()
    included_vars = fill(false, length(operations(ls)))
    nested_loop_syms = Symbol[]#Set{Symbol}()
    total_cost = 0.0
    iter = 1.0
    # Need to check if fusion is possible
    size_T = biggest_type_size(ls)
    W, Wshift = VectorizationBase.pick_vector_width_shift(length(ls, vectorized), size_T)::Tuple{Int,Int}
    for itersym ∈ order
        # Add to set of defined symbles
        push!(nested_loop_syms, itersym)
        liter = Float64(length(ls, itersym))
        if itersym === vectorized
            liter /= W
        end
        iter *= liter
        # check which vars we can define at this level of loop nest
        for (id,op) ∈ enumerate(operations(ls))
            # won't define if already defined...
            # id = identifier(op)
            included_vars[id] && continue
            # it must also be a subset of defined symbols
            loopdependencies(op) ⊆ nested_loop_syms || continue
            # hasintersection(reduceddependencies(op), nested_loop_syms) && return Inf
            rd = reduceddependencies(op)
            hasintersection(rd, nested_loop_syms[1:end-length(rd)]) && return Inf
            included_vars[id] = true

            total_cost += iter * first(cost(op, vectorized, Wshift, size_T))
            total_cost > max_cost && return total_cost # abort if more expensive; we only want to know the cheapest
        end
    end
    total_cost
end

# only covers vectorized ops; everything else considered lifted?
function depchain_cost!(
    skip::Vector{Bool}, op::Operation, vectorized::Symbol, Wshift::Int, size_T::Int, rt::Float64 = 0.0, sl::Int = 0
)
    skip[identifier(op)] = true
    # depth first search
    for opp ∈ parents(op)
        skip[identifier(opp)] && continue
        rt, sl = depchain_cost!(skip, opp, vectorized, Wshift, size_T, rt, sl)
    end
    # Basically assuming memory and compute don't conflict, but everything else does
    # Ie, ignoring the fact that integer and floating point operations likely don't either
    if iscompute(op)
        rtᵢ, slᵢ = cost(op, vectorized, Wshift, size_T)
        rt += rtᵢ; sl += slᵢ
    end
    rt, sl
end
function parentsnotreduction(op::Operation)
    for opp ∈ parents(op)
        isreduction(opp) && return false
    end
    return true
end
function unroll_no_reductions(ls, order, vectorized, Wshift, size_T)
    innermost = last(order)
    compute_rt = 0.0
    load_rt = 0.0
    # latency not a concern, because no depchains
    for op ∈ operations(ls)
        dependson(op, innermost) || continue
        if iscompute(op)
            compute_rt += first(cost(op, vectorized, Wshift, size_T))
        elseif isload(op)
            load_rt += first(cost(op, vectorized, Wshift, size_T))
        end
    end
    # heuristic guess
    # @show compute_rt, load_rt
    min(4, round(Int, (compute_rt + load_rt + 1) / compute_rt))
end
function determine_unroll_factor(
    ls::LoopSet, order::Vector{Symbol}, unrolled::Symbol, vectorized::Symbol = first(order)
)
    size_T = biggest_type_size(ls)
    W, Wshift = VectorizationBase.pick_vector_width_shift(length(ls, vectorized), size_T)::Tuple{Int,Int}

    # The strategy is to use an unroll factor of 1, unless there appears to be loop carried dependencies (ie, num_reductions > 0)
    # The assumption here is that unrolling provides no real benefit, unless it is needed to enable OOO execution by breaking up these dependency chains
    num_reductions = 0#sum(isreduction, operations(ls))
    for op ∈ operations(ls)
        if isreduction(op) & iscompute(op) && parentsnotreduction(op)
            num_reductions += 1
        end
    end
    if iszero(num_reductions)
        # if only 1 loop, no need to unroll
        # if more than 1 loop, there is some cost. Picking 2 here as a heuristic.
        return length(order) == 1 ? 1 : unroll_no_reductions(ls, order, vectorized, Wshift, size_T)
    end
    # So if num_reductions > 0, we set the unroll factor to be high enough so that the CPU can be kept busy
    # if there are, U = max(1, round(Int, max(latency) * throughput / num_reductions)) = max(1, round(Int, latency / (recip_throughput * num_reductions)))
    # We also make sure register pressure is not too high.
    latency = 0
    compute_recip_throughput = 0.0
    visited_nodes = fill(false, length(operations(ls)))
    load_recip_throughput = 0.0
    store_recip_throughput = 0.0
    for op ∈ operations(ls)
        dependson(op, unrolled) || continue
        if isreduction(op)
            rt, sl = depchain_cost!(visited_nodes, op, vectorized, Wshift, size_T)
            latency = max(sl, latency)
            compute_recip_throughput += rt
        elseif isload(op)
            load_recip_throughput += first(cost(op, vectorized, Wshift, size_T))
        elseif isstore(op)
            store_recip_throughput += first(cost(op, vectorized, Wshift, size_T))
        end
    end
    recip_throughput = max(
        compute_recip_throughput,
        load_recip_throughput,
        store_recip_throughput
    )
    max(1, round(Int, latency / (recip_throughput * num_reductions) ) )
end

function tile_cost(X, U, T)
    X[1] + X[4] + X[2] / T + X[3] / U
end
function solve_tilesize(X, R)
    @inbounds any(iszero, (R[1],R[2],R[3])) && return -1,-1,Inf #solve_smalltilesize(X, R, Umax, Tmax)
    # @inbounds any(iszero, (R[1],R[2],R[3])) && return -1,-1,Inf #solve_smalltilesize(X, R, Umax, Tmax)
    # We use a lagrange multiplier to find floating point values for U and T
    # first solving for U via quadratic formula
    # X is vector of costs, and R is of register pressures
    RR = REGISTER_COUNT - R[3] - R[4]
    a = (R[1])^2*X[2] - (R[2])^2*R[1]*X[3]/RR
    b = 2*R[1]*R[2]*X[3]
    c = -RR*R[1]*X[3]
    Ufloat = (sqrt(b^2 - 4a*c) - b) / (2a)
    Tfloat = (RR - Ufloat*R[2])/(Ufloat*R[1])
    # @show Ufloat, Tfloat
    (isfinite(Tfloat) && isfinite(Ufloat)) || return -1,-1,Inf
    Ulow = max(1, floor(Int, Ufloat)) # must be at least 1
    Tlow = max(1, floor(Int, Tfloat)) # must be at least 1
    Uhigh = Ulow + 1 #ceil(Int, Ufloat)
    Thigh = Tlow + 1 #ceil(Int, Tfloat)

    RR = REGISTER_COUNT - R[3] - R[4]
    U, T = Ulow, Tlow
    tcost = tile_cost(X, Ulow, Tlow)
    # @show Ulow*Thigh*R[1] + Ulow*R[2]
    if RR ≥ Ulow*Thigh*R[1] + Ulow*R[2]
        tcost_temp = tile_cost(X, Ulow, Thigh)
        # @show tcost_temp, tcost
        if tcost_temp < tcost
            tcost = tcost_temp
            U, T = Ulow, Thigh
        end
    end
    # The RR + 1 is a hack to get it to favor Uhigh in more scenarios
    Tl = Tlow
    while RR < Uhigh*Tl*R[1] + Uhigh*R[2]
        Tl -= 1
    end
    tcost_temp = tile_cost(X, Uhigh, Tl)
    if tcost_temp < tcost
        tcost = tcost_temp
        U, T = Uhigh, Tl
    end
    if RR > Uhigh*Thigh*R[1] + Uhigh*R[2]
        throw("Something went wrong when solving for Tfloat and Ufloat.")
    end
    U, T, tcost
end
function solve_tilesize_constU(X, R, U)
    floor(Int, (REGISTER_COUNT - R[3] - R[4] - U*R[2]) / (U * R[1]))
end
function solve_tilesize_constT(X, R, T)
    floor(Int, (REGISTER_COUNT - R[3] - R[4]) / (T * R[1] + R[2]))
end
function solve_tilesize_constT(ls, T)
    R = @view ls.reg_pres[:,1]
    floor(Int, (REGISTER_COUNT - R[3] - R[4]) / (T * R[1] + R[2]))
end
# Tiling here is about alleviating register pressure for the UxT
function solve_tilesize(X, R, Umax, Tmax)
    first(R) == 0 && return -1,-1,Inf #solve_smalltilesize(X, R, Umax, Tmax)
    U, T, cost = solve_tilesize(X, R)
    # T -= T & 1
    # U = min(U, T)
    U_too_large = U > Umax
    T_too_large = T > Tmax
    if U_too_large
        if T_too_large
            U = Umax
            T = Tmax
        else # U too large, resolve T
            U = Umax
            T = solve_tilesize_constU(X, R, U)
        end
    elseif T_too_large
        T = Tmax
        U = solve_tilesize_constT(X, R, T)
    end
    U, T, cost
end
function solve_tilesize(
    ls::LoopSet, unrolled::Symbol, tiled::Symbol,
    cost_vec::AbstractVector{Float64} = @view(ls.cost_vec[:,1]),
    reg_pressure::AbstractVector{Int} = @view(ls.reg_pres[:,1])
)
    maxT = 4#8
    maxU = 4#8
    if isstaticloop(ls, tiled)
        maxT = min(2maxT, looprangehint(ls, tiled))
    end
    if isstaticloop(ls, unrolled)
        maxU = min(2maxU, looprangehint(ls, unrolled))
    end
    solve_tilesize(cost_vec, reg_pressure, maxU, maxT)
end

function set_upstream_family!(adal::Vector{T}, op::Operation, val::T) where {T}
    adal[identifier(op)] = val
    for opp ∈ parents(op)
        set_upstream_family!(adal, opp, val)
    end
end

# Just tile outer two loops?
# But optimal order within tile must still be determined
# as well as size of the tiles.
function evaluate_cost_tile(
    ls::LoopSet, order::Vector{Symbol}, vectorized::Symbol
)
    N = length(order)
    @assert N ≥ 2 "Cannot tile merely $N loops!"
    tiled = order[1]
    unrolled = order[2]
    ops = operations(ls)
    nops = length(ops)
    included_vars = fill(false, nops)
    unrolledtiled = fill(false, 2, nops)
    descendentsininnerloop = fill(false, nops)
    innerloop = last(order)
    iters = fill(-99.9, nops)
    nested_loop_syms = Symbol[]# Set{Symbol}()
    iter = 1.0
    # Need to check if fusion is possible
    size_T = biggest_type_size(ls)
    W, Wshift = VectorizationBase.pick_vector_width_shift(length(ls, vectorized), size_T)::Tuple{Int,Int}
    # costs =
    # cost_mat[1] / ( unrolled * tiled)
    # cost_mat[2] / ( tiled)
    # cost_mat[3] / ( unrolled)
    # cost_mat[4]
    # @show order
    cost_vec = cost_vec_buf(ls)
    reg_pressure = reg_pres_buf(ls)
    # @inbounds reg_pressure[2] = 1
    # @inbounds reg_pressure[3] = 1
    for n ∈ 1:N
        itersym = order[n]
        # Add to set of defined symbles
        push!(nested_loop_syms, itersym)
        if n == 1
            iter = length(ls, itersym) * length(ls, order[2]) / W
        elseif n > 2
            iter *= Float64(length(ls, itersym))
        end
        # check which vars we can define at this level of loop nest
        for (id, op) ∈ enumerate(ops)
            # isconstant(op) && continue
            # @assert id == identifier(op)+1 # testing, for now
            # won't define if already defined...
            included_vars[id] && continue
            # it must also be a subset of defined symbols
            loopdependencies(op) ⊆ nested_loop_syms || continue
            # # @show nested_loop_syms
            # # @show reduceddependencies(op)
            rd = reduceddependencies(op)
            hasintersection(rd, nested_loop_syms[1:end-length(rd)]) && return 0,0,Inf
            included_vars[id] = true
            unrolledtiled[1,id] = unrolled ∈ loopdependencies(op)
            unrolledtiled[2,id] = tiled ∈ loopdependencies(op)
            iters[id] = iter
            innerloop ∈ loopdependencies(op) && set_upstream_family!(descendentsininnerloop, op, true)
        end
    end
    for (id, op) ∈ enumerate(ops)
        iters[id] == -99.9 && continue
        descendentsininnerloop[id] || continue
        isunrolled = unrolledtiled[1,id]
        istiled = unrolledtiled[2,id]
        rt, lat, rp = cost(op, vectorized, Wshift, size_T)
            # @show instruction(op), rt, lat, rp, iter
        rt *= iters[id]
            # @show isunrolled, istiled
        if isunrolled && istiled # no cost decrease; cost must be repeated
            cost_vec[1] += rt
            reg_pressure[1] += rp
        elseif isunrolled # cost decreased by tiling
            cost_vec[2] += rt
            reg_pressure[2] += rp
        elseif istiled # cost decreased by unrolling
            cost_vec[3] += rt
            reg_pressure[3] += rp
        else# neither unrolled or tiled
            cost_vec[4] += rt
            reg_pressure[4] += rp
        end
    end
    # @show order, vectorized cost_vec reg_pressure
    # @show solve_tilesize(ls, unrolled, tiled, cost_vec, reg_pressure)
    solve_tilesize(ls, unrolled, tiled, cost_vec, reg_pressure)
end

struct LoopOrders
    syms::Vector{Symbol}
    buff::Vector{Symbol}
end
function LoopOrders(ls::LoopSet)
    syms = [s for s ∈ keys(ls.loops)]
    LoopOrders(syms, similar(syms))
end
function Base.iterate(lo::LoopOrders)
    lo.syms, zeros(Int, length(lo.syms))# - 1)
end

function swap!(x, i, j)
    xᵢ, xⱼ = x[i], x[j]
    x[j], x[i] = xᵢ, xⱼ
end
function advance_state!(state)
    N = length(state)
    for n ∈ 1:N
        sₙ = state[n]
        if sₙ == N - n
            if n == N
                return false
            else
                state[n] = 0
            end
        else
            state[n] = sₙ + 1
            break
        end
    end
    true
end
# I doubt this is the most efficient algorithm, but it's the simplest thing
# that I could come up with.
function Base.iterate(lo::LoopOrders, state)
    advance_state!(state) || return nothing
    # # @show state
    syms = copyto!(lo.buff, lo.syms)
    for i ∈ eachindex(state)
        sᵢ = state[i]
        sᵢ == 0 || swap!(syms, i, i + sᵢ)
    end
    syms, state
end
function choose_unroll_order(ls::LoopSet, lowest_cost::Float64 = Inf)
    lo = LoopOrders(ls)
    best_order = lo.syms
    new_order, state = iterate(lo) # right now, new_order === best_order
    best_vec = first(new_order)
    while true
        for new_vec ∈ new_order
            cost_temp = evaluate_cost_unroll(ls, new_order, lowest_cost, new_vec)
            if cost_temp < lowest_cost
                lowest_cost = cost_temp
                best_order = new_order
                best_vec = new_vec
            end
        end
        iter = iterate(lo, state)
        iter === nothing && return best_order, best_vec, lowest_cost
        new_order, state = iter
    end
end
function choose_tile(ls::LoopSet)
    lo = LoopOrders(ls)
    best_order = copyto!(ls.loop_order.bestorder, lo.syms)
    best_vec = first(best_order) # filler
    new_order, state = iterate(lo) # right now, new_order === best_order
    U, T, lowest_cost = 0, 0, Inf
    while true
        for new_vec ∈ @view(new_order[2:end]) # view to skip first
            U_temp, T_temp, cost_temp = evaluate_cost_tile(ls, new_order, new_vec)
            if cost_temp < lowest_cost
                lowest_cost = cost_temp
                U, T = U_temp, T_temp
                best_vec = new_vec
                copyto!(best_order, new_order)
                save_tilecost!(ls)
            end
        end
        iter = iterate(lo, state)
        iter === nothing && return best_order, best_vec, U, T, lowest_cost
        new_order, state = iter
    end
end
# Last in order is the inner most loop
function choose_order(ls::LoopSet)
    if num_loops(ls) > 1
        torder, tvec, tU, tT, tc = choose_tile(ls)
    else
        tc = Inf
    end
    uorder, uvec, uc = choose_unroll_order(ls, tc)
    if num_loops(ls) > 1 && tc ≤ uc
        return torder, tvec, min(tU, tT), tT
        # return torder, tvec, 4, 4#5, 5
    else
        return uorder, uvec, determine_unroll_factor(ls, uorder, first(uorder), uvec), -1
    end
end