text:
#Kacper Kulczak 279079
#Program do zadania 2.12 na pracownie z przedmiotu Analiza Numeryczna M

#iloraz roznicowy
function il_roz(xs,ys,n)
          b = copy(ys)
        for  j = 1:n
            for k = n : -1 : j
                b[k+1] = (b[k+1]- b[k-1+1])/(xs[k+1]- xs[k-j+1])
            end
        end
        return b
end
#Obliczenia wartosci wielomianu interpolowanego w wezlach xs
function inter_wiel(x,xs, iloraz, n)
    temp = 1
    res = iloraz[0+1]
    for i=1:n
        temp *= (x - xs[i-1+1])
         res += iloraz[i+1] *temp
    end
    return res
end
#Automat do interpolacji Newtona
function interpNewton(xs,ys)
        n = length(xs)-1
        iloraz = il_roz(xs,ys,n)
        return x -> inter_wiel(x,xs,iloraz,n)
end

#naturalna funkcja sklejana obliczanie macierzy
function spline_matrix(xs, ys,n)
    h = zeros(n+1)
    lambda = zeros(n+1)
    q=zeros(n+1)
    p=zeros(n+1)
    u=zeros(n+1)
    d=zeros(n+1)
    M=zeros(n+1)
    #parametry h[i]
    for k = 1:n
        h[k+1] = xs[k+1]-xs[k-1+1]
    end
    #parametry lambda[i]
    for k = 0:(n-1)
        lambda[k+1] = h[k+1] / (h[k+1]+h[k+1+1])
    end
    #parametry d[k] (ilorazy różnicowe)

    for k = 1 : (n-1)
        #println("it:",k)
        temp =  il_roz( [ xs[k-1+1], xs[k+1], xs[k+1+1] ],   [ ys[k-1+1], ys[k+1], ys[k+1+1]  ]   ,   2)
        d[k+1] = 6 * temp[3]
    end
    #Wyliczanie M[k] drugiej pochodnej
    for k = 1:(n-1)
        p[k+1] = lambda[k+1] * q[k-1+1] + 2
        q[k+1] = (lambda[k+1] - 1) / p[k+1]
        u[k+1] = (d[k+1] - lambda[k+1] * u[k-1+1]) / p[k+1]
    end
        M[n-1+1]=u[n-1+1]
        for k = (n-2): -1 : 1
            M[k+1]=u[k+1] + q[k+1] *M[k+1+1]
        end
        #=
        println("h: ", h)
        println("l\: ",lambda )
        println("p: ",p)
        println("q: ",q)
        println("u: ",u)
        println("d: ",d)
        println("M: ",M)
        =#
        #@show M
        return  (h,M)
end

#wyliczanie macierzy dla funkcji sklejanych przy krzywej parametrycznej na raz dla x oraz y
function spline_matrix_2(xs,ys,ts,n)
     #  @show n
    h = zeros(n+1)
    lambda = zeros(n+1)
    q=zeros(n+1)
    p=zeros(n+1)
    ux=zeros(n+1)
    uy = zeros(n+1)
    dx=zeros(n+1)
    dy = zeros(n+1)
    Mx=zeros(n+1)
    My=zeros(n+1)
    #parametry h[i]
    for k = 1:n
        h[k+1] = ts[k+1]-ts[k-1+1]
    end
    #parametry lambda[i]
    for k = 0:(n-1)
        lambda[k+1] = h[k+1] / (h[k+1]+h[k+1+1])
    end
    #parametry dy[k] (ilorazy różnicowe)
    for k = 1 : (n-1)
        temp =  il_roz( [ ts[k-1+1], ts[k+1], ts[k+1+1] ],   [ ys[k-1+1], ys[k+1], ys[k+1+1]  ]   ,   2)
        dy[k+1] = 6 * temp[3]
    end
    #parametry dx[k] (ilorazy różnicowe )x
    for k = 1 : (n-1)
        temp =  il_roz( [ ts[k-1+1], ts[k+1], ts[k+1+1] ],   [ xs[k-1+1], xs[k+1], xs[k+1+1]  ]   ,   2)
        dx[k+1] = 6 * temp[3]
    end
    #Wyliczanie M
   for k = 1:(n-1)
        p[k+1] = lambda[k+1] * q[k-1+1] + 2
        q[k+1] = (lambda[k+1] - 1) / p[k+1]
        ux[k+1] = (dx[k+1] - lambda[k+1] * ux[k-1+1]) / p[k+1]
        uy[k+1] = (dy[k+1] - lambda[k+1] * uy[k-1+1]) / p[k+1]
    end
        Mx[n-1+1]=ux[n-1+1]
        My[n-1+1]=uy[n-1+1]
        for k = (n-2): -1 : 1
            Mx[k+1]=ux[k+1] + q[k+1] *Mx[k+1+1]
            My[k+1]=uy[k+1] + q[k+1] *My[k+1+1]
        end
    return  (h,Mx,My)
end

 #Wzór naturalnej funkcji sklejanej
function spline_equation(x,xs,ys,h,M,n)
    if x < xs[1] || x >xs[end]
        error("funkcja sklejana nie zdefiniowana dla tego argumentu")
    end
    k=1
    for i=1:n
        if x <= xs[i+1]
            k = i
            break
        end
    end
    res =  (M[k-1+1] * (xs[k+1]-x)^3 + M[k+1]*((x-xs[k-1+1])^3))/6
    res += (ys[k-1+1]-M[k-1+1]*(h[k+1]^2)/6)*(xs[k+1]-x)
    res += (ys[k+1]   -M[k+1    ]*(h[k+1]^2)/6)*(x-xs[k-1+1])
    res *= (1 / h[k+1])
    return res
end
#Automat do budowania funkcji interpolacyjnej
function interpspline(xs,ys)
    n =length(xs)-1
    h,m = spline_matrix(xs,ys,n)
    return x-> spline_equation(x,xs,ys,h,m,n)
end
#Oszczedne tworzenie naturalnej funkcji sklejanej dla x i y na raz przy krzywej parametrcznej
function interpSplineParam(xs,ys,ts)
    n = length(xs) - 1
    h,mx,my = spline_matrix_2(xs,ys,ts,n)
return (x->spline_equation(x,ts,xs,h,mx,n), x->spline_equation(x, ts, ys, h, my, n)  )
end

#Automat do interpolacji krzywej parametrycznej
function interpparametric(xs,ys)
    n = length(xs)-1
    ns = [1:(n+1);]
    Xfun = interpspline(ns,xs)
    Yfun = interpspline(ns,ys)
    points =  linspace(1, n+1, 1000)
    return (map(Xfun,points),map(Yfun,points))
end

#Automat do tworzenia wykresow krzywej parametrycznej interpolowanej sklejana funckja naturalna.
#Wyswietla wykres, oraz wykres bledow, i norme
function make_param_curve(fx,fy,ts,xss, mytitle)
    xs = map(fx,ts)
    ys = map(fy,ts)
    Splinex,Spliney = interpSplineParam(xs,ys,ts)
   #Krzywa parametryczna
    trace1= scatter(;x = map(Splinex, xss) , y =map(Spliney, xss), mode = "lines", name = "Skejana")
    trace2 = scatter(;x = map(fx,xss), y=map(fy,xss), mode= "lines", name="Krzywa Parametryczna" )
    layout1 = Layout(;title=mytitle)
    #Funkcje błedow i ich wykres
    trace3 = scatter(;x = xss , y =map(t-> abs(Splinex(t) - fx(t)), xss), mode = "lines", name = "|Sx(t) - fx(t)|")
    trace4 = scatter(;x = xss , y =map( t-> abs(Spliney(t) - fy(t)), xss), mode = "lines", name = "|Sy(t) - fy(t)|")
    layout2 = Layout(;title="Funkcje błędu")
    #norma w liczbie
    norm =  max(map(t-> sqrt((fx(t)-Splinex(t))^2+(fy(t)-Spliney(t))^2), xss   ) )
    return (plot([trace2, trace1],layout1),    plot([trace3,trace4],layout2)      , norm)
end

max_stars_repo_name: Magikis/Uniwersity

max_issues_repo_name: Magikis/Uniwersity

max_forks_repo_name: Magikis/Uniwersity


text:
include("bBalancedCut.jl")

##
mat = matread("Netscience_xy.mat")
A = mat["A"]

n = size(A, 1)
T = Int(ceil(10 * log(n)^2))
b = 1/3

min_S, min_expansion = pseudo_balanced_cut_matching(A, T, b)
@show min_S
@show length(min_S) / n
@show min_expansion
@show compute_expansion(A, min_S)

##
using MatrixMarket

A = MatrixMarket.mmread("ca-AstroPh.mtx")
n = size(A, 1)
T = Int(ceil(10 * log(n)^2))

min_S, min_expansion = cut_matching(A, T)
@show min_S
@show length(min_S) / n
@show min_expansion
@show compute_expansion(A, min_S)
##

ext: jl

max_stars_repo_name: luotuoqingshan/Sparsest-Cut

max_issues_repo_name: luotuoqingshan/Sparsest-Cut

max_forks_repo_name: luotuoqingshan/Sparsest-Cut


text:
using Test, MyToolbox
using Base.Iterators

@testset "ziplongest" begin
    @test collect(ziplongest([1,2,3], [1,2])) == [(1,1), (2,2), (3,nothing)]
    @test collect(ziplongest([1,2], [1,2,3])) == collect(zip([1,2,nothing], [1,2,3]))
    @test collect(ziplongest(repeated(1), [1,2,3])) == [(1,1), (1,2), (1,3)]
    @test ziplongest(repeated(1), repeated(2)) == zip(repeated(1), repeated(2))
    @test collect(ziplongest([1,2,3], [1,2], [1]; padval=9)) == [(1,1,1), (2,2,9), (3,9,9)]
end


text:
# Epema propeller
nblades_list = [6, 6, 6]
D = EpemaData.rotor["diameter"]
rtip_list = ones(3) * D / 2
radii_list = fill(rs_desired .* rtip_list[1], 3)
rhub_list = fill(radii_list[1][1], 3)
nJs = 27
Js = ones(3,nJs) .* range(0.01, stop=2.0, length=nJs)'
# Vinf = 19 # m/s
# ns = Vinf ./ Js / 2 / rtip_list[1]
# omegas = ns .* 2 * pi / 60
omegas = ones(3,nJs) .* 5000 * 2 * pi / 60 # rad/s, arbitrary RPM
vinf = 10.0 # or set vinf

chords_raw = EpemaData.rotor["chords_raw"]
chords_list = fill(AS.FM.linear(chords_raw[:,1], chords_raw[:,2], radii_list[1]), 3)

twists_raw = EpemaData.rotor["twists_raw"]
twist = AS.FM.linear(twists_raw[:,1] .* rtip_list[1], twists_raw[:,2], radii_list[1]) .* pi/180 # radians
# set collective at 70\% radius
twists_07_desired = [25.0, 30.0, 35.0] .* pi/180
twists_07 = AS.FM.linear(radii_list[1], twist, 0.7 * rtip_list[1])
twists_list = Vector{Vector{typeof(twist[1])}}(undef, 3)
twist_extra = 3.0 * pi/180
for (i, twistval) in enumerate(twists_07_desired)
    twists_list[i] = twist .+ twistval .- twists_07 .+ twist_extra
end

# extract scanned contour data, smooth, and interpolate to rs_desired
locations = string.(Int.(round.(rs_desired .* 1000, digits=0)))
locations[1:end-1] = "0" .* locations[1:end-1]

##########################################################
# uncomment the next 2 lines to use all the distinct polars:
# contourfilenames = "epema_interpolated_bspline_n30_" .* locations
# repeattag = ""
# or uncomment these lines to try just repeating the r/R = 0.7 airfoil:
contourfilenames = fill("epema_interpolated_bspline_n30_" * locations[8],length(locations))
repeattag = "_repeat07af"
##########################################################
contour_paths_list = fill(joinpath(contour_directory, contourfilenames[1] * ".dat"),3)
# contour_paths_list = fill(joinpath.(contour_directory, contourfilenames .* ".dat"),3)
airfoilnames = fill(contourfilenames, 3)

plotstepi = 1:length(Js)

mach_correction = nothing
# mach_correction = CC.PrandtlGlauert()

Res = [5e4, 1e5, 5e5, 1e6, 1e7]
Machs = [0.0, 0.1, 0.2, 0.3]

Res_lists = fill(Res, 3)
Ms_lists = fill(Machs, 3)
# Res_lists = fill(fill(Res, length(radii_list[1])), 3)
# Ms_lists = fill(fill(Machs, length(radii_list[1])), 3)

index = [1,2,3]
positions = [[0.0, 0.0, -10.0], [0.0, 0.0, 0.0], [0.0, 0.0, 10.0]]
orientations = fill([-1.0, 0.0, 0.0],3)
spin_directions = fill(true,3)

plot_base_name = "epema"
# data_epema_rotor_sweep = AS.rotor_sweep_template(Js, omegas, nblades, rhub_list, rtip_list, radii, chords, twists, airfoilcontours, airfoilnames, index, positions, orientations, spin_directions, Res_list, Ms_list;
#     polar_directory = polar_directory,
#     closefigure = true,
#     useoldfiles = true,
#     rotor_names = [LS.L"\theta_{r/R=0.7} = 25^\circ", LS.L"\theta_{r/R=0.7} = 30^\circ", LS.L"\theta_{r/R=0.7} = 35^\circ"],
#     plot_directory = joinpath(AS.topdirectory, "data", "plots", AS.TODAY),
#     plot_base_name = plot_base_name,
#     plot_extension = ".png"
#     )
data_epema_rotor_sweep = AS.rotor_sweep_template(Js, omegas, nblades_list, rhub_list, rtip_list, mach_correction,
    radii_list, chords_list, twists_list,
    Res_lists, Ms_lists, contour_paths_list,
    index, positions, orientations, spin_directions;
        polar_directory = joinpath(AS.topdirectory, "test", "data", "airfoil", "polars"),
        rotor_names = [LS.L"\theta_{r/R=0.7} = 25^\circ", LS.L"\theta_{r/R=0.7} = 30^\circ", LS.L"\theta_{r/R=0.7} = 35^\circ"],
        plot_polars = true,
        plot_directory = joinpath(AS.topdirectory, "data", "plots", AS.TODAY),
        plot_base_name = plot_base_name,
        plot_extension = ".png",
        step_symbol = LS.L"J",
        close_figure = false,
        use_old_files = true,
        use_old_corrected_files = true
)

step_range = data_epema_rotor_sweep[end-1]
objective = AS.runsimulation!(data_epema_rotor_sweep...)

epema_eta = EpemaData.epema_eta

epema_ct = EpemaData.epema_ct

# Epema xrotor data
ct_xrotor = EpemaData.ct_xrotor

fig = plt.figure(plot_base_name * "_rotor_sweep")
axs = fig.get_axes()

# plot CTs
for (i,data) in enumerate(epema_ct)
    cratio = i / length(epema_ct)
    axs[1].scatter(data[:,1], data[:,2], marker = "o", s = 100, facecolors="none", edgecolors = (0.05, 0.85-cratio*0.7, 0.15 + 0.75 * cratio), label="Experiment")
    axs[1].scatter(ct_xrotor[i][:,1], marker = "x", ct_xrotor[i][:,2], color = (0.05, 0.85-cratio*0.7, 0.15 + 0.75 * cratio), label="XROTOR")
end
axs[1].legend(loc="upper left", bbox_to_anchor=(1.01,1))

# plot efficiencies
for (i,data) in enumerate(epema_eta)
    cratio = i / length(epema_eta)
    axs[3].scatter(data[:,1], data[:,2], marker = "o", s = 100, facecolors="none", edgecolors = (0.05, 0.85-cratio*0.7, 0.15 + 0.75 * cratio), label="Experiment")
end

# save figure
savename = plot_base_name * "_rotor_sweep_extratwist_$(round(twist_extra * 180/pi; digits = 0))" * repeattag * ".pdf"
fig.savefig(joinpath(plot_directory, savename), bbox_inches="tight")

max_stars_repo_name: byuflowlab/AircraftSystems

max_issues_repo_name: byuflowlab/AircraftSystems

max_forks_repo_name: byuflowlab/AircraftSystems


text:
include("sec_6_2.jl")

using PyPlot

r = 3
M = 2r
ppI = 40
store = Store(α, r, M, ppI)

N = 3
t, U = FODEdG!(λ, tmax, f, u0, N, r, M, store)

figure(3)
tt = range(0, tmax, length=201)
exact_u = [ u(t, λ, u0, f, M, store) for t in tt ]
pcwise_t, pcwise_U = evaluate_pcwise_poly!(U, t, ppI, store)
line1 = plot(tt, exact_u, "k:")
line2 = [ plot(pcwise_t[:,n], pcwise_U[:,n], "c") for n = 1:N ]
legend((line1[1], line2[1][1]), (L"$u$", L"$U$"))
grid(true)
xlabel(L"$t$", fontsize=12)
xticks((0.0, 2/3, 4/3, 2), ("0", "2/3", "4/3", "2"))
savefig("fig3.eps")


text:
# Modified from Diercxk.jl:
# https://github.com/kbarbary/Dierckx.jl/blob/master/deps/build.jl  (3-clause BSD license)

# do it by hand since problem with BinDeps

suffix = ""

@unix_only begin
    my_download = @osx ? `curl -O` : `wget`
    suffix = @osx ? "dylib" : "so"
end

@windows_only begin
    error("Package not currently available on Windows")
end

lib_name = "crlibm-1.0beta4"

src_dir = "src"
cd(src_dir)

file = "http://lipforge.ens-lyon.fr/frs/download.php/162/$(lib_name).tar.gz"

println("Downloading the library files from $file")
println("Working in ", pwd())

run(`$(my_download) $file`)
#download(file)
run(`tar xzf $(lib_name).tar.gz`)

#srcdir = "$(src_dir)/$(lib_name)"

cd(lib_name)
println("Working in ", pwd())

suffix = @osx? "dylib" : "so"
run(`./configure CFLAGS=-fpic --silent`)
println("Working in ", pwd())

run(`make -s V=0`)
run(`make -s -f ../shared.mk SUFFIX=$suffix V=0`)


text:
# backproject.jl

export backproject, backproject!

"""
    backproject!(image, view, plan, viewidx)
Backproject a single view.
"""
function backproject!(
    image::AbstractArray{<:RealU, 3},
    view::AbstractMatrix{<:RealU},
    plan::SPECTplan,
    viewidx::Int,
)

    Threads.@threads for z = 1:plan.imgsize[3] # 1:nz
        thid = Threads.threadid() # thread id

        # rotate mumap
        imrotate!((@view plan.mumapr[:, :, z]),
                  (@view plan.mumap[:, :, z]),
                  plan.viewangle[viewidx],
                  plan.planrot[thid],
                  )

    end

    # adjoint of convolving img with psf and applying attenuation map
    Threads.@threads for y = 1:plan.imgsize[2] # 1:ny
        thid = Threads.threadid() # thread id
        # account for half of the final slice thickness
        scale3dj!(plan.exp_mumapr[thid], plan.mumapr, y, -0.5)
        for j = 1:y
            plus3dj!(plan.exp_mumapr[thid], plan.mumapr, j)
        end

        broadcast!(*, plan.exp_mumapr[thid], plan.exp_mumapr[thid], - plan.dy)

        broadcast!(exp, plan.exp_mumapr[thid], plan.exp_mumapr[thid])

        fft_conv_adj!((@view plan.imgr[:, y, :]),
                       view,
                       (@view plan.psfs[:, :, y, viewidx]),
                       plan.planpsf[thid],
                       )

        mul3dj!(plan.imgr, plan.exp_mumapr[thid], y)
    end

    # adjoint of rotating image
    Threads.@threads for z = 1:plan.imgsize[3] # 1:nz
        thid = Threads.threadid()

        imrotate_adj!((@view image[:, :, z]),
                      (@view plan.imgr[:, :, z]),
                      plan.viewangle[viewidx],
                      plan.planrot[thid],
                      )
    end

    return image
end

"""
    backproject!(image, view, plan, thid, viewidx)
Backproject a single view.
"""
function backproject!(
    image::AbstractArray{<:RealU, 3},
    view::AbstractMatrix{<:RealU},
    plan::SPECTplan,
    thid::Int,
    viewidx::Int,
)

    for z = 1:plan.imgsize[3] # 1:nz
        # thid = Threads.threadid() # thread id

        # rotate mumap
        imrotate!((@view plan.mumapr[thid][:, :, z]),
                  (@view plan.mumap[:, :, z]),
                  plan.viewangle[viewidx],
                  plan.planrot[thid],
                  )

    end

    # adjoint of convolving img with psf and applying attenuation map
    for y = 1:plan.imgsize[2] # 1:ny
        thid = Threads.threadid() # thread id
        # account for half of the final slice thickness
        scale3dj!(plan.exp_mumapr[thid], plan.mumapr[thid], y, -0.5)
        for j = 1:y
            plus3dj!(plan.exp_mumapr[thid], plan.mumapr[thid], j)
        end

        broadcast!(*, plan.exp_mumapr[thid], plan.exp_mumapr[thid], - plan.dy)

        broadcast!(exp, plan.exp_mumapr[thid], plan.exp_mumapr[thid])

        fft_conv_adj!((@view plan.imgr[thid][:, y, :]),
                       view,
                       (@view plan.psfs[:, :, y, viewidx]),
                       plan.planpsf[thid],
                       )

        mul3dj!(plan.imgr[thid], plan.exp_mumapr[thid], y)
    end

    # adjoint of rotating image
    for z = 1:plan.imgsize[3] # 1:nz
        imrotate_adj!((@view plan.imgr[thid][:, :, z]),
                      (@view plan.imgr[thid][:, :, z]),
                      plan.viewangle[viewidx],
                      plan.planrot[thid],
                      )
    end

    broadcast!(+, image, image, plan.imgr[thid])

    return image
end

"""
    backproject!(image, views, plan ; index)
Backproject multiple views into `image`.
Array `image` is not initialized to zero; caller must do that.
"""
function backproject!(
    image::AbstractArray{<:RealU, 3},
    views::AbstractArray{<:RealU, 3},
    plan::SPECTplan;
    index::AbstractVector{<:Int} = 1:plan.nview, # all views
)

    # loop over each view index
    image .= zero(plan.T) # must be initialized as zero
    if plan.mode === :fast
        [plan.add_img[i] .= zero(plan.T) for i = 1:plan.nthread]
        Threads.@threads for i in index
            thid = Threads.threadid()
            backproject!(plan.add_img[thid], (@view views[:, :, i]), plan, thid, i)
        end

        for i = 1:plan.nthread
            broadcast!(+, image, image, plan.add_img[i])
        end
    else
        for i in index
            backproject!(plan.add_img, (@view views[:, :, i]), plan, i)
            broadcast!(+, image, image, plan.add_img)
        end
    end

    return image
end

"""
    image = backproject(views, plan ; kwargs...)
SPECT backproject `views`; this allocates the returned 3D array.
"""
function backproject(
    views::AbstractArray{<:RealU, 3},
    plan::SPECTplan;
    kwargs...,
)
    image = zeros(plan.T, plan.imgsize)
    backproject!(image, views, plan; kwargs...)
    return image
end

"""
    image = backproject(views, mumap, psfs, dy; interpmeth, kwargs...)
SPECT backproject `views` using attenuation map `mumap` and PSF array `psfs` for pixel size `dy`.
This method initializes the `plan` as a convenience.
Most users should use `backproject!` instead after initializing those, for better efficiency.
"""
function backproject(
    views::AbstractArray{<:RealU, 3}, # [nx,nz,nview]
    mumap::AbstractArray{<:RealU, 3}, # [nx,ny,nz] attenuation map, must be 3D, possibly zeros()
    psfs::AbstractArray{<:RealU, 4},
    dy::RealU;
    interpmeth::Symbol = :two,
    mode::Symbol = :fast,
    kwargs...,
)

    size(mumap,1) == size(mumap,1) == size(views,1) ||
        throw(DimensionMismatch("nx"))
    size(mumap,3) == size(views,2) || throw(DimensionMismatch("nz"))
    plan = SPECTplan(mumap, psfs, dy; interpmeth, mode, kwargs...)
    return backproject(views, plan; kwargs...)
end


text:
module DataUtils
#
# FIXME messing with Generics and Show() in particular
#
import DataFrames: DataFrame
import Statistics: mean, median, std, quantile
import Parameters: @with_kw
using Definitions

export summarise_over_positive, add_to!
export initialise, MinMaxes, show

"""
FIXME passing a type like this can't be idiomatic ...
d needs to be some type of enum (something which instances can work with)
"""
function initialise( x :: Real, d :: DataType )
    dict = Dict{d,Real}()
    for i in instances( d )
        dict[i] = x
    end
    dict
end

@with_kw mutable struct MinMaxes{T}
    ## FIXME I don't think passing T is idiomatic ..
    max ::Dict{T,Real}=initialise( -99999999999999999.9, T )
    min ::Dict{T,Real}=initialise( 99999999999999.999, T )
    sum ::Dict{T,Real}=initialise( 0.0, T )
    poscounts ::Dict{T,Real} = initialise( 0.0, T )
    n :: Integer = 0
end

function add_to!( mms :: MinMaxes, addn :: Dict )
    kys = keys( addn )
    mms.n += 1
    for k in kys
        mms.max[k] = max( mms.max[k], addn[k] )
        mms.min[k] = min( mms.min[k], addn[k] )
        mms.sum[k] += addn[k]
        if addn[k]>0
            mms.poscounts[k]+= 1.0
        end
    end
end

import Base.show

function show( io::IO, mms :: MinMaxes{T} ) where T
    s = ""
    ks = sort(collect(keys(mms.poscounts)))
    for k in ks
        maxx = mms.max[k]
        minx = mms.min[k]
        pc = mms.poscounts[k]
        if pc>0
            mean = mms.sum[k] / pc
        else
            mean = 0.0
        end
        #show( io, "$k = (max=$(maxx), min=$(minx), mean=$mean poscounts=$pc) \n");
        println( io, "$k = (max=$(maxx), min=$(minx), mean=$mean poscounts=$pc)");
    end
end

function show( mms :: MinMaxes{T} ) where T
    show( stdout, mms )
end

"""
Means, etc. over just the positive and non-missing elements of a df column
"""
function summarise_over_positive( df::DataFrame, col::Symbol )::NamedTuple
   m = (df[!,col].!== missing) .&
       (df[!,col] .> 0.0)
   sz = size(df[m,col])[1]
   if sz > 0
       (
       num_positive = sz,
       frame_size = size(df[!,col])[1],
       max=maximum(df[m,col]),
       min=minimum(df[m,col]),
       mean=mean(df[m,col]),
       stdev=std(df[m,col]),
       deciles = quantile( df[m,col], [0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9] )
      )
   else
      ( num_positive = 0,
        frame_size = size(df[!,col])[1],
        max=0,
        min=0,
        mean=0,
        stdev=0,
        deciles=[] )
   end
end

end


text:
using EchelleInstruments
using Test

@testset "EchelleInstruments.jl" begin
    @warn "We will need public data files before incorporating tests here."
    # Write your tests here.
end

max_stars_repo_name: RvSpectML/EchelleInstruments

max_issues_repo_name: RvSpectML/EchelleInstruments

max_forks_repo_name: RvSpectML/EchelleInstruments


text:
# fallbacks for inbands_getindex and inbands_setindex!
@inline function inbands_getindex(x::AbstractArray, i::Integer, j::Integer)
    @inbounds r = getindex(x, i, j)
    r
end
@inline function inbands_setindex!(x::AbstractArray, v, i::Integer, j::Integer)
    @inbounds r = setindex!(x, v, i, j)
    r
end

inbands_getindex(x::Adjoint, i::Integer, j::Integer) =
    inbands_getindex(parent(x), j, i)'
inbands_getindex(x::Transpose, i::Integer, j::Integer) =
    transpose(inbands_getindex(parent(x), j, i))
inbands_setindex!(x::Adjoint, v, i::Integer, j::Integer) =
    inbands_setindex!(parent(x), v', j, i)
inbands_setindex!(x::Transpose, v, i::Integer, j::Integer) =
    inbands_setindex!(parent(x), transpose(v), j, i)

###
# Lazy getindex
# this uses a lazy-materialize idiom to construct a matrix based
# on the memory layout
###

@inline _sub_materialize(_, V) = Array(V)
@inline _sub_materialize(::AbstractBandedLayout, V) = BandedMatrix(V)
@inline _materialize(V::SubArray) = _sub_materialize(MemoryLayout(V), V)

@inline _lazy_getindex(A, I...) = _materialize(view(A, I...))
@inline getindex(A::AbstractBandedMatrix, kr::Colon, jr::Colon) = _lazy_getindex(A, kr, jr)
@inline getindex(A::AbstractBandedMatrix, kr::Colon, jr::AbstractUnitRange) = _lazy_getindex(A, kr, jr)
@inline getindex(A::AbstractBandedMatrix, kr::AbstractUnitRange, jr::Colon) = _lazy_getindex(A, kr, jr)
@inline getindex(A::AbstractBandedMatrix, kr::AbstractUnitRange, jr::AbstractUnitRange) = _lazy_getindex(A, kr, jr)

@inline getindex(A::AbstractMatrix, b::Band) = _lazy_getindex(A, b)
@inline getindex(A::AbstractMatrix, kr::BandRangeType, j::Integer) = _lazy_getindex(A, kr, j)
@inline getindex(A::AbstractMatrix, k::Integer, jr::BandRangeType) = _lazy_getindex(A, k, jr)


text:
# takes in both legendre coefficients and rescaling,
# matches meta-data and writes to single file
using SARA
using JSON
using HDF5

# returns true if dwell time and peak temperature are the same
same_meta(f) = g->same_meta(f, g)
function same_meta(f, g)
    fm = f["meta"]
    gm = g["meta_data"]
    fm["temp"] == gm["Tpeak"] && fm["dwell"] == gm["dwell"]
end

# put path to legendre_coefficients.json and bias.json here
path = "/Users/sebastianament/Documents/SEA/XRD Analysis/SARA/Bi2O3_19F44_01/"
file = "legendre_coefficients.json"
f = JSON.parsefile(path * file)
f = f["spectra"] # don't need wavelengths here

file = "bias.json"
g = JSON.parsefile(path * file)

const broken_index = 210 # not usable data
deleteat!(f, broken_index)
indices = zeros(Int, length(f))
for i in 1:length(f)
    fi = f[i]
    j = findfirst(same_meta(fi), g)
    if isnothing(j) # used to find broken index
        println(i)
        continue
    end
    indices[i] = j
    fi["rescaling_parameters"] = g[j]["rescaling_parameters"]
    f[i] = fi # technically not necessary since dictionaries are references
end

h = open(path * "Bi2O3_19F44_01_inner_loop_data.json", "w")
JSON.print(h, f)
close(h)


text:
# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.

function _add_transmission_line!(
    model::JuMP.Model,
    lm::TransmissionLine,
    f::ShiftFactorsFormulation,
)::Nothing
    overflow = _init(model, :overflow)
    for t in 1:model[:instance].time
        overflow[lm.name, t] = @variable(model, lower_bound = 0)
        add_to_expression!(
            model[:obj],
            overflow[lm.name, t],
            lm.flow_limit_penalty[t],
        )
    end
    return
end

function _setup_transmission(
    model::JuMP.Model,
    formulation::ShiftFactorsFormulation,
)::Nothing
    instance = model[:instance]
    isf = formulation.precomputed_isf
    lodf = formulation.precomputed_lodf
    if length(instance.buses) == 1
        isf = zeros(0, 0)
        lodf = zeros(0, 0)
    elseif isf === nothing
        @info "Computing injection shift factors..."
        time_isf = @elapsed begin
            isf = UnitCommitment._injection_shift_factors(
                lines = instance.lines,
                buses = instance.buses,
            )
        end
        @info @sprintf("Computed ISF in %.2f seconds", time_isf)
        @info "Computing line outage factors..."
        time_lodf = @elapsed begin
            lodf = UnitCommitment._line_outage_factors(
                lines = instance.lines,
                buses = instance.buses,
                isf = isf,
            )
        end
        @info @sprintf("Computed LODF in %.2f seconds", time_lodf)
        @info @sprintf(
            "Applying PTDF and LODF cutoffs (%.5f, %.5f)",
            formulation.isf_cutoff,
            formulation.lodf_cutoff
        )
        isf[abs.(isf).<formulation.isf_cutoff] .= 0
        lodf[abs.(lodf).<formulation.lodf_cutoff] .= 0
    end
    model[:isf] = isf
    model[:lodf] = lodf
    return
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Test of key bindings.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

@testset "Key bindings" begin
    pagerd = _create_pagerd("")

    TerminalPager.set_keybinding('L', :left)
    k = TerminalPager.Keystroke(value = "L")

    pagerd.redraw          = false
    pagerd.start_row       = 10
    pagerd.start_col       = 10
    pagerd.lines_cropped   = 20
    pagerd.columns_cropped = 20
    pagerd.display_size    = (20,20)

    TerminalPager._pager_key_process!(pagerd, k)

    @test pagerd.start_row == 10
    @test pagerd.start_col == 9
    @test pagerd.redraw    == true
    @test pagerd.event     == nothing

    TerminalPager.delete_keybinding('L')

    pagerd.redraw          = false
    pagerd.start_row       = 10
    pagerd.start_col       = 10
    pagerd.lines_cropped   = 20
    pagerd.columns_cropped = 20
    pagerd.display_size    = (20,20)

    TerminalPager._pager_key_process!(pagerd, k)

    @test pagerd.start_row == 10
    @test pagerd.start_col == 10
    @test pagerd.redraw    == false
    @test pagerd.event     == nothing

    TerminalPager.set_keybinding(:left, :bol, shift = true)
    k = TerminalPager.Keystroke(value = :left, shift = true)

    pagerd.redraw          = false
    pagerd.start_row       = 10
    pagerd.start_col       = 10
    pagerd.lines_cropped   = 20
    pagerd.columns_cropped = 20
    pagerd.display_size    = (20,20)

    TerminalPager._pager_key_process!(pagerd, k)

    @test pagerd.start_row == 10
    @test pagerd.start_col == 1
    @test pagerd.redraw    == true
    @test pagerd.event     == nothing

    TerminalPager.delete_keybinding(:left, shift = true)

    pagerd.redraw          = false
    pagerd.start_row       = 10
    pagerd.start_col       = 10
    pagerd.lines_cropped   = 20
    pagerd.columns_cropped = 20
    pagerd.display_size    = (20,20)

    TerminalPager._pager_key_process!(pagerd, k)

    @test pagerd.start_row == 10
    @test pagerd.start_col == 10
    @test pagerd.redraw    == false
    @test pagerd.event     == nothing

    TerminalPager.reset_keybindings()

    pagerd.redraw          = false
    pagerd.start_row       = 10
    pagerd.start_col       = 15
    pagerd.lines_cropped   = 20
    pagerd.columns_cropped = 20
    pagerd.display_size    = (20,20)

    TerminalPager._pager_key_process!(pagerd, k)

    @test pagerd.start_row == 10
    @test pagerd.start_col == 5
    @test pagerd.redraw    == true
    @test pagerd.event     == nothing
end


text:
include("ReflectBlock.jl")

num_grover_step(prob::Real) = Int(round(pi/4/sqrt(prob))) - 1

function GroverSearch(oracle, num_bit::Int; psi = uniform_state(num_bit))
    uni_reflect = Reflect(uniform_state(num_bit))
    # solve a real search problem
    num_iter = num_grover_step(1.0/(1<<num_bit))

    for i in 1:num_iter
        psi = uni_reflect(oracle(psi))
    end
    psi
end

max_stars_repo_name: Roger-luo/QCBM

max_issues_repo_name: Roger-luo/QCBM

max_forks_repo_name: Roger-luo/QCBM


text:
mutable struct DelayDirectJumpAggregation{T,S,F1,F2,RNG,IType} <: AbstractDSSAJumpAggregator
    next_jump::Int
    prev_jump::Int
    next_jump_time::T
    end_time::T
    cur_rates::Vector{T} # cur_rates here is the cumsum of rates
    sum_rate::T
    ma_jumps::S
    rates::F1
    affects!::F2
    save_positions::Tuple{Bool,Bool}
    rng::RNG
    next_delay::Union{Nothing,Vector{Int}}
    num_next_delay::Union{Nothing,Vector{Int}}
    time_to_next_jump::T
    next_delay_time::T
    shadow_integrator::IType
    copied::Bool
end
mutable struct ShadowIntegrator{uType,chanType,T}
    u::uType
    de_chan::chanType
    delayjumpsets::DelayJumpSet
    cur_rates::Vector{T}
end

function DelayDirectJumpAggregation(nj::Int, njt::T, et::T, crs::Vector{T}, sr::T, maj::S, rs::F1, affs!::F2, sps::Tuple{Bool,Bool}, rng::RNG; u0, kwargs...) where {T,S,F1,F2,RNG}
    ttnj = zero(et)
    ndt = zero(et)
    shadow_integrator = ShadowIntegrator{typeof(u0),Vector{Vector{T}}, T}(copy(u0), [Vector{T}()], DelayJumpSet(Dict(), Dict(), Dict()), copy(crs))
    nd = nothing
    nnd = [] # in the Direct Method the number of next delay equals always 1
    DelayDirectJumpAggregation{T,S,F1,F2,RNG,typeof(shadow_integrator)}(nj, nj, njt, et, crs, sr, maj, rs, affs!, sps, rng, nd, nnd, ttnj, ndt, shadow_integrator, false)
end

function aggregate(aggregator::DelayDirect, u, p, t, end_time, constant_jumps, ma_jumps, save_positions, rng; kwargs...)

    # handle constant jumps using function wrappers
    rates, affects! = get_jump_info_fwrappers(u, p, t, constant_jumps)
    build_jump_aggregation(DelayDirectJumpAggregation, u, p, t, end_time, ma_jumps, rates, affects!, save_positions, rng; u0=copy(u), kwargs...)
end

function initialize!(p::DelayDirectJumpAggregation, integrator, u, params, t)
    p.shadow_integrator.delayjumpsets = integrator.delayjumpsets
    generate_jumps!(p, integrator, u, params, t)
    nothing
end

@inline function execute_jumps!(p::DelayDirectJumpAggregation, integrator, u, params, t)
    update_state_delay_direct!(p, integrator, integrator.u, t)
    nothing
end

function generate_jumps!(p::DelayDirectJumpAggregation, integrator, u, params, t)
    generate_time_to_next_jump!(p, integrator, params, t)
    @fastmath p.next_jump_time = t + p.time_to_next_jump
    @inbounds p.next_jump = searchsortedfirst(p.cur_rates, rand(p.rng) * p.sum_rate) #
    nothing
end
"""
    Create time_to_next_jump based on the shawdow integrator, the goal is to use p.copied to avoid the case when no changes happened wrt integrator.u such that the allocation can be optimised
"""
@inline function generate_time_to_next_jump!(p::DelayDirectJumpAggregation, integrator, params, t)

    fill_cum_rates_and_sum!(p, integrator.u, params, t)
    r1 = rand(p.rng)
    if isempty(reduce(vcat, integrator.de_chan))
        ttnj = -log(r1) / p.cur_rates[end]
        ttnj_last = ttnj
        p.copied = false
    else
        find_next_delay_num!(p, integrator.de_chan)
        aₜ = p.cur_rates[end] * p.next_delay_time
        F = one(t) - exp(-aₜ)
        aₜ_ = zero(aₜ)
        prev_T1 = zero(t)
        cur_T1 = zero(t)
        cur_T1 += p.next_delay_time
        if F < r1 # if enter the loop
            shadow_integrator = p.shadow_integrator
            shadow_integrator.u = copy(integrator.u)
            shadow_integrator.de_chan = deepcopy(integrator.de_chan)
            p.copied = true
        else
            p.copied = false
        end
        while F < r1
            shift_delay_channel!(shadow_integrator.de_chan, p.next_delay_time)
            update_delay_channel!(shadow_integrator.de_chan)
            update_delay_complete!(p, shadow_integrator)

            calculate_sum_rate!(p, shadow_integrator, shadow_integrator.u, params, t + cur_T1)

            find_next_delay_num!(p, shadow_integrator.de_chan)
            prev_T1 = cur_T1 # to avoid cur_T1 = Inf
            cur_T1 += p.next_delay_time

            aₜ_ = aₜ # backup aₜ
            aₜ += shadow_integrator.cur_rates[end] * (p.next_delay_time)
            F = one(t) - exp(-aₜ)
        end
        sum_ = p.copied ? shadow_integrator.cur_rates[end] : p.sum_rate
        ttnj_last = (-log(one(t) - r1) - aₜ_) / sum_
        ttnj = prev_T1 + ttnj_last
    end

    if p.copied
        shift_delay_channel!(shadow_integrator.de_chan, ttnj_last)
        update_delay_channel!(shadow_integrator.de_chan)
        fill_cum_rates_and_sum!(p, shadow_integrator.u, params, t + ttnj)
    end
    p.time_to_next_jump = ttnj
    nothing
end

"""
    update_delay_at_tstop_test!

this function will change
- p.next_delay, p.num_next_delay, p.next_delay_time
- integrator.u, integrator.de_chan
"""
function update_delay_at_tstop_test!(p, integrator, params, t, tgap)
    # direct_algo!(p, integrator, params, t; tgap = tgap)
    cur_T1 = zero(t)
    prev_T1 = zero(t)
    find_next_delay_num!(p, integrator.de_chan)
    cur_T1 += p.next_delay_time
    while cur_T1 <= tgap
        shift_delay_channel!(integrator.de_chan, p.next_delay_time)
        update_delay_channel!(integrator.de_chan)
        update_delay_complete!(p, integrator)
        find_next_delay_num!(p, integrator.de_chan)
        prev_T1 = cur_T1 # to avoid cur_T1 = Inf
        cur_T1 += p.next_delay_time
    end
    ttnj_last = tgap - prev_T1
    shift_delay_channel!(integrator.de_chan, ttnj_last)
    update_delay_channel!(integrator.de_chan)
    nothing
end

"""
    function fill_cum_rates_and_sum!(p::DelayDirectJumpAggregation, u, params, t)

This function modifies `p.cur_rates` and `p.sum_rate`
"""
function fill_cum_rates_and_sum!(p::DelayDirectJumpAggregation, u, params, t)
    prev_rate = zero(t)
    new_rate = zero(t)
    cur_rates = p.cur_rates
    # mass action rates
    majumps = p.ma_jumps
    idx = get_num_majumps(majumps)
    @inbounds for i in 1:idx
        new_rate = evalrxrate(u, i, majumps)
        cur_rates[i] = new_rate + prev_rate
        prev_rate = cur_rates[i]
    end
    # constant jump rates
    idx += 1
    rates = p.rates
    @inbounds for i in 1:length(p.rates)
        new_rate = rates[i](u, params, t)
        cur_rates[idx] = new_rate + prev_rate
        prev_rate = cur_rates[idx]
        idx += 1
    end

    @inbounds sum_rate = cur_rates[end]
    p.sum_rate = sum_rate
    nothing
end
"""
    Only changes s.sum_rate
"""
function calculate_sum_rate!(p, s::ShadowIntegrator, u, params, t)
    prev_rate = zero(t)
    new_rate = zero(t)
    cur_rates = s.cur_rates
    majumps = p.ma_jumps
    idx = get_num_majumps(majumps)
    @inbounds for i in 1:idx
        new_rate = evalrxrate(u, i, majumps)
        cur_rates[i] = new_rate + prev_rate
        prev_rate = cur_rates[i]
    end
    # constant jump rates
    idx += 1
    rates = p.rates
    @inbounds for i in 1:length(rates)
        new_rate = rates[i](u, params, t)
        cur_rates[idx] = new_rate + prev_rate
        prev_rate = cur_rates[idx]
        idx += 1
    end
end

@inline function update_state_delay_direct!(p::DelayDirectJumpAggregation, integrator, u, t)
    @unpack ma_jumps, next_jump, time_to_next_jump = p
    @unpack delay_trigger_set, delay_interrupt_set = integrator.delayjumpsets
    if p.copied
        integrator.u = copy(p.shadow_integrator.u)
        integrator.de_chan = deepcopy(p.shadow_integrator.de_chan)
    else
        # update_delay_at_tstop_test!(p, integrator, p, t, time_to_next_jump)
        shift_delay_channel!(integrator.de_chan, time_to_next_jump)
        update_delay_channel!(integrator.de_chan)
    end

    num_ma_rates = get_num_majumps(ma_jumps)
    if next_jump <= num_ma_rates # if the next jump is a mass action jump
        if u isa SVector
            integrator.u = executerx(u, next_jump, ma_jumps)
        else
            @inbounds executerx!(integrator.u, next_jump, ma_jumps)
        end
    else
        idx = next_jump - num_ma_rates
        @inbounds p.affects![idx](integrator)
    end

    if next_jump in delay_interrupt_set
        update_delay_interrupt!(p, integrator)
    elseif next_jump in delay_trigger_set
        update_delay_trigger!(p, integrator)
    end
    # save jump that was just executed
    p.prev_jump = next_jump
    nothing
end


text:
using Plots, Plots.PlotMeasures, LaTeXStrings

# We don't use BenchmarkTools for this model

model = "HEAT3D"
cases = ["HEAT01-discrete", "HEAT02-discrete"]
if TEST_LONG
    push!(cases, "HEAT03-discrete")
    push!(cases, "HEAT04-discrete")
end
push!(cases, "HEAT01-continuous")
push!(cases, "HEAT02-continuous")

results = Dict(model => Dict(c => -1.0 for c in cases))
max_temp = Dict()
validation = []
Tmax = [0.10369, 0.02966, 0.01716, 0.01161, 0.01005]

include("heat3d.jl")

Δ = 1e-4
NSTEPS = 2_000
idx = 0 # index for the case considered

# ----------------------------------------
#  HEAT01 (discrete time)
# ----------------------------------------
idx += 1
A, _, Ω₀, ℓ = heat01(δ=0.02)
ivp = @ivp(x' = A * x, x ∈ Universe(size(A, 1)), x(0) ∈ Ω₀);
sol = solve(ivp, T=40.0, alg=LGG09(δ=0.02, template=[ℓ], approx_model=NoBloating())) # warm-up run
results[model][cases[idx]] = @elapsed sol = solve(ivp, T=40.0, alg=LGG09(δ=0.02, template=[ℓ], approx_model=NoBloating()))
max_out = ρ(ℓ, sol)
max_temp[cases[idx]] = max_out
property = max_out ∈ Tmax[1] .. Tmax[1] + Δ
push!(validation, Int(property))

sol = nothing
GC.gc()

# ----------------------------------------
#  HEAT02 (discrete time)
# ----------------------------------------
idx += 1
A, _, Ω₀, ℓ = heat02(δ=0.02)
ivp = @ivp(x' = A * x, x ∈ Universe(size(A, 1)), x(0) ∈ Ω₀);
sol = solve(ivp, T=40.0, alg=LGG09(δ=0.02, template=[ℓ], approx_model=NoBloating())) # warm-up run
results[model][cases[idx]] = @elapsed sol = solve(ivp, T=40.0, alg=LGG09(δ=0.02, template=[ℓ], approx_model=NoBloating()))
max_out = ρ(ℓ, sol)
max_temp[cases[idx]] = max_out
property = max_out ∈ Tmax[2] .. Tmax[2] + Δ
push!(validation, Int(property))

sol = nothing
GC.gc()

# ----------------------------------------
#  HEAT03 (discrete time)
# ----------------------------------------
if TEST_LONG
    idx += 1
    A, Aᵀδ, Ω₀, ℓ = heat03(δ=0.02)

    # warm-up run
    out = Vector{Float64}(undef, 1)
    reach_homog_krylov_LGG09!(out, Ω₀, Aᵀδ, sparse(ℓ), 1, m=94, tol=1e-8, hermitian=true)

    out = Vector{Float64}(undef, NSTEPS)
    results[model][cases[idx]] = @elapsed reach_homog_krylov_LGG09!(out, Ω₀, Aᵀδ,
                            sparse(ℓ), NSTEPS, m=94, tol=1e-8, hermitian=true)
    max_out = maximum(out)
    max_temp[cases[idx]] = max_out
    property = max_out ∈ Tmax[3] .. Tmax[3] + Δ
    push!(validation, Int(property))

    out = nothing
    GC.gc()
end

# ----------------------------------------
#  HEAT04 (discrete time)
# ----------------------------------------
if TEST_LONG
    idx += 1
    A, Aᵀδ, Ω₀, ℓ = heat04(δ=0.02)
    out = Vector{Float64}(undef, NSTEPS)
    results[model][cases[idx]] = @elapsed reach_homog_krylov_LGG09!(out, Ω₀, Aᵀδ,
                            sparse(ℓ), NSTEPS, m=211, tol=1e-8, hermitian=true)
    max_out = maximum(out)
    max_temp[cases[idx]] = max_out
    property = max_out ∈ Tmax[4] .. Tmax[4] + Δ
    push!(validation, Int(property))

    out = nothing
    GC.gc()
end

# ----------------------------------------
#  HEAT01 (continuous time)
# ----------------------------------------
idx += 1
A, _, Ω₀, ℓ = heat01(δ=0.02)
ivp = @ivp(x' = A * x, x ∈ Universe(size(A, 1)), x(0) ∈ Ω₀);
sol = solve(ivp, T=1.0, alg=LGG09(δ=0.02, template=[ℓ])) # warm-up run
results[model][cases[idx]] = @elapsed sol = solve(ivp, T=40.0, alg=LGG09(δ=0.02, template=[ℓ]))
max_out = ρ(ℓ, sol)
max_temp[cases[idx]] = max_out
property = Tmax[1] ≤ max_out
push!(validation, Int(property))

sol = nothing
GC.gc()

# ----------------------------------------
#  HEAT02 (continuous time)
# ----------------------------------------
idx += 1
A, _, Ω₀, ℓ = heat02(δ=0.02)
ivp = @ivp(x' = A * x, x ∈ Universe(size(A, 1)), x(0) ∈ Ω₀)
results[model][cases[idx]] = @elapsed sol = solve(ivp, T=40.0, alg=LGG09(δ=0.02, template=[ℓ]))
max_out = ρ(ℓ, sol)
max_temp[cases[idx]] = max_out
property = Tmax[2] ≤ max_out
push!(validation, Int(property))

sol = nothing
GC.gc()

# ==============================================================================
# Save benchmark results
# ==============================================================================

# export runtimes
runtimes = Dict()
for (i, c) in enumerate(cases)
    t = results[model][c] # in seconds
    runtimes[c] = round(t, digits=4)
end

for (i, c) in enumerate(cases)
    max_temp_c = round(max_temp[c], digits=6)
    print(io, "JuliaReach, $model, $c, $(validation[i]), $(runtimes[c]), $max_temp_c\n")
end

# ==============================================================================
# Plot
# ==============================================================================

function flowpipe_HEAT01()
    δ = 0.02
    NSTEPS = 2000
    A, Aᵀδ, Ω₀, ℓ = heat01(δ=δ)

    out_plus = Vector{Float64}(undef, NSTEPS)
    reach_homog_krylov_LGG09!(out_plus, Ω₀, Aᵀδ, sparse(ℓ), NSTEPS, hermitian=true)

    out_minus = Vector{Float64}(undef, NSTEPS)
    reach_homog_krylov_LGG09!(out_minus, Ω₀, Aᵀδ, sparse(-ℓ), NSTEPS, hermitian=true)

    [Interval(i*δ, (i+1)*δ) × Interval(-out_minus[i+1], out_plus[i+1]) for i in 0:NSTEPS-1]
end

sol = flowpipe_HEAT01()
fig = Plots.plot()

Plots.plot!(fig, sol, linecolor=:blue, color=:blue, alpha=0.8,
    tickfont=font(30, "Times"), guidefontsize=45,
    xlab=L"t",
    ylab=L"x_{62}", # in the original model this is called x62
    xtick=[0., 10., 20., 30., 40.], ytick=[0, 0.025, 0.05, 0.075, 0.1],
    xlims=(0., 40.), ylims=(0.0, 0.11),
    bottom_margin=6mm, left_margin=2mm, right_margin=4mm, top_margin=3mm,
    size=(1000, 1000))

savefig(fig, joinpath(@__DIR__, "ARCH-COMP21-JuliaReach-Heat3D.pdf"))

sol = nothing
GC.gc()


text:
using Test

@testset "TESTSETS: Test read TCX dir" begin
    @testset "CASE: Test dir with NO TCX file" begin
	err, _ =  TCX.parse_tcx_dir(tempdir())
        @test err == 404
    end

    @testset "CASE: Test dir with TCX file" begin
        err, _ =  TCX.parse_tcx_dir(@__DIR__)
        @test err == 200
    end

    @testset "CASE: Test relative dir with TCX file" begin
        err, _ =  TCX.parse_tcx_dir(".")
        @test err == 200
    end
    @testset "CASE: Test DataFrame after process dir" begin
        err, ta =  TCX.parse_tcx_dir(".")
        @test (err == 200) & (size(getDataFrame(ta), 1) == 9273)
    end
end


text:
cd(@__DIR__); include("setups/path128.jl")
gr(dpi = 200)

##
f = 1.5 * 𝚽[:, 11] + 𝚽[:, 31] .* characteristic(1:32, N) +
    𝚽[:, 61] .* characteristic(1:64, N) + 0.5 * 𝚽[:, 111]

## (a)
plt = path_spectrogram(f, distDCT, 𝚽; c = 0.01)
savefig(plt, "../figs/Path128_Spectrogram_CoeffEnergy_sig001dmax.png")

## (b)
plt = path_spectrogram(f, distDCT, 𝚽; c = 0.02)
savefig(plt, "../figs/Path128_Spectrogram_CoeffEnergy_sig002dmax.png")

## (c)
plt = path_spectrogram(f, distDCT, 𝚽; c = 0.05)
savefig(plt, "../figs/Path128_Spectrogram_CoeffEnergy_sig005dmax.png")

## (d)
plt = path_spectrogram(f, distDCT, 𝚽; c = 0.1)
savefig(plt, "../figs/Path128_Spectrogram_CoeffEnergy_sig01dmax.png")


text:
using DataDeps

register(DataDep(
  "cpu-act",
  "https://www.openml.org/d/197",
  "https://www.openml.org/data/get_csv/3634/dataset_2183_cpu_act.arff",
  "316cdd36b28d1af79584b043b1ffc16bc492fe03f49b2d653eb09468f55392ed",
  post_fetch_method=(path -> begin
    UCIData.process_dataset(path,
        target_index=22,
        feature_indices=1:21,
        header=true,
    )
  end),
))

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


text:
isdefined(Base, :__precompile__) && __precompile__()

module LibCURL

import Compat: is_windows

const time_t = Int
const size_t = Csize_t
const curl_off_t = Int64

include("lC_exports_h.jl")
include("lC_common_h.jl")

const libcurl = if is_windows()
    Pkg.dir("WinRPM","deps","usr","$(Sys.ARCH)-w64-mingw32","sys-root","mingw","bin","libcurl-4")
else
    "libcurl"
end

include("lC_curl_h.jl")

curl_easy_setopt(handle, opt, ptrval::Array) = ccall((:curl_easy_setopt, libcurl), CURLcode, (Ptr{CURL}, CURLoption, Ptr{Void}), handle, opt, pointer(ptrval))
curl_easy_setopt(handle, opt, ptrval::Integer) = ccall((:curl_easy_setopt, libcurl), CURLcode, (Ptr{CURL}, CURLoption, Clong), handle, opt, ptrval)
curl_easy_setopt{T}(handle, opt, ptrval::Ptr{T}) = ccall((:curl_easy_setopt, libcurl), CURLcode, (Ptr{CURL}, CURLoption, Ptr{T}), handle, opt, ptrval)
curl_easy_setopt(handle, opt, ptrval::AbstractString) = ccall((:curl_easy_setopt, libcurl), CURLcode, (Ptr{CURL}, CURLoption, Ptr{UInt8}), handle, opt, ptrval)

curl_multi_setopt(handle, opt, ptrval::Array) = ccall((:curl_multi_setopt, libcurl), CURLMcode, (Ptr{CURLM}, CURLMoption, Ptr{Void}), handle, opt, pointer(ptrval))
curl_multi_setopt(handle, opt, ptrval::Integer) = ccall((:curl_multi_setopt, libcurl), CURLMcode, (Ptr{CURLM}, CURLMoption, Clong), handle, opt, ptrval)
curl_multi_setopt{T}(handle, opt, ptrval::Ptr{T}) = ccall((:curl_multi_setopt, libcurl), CURLMcode, (Ptr{CURLM}, CURLMoption, Ptr{T}), handle, opt, ptrval)
curl_multi_setopt(handle, opt, ptrval::AbstractString) = ccall((:curl_multi_setopt, libcurl), CURLMcode, (Ptr{CURLM}, CURLMoption, Ptr{UInt8}), handle, opt, ptrval)

curl_easy_getinfo(handle, opt, ptrval::Array) = ccall((:curl_easy_getinfo, libcurl), CURLcode, (Ptr{CURL}, CURLoption, Ptr{Void}), handle, opt, pointer(ptrval))
curl_easy_getinfo(handle, opt, ptrval::AbstractString) = ccall((:curl_easy_getinfo, libcurl), CURLcode, (Ptr{CURL}, CURLoption, Ptr{UInt8}), handle, opt, ptrval)

include("lC_defines_h.jl")

include("Mime_ext.jl")
export Mime_ext

end


text:
module MechGluecode
using Requires
export value, ODE_DEFAULT_NORM, UNITLESS_ABS2, Unitfu, norm

function __init__()
    @require MechanicalUnits = "e6be9192-89dc-11e9-36e6-5dbcb28f419e" begin
        @require Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" begin
            @info "Plots => using MechGluePlots"
            @eval using MechGluePlots
        end
        @require Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" begin
            @info "Interpolations => using MechGlueInterpolations"
            @eval using MechGlueInterpolations
        end
        @require ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78" begin
            #@info "ModelingToolkit => using MechGlueModelingToolkit"
            #@eval using MechGlueModelingToolkit
        end
        @require RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd" begin
           # @info "RecursiveArrayTools => using MechGlueRecursiveArrayTools"
           # @eval using MechGlueDiffEqBase
           #  @info "RecursiveArrayTools => using MechGlueRecursiveArrayTools"
           #  @eval using MechGlueRecursiveArrayTools
        end
        @require DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e" begin
            @info "DiffEqBase => using MechGlueDiffEqBase"
            @eval using MechGlueDiffEqBase
        end
    end
    @info "MechGluecode init"
end
end


text:
@testset "kernels" begin
    rng = StableRNG(25)
    v = [
        [3.0, 4.0, 5.0],
        [2.0, 5.0, 5.0],
        [1.0, 4.0, 5.0],
        [2.0, 2.0, 5.0],
        [1.0, 1.0, 3.0],
    ]

    v1 = [sort!(rand(rng, 1.0:5.0, 3)) for _ in 1:5]
    dists_matrix = NeighDistsMatrix(v)

    # User Defined Kernel (UDK)
    kernel  = UDK(func=identity, sort=true)
    @test sort_idxs(kernel) == true
    @test get_weights(kernel, dists_matrix) == dists_matrix
    kernel  = UDK(func=x->rand(rng, 3, 5), sort=true)
    @test_throws DimensionMismatch get_weights(kernel, dists_matrix)
    kernel  = UDK(func=x->rand(rng, "String", 5, 3), sort=true)
    @test_throws Exception  get_weights(kernel, dists_matrix)

    # Uniform Kernel (Uniform)
    kernel = Uniform()
    @test sort_idxs(kernel) == false
    @test get_weights(kernel, dists_matrix) === nothing

    # Dudani
    kernel = Dudani()
    @test sort_idxs(kernel) == true
    @test get_weights(
        kernel, dists_matrix
    ) == [1.0 0.5 0.0;1.0 0.0 0.0;1.0 0.25 0.0;1.0 1.0 0.0;1.0 1.0 0.0]

    # Inverse
    kernel = Inverse()
    @test sort_idxs(kernel) == false
    dists_matrix2 = NeighDistsMatrix(deepcopy(v))
    dists_matrix2[1,1] = 0
    dists_matrix2[3,3] = 0
    @test get_weights(
        kernel, dists_matrix2
    ) == [1.0 0.0 0.0;0.5 0.2 0.2;0.0 0.0 1.0;0.5 0.5 0.2;1.0 1.0 1/3]

    # Inverse Squared (ISquared)
    kernel = ISquared()
    @test sort_idxs(kernel) == false
    dists_matrix3 = NeighDistsMatrix(deepcopy(v))
    dists_matrix3[1,1] = 0
    dists_matrix3[3,3] = 0
    @test get_weights(
        kernel, dists_matrix3
    ) == [1.0 0.0 0.0;0.25 0.04 0.04;0.0 0.0 1.0;0.25 0.25 0.04;1.0 1.0 1/9]

    # Rank
    kernel = Rank()
    @test sort_idxs(kernel) == true
    @test get_weights(
        kernel, dists_matrix
    ) == repeat(collect(size(dists_matrix, 2):-1.0:1.0)', size(dists_matrix, 1))

    # Macleod
    kernel = Macleod(;a=1.5)
    @test sort_idxs(kernel) == true
    @test get_weights(
        kernel, dists_matrix
    ) == [4.0 3.2 2.4;9.0   5.4  5.4;16.0 11.2 9.6;9.0 9.0 5.4;4.0 4.0 2.4]

    # Zavreal
    kernel = Zavreal(a=0.5, b=1.5)
    @test sort_idxs(kernel) == false
    @test get_weights(
        kernel, dists_matrix
    ) == exp.(kernel.a .* dists_matrix .^ kernel.b)

    # ReciprocalRank
    kernel = ReciprocalRank()
    @test sort_idxs(kernel) == true
    @test get_weights(
        kernel, dists_matrix
    ) == repeat(1 ./ collect(1:size(dists_matrix, 2))', size(dists_matrix, 1))

    # DualU
    kernel = DualU()
    @test sort_idxs(kernel) == true
    @test get_weights(
        kernel, dists_matrix
    ) == [1.0 0.5 0.0;1.0 0.0 0.0;1.0 0.25 0.0;1.0 1.0 0.0;1.0 1.0 0.0]
    # DualD
    kernel = DualD()
    @test sort_idxs(kernel) == true
    @test get_weights(
        kernel, dists_matrix
    ) == [1.0 0.5 0.0;1.0 0.0 0.0;1.0 0.25 0.0;1.0 1.0 0.0;1.0 1.0 0.0]
    # Fibonacci
    kernel = Fibonacci()
    @test sort_idxs(kernel) == true
    @test get_weights(
        kernel, dists_matrix
    ) == repeat([2 1 1], size(dists_matrix, 1))
end


text:
using Oceananigans.Advection: AbstractAdvectionScheme
using Oceananigans.Operators: ℑxᶠᵃᵃ, ℑxᶜᵃᵃ, ℑyᵃᶠᵃ, ℑyᵃᶜᵃ, ℑzᵃᵃᶠ, ℑzᵃᵃᶜ
using Oceananigans.TurbulenceClosures: AbstractTurbulenceClosure, AbstractTimeDiscretization

const ATC = AbstractTurbulenceClosure
const ATD = AbstractTimeDiscretization

"""
    conditional_flux(i, j, k, ibg::IBG, ℓx, ℓy, ℓz, qᴮ, qᴵ, nc)

Return either

    i) The boundary flux `qᴮ` if the node condition `nc` is true (default: `nc = peripheral_node`), or
    ii) The interior flux `qᴵ` otherwise.

This can be used either to condition intrinsic flux functions, or immersed boundary flux functions.
"""
@inline conditional_flux(i, j, k, ibg, ℓx, ℓy, ℓz, qᴮ, qᴵ) =
    ifelse(immersed_peripheral_node(ℓx, ℓy, ℓz, i, j, k, ibg), qᴮ, qᴵ)

# Conveniences
@inline conditional_flux_ccc(i, j, k, ibg::IBG, qᴮ, qᴵ) = conditional_flux(i, j, k, ibg, c, c, c, qᴮ, qᴵ)
@inline conditional_flux_ffc(i, j, k, ibg::IBG, qᴮ, qᴵ) = conditional_flux(i, j, k, ibg, f, f, c, qᴮ, qᴵ)
@inline conditional_flux_fcf(i, j, k, ibg::IBG, qᴮ, qᴵ) = conditional_flux(i, j, k, ibg, f, c, f, qᴮ, qᴵ)
@inline conditional_flux_cff(i, j, k, ibg::IBG, qᴮ, qᴵ) = conditional_flux(i, j, k, ibg, c, f, f, qᴮ, qᴵ)

@inline conditional_flux_fcc(i, j, k, ibg::IBG, qᴮ, qᴵ) = conditional_flux(i, j, k, ibg, f, c, c, qᴮ, qᴵ)
@inline conditional_flux_cfc(i, j, k, ibg::IBG, qᴮ, qᴵ) = conditional_flux(i, j, k, ibg, c, f, c, qᴮ, qᴵ)
@inline conditional_flux_ccf(i, j, k, ibg::IBG, qᴮ, qᴵ) = conditional_flux(i, j, k, ibg, c, c, f, qᴮ, qᴵ)

#####
##### Diffusive fluxes
#####

# ccc, ffc, fcf
@inline _viscous_flux_ux(i, j, k, ibg::GFIBG, args...) = conditional_flux_ccc(i, j, k, ibg, zero(eltype(ibg)), viscous_flux_ux(i, j, k, ibg, args...))
@inline _viscous_flux_uy(i, j, k, ibg::GFIBG, args...) = conditional_flux_ffc(i, j, k, ibg, zero(eltype(ibg)), viscous_flux_uy(i, j, k, ibg, args...))
@inline _viscous_flux_uz(i, j, k, ibg::GFIBG, args...) = conditional_flux_fcf(i, j, k, ibg, zero(eltype(ibg)), viscous_flux_uz(i, j, k, ibg, args...))

 # ffc, ccc, cff
@inline _viscous_flux_vx(i, j, k, ibg::GFIBG, args...) = conditional_flux_ffc(i, j, k, ibg, zero(eltype(ibg)), viscous_flux_vx(i, j, k, ibg, args...))
@inline _viscous_flux_vy(i, j, k, ibg::GFIBG, args...) = conditional_flux_ccc(i, j, k, ibg, zero(eltype(ibg)), viscous_flux_vy(i, j, k, ibg, args...))
@inline _viscous_flux_vz(i, j, k, ibg::GFIBG, args...) = conditional_flux_cff(i, j, k, ibg, zero(eltype(ibg)), viscous_flux_vz(i, j, k, ibg, args...))

 # fcf, cff, ccc
@inline _viscous_flux_wx(i, j, k, ibg::GFIBG, args...) = conditional_flux_fcf(i, j, k, ibg, zero(eltype(ibg)), viscous_flux_wx(i, j, k, ibg, args...))
@inline _viscous_flux_wy(i, j, k, ibg::GFIBG, args...) = conditional_flux_cff(i, j, k, ibg, zero(eltype(ibg)), viscous_flux_wy(i, j, k, ibg, args...))
@inline _viscous_flux_wz(i, j, k, ibg::GFIBG, args...) = conditional_flux_ccc(i, j, k, ibg, zero(eltype(ibg)), viscous_flux_wz(i, j, k, ibg, args...))

# fcc, cfc, ccf
@inline _diffusive_flux_x(i, j, k, ibg::GFIBG, args...) = conditional_flux_fcc(i, j, k, ibg, zero(eltype(ibg)), diffusive_flux_x(i, j, k, ibg, args...))
@inline _diffusive_flux_y(i, j, k, ibg::GFIBG, args...) = conditional_flux_cfc(i, j, k, ibg, zero(eltype(ibg)), diffusive_flux_y(i, j, k, ibg, args...))
@inline _diffusive_flux_z(i, j, k, ibg::GFIBG, args...) = conditional_flux_ccf(i, j, k, ibg, zero(eltype(ibg)), diffusive_flux_z(i, j, k, ibg, args...))

#####
##### Advective fluxes
#####

# dx(uu), dy(vu), dz(wu)
# ccc,    ffc,    fcf
@inline _advective_momentum_flux_Uu(i, j, k, ibg::GFIBG, args...) = conditional_flux_ccc(i, j, k, ibg, zero(eltype(ibg)), advective_momentum_flux_Uu(i, j, k, ibg, args...))
@inline _advective_momentum_flux_Vu(i, j, k, ibg::GFIBG, args...) = conditional_flux_ffc(i, j, k, ibg, zero(eltype(ibg)), advective_momentum_flux_Vu(i, j, k, ibg, args...))
@inline _advective_momentum_flux_Wu(i, j, k, ibg::GFIBG, args...) = conditional_flux_fcf(i, j, k, ibg, zero(eltype(ibg)), advective_momentum_flux_Wu(i, j, k, ibg, args...))

# dx(uv), dy(vv), dz(wv)
# ffc,    ccc,    cff
@inline _advective_momentum_flux_Uv(i, j, k, ibg::GFIBG, args...) = conditional_flux_ffc(i, j, k, ibg, zero(eltype(ibg)), advective_momentum_flux_Uv(i, j, k, ibg, args...))
@inline _advective_momentum_flux_Vv(i, j, k, ibg::GFIBG, args...) = conditional_flux_ccc(i, j, k, ibg, zero(eltype(ibg)), advective_momentum_flux_Vv(i, j, k, ibg, args...))
@inline _advective_momentum_flux_Wv(i, j, k, ibg::GFIBG, args...) = conditional_flux_cff(i, j, k, ibg, zero(eltype(ibg)), advective_momentum_flux_Wv(i, j, k, ibg, args...))

# dx(uw), dy(vw), dz(ww)
# fcf,    cff,    ccc
@inline _advective_momentum_flux_Uw(i, j, k, ibg::GFIBG, args...) = conditional_flux_fcf(i, j, k, ibg, zero(eltype(ibg)), advective_momentum_flux_Uw(i, j, k, ibg, args...))
@inline _advective_momentum_flux_Vw(i, j, k, ibg::GFIBG, args...) = conditional_flux_cff(i, j, k, ibg, zero(eltype(ibg)), advective_momentum_flux_Vw(i, j, k, ibg, args...))
@inline _advective_momentum_flux_Ww(i, j, k, ibg::GFIBG, args...) = conditional_flux_ccc(i, j, k, ibg, zero(eltype(ibg)), advective_momentum_flux_Ww(i, j, k, ibg, args...))

@inline _advective_tracer_flux_x(i, j, k, ibg::GFIBG, args...) = conditional_flux_fcc(i, j, k, ibg, zero(eltype(ibg)), advective_tracer_flux_x(i, j, ibg, args...))
@inline _advective_tracer_flux_y(i, j, k, ibg::GFIBG, args...) = conditional_flux_cfc(i, j, k, ibg, zero(eltype(ibg)), advective_tracer_flux_y(i, j, ibg, args...))
@inline _advective_tracer_flux_z(i, j, k, ibg::GFIBG, args...) = conditional_flux_ccf(i, j, k, ibg, zero(eltype(ibg)), advective_tracer_flux_z(i, j, ibg, args...))

#####
##### "Boundary-aware" reconstruct
#####
##### Don't reconstruct with immersed cells!
#####

@inline near_x_boundary(i, j, k, ibg, ::AbstractAdvectionScheme{0}) = false
@inline near_y_boundary(i, j, k, ibg, ::AbstractAdvectionScheme{0}) = false
@inline near_z_boundary(i, j, k, ibg, ::AbstractAdvectionScheme{0}) = false

@inline near_x_boundary(i, j, k, ibg, ::AbstractAdvectionScheme{1}) = inactive_cell(i-1, j, k, ibg) | inactive_cell(i, j, k, ibg) | inactive_cell(i+1, j, k, ibg)
@inline near_y_boundary(i, j, k, ibg, ::AbstractAdvectionScheme{1}) = inactive_cell(i, j-1, k, ibg) | inactive_cell(i, j, k, ibg) | inactive_cell(i, j+1, k, ibg)
@inline near_z_boundary(i, j, k, ibg, ::AbstractAdvectionScheme{1}) = inactive_cell(i, j, k-1, ibg) | inactive_cell(i, j, k, ibg) | inactive_cell(i, j, k+1, ibg)

@inline near_x_boundary(i, j, k, ibg, ::AbstractAdvectionScheme{2}) = inactive_cell(i-2, j, k, ibg) | inactive_cell(i-1, j, k, ibg) | inactive_cell(i, j, k, ibg) | inactive_cell(i+1, j, k, ibg) | inactive_cell(i+2, j, k, ibg)
@inline near_y_boundary(i, j, k, ibg, ::AbstractAdvectionScheme{2}) = inactive_cell(i, j-2, k, ibg) | inactive_cell(i, j-1, k, ibg) | inactive_cell(i, j, k, ibg) | inactive_cell(i, j+1, k, ibg) | inactive_cell(i, j+2, k, ibg)
@inline near_z_boundary(i, j, k, ibg, ::AbstractAdvectionScheme{2}) = inactive_cell(i, j, k-2, ibg) | inactive_cell(i, j, k-1, ibg) | inactive_cell(i, j, k, ibg) | inactive_cell(i, j, k+1, ibg) | inactive_cell(i, j, k+2, ibg)

using Oceananigans.Advection: WENOVectorInvariantVel, VorticityStencil, VelocityStencil

@inline function near_horizontal_boundary_x(i, j, k, ibg, scheme::WENOVectorInvariantVel)
    return inactive_node(c, f, c, i, j, k, ibg)     |
           inactive_node(c, f, c, i-3, j, k, ibg)   | inactive_node(c, f, c, i+3, j, k, ibg) |
           inactive_node(c, f, c, i-2, j, k, ibg)   | inactive_node(c, f, c, i+2, j, k, ibg) |
           inactive_node(c, f, c, i-1, j, k, ibg)   | inactive_node(c, f, c, i+1, j, k, ibg) |
           inactive_node(f, c, c, i, j, k, ibg)     |
           inactive_node(f, c, c, i-2, j, k, ibg)   | inactive_node(f, c, c, i+2, j, k, ibg) |
           inactive_node(f, c, c, i-1, j, k, ibg)   | inactive_node(f, c, c, i+1, j, k, ibg) |
           inactive_node(f, c, c, i-2, j+1, k, ibg) | inactive_node(f, c, c, i+2, j+1, k, ibg) |
           inactive_node(f, c, c, i-1, j+1, k, ibg) | inactive_node(f, c, c, i+1, j+1, k, ibg) |
           inactive_node(f, c, c, i, j+1, k, ibg)
end

@inline function near_horizontal_boundary_y(i, j, k, ibg, scheme::WENOVectorInvariantVel)
    return inactive_node(f, c, c, i, j, k, ibg)     |
           inactive_node(f, c, c, i, j+3, k, ibg)   | inactive_node(f, c, c, i, j+3, k, ibg) |
           inactive_node(f, c, c, i, j+2, k, ibg)   | inactive_node(f, c, c, i, j+2, k, ibg) |
           inactive_node(f, c, c, i, j+1, k, ibg)   | inactive_node(f, c, c, i, j+1, k, ibg) |
           inactive_node(c, f, c, i, j, k, ibg)     |
           inactive_node(c, f, c, i, j-2, k, ibg)   | inactive_node(c, f, c, i, j+2, k, ibg) |
           inactive_node(c, f, c, i, j-1, k, ibg)   | inactive_node(c, f, c, i, j+1, k, ibg) |
           inactive_node(c, f, c, i+1, j-2, k, ibg) | inactive_node(c, f, c, i+1, j+2, k, ibg) |
           inactive_node(c, f, c, i+1, j-1, k, ibg) | inactive_node(c, f, c, i+1, j+1, k, ibg) |
           inactive_node(c, f, c, i+1, j, k, ibg)
end

# Takes forever to compile, but works.
# @inline near_x_boundary(i, j, k, ibg, ::AbstractAdvectionScheme{buffer}) where buffer = any(ntuple(δ -> inactive_node(i - buffer - 1 + δ, j, k, ibg), Val(2buffer + 1)))
# @inline near_y_boundary(i, j, k, ibg, ::AbstractAdvectionScheme{buffer}) where buffer = any(ntuple(δ -> inactive_node(i, j - buffer - 1 + δ, k, ibg), Val(2buffer + 1)))
# @inline near_z_boundary(i, j, k, ibg, ::AbstractAdvectionScheme{buffer}) where buffer = any(ntuple(δ -> inactive_node(i, j, k - buffer - 1 + δ, ibg), Val(2buffer + 1)))

for bias in (:symmetric, :left_biased, :right_biased)
    for (d, ξ) in enumerate((:x, :y, :z))

        code = [:ᵃ, :ᵃ, :ᵃ]

        for loc in (:ᶜ, :ᶠ)
            code[d] = loc
            second_order_interp = Symbol(:ℑ, ξ, code...)
            interp = Symbol(bias, :_interpolate_, ξ, code...)
            alt_interp = Symbol(:_, interp)

            near_boundary = Symbol(:near_, ξ, :_boundary)

            # Conditional high-order interpolation in Bounded directions
            @eval begin
                import Oceananigans.Advection: $alt_interp
                using Oceananigans.Advection: $interp

                @inline $alt_interp(i, j, k, ibg::ImmersedBoundaryGrid, scheme, ψ) =
                    ifelse($near_boundary(i, j, k, ibg, scheme),
                           $second_order_interp(i, j, k, ibg.underlying_grid, ψ),
                           $interp(i, j, k, ibg.underlying_grid, scheme, ψ))
            end
            if ξ == :z
                @eval begin
                    import Oceananigans.Advection: $alt_interp
                    using Oceananigans.Advection: $interp

                    @inline $alt_interp(i, j, k, ibg::ImmersedBoundaryGrid, scheme::WENOVectorInvariant, ∂z, VI, u) =
                        ifelse($near_boundary(i, j, k, ibg, scheme),
                            $second_order_interp(i, j, k, ibg.underlying_grid, ∂z, u),
                            $interp(i, j, k, ibg.underlying_grid, scheme, ∂z, VI, u))
                end
            else
                @eval begin
                    import Oceananigans.Advection: $alt_interp
                    using Oceananigans.Advection: $interp

                    @inline $alt_interp(i, j, k, ibg::ImmersedBoundaryGrid, scheme::WENOVectorInvariant, ζ, VI, u, v) =
                        ifelse($near_boundary(i, j, k, ibg, scheme),
                                $second_order_interp(i, j, k, ibg.underlying_grid, ζ, u, v),
                                $interp(i, j, k, ibg.underlying_grid, scheme, ζ, VI, u, v))
                end
            end
        end
    end
end

for bias in (:left_biased, :right_biased)
    for (d, dir) in zip((:x, :y), (:xᶜᵃᵃ, :yᵃᶜᵃ))
        interp     = Symbol(bias, :_interpolate_, dir)
        alt_interp = Symbol(:_, interp)

        near_horizontal_boundary = Symbol(:near_horizontal_boundary_, d)
        @eval begin
            @inline $alt_interp(i, j, k, ibg::ImmersedBoundaryGrid, scheme::WENOVectorInvariantVel, ζ, ::Type{VelocityStencil}, u, v) =
            ifelse($near_horizontal_boundary(i, j, k, ibg, scheme),
               $alt_interp(i, j, k, ibg, scheme, ζ, VorticityStencil, u, v),
               $interp(i, j, k, ibg.underlying_grid, scheme, ζ, VelocityStencil, u, v))
        end
    end
end


text:
function construct_ohp_curve(nturn,pitch,height,gap,ds,x0,y0,flipx,flipy,angle)
    rad = 0.5*pitch
    len = height
    x, y = Float64[], Float64[]
    xf, yf = 0.0, -gap-rad
    for i in 1:nturn
        xi, yi, xf, yf = _channel_pair(pitch,len,ds,xf,yf,false,false,0.0)
        append!(x,xi)
        append!(y,yi)
    end
    xi, yi, xf, yf = _closure(pitch,len,gap,2*nturn*pitch,ds,xf,yf,false,false,0.0)
    append!(x,xi)
    append!(y,yi)
    push!(x,xf)
    push!(y,yf)
    _transform!(x,y,x0,y0,flipx,flipy,angle)
    x[1:end-1], y[1:end-1], x[end], y[end]
end

function _transform!(x,y,x0,y0,flipx,flipy,angle)
    xtmp, ytmp = copy(x), copy(y)
    if flipx
        ytmp .= -ytmp
    end
    if flipy
        xtmp .= -xtmp
    end
    ca, sa = cos(angle), sin(angle)
    x .= xtmp.*ca .- ytmp.*sa
    y .= xtmp.*sa .+ ytmp.*ca
    x .+= x0
    y .+= y0
    x, y
 end

 function _channel_pair(pitch,len,ds,x0,y0,flipx,flipy,angle)
     rad = 0.5*pitch
     x, y = Float64[], Float64[]
     xi, yi, xf, yf = _line(len,ds,0.0,0.0,false,false,-π/2)
     append!(x,xi)
     append!(y,yi)
     xi, yi, xf, yf = _halfturn(rad,ds,xf,yf,true,true,0.0)
     append!(x,xi)
     append!(y,yi)
     xi, yi, xf, yf = _line(len,ds,xf,yf,false,false,π/2)
     append!(x,xi)
     append!(y,yi)
     xi, yi, xf, yf = _halfturn(rad,ds,xf,yf,false,true,0.0)
     append!(x,xi)
     append!(y,yi)
     push!(x,xf)
     push!(y,yf)
     _transform!(x,y,x0,y0,flipx,flipy,angle)
     x[1:end-1], y[1:end-1], x[end], y[end]
 end

 function _closure(pitch,channellen,closuregap,closurelen,ds,x0,y0,flipx,flipy,angle)
    rad = 0.5*pitch
    x, y = Float64[], Float64[]
    xi, yi, xf, yf = _line(channellen,ds,0.0,0.0,false,false,-π/2)
    append!(x,xi)
    append!(y,yi)
    xi, yi, xf, yf = _halfturn(rad,ds,xf,yf,true,true,0.0)
    append!(x,xi)
    append!(y,yi)
    xi, yi, xf, yf = _line(channellen+closuregap,ds,xf,yf,false,false,π/2)
    append!(x,xi)
    append!(y,yi)
    xi, yi, xf, yf = _quarterturn(rad,ds,xf,yf,false,false,0.0)
    append!(x,xi)
    append!(y,yi)
    xi, yi, xf, yf = _line(closurelen,ds,xf,yf,false,false,π)
    append!(x,xi)
    append!(y,yi)
    xi, yi, xf, yf = _quarterturn(rad,ds,xf,yf,false,false,π/2)
    append!(x,xi)
    append!(y,yi)
    xi, yi, xf, yf = _line(closuregap,ds,xf,yf,false,false,-π/2)
    append!(x,xi)
    append!(y,yi)
    push!(x,xf)
    push!(y,yf)
    _transform!(x,y,x0,y0,flipx,flipy,angle)
    x[1:end-1], y[1:end-1], x[end], y[end]
end

####  Basic elements  ####

# in each element, adjust the spacing so that there is an integral number of points, uniformly spaced
# along entire shape. Sets the first point at the origin.

#=
create a set of points in an arc of radius `radius` and turning angle `turnangle`,
with a nominal spacing `ds` at (0,0) and proceeding counterclockwise about the center
at (-radius,0). Then the initial point is moved to (x0,y0). If either `flipx` and
`flipy` are set to true, then it flips the curve about that axis. The curve is
rotated by angle `angle`.
=#
function _turn(radius,turnangle,ds,x0,y0,flipx,flipy,angle)
    np1 = ceil(Int,turnangle*radius/ds)
    θ = range(0.0,turnangle,length=np1)
    x, y = Float64[], Float64[]
    for θi in θ
        push!(x,radius*cos(θi)-radius)
        push!(y,radius*sin(θi))
    end
    _transform!(x,y,x0,y0,flipx,flipy,angle)
    x[1:end-1], y[1:end-1], x[end], y[end]
end

_halfturn(radius,ds,x0,y0,flipx,flipy,angle) = _turn(radius,π,ds,x0,y0,flipx,flipy,angle)

_quarterturn(radius,ds,x0,y0,flipx,flipy,angle) = _turn(radius,π/2,ds,x0,y0,flipx,flipy,angle)

#=
create a set of points in a straight line of length `length` with nominal spacing `ds`
at (0,0) and proceeding in the x direction.
The initial point is moved to (x0,y0). If either `flipx` and `flipy` are set to
true, then it flips the curve about that axis. The curve is rotated by angle `angle`.
=#

function _line(length,ds,x0,y0,flipx,flipy,angle)
    np1 = ceil(Int,length/ds)
    s = range(0.0,length,length=np1)
    x, y = Float64[], Float64[]
    for si in s
        push!(x,si)
        push!(y,0.0)
    end
    _transform!(x,y,x0,y0,flipx,flipy,angle)
    x[1:end-1], y[1:end-1], x[end], y[end]
end


text:
# =====================
# Utility functions
# =====================

# given an nxn matrix, return an (n+m)×(n+m) matrix with m new zero rows
# and m new zero columns
function _add_dimension(A::AbstractMatrix, m=1)
    n = size(A, 1)
    return vcat(hcat(A, zeros(n, m)), zeros(m, n+m))
end


text:
using HTTP
using Genie.Renderer.Json
import JSON
using Test

host = "localhost"
port = 8001

# Test geniews.jl
function test_learn_interactive_environment()
    r = HTTP.request("GET", "http://localhost:$port/hello/world")
    @test "Hello World!" == String(r.body)

    r = HTTP.request("GET", "http://localhost:$port/echo/ciao")
    @test "ciao" == String(r.body)

    r = HTTP.request("GET", "http://localhost:$port/sum/3/2")
    @test "5" == String(r.body)

    r = HTTP.request("GET", "http://localhost:$port/sum/3/2?initial_value=10")
    @test "15" == String(r.body)

    r = HTTP.request("GET", "http://localhost:$port/hello.html")
    @test "Hello World" == String(r.body)

    r = HTTP.request("GET", "http://localhost:$port/hello.json")
    @test "\"Hello World\"" == String(r.body)

    r = HTTP.request("GET", "http://localhost:$port/hello.txt")
    @test "Hello World" == String(r.body)
end

function test_learn_simple_api_backend()
    r = HTTP.request("GET", "http://$host:$port/")
    j = JSON.parse(IOBuffer(r.body))
    #JSON.print(stdout, j, 4)
    @test "Hi there!" == j["message"]

    #response = HTTP.request("POST", "http://$host:$port/echo", [("Content-Type", "application/json")], """{"message":"hello", "repeat":3}""")

    payload = Dict(
        "message" => "hello",
        "repeat" => 3,
    )
    response = HTTP.request(
        "POST",
        "http://$host:$port/echo",
        [("Content-Type", "application/json")],
        JSON.json(payload)
    )

    j = JSON.parse(IOBuffer(response.body))
    JSON.print(stdout, j, 4)
end

#test_learn_interactive_environment()
test_learn_simple_api_backend()

nothing


text:
"""
Calculate initial vortex condition

κ: vortex strength
μ: decay rate
rc: critical radius

"""
function vortex_ic!(
    prim,
    γ,
    x,
    y,
    κ = 0.25,
    μ = 0.204,
    rc = 0.05;
    x0 = 0.8,
    y0 = 0.5,
)
    T0 = 1 / prim[end]
    s = prim[1]^(1-γ) / (2 * prim[end])

    r = sqrt((x - x0)^2 + (y - y0)^2)
    η = r / rc

    δu = κ * η * exp(μ * (1-η^2)) * (y - y0) / r
    δv = -κ * η * exp(μ * (1-η^2)) * (x - x0) / r
    δT = -(γ-1)*κ^2/(8*μ*γ)*exp(2*μ*(1-η^2))

    ρ = prim[1]^(γ-1) * (T0+δT) / T0^(1/(γ-1))
    prim .= [ρ, prim[2]+δu, prim[3]+δv, 1/(1/prim[4]+δT)]

    return nothing
end

max_stars_repo_name: CSMMLab/FRSG

max_issues_repo_name: CSMMLab/FRSG

max_forks_repo_name: CSMMLab/FRSG


text:
module MaximumLikelihoodProblems

include("public.jl")

include("types.jl")

include("fit.jl")
include("linear_algebra.jl")
include("loglikelihood.jl")

include("Internal.jl")

end # module


text:
module LoadCell
export HX711, __init__, read, read_avg, get_grams, tare, power_down, power_up
using Pkg
using PyCall
GPIO = pyimport("RPi.GPIO")

#ENV["PYTHON"] = raw"/usr/bin/python3.6"
#Pkg.build("PyCall")

mutable struct HX711
	output::Int64
	pwr_clock::Int64

	gain::Int64
	offset::Int64
	scale::Int64
	power::Bool
end

# output, pwr_clock, gain, offset, scale

function __init__(output::Int64, pwr_clock::Int64, gain::Int64 = 128)
	#output: Serial Data Output pin
	#pwr_clock: Power Down and Serial Clock Input pin
	#gain: set gain 128, 64, 32
	@assert (sum(gain .== [32,64,128]) == 1) "gain must be 32,64, or 128"

	offset = 0
    scale = 1
	gain = log2(gain)
	power = true

	GPIO.setmode(GPIO.BCM)
	GPIO.setup(pwr_clock, GPIO.OUT)
	GPIO.setup(output, GPIO.IN)

	self = HX711(output, pwr_clock, gain, offset, scale, power)
	power_up(self)
	return self
end

function read(self::HX711)

	if GPIO.input(self.output) != 0
		error("NO INPUT"); end

	count = 0

	for i in 1:24
		GPIO.output(self.pwr_clock, GPIO.HIGH)
		count = count << 1
		GPIO.output(self.pwr_clock, GPIO.LOW)

		if GPIO.input(self.output)
			count += 1
		end
	end

	GPIO.output(self.pwr_clock, GPIO.HIGH)
	count = count ^ 0x800000
	GPIO.output(self.pwr_clock, GPIO.LOW)

	# set channel and gain factor for next reading
	for i in 1:self.gain
		GPIO.output(self.pwr_clock, GPIO.HIGH)
		GPIO.output(self.pwr_clock, GPIO.LOW)
	end

	return count
end

function read_avg(self::HX711, times::Int = 16)
	sum = 0
	for i in times; sum += read(self); end
	return sum / times
end

function get_grams(self::HX711, times::Int = 16)
	#high number of times will have a slower runtime speed
	value = read_avg(self, times) - self.offset
	return  value / self.scale
end

function tare(self::HX711, times::Int = 16)
	self.offset = read_avg(self, times)
end

function power_down(self::HX711)
	self.power = false
	GPIO.output(self.pwr_clock, GPIO.LOW)
	GPIO.output(self.pwr_clock, GPIO.HIGH)
end

function power_up(self::HX711)
	self.power = true
	GPIO.output(self.pwr_clock, GPIO.LOW)
end
end # module


text:
function calc_evals( LF, Potentials, evecs )
    Ncols   = size(evecs)[2]
    Npoints = size(evecs)[1]
    ΔV = LF.LFx.h * LF.LFy.h * LF.LFz.h
    Hred = zeros(Float64,Ncols,Ncols)
    Hv = op_H( LF, Potentials, evecs )
    # Hred = evecs' * Hv
    BLAS.gemm!('T','N',1.0,evecs,Hv,0.0,Hred)
    #
    D, V = eig(Hred*ΔV)  # XXX deltaV is applied here
    return D
end

function calc_evals( LF, ∇2::SparseMatrixCSC{Float64,Int64}, Potentials, evecs )
    Ncols   = size(evecs)[2]
    Npoints = size(evecs)[1]
    ΔV = LF.LFx.h * LF.LFy.h * LF.LFz.h
    Hv = zeros(Float64,Npoints,Ncols)
    Hred = zeros(Float64,Ncols,Ncols)
    Hv = op_H( LF, ∇2, Potentials, evecs )
    # Hred = evecs' * Hv
    BLAS.gemm!('T','N',1.0,evecs,Hv,0.0,Hred)
    #
    D, V = eig(Hred*ΔV)
    return D
end


text:
function FindMarker(r)
    counter = 1
    while !ismarker(r)
        movementsbyfixmovesandcheckmarker(r, Sud, counter)
        if ismarker(r)
            break
        end
        movementsbyfixmovesandcheckmarker(r, Ost, counter)
        if ismarker(r)
            break
        end
        counter += 1
        movementsbyfixmovesandcheckmarker(r, Nord, counter)
        if ismarker(r)
            break
        end
        movementsbyfixmovesandcheckmarker(r, West, counter)
        if ismarker(r)
            break
        end
        counter += 1
    end
end

function movementsbyfixmovesandcheckmarker(r, side, moves)
    while moves != 0
        move!(r, side)
        if ismarker(r)
            break
        end
        moves -= 1
    end
end


text:
module TypeVarExtraction

using MLStyle
using MLStyle.Toolz.List: cons, nil
using MLStyle.Infras
using MLStyle.MatchCore
using MLStyle.Qualification

export TypeVar, Relation, ChainRelation
export any_constraint, to_tvar, extract_tvars

struct TypeVar
    t :: Symbol
end

struct Relation
    l :: Symbol
    op :: Symbol
    r
end

struct ChainRelation
    var :: Symbol
    lower
    super
end

function any_constraint(t, forall)

    function is_rel(::Relation)
        true
    end

    function is_rel(::ChainRelation)
        true
    end

    function is_rel(::TypeVar)
        false
    end

    !(t isa Symbol) || any(is_rel, collect(extract_tvars(forall)))
end

function extract_tvars(t :: AbstractArray)
    @match t begin
        [] => nil()
        [hd && if hd isa Symbol end, tl...] => cons(TypeVar(hd), extract_tvars(tl))
        [:($hd <: $r), tl...] =>  cons(Relation(hd, :<:, r), extract_tvars(tl))
        [:($hd >: $(r)), tl...] =>  cons(Relation(hd, Symbol(">:"), r), extract_tvars(tl))
        [:($lower <: $hd <: $super), tl...] ||
        [:($super >: $hd >: $lower), tl...] => cons(ChainRelation(hd, lower, super))
        _ => @syntax_err "invalid tvars($t)"
    end
end

to_tvar(t::TypeVar) = t.t
to_tvar(t::Relation) = t.l
to_tvar(t::ChainRelation) = t.var

end


text:
set_theme!(
    rowgap = 5,
    colgap = 5,
    fontsize = 16,
    backgroundcolor = "#212946",
    Axis = (
        xgridcolor = :white,
        ygridcolor = :white,
        backgroundcolor = "#212946",
        xgridstyle=:dash,
        ygridstyle=:dash,
        xgridwidth=0.35,
        ygridwidth=0.35,
        xtickalign=0,
        ytickalign=0,
        xtickcolor = :white,
        ytickcolor = :white,
        bottomspinecolor = :white,
        topspinecolor =:white,
        leftspinecolor =:white,
        rightspinecolor =:white,
        xlabelcolor=:white,
        ylabelcolor = :white,
        xticklabelcolor = :white,
        yticklabelcolor = :white,
        xgridvisible = true,
        ygridvisible = true),
    Colorbar = (
        topspinevisible = false,
        rightspinevisible = false,
        bottomspinevisible = false,
        leftspinevisible = false,
        width = 12,
        height = Relative(4/4),
        labelcolor = :white,
        ticklabelcolor = :white,
        tickalign = 1,
        tickcolor = :white,
        labelpadding = -5,
        ticklabelsize = 16),
    Legend = (
        bgcolor = ("#748AA6FF", 0.25),
        labelcolor = :white,
        titlecolor = :white,
        framecolor=:white)
)

max_stars_repo_name: aivodji/MakieNotebooks

max_issues_repo_name: aivodji/MakieNotebooks

max_forks_repo_name: aivodji/MakieNotebooks


text:
## Adapted from Protein3DAlign.jl (Tim Holy)
@enum NWTrace::UInt8 UNDET LEFT TOP DIAG

Base.show(io::IO, t::NWTrace) = show(io, t == UNDET ? 'U' :
                                         t == LEFT ? 'L' :
                                         t == TOP ? 'T' : 'D')
Base.show(io::IO, ::typeof(MIME("text/plain")), t::NWTrace) = show(io, t)

"""
    PairIterator(P, start)

Construct an iterator over the aligned pairs, visiting only `DIAG` elements of `P`.
`start` is always visited regardless of the value of `P`, and it is implicitly treated as
`DIAG`.

# Examples

```jldoctest
julia> using Protein3DAlign: UNDET, LEFT, TOP, DIAG

julia> P = [UNDET LEFT LEFT LEFT;
            TOP   LEFT LEFT LEFT;
            TOP   DIAG LEFT LEFT;
            TOP   TOP  LEFT UNDET];

julia> for pair in Protein3DAlign.PairIterator(P, (4, 4))
           @show pair
       end
pair = (4, 4)
pair = (3, 2)
```

This returned sequence arises from the following:

- the iterator starts at (4, 4) and returns it regardless of the value of `P`
- the iterator moves to (3, 3) (treating the starting point as if it were `DIAG`).
  Since the value of `P` here is `LEFT`, it moves on to (3, 2).
  This is `DIAG`, so it returns this pair.
- The iterator then progresses to (2, 1) `TOP` followed by (1, 1) `UNDET`, which causes termination.
  None of these items are returned.
"""
struct PairIterator{M<:AbstractMatrix{NWTrace}}
    P::M
    start::Tuple{Int,Int}
end

Base.IteratorSize(::Type{<:PairIterator}) = Base.SizeUnknown()
Base.eltype(::Type{<:PairIterator}) = Tuple{Int,Int}

Base.iterate(iter::PairIterator) = iterate(iter, iter.start)

function Base.iterate(iter::PairIterator, (i,j)::Tuple{Int,Int})
    item = (i, j)
    P = iter.P
    ax1f, ax2f = map(first, axes(P))
    i -= 1; j -= 1    # treat the current item as DIAG (on entry, it might not have been set yet)
    while i >= ax1f && j >= ax2f
        p = P[i, j]
        p == UNDET && break
        p == DIAG && return item, (i, j)
        if p == TOP
            i -= 1
        else # p == LEFT
            j -= 1
        end
    end
    return item, nothing
end
Base.iterate(iter::PairIterator, ::Nothing) = nothing

mutable struct DistanceNeedlemanWunsch{T}
    similarity::Matrix{T}
    score::Matrix{T}
    trace::Matrix{NWTrace}
    gapopen::T
    gapextend::T
end

DistanceNeedlemanWunsch(similarity, gapopen=-0.6, gapextend=0.0
    ) = DistanceNeedlemanWunsch(similarity,
                                 zeros(size(similarity).+1...),
                                 similar(zeros(size(similarity).+1...), NWTrace),
                                 gapopen,
                                 gapextend)

function fill!(nw::DistanceNeedlemanWunsch)
    similarity = nw.similarity
    score = nw.score
    trace = nw.trace
    gapopen, gapextend = nw.gapopen, nw.gapextend

    m,n = size(similarity)

    # First column & row
    for i in 1:m
        score[i+1,1] = score[i,1] + (i == 1 ? gapopen : gapextend)
        trace[i+1,1] = TOP
    end

    for j in 1:n
        score[1,j+1] = score[1,j] + (j == 1 ? gapopen : gapextend)
        trace[1,j+1] = LEFT
    end

    biasleft = m > n
    for j in 1:n, i in 1:m
        scorediag = score[i,j] + similarity[i,j]     # only distance between points is considered (penalized)
        scoreleft = score[i+1,j] + (i == m ? gapextend : gappenalty(gapopen, gapextend, trace[i+1,j], LEFT))
        scoretop  = score[i,j+1] + (j == n ? gapextend : gappenalty(gapopen, gapextend, trace[i,j+1], TOP))

        if scorediag > max(scoreleft, scoretop)
            score[i+1,j+1] = scorediag
            trace[i+1,j+1] = DIAG
        elseif scoreleft > scoretop || (scoreleft == scoretop && biasleft)
            score[i+1,j+1] = scoreleft
            trace[i+1,j+1] = LEFT
        else
            score[i+1,j+1] = scoretop
            trace[i+1,j+1] = TOP
        end
    end
end

gappenalty(gap_open::Real, gap_extend::Real, p, pnew) = p == pnew ? gap_extend : gap_open

function traceback(F::AbstractMatrix, P::AbstractMatrix, i=nothing, j=nothing)
    ax1, ax2 = axes(P)
    i === nothing && (i = last(ax1))
    j === nothing && (j = last(ax2))
    score = F[i,j]
    # Trace back
    matchedpairs = Tuple{Int,Int}[]
    while i > first(ax1) || j > first(ax2)
        p = P[i,j]
        if p == DIAG
            i -= 1
            j -= 1
            push!(matchedpairs, (i,j))
        elseif p == LEFT
            j -= 1
        else
            i -= 1
        end
    end
    reverse!(matchedpairs)
    return score, matchedpairs
end

traceback(nw::DistanceNeedlemanWunsch, i=nothing, j=nothing) = traceback(nw.score, nw.trace, i, j)


text:
#!/usr/bin/env julia --project
using ClimateMachine
ClimateMachine.init(parse_clargs = true)

using ClimateMachine.Atmos
using ClimateMachine.Orientations
using ClimateMachine.ConfigTypes
using ClimateMachine.NumericalFluxes
using ClimateMachine.Diagnostics
using ClimateMachine.GenericCallbacks
using ClimateMachine.ODESolvers
using ClimateMachine.TurbulenceClosures
using ClimateMachine.SystemSolvers: ManyColumnLU
using ClimateMachine.Mesh.Filters
using ClimateMachine.Mesh.Grids
using ClimateMachine.Mesh.Topologies
using ClimateMachine.Mesh.Interpolation
using Thermodynamics.TemperatureProfiles
using ClimateMachine.VariableTemplates
using Thermodynamics: air_density, total_energy

using LinearAlgebra
using StaticArrays

using CLIMAParameters
using CLIMAParameters.Planet: day, planet_radius
struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

function set_topofun(f, r_inner, r_outer, topography)
    function wrapper_topo(a, b, c)
        return f(
            a,
            b,
            c,
            max(abs(a), abs(b), abs(c));
            r_inner = r_inner,
            r_outer = r_outer,
            topography = topography,
        )
    end
    return wrapper_topo
end

function init_solid_body_rotation!(problem, bl, state, aux, localgeo, t)
    FT = eltype(state)

    # initial velocity profile (we need to transform the vector into the Cartesian
    # coordinate system)
    u_0::FT = 0
    u_sphere = SVector{3, FT}(u_0, 0, 0)
    u_init = sphr_to_cart_vec(bl.orientation, u_sphere, aux)
    e_kin::FT = 0.5 * sum(abs2.(u_init))

    # Assign state variables
    state.ρ = aux.ref_state.ρ
    state.ρu = u_init
    state.energy.ρe = aux.ref_state.ρe + state.ρ * e_kin

    nothing
end

function config_solid_body_rotation(FT, poly_order, resolution, ref_state)

    # Set up the atmosphere model
    exp_name = "SolidBodyRotation"
    domain_height::FT = 30e3 # distance between surface and top of atmosphere (m)

    _planet_radius = FT(planet_radius(param_set))

    physics = AtmosPhysics{FT}(
        param_set;
        ref_state = ref_state,
        turbulence = ConstantKinematicViscosity(FT(0)),
        #hyperdiffusion = DryBiharmonic(FT(8 * 3600)),
        moisture = DryModel(),
    )
    model = AtmosModel{FT}(
        AtmosGCMConfigType,
        physics;
        init_state_prognostic = init_solid_body_rotation!,
        source = (Gravity(), Coriolis()),
    )

    config = ClimateMachine.AtmosGCMConfiguration(
        exp_name,
        poly_order,
        resolution,
        domain_height,
        param_set,
        init_solid_body_rotation!;
        model = model,
        numerical_flux_first_order = RoeNumericalFlux(),
        meshwarp = set_topofun(
            cubed_sphere_topo_warp,                   ## Topography warp function
            _planet_radius,                       ## Domain inner radius
            _planet_radius + domain_height,       ## Domain outer radius
            DCMIPMountain(),                      ## Problem specific dispatch
        ),
    )

    return config
end

function main()
    # Driver configuration parameters
    FT = Float64                             # floating type precision
    poly_order = 5                           # discontinuous Galerkin polynomial order
    n_horz = 8                               # horizontal element number
    n_vert = 4                               # vertical element number
    timestart::FT = 0                        # start time (s)
    timeend::FT = 7200

    # Set up a reference state for linearization of equations
    temp_profile_ref =
        DecayingTemperatureProfile{FT}(param_set, FT(290), FT(220), FT(8e3))
    ref_state = HydrostaticState(temp_profile_ref)

    # Set up driver configuration
    driver_config =
        config_solid_body_rotation(FT, poly_order, (n_horz, n_vert), ref_state)

    # Set up experiment
    ode_solver_type = ClimateMachine.IMEXSolverType(
        implicit_model = AtmosAcousticGravityLinearModel,
        implicit_solver = ManyColumnLU,
        solver_method = ARK2GiraldoKellyConstantinescu,
        split_explicit_implicit = false,
        discrete_splitting = true,
    )

    CFL = FT(0.2) # target acoustic CFL number

    # time step is computed such that the horizontal acoustic Courant number is CFL
    solver_config = ClimateMachine.SolverConfiguration(
        timestart,
        timeend,
        driver_config,
        Courant_number = CFL,
        ode_solver_type = ode_solver_type,
        CFL_direction = HorizontalDirection(),
        diffdir = HorizontalDirection(),
    )

    # initialize using a different ref state (mega-hack)
    temp_profile_init =
        DecayingTemperatureProfile{FT}(param_set, FT(280), FT(230), FT(9e3))
    init_ref_state = HydrostaticState(temp_profile_init)

    init_driver_config = config_solid_body_rotation(
        FT,
        poly_order,
        (n_horz, n_vert),
        init_ref_state,
    )
    init_solver_config = ClimateMachine.SolverConfiguration(
        timestart,
        timeend,
        init_driver_config,
        Courant_number = CFL,
        ode_solver_type = ode_solver_type,
        CFL_direction = HorizontalDirection(),
        diffdir = HorizontalDirection(),
    )

    # initialization
    solver_config.Q .= init_solver_config.Q

    # Set up diagnostics
    dgn_config = config_diagnostics(FT, driver_config)

    # Set up user-defined callbacks
    filterorder = 20
    filter = ExponentialFilter(solver_config.dg.grid, 0, filterorder)
    cbfilter = GenericCallbacks.EveryXSimulationSteps(1) do
        Filters.apply!(
            solver_config.Q,
            AtmosFilterPerturbations(driver_config.bl),
            solver_config.dg.grid,
            filter,
            # filter perturbations from the initial state
            state_auxiliary = init_solver_config.dg.state_auxiliary,
        )
        nothing
    end

    # Run the model
    result = ClimateMachine.invoke!(
        solver_config;
        diagnostics_config = dgn_config,
        user_callbacks = (cbfilter,),
        check_euclidean_distance = false,
    )

    relative_error =
        norm(solver_config.Q .- init_solver_config.Q) /
        norm(init_solver_config.Q)
    @info "Relative error = $relative_error"
end

function config_diagnostics(FT, driver_config)
    interval = "0.5shours" # chosen to allow diagnostics every 30 simulated minutes

    _planet_radius = FT(planet_radius(param_set))

    info = driver_config.config_info
    boundaries = [
        FT(-90.0) FT(-180.0) _planet_radius
        FT(90.0) FT(180.0) FT(_planet_radius + info.domain_height)
    ]
    resolution = (FT(2), FT(2), FT(1000)) # in (deg, deg, m)
    interpol = ClimateMachine.InterpolationConfiguration(
        driver_config,
        boundaries,
        resolution,
    )

    dgngrp = setup_atmos_default_diagnostics(
        AtmosGCMConfigType(),
        interval,
        driver_config.name,
        interpol = interpol,
    )

    return ClimateMachine.DiagnosticsConfiguration([dgngrp])
end

main()


text:
using CSV
using SparseArrays

function load_points(filename)
    # ripser uses Float32
    table = CSV.File(filename; header=0, type=Float32, delim=' ')
    nrow = length(table)
    dim = length(table[1])
    result = Vector{NTuple{dim,Float32}}(undef, nrow)
    for i in 1:nrow
        result[i] = tuple(table[i]...)
    end
    return result
end

function load_dist(filename)
    table = CSV.File(filename; header=0, type=Float32, delim=' ')
    result = Matrix{Float32}(undef, (length(table), length(table)))
    for (i, r) in enumerate(table)
        result[i, :] .= r
    end
    return result
end

function load_sparse(filename)
    table = CSV.Rows(filename; header=0, type=Float32, delim=' ')
    I = getindex.(table, 1) .+ 1
    J = getindex.(table, 2) .+ 1
    V = getindex.(table, 3) .+ 1
    return sparse([I; J], [J; I], [V; V])
end

function load_dipha(filename)
    open(filename, "r") do f
        magic = read(f, Int)
        @assert magic == 8067171840
        type = read(f, Int)
        @assert type == 1
        len = read(f, Int)
        dim = read(f, Int)
        if dim == 2
            m = read(f, Int)
            n = read(f, Int)
            result = Array{Float64}(undef, (m, n))
            read!(f, result)

            return result
        elseif dim == 3
            m = read(f, Int)
            n = read(f, Int)
            o = read(f, Int)
            result = Array{Float64}(undef, (m, n, o))
            read!(f, result)

            return result
        else
            error("not implemented for dim=$dim")
        end
    end
end

function load_data(filename)
    ext = splitext(filename)[2]
    if ext == ".pts"
        return load_points(filename)
    elseif ext == ".dist"
        return load_dist(filename)
    elseif ext == ".spdist"
        return load_sparse(filename)
    elseif ext == ".dipha"
        return load_dipha(filename)
    else
        error("unsupported filetype \"$(ext)\"")
    end
end


text:
"""
	struct ExtractArray{T}
		item::T
	end

Convert array of values to a `Mill.BagNode` with items converted
by `item`. The entire array is assumed to be a single bag.

# Examples

```jldoctest
julia> ec = ExtractArray(ExtractCategorical(2:4));

julia> ec([2,3,1,4]).data
4×4 Mill.ArrayNode{Mill.MaybeHotMatrix{Union{Missing, UInt32}, UInt32, Union{Missing, Bool}}, Nothing}:
  true    ⋅      ⋅      ⋅
   ⋅     true    ⋅      ⋅
   ⋅      ⋅      ⋅     true
   ⋅      ⋅     true    ⋅

julia> es = ExtractArray(ExtractScalar());

julia> es([2,3,4])
BagNode \t# 1 obs, 80 bytes
  └── ArrayNode(1×3 Array with Union{Missing, Float32} elements) \t# 3 obs, 63 bytes

julia> es([2,3,4]).data
1×3 Mill.ArrayNode{Matrix{Union{Missing, Float32}}, Nothing}:
 2.0f0  3.0f0  4.0f0
```
"""
struct ExtractArray{T} <: BagExtractor
	item::T
end

extract_empty_bag_item(::BagExtractor, store_input) = @error "Abstract method, please, implement it for your extractor"
extract_empty_bag_item(s::ExtractArray, store_input) = s.item(extractempty; store_input)

"""
returns missing bag of 1 observation
"""
function extract_missing_bag(s::BagExtractor, v; store_input=false)
	Mill._emptyismissing[] && return _make_bag_node(missing, [0:-1], [v], store_input)
	ds = extract_empty_bag_item(s, store_input)
	_make_bag_node(ds, [0:-1], [v], store_input)
end

(s::ExtractArray)(v::MissingOrNothing; store_input=false) = extract_missing_bag(s, v; store_input)

(s::ExtractArray)(v::AbstractVector; store_input=false) =
    isempty(v) ?
	extract_missing_bag(s, v; store_input) :
	BagNode(mapreduce(x->s.item(x; store_input), catobs, v),[1:length(v)])
(s::ExtractArray)(v::HierarchicType; store_input=false) = extract_missing_bag(s, v; store_input)

(s::ExtractArray)(v::ExtractEmpty; store_input=false) = make_empty_bag(s.item(extractempty; store_input), v)

Base.hash(e::ExtractArray, h::UInt) = hash(e.item, h)
Base.:(==)(e1::ExtractArray, e2::ExtractArray) = e1.item == e2.item


text:
function johnsontrottermove!(ints, isleft)
    len = length(ints)
    function ismobile(pos)
        if isleft[pos] && (pos > 1) && (ints[pos-1] < ints[pos])
            return true
        elseif !isleft[pos] && (pos < len) && (ints[pos+1] < ints[pos])
            return true
        end
        false
    end
    function maxmobile()
        arr = [ints[pos] for pos in 1:len if ismobile(pos)]
        if isempty(arr)
            0, 0
        else
            maxmob = maximum(arr)
            maxmob, findfirst(x -> x == maxmob, ints)
        end
    end
    function directedswap(pos)
        tmp = ints[pos]
        tmpisleft = isleft[pos]
        if isleft[pos]
            ints[pos] = ints[pos-1]; ints[pos-1] = tmp
            isleft[pos] = isleft[pos-1]; isleft[pos-1] = tmpisleft
        else
            ints[pos] = ints[pos+1]; ints[pos+1] = tmp
            isleft[pos] = isleft[pos+1]; isleft[pos+1] = tmpisleft
        end
    end
    (moveint, movepos) = maxmobile()
    if movepos > 0
        directedswap(movepos)
        for (i, val) in enumerate(ints)
            if val > moveint
                isleft[i] = !isleft[i]
            end
        end
        ints, isleft, true
    else
        ints, isleft, false
    end
end
function johnsontrotter(low, high)
    ints = collect(low:high)
    isleft = [true for i in ints]
    firstconfig = copy(ints)
    iters = 0
    while true
        iters += 1
        println("$ints $(iters & 1 == 1 ? "+1" : "-1")")
        if johnsontrottermove!(ints, isleft)[3] == false
            break
        end
    end
    println("There were $iters iterations.")
end
johnsontrotter(1,4)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# ==================================
# Beam
#
# system type: continuous LTI system
# state dimension: 348
# input dimension: 1
# ==================================
using ReachabilityBenchmarks, MathematicalSystems, LazySets, MAT

function beam_model()
    file = matopen(@relpath "beam.mat")

    # system matrix
    A = read(file, "A")

    # input matrix
    B = read(file, "B")

    # state domain
    X = Universe(348)

    # input domain
    U = BallInf([0.5], 0.3)

    # continuous LTI system
    S = @system(x' = Ax + Bu, x ∈ X, u ∈ U)

    return S
end

max_stars_repo_name: JuliaReach/ReachabilityBenchmarks

max_issues_repo_name: JuliaReach/ReachabilityBenchmarks

max_forks_repo_name: JuliaReach/ReachabilityBenchmarks


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    for T in (Bool, Int, Float32, Float64)
        for A in (Vector, Matrix)
            precompile(stridedpointer, (A{T},))
        end
    end
    function precompile_nt(@nospecialize(T))
        for I ∈ (Int8, UInt8, Int16, UInt16, Int32, UInt32, Int64, UInt64)
            precompile(vload_quote, (Type{T}, Type{I}, Symbol, Int, Int, Int, Int, Bool, Bool))
            precompile(vload_quote, (Type{T}, Type{I}, Symbol, Int, Int, Int, Int, Bool, Bool, Bool, Bool))
        end
    end
    U = NativeTypes
    while isa(U, Union)
        T, U = U.a, U.b
        precompile_nt(T)
    end
    precompile_nt(U)
    precompile(_pick_vector_width, (Type, Vararg{Type,100}))
    if VERSION ≥ v"1.7.0-DEV.346"
        for T ∈ (Float32, Float64)
            W = pick_vector_width(T)
            precompile(>=, (Int, MM{W, 1, Int}))
            for op ∈ (-, Base.FastMath.sub_fast)
                precompile(op, (Vec{W, T}, ))
            end
            for op ∈ (+, -, *, Base.FastMath.add_fast, Base.FastMath.sub_fast, Base.FastMath.mul_fast)
                precompile(op, (Vec{W, T}, Vec{W, T}))
            end
            for op ∈ (VectorizationBase.vfmadd, VectorizationBase.vfmadd_fast)
                precompile(op, (Vec{W, T}, Vec{W, T}, Vec{W, T}))
            end
        end
    end
end


text:
using Test

@test true


text:
using Test
using ProBeTune
using Statistics
using DiffEqFlux

# test samplers
@testset "Samplers" begin
    @testset "Continuos sampler" begin
        f(rng, t_span, n) = t -> rand(rng, 1)

        s = ContinousInputSampler(f, (0, 10))
        @test s.tspan == (0, 10)

        A = [f(0.0) for f in s[1:10]]
        resample!(s)
        B = [f(0.0) for f in s[1:10]]
        @test A !== B

        reset!(s)
        C = [f(0.0) for f in s[1:10]]
        @test A == C
    end

    @testset "Constant sampler" begin
        f(rng, n) =rand(rng, 5)

        s = ConstantInputSampler(f)

        A = s[1:10]
        resample!(s)
        B = s[1:10]
        @test A !== B

        reset!(s)
        C = s[1:10]
        @test A == C
    end
end

# Test a simple example
@testset "simple_example" begin
    dim_sys = 10
    diffusive_nl_example = PBTExample.create_diffusive_nl_example(dim_sys, 2, 0.:0.1:10., 10);

    p_initial = ones(2 * 10 + dim_sys)

    d_init, p_tuned, dists = behavioural_distance(diffusive_nl_example, p_initial; abstol=1e-2, reltol=1e-2)

    @test typeof(p_tuned) == typeof(p_initial)
    @test p_initial != p_tuned # test if the optimizer actually changes the parameters, this way we can also indirectly check if AD works

    res_10 = pbt_tuning(diffusive_nl_example, p_tuned; abstol=1e-4, reltol=1e-4, optimizer = ADAM(0.5), optimizer_options = (:maxiters => 50,))

    @test res_10.retcode == :Default || res_10.retcode == :Success

    d_end, p_tuned, dists = behavioural_distance(diffusive_nl_example, res_10.u; abstol=1e-2, reltol=1e-2)

    @test d_init > d_end # Testing if the tuning results in a smaller behavioural distance
end