text:
@recipe function heatmap(mod::Medium, field=:vp)
	mx = mod.mgrid[2]
	mz = mod.mgrid[1]
	xlabel --> "x"
	ylabel --> "z"
	color --> :pu_or
	xlim --> (mx[1], mx[end])
	zlim --> (mz[1], mz[end])
	title --> string(field)
	clim --> (mod.bounds[field][1], mod.bounds[field][2])
	yflip := true
	mx, mz, mod[field]
end

@recipe function scatter(ageom::AGeom, ::SSrcs)
	sx = vcat([ag.sx for ag in ageom]...)
	sz = vcat([ag.sz for ag in ageom]...)
	legend --> false
	markersize --> 7
	markercolor := :red
	markershape := :xcross
	sx, sz
end

@recipe function scatter(ageom::AGeom, ::Recs)
	rx = vcat([ag.rx for ag in ageom]...)
	rz = vcat([ag.rz for ag in ageom]...)
	markersize --> 7
	legend --> false
	markercolor := :blue
	markershape := :utriangle
	rx, rz
end

@recipe function heatmap(dat::NamedD, field::Symbol=:P, wclip_perc=0.0, bclip_perc=0.0)

	#wclip::Vector{Float64}=[maximum(broadcast(maximum, td[id].d)) for id in 1:length(td)],
	#bclip::Vector{Float64}=[minimum(broadcast(minimum, td[id].d)) for id in 1:length(td)],

	dp=dat.d[field]
	dz = dat.grid
	dx = 1:size(dp,2)
	dmin=minimum(dp)
	dmax=maximum(dp)
	dmin=dmin+bclip_perc*inv(100)*abs(dmin)
	dmax=dmax-wclip_perc*inv(100)*abs(dmax)
	legend --> true
	xlabel --> "channel index"
	ylabel --> "time"
	color --> :grays
	clim --> (dmin, dmax)
	yflip := true
	dx, dz, dp
end

#=
Plot time-domain data of type `Data.TD`

# Arguments
* `td::Vector{Data.TD}` : time-domain data to be compared

# Keyword Arguments
* `ssvec::Vector{Vector{Int64}}=fill([1], length(td))` : supersource vector to be plotted
* `fields::Vector{Int64}=[1]` : field vector to be plotted
* `tr_flag::Bool=false` : plot time-reversed data when true
* `attrib::Symbol=:wav` : specify type of plot
"""

	if(ssvec===nothing)
		if(:r ∈ fields)
			urpos = AGeom_get([ageom],:urpos)
			rx=urpos[2]
			rz=urpos[1]
		else
			b=nothing
		end
		if(:s ∈ fields)
			uspos = AGeom_get([ageom],:uspos)
			sx=uspos[2]
			sz=uspos[1]
		else
			a=nothing
		end
	else
		if(:r ∈ fields)
			rxpos = [ageom.rx[iss] for iss in ssvec]
			rzpos = [ageom.rz[iss] for iss in ssvec]
			rx=vcat(rxpos...); rz=vcat(rzpos...)
		else
			b=nothing
		end
		if(:s ∈ fields)
			sxpos = [ageom.sx[iss] for iss in ssvec]
			szpos = [ageom.sz[iss] for iss in ssvec]
			sx=vcat(sxpos...); sz=vcat(szpos...)
		else
			a=nothing
		end
	end

	if(:r ∈ fields)
		@series begin
		end
	end
	if(:s ∈ fields)
		@series begin
		end
	end

end
=#

#=
@userplot Seismic
"""
Plot the velocity and density seismic models.

# Arguments

* `model::Medium` : model that should be plotted

# Keyword Arguments

* `xlim::Vector{Float64}=[model.mgrid[2][1],model.mgrid[2][end]]` : minimum and maximum limits of the second dimension while plotting
* `zlim::Vector{Float64}=[model.mgrid[1][1],model.mgrid[1][end]]` : minimum and maximum limits of the first dimension while plotting
* `fields::Vector{Symbol}=[:vp, :ρ]` : fields that are to be plotted
* `contrast_flag=false` : plot only the edges of the model
* `use_bounds=false` : adjust `clim` to the bounds in the seismic model

"""
@recipe function fseismic(p::Seismic;
		   fields=[:vp, :ρ],
		   contrast_flag=false,
		   use_bounds=false,
		  )
	if(contrast_flag)
		warn("ImageFiltering bug needs to be fixed")
	end
	model=p.args[1]
	nz,nx=length.(model.mgrid)

	nrow = (nx <= nz) ? 1 : length(fields)
	ncol = (nx <= nz) ? length(fields) : 1

	layout --> (nrow,ncol)
	for (i,iff) in enumerate(fields)
		name=string(iff)

		m = model[iff]
		mx = model.mgrid[2]
		mz = model.mgrid[1]
		if(contrast_flag)
			mmin=minimum(m)
			mmax=maximum(m)
			#m=imfilter(m, Kernel.Laplacian())
			mmmin=minimum(m)
			mmmax=maximum(m)
			for j in eachindex(m)
				m[j]=mmin+((m[j]-mmmin)*inv(mmmax-mmmin))*(mmax-mmin)
			end
		end
		mmin = use_bounds ? model.bounds[iff][1] : minimum(m)
		mmax = use_bounds ? model.bounds[iff][2] : maximum(m)
		@series begin
			subplot := i
			aspect_ratio := :equal
			seriestype := :heatmap
			legend --> true
			xlabel --> "x [m]"
			ylabel --> "z [m]"
			color --> :grays
			#xlim --> (mx[1], mx[end])
			#zlim --> (mz[1], mz[end])
			title --> name
			clim --> (mmin, mmax)
			yflip := true
			mx, mz, m
		end
	end

end
#

#"""
#save current fig using matlab2tikz
#"""
#function printfig(fig; filename::String = "FIG")
#	if filename == "FIG"
#		temp = 1;
#		file = join([filename string(temp) ".tex"])
#		while isfile(file)
#			temp += 1;
#			file = join([filename string(temp) ".tex"])
#		end
#		tikzfile = join([file[1:end-4] ".tikz"]);
#	else
#		file = filename;
#		tikzfile = join([file[1:end-4] ".tikz"]);
#	end
#	#matplotlib2tikz.save(tikzfile,fig);
#	fin = open(tikzfile,"r");
#	fout = open(file,"w");
#	lines = readlines(fin);
#	close(fin);
#	rm(tikzfile);
#	lines_old = ["\\documentclass[tikz]{standalone}\n";
#	"\\usepackage{pgfplots}\n";
#	"\\pgfplotsset{compat=newest}\n";
#	"\\usepackage{amsmath}\n";
#	"\\usepackage{siunitx}\n";
#	"\\begin{document}\n";lines;"\n\\end{document}\n"];
#
#	lines = ["% Generated from Julia\n";"\\begin{center}\n";lines;"\\end{center}"]
#	write(fout,lines)
#	close(fout)
#end

@userplot Spectrum

@recipe function pspectrum(p::Spectrum; tgrid=nothing)
	wav=p.args[1]
	if(tgrid===nothing)
		x=0:length(wav)
		tgrid=range(0.0, stop=Float64(length(wav)-1), step=1.0)
	end

	fgrid= DSP.rfftfreq(length(tgrid), inv(step(tgrid)))
	powwav = (abs.(FFTW.rfft(wav, [1])).^2)
	powwavdb = 10. * log10.(powwav./maximum(powwav)) # power in decibel after normalizing

	@series begin
		subplot := 1
		legend := false
		ylabel := "power (dB)"
		xlabel := "frequency (Hz)"
		fgrid, powwavdb
	end
end

@userplot Src

"""
Plot the source wavelet used for acquisition.

# Arguments

* `srcwav::Src` : source acquisition parameters
"""
@recipe function psrc(p::Src)
	srcwav=p.args[1]
	tgrid=srcwav.tgrid

	wav=srcwav.wav[1,1]
	layout := (2,1)
	@series begin
		subplot := 1
		legend := false
		xlabel := "time [s]"
		ylabel := "amplitude"
		tgrid, wav
	end

	fgrid= DSP.rfftfreq(length(tgrid), inv(step(tgrid)))
	powwav = (abs.(FFTW.rfft(wav, [1])).^2)
	powwavdb = 10. * log10.(powwav./maximum(powwav)) # power in decibel after normalizing
	@series begin
		subplot := 2
		ylabel := "power [dB]"
		xlabel := "frequency [Hz]"
		legend := false
		fgrid, powwavdb
	end

end

#=
function mscatter(x, y, titname="", axla=["",""])
    fact=1
    x=x[1:fact:end]
    x[:] /= norm(x) # normalize inputs before scatter
    y=y[1:fact:end]
    y[:] /= norm(y) # normalize inputs before scatter
    scatter(x,y,marker="x",s=0.5,color="gray")
    ax=gca()
    ax[:set_aspect]("equal")
    ax[:spines]["top"][:set_visible](false)
    ax[:spines]["right"][:set_visible](false)
    ax[:spines]["left"][:set_position]("center")
    ax[:spines]["bottom"][:set_position]("center")
    setp(ax[:get_yticklabels](false),color="none")
    setp(ax[:get_xticklabels](false),color="none")
    xmax = maximum(abs, x)
    xmax = iszero(xmax) ? 1.0 : xmax
    ymax = maximum(abs, y)
    ymax = iszero(ymax) ? 1.0 : ymax
    li = 1.2*max(xmax, ymax)
    #xlim(-li, li)
    #ylim(-li, li)
    xlabel(axla[1])
    ylabel(axla[2])
    title(titname)
    return nothing
end

function myhist(s, dist)
    fact=1
    x=s[1,1:fact:end]
    y=s[2,1:fact:end]
    nbins=div(length(x), 10)
    h1 = normalize(fit(Histogram, x, nbins=200, closed=:right))
    h2 = normalize(fit(Histogram, y, nbins=200, closed=:right))

#    h = PyPlot.plt[:hist](x, 200, color="red", histtype="barstacked", normed=true)
 #   h = PyPlot.plt[:hist](y, 200, color="blue",histtype="barstacked", normed=true)
    plot(h1.edges[1][2:end], h1.weights,color="red", linewidth=0.7)
    plot(h2.edges[1][2:end], h2.weights,color="blue", linewidth=0.7)
    ax=gca()
    xmax = 1.2*maximum(abs, x)
    xx=linspace(-xmax, xmax, 300)
    xlim(-xmax, xmax)
    yymax=0.
    for idist in 1:length(dist)
        yy = pdf(dist[idist], xx)
        yymax=maximum([maximum(yy), yymax])
        ax[:plot](xx, yy, color="black",linestyle="--", linewidth=1)
    end
    (yymax ≠ 0.) &&    ylim([0, 1.2*yymax])
    setp(ax[:get_yticklabels](),color="none")
    setp(ax[:get_xticklabels](),color="none")
    #setp(ax[:tick_params](false), labelbottom="off")
    title("Marginal Densities")
    return nothing
end
=#

end # module

=#


text:
module ExpressCommands

using AbInitioSoftwareBase: load, extension
using Comonicon: @cast, @main
using EquationsOfStateOfSolids: Parameters, EquationOfStateOfSolids
using Pkg: add, rm, gc, @pkg_str
using Preferences: @set_preferences!, @load_preference
using PrettyPrint: pprint
using Serialization: deserialize
using SimpleWorkflows: Workflow
import TikzGraphs
import TikzPictures

"""
Print the `file` in a pretty format.

# Arguments
- `file`: the file to be printed. Available extensions are `.jls`, `.json`, `.yaml`, `.yml` or `.toml`.
"""
@cast function print(file)
    ext = lowercase(extension(file))
    if ext == "jls"
        data = deserialize(file)
        if data isa Union{Parameters,EquationOfStateOfSolids}
            display(data)
        else
            pprint(data)
        end
    elseif ext in ("json", "yaml", "yml", "toml")
        data = load(file)
        pprint(data)
    else
        error("unsupported file extension `$ext`!")
    end
end

# @cast function install(plugin)
#     name = lowercase(plugin)
#     if name == "qe"
#         pkg"add https://github.com/MineralsCloud/QuantumESPRESSOExpress.jl.git"
#         # add("QuantumESPRESSOExpress")
#     else
#         error("unsupported plugin `$name`!")
#     end
# end

# @cast function uninstall(plugin)
#     name = lowercase(plugin)
#     if name == "qe"
#         rm("QuantumESPRESSOExpress")
#         gc()
#     else
#         error("unsupported plugin `$name`!")
#     end
# end

# @cast function use(plugin)
#     name = lowercase(plugin)
#     if name == "qe"
#         @set_preferences!("plugin" => "Quantum ESPRESSO")
#     else
#         error("unsupported plugin `$name`!")
#     end
# end

# macro load_plugin()
#     name = @load_preference("plugin")
#     if name == "Quantum ESPRESSO"
#         return :(@eval using QuantumESPRESSOExpress)
#     else
#         error("unsupported plugin `$name`!")
#     end
# end

@cast function graph(file)
    ext = extension(file)
    workflow = deserialize(file)
    typeassert(workflow, Workflow)
    if ext == "jls"
        TikzGraphs.save(
            TikzPictures.PDF(replace(file, ".jls" => ".pdf")),
            TikzGraphs.plot(workflow.graph),
        )
    else
        error("unsupported extension `$ext`!")
    end
end

include("EOS.jl")
@cast EOS

include("Phonon.jl")
@cast Phonon

include("QHA.jl")
@cast QHA

function whichmodule(recipe)
    recipe = lowercase(recipe)
    return if recipe == "eos"
        EOS
    elseif recipe in ("phonon dispersion", "vdos")
        Phonon
    elseif recipe in ("qha single", "multi qha")
        QHA
    else
        error("workflow `$recipe` is not recognized!")
    end
end

@cast function run(cfgfile)
    dict = load(cfgfile)
    mod = whichmodule(dict["recipe"])
    return mod.run(cfgfile)
end

"""
The main command `xps`.
"""
@main

end

ext: jl


text:
import .DataFrames

function build_default_getter(whitelist, blacklist, replace)
    function default_data(mc)
        d = Dict{Symbol, Any}()
        param = parameters(mc)
        for (k, v) in pairs(param)
            if !(k in blacklist) && (isempty(whitelist) || (k in whitelist))
                push!(d, replace(k => v))
            end
        end
        # push!(d, :mc => mc)
        d
    end
end

function DataFrames.DataFrame(mcs::Vector{Union{Nothing, T}}; kwargs...) where {T <: MonteCarloFlavor}
    clean = [mc for mc in mcs if mc !== nothing]
    DataFrame(clean; kwargs...)
end
function DataFrames.DataFrame(
        mcs::Vector{<: MonteCarloFlavor};
        whitelist = [], blacklist = [], replace = identity,
        getter = build_default_getter(whitelist, blacklist, replace)
    )

    df = mcs |> first |> getter |> DataFrames.DataFrame

    for i in 2:length(mcs)
        data = mcs[i] |> getter
        push!(df, data)
    end

    df.mc = mcs

    df
end


text:
# This file is a part of ValueShapes.jl, licensed under the MIT License (MIT).

using ValueShapes
using Test

using Distributions
using Random

@testset "const_value_dist" begin
    univariate_cvd = ConstValueDist(Int64(42))
    v = broadcast(Int64, [1,2, 3])
    multivariate_cvd = ConstValueDist(v)
    @test @inferred(length(multivariate_cvd)) == length(v)

    shape = varshape(univariate_cvd)

    @test @inferred(insupport(univariate_cvd, 42)) == true
    @test @inferred(insupport(univariate_cvd, 41.999)) == false
    @test @inferred(insupport(univariate_cvd, 42.999)) == false
    @test @inferred(totalndof(shape)) == 0
    @test @inferred(size(univariate_cvd)) == ()
    @test @inferred(length(univariate_cvd)) == 1

    @test @inferred(minimum(univariate_cvd)) == 42
    @test @inferred(maximum(univariate_cvd)) == 42
    @test @inferred(maximum(@inferred(pdf(univariate_cvd, -50:50)))) == Inf

    @test @inferred(pdf(univariate_cvd, 42)) == Inf

    @test @inferred(cdf(univariate_cvd, 41.999)) == 0
    @test @inferred(cdf(univariate_cvd, 42)) == 1
    @test @inferred(cdf(univariate_cvd, Inf)) == 1

    @test @inferred(mode(univariate_cvd)) == 42

    @test @inferred(rand(univariate_cvd, 10)) == fill(42, 10)

    @test @inferred(logpdf(univariate_cvd, -1)) == -Inf
    @test @inferred(logpdf(univariate_cvd, 41.999)) == -Inf
    @test @inferred(logpdf(univariate_cvd, 42)) == Inf
    @test @inferred(logpdf(univariate_cvd, 43.999)) == -Inf

    @test @inferred(eltype(univariate_cvd)) == Int64
end


text:
using Test
using QuantumOptics, CollectiveSpins

@testset "geometry" begin

cs = CollectiveSpins

@test cs.geometry.chain(2., 2) isa Vector{Vector{Float64}}
@test length(cs.geometry.chain(1., 8)) == 8
@test cs.geometry.chain(1., 3) == Vector{Float64}[[0.,0.,0.], [1.,0.,0.], [2.,0.,0.]]
@test cs.geometry.chain(2, 2) == Vector{Float64}[[0.,0.,0.], [2.,0.,0.]]

@test cs.geometry.triangle(1.0) isa Vector{Vector{Float64}}
@test length(cs.geometry.triangle(1)) == 3
@test cs.geometry.triangle(1.0) == Vector{Float64}[[0.,0.,0.], [1.,0.,0.], [0.5,0.5*sqrt(3),0.]]

@test cs.geometry.rectangle(2., 1.) isa Vector{Vector{Float64}}
@test length(cs.geometry.rectangle(2., 1.; Nx=4, Ny=5)) == 20
@test cs.geometry.rectangle(2., 1.) == Vector{Float64}[[0.,0.,0.], [2.,0.,0.], [0.,1.,0.], [2.,1.,0.]]
@test cs.geometry.rectangle(2, 1; Nx=2, Ny=1) == Vector{Float64}[[0.,0.,0.], [2.,0.,0.]]

@test cs.geometry.square(2.) isa Vector{Vector{Float64}}
@test length(cs.geometry.square(2; Nx=3, Ny=4)) == 12
@test cs.geometry.square(2) == Vector{Float64}[[0.,0.,0.], [2.,0.,0.], [0.,2.,0.], [2.,2.,0.]]
@test cs.geometry.square(2; Nx=1, Ny=1) == Vector{Float64}[[0.,0.,0.]]

@test cs.geometry.hexagonal(1) isa Vector{Vector{Float64}}
@test length(cs.geometry.hexagonal(1; Nr=1)) == 7
@test length(cs.geometry.hexagonal(1; Nr=2)) == 19
@test cs.geometry.hexagonal(2; Nr=0) == Vector{Float64}[[0.,0.,0.]]
a = sqrt(3/4)
@test cs.geometry.hexagonal(1; Nr=1) == Vector{Float64}[[-a,-1/2,0.], [-a,1/2,0.], [0,-1.,0.], [0.,0.,0.], [0.,1.,0.], [a,-1/2,0.], [a,1/2,0.]]

@test cs.geometry.box(2., 1, 3) isa Vector{Vector{Float64}}
@test length(cs.geometry.box(2, 1, 3; Nx=2, Ny=3, Nz=5)) == 30
@test cs.geometry.box(2, 1, 3; Nx=2, Ny=3, Nz=1) == cs.geometry.rectangle(2, 1; Nx=2, Ny=3)
rect = cs.geometry.rectangle(2, 1; Nx=2, Ny=3)
box = cs.geometry.box(2, 1, 3; Nx=2, Ny=3, Nz=4)
for i=0:3
    plane = box[i*6+1:i*6+6]
    plane_shifted = Vector{Float64}[x-Float64[0.,0.,i*3] for x=plane]
    @test plane_shifted == rect
end

@test cs.geometry.cube(2.) isa Vector{Vector{Float64}}
@test length(cs.geometry.cube(2; Nx=2, Ny=3, Nz=4)) == 24
rect = cs.geometry.rectangle(2, 2; Nx=2, Ny=3)
cube = cs.geometry.cube(2; Nx=2, Ny=3, Nz=4)
for i=0:3
    plane = cube[i*6+1:i*6+6]
    plane_shifted = Vector{Float64}[x-Float64[0.,0.,i*2] for x=plane]
    @test plane_shifted == rect
end

end # testset


text:
module JuliaDot

include("_dot.jl")

end


text:
using LowRank
using Plots

function plot_changing_n(method, method_name, d, k, std, num_times, ns)
    mean_errors = []
    mean_indices_correctnesses = []
    mean_times = []
    baseline_errors = []
    for n in ns
        N = 10*n
        A = LowRank.matrix_with_normal_noise(N, d, n, k, std)
        errors, indicess, times = LowRank.run_method(B -> method(B, n, k), A, num_times)
        baseline_error, _ = LowRank.truncated_svd(A[:,1:n], k)
        baseline_indices = Set(1:n)
        mean_error = sum(errors)/num_times
        mean_indices_correctness = sum([length(intersect(inds,baseline_indices))/n for inds in indicess])/num_times
        mean_time = sum(times)/num_times
        push!(mean_errors, mean_error)
        push!(mean_indices_correctnesses, mean_indices_correctness)
        push!(mean_times, mean_time)
        push!(baseline_errors, baseline_error)
    end

    fontsize = 10
    texts = ["method: $method_name", "number of points N: 10n", "dimension d: $d", "target rank k: $k", "noise std: $std", "number of runs: $num_times"]
    y_locations = 0.9:-0.15:0
    annotations = [(0.2,y_locations[i], text(t, :left, fontsize)) for (i,t) in enumerate(texts)]
    p_text = plot([],[], ann=annotations, border =:none, label=:none)
    p_error = plot()
    plot!(p_error, ns, baseline_errors, label = "baseline errors")
    plot!(p_error, ns, mean_errors, label = "mean errors")
    p_indices = plot(ns, mean_indices_correctnesses, label = "mean indices correctness")
    ylims!(p_indices, (0,1))
    p_time = plot(ns, mean_times, label = "mean time")
    xaxis!(p_error, "size of target set, n")
    xaxis!(p_indices, "size of target set, n")
    xaxis!(p_time, "size of target set, n")
    p = plot(p_text,p_error, p_indices, p_time)
end

method_name = "Global FastLowRank"
method = FastLowRank
d = 10
k = 4
std = 0.1
num_times = 10
ns = 5:5:50
p = plot_changing_n(method, method_name, d, k, std, num_times, ns)
savefig(p, "GlobalFastLowRank.png")
method_name = "Local FastLowRank"
method = (B, n, k) -> FastLowRank(B, n, k, step = LowRank.L2_step_local)
p = plot_changing_n(method, method_name, d, k, std, num_times, ns)
savefig(p, "LocalFastLowRank.png")


text:
import CUDA.CUDNN:
    cudnnSetTensor!,
    cudnnScaleTensor!,
    cudnnScaleTensor,
    cudnnAddTensor!,
    cudnnAddTensor,
    CUDNN_TENSOR_NHWC

@testset "cudnn/inplace" begin
    x = CUDA.rand(10)
    cudnnSetTensor!(x, 7)
    @test all(isequal(7), Array(x))
    ax = rand(10)
    cx = CuArray(ax)
    @test 7*ax ≈ cudnnScaleTensor(cx, 7) |> Array
    @test 7*ax ≈ cudnnScaleTensor!(similar(cx), cx, 7) |> Array
    ax,ab = rand(5,4,3,2),rand(1,1,3,1)
    cx,cb = CuArray.((ax,ab))
    @test ax .+ ab ≈ cudnnAddTensor(cx, cb) |> Array
    @test ax .+ 7*ab ≈ cudnnAddTensor(cx, cb, alpha=7) |> Array
    @test 7*ax .+ ab ≈ cudnnAddTensor(cx, cb, beta=7) |> Array
    @test ax .+ ab ≈ cudnnAddTensor!(similar(cx), cx, cb) |> Array
    @test ax .+ 7*ab ≈ cudnnAddTensor!(similar(cx), cx, cb, alpha=7) |> Array
    @test 7*ax .+ ab ≈ cudnnAddTensor!(similar(cx), cx, cb, beta=7) |> Array
    @test ax .+ ab ≈ cudnnAddTensor!(cx, cx, cb) |> Array
    @test ax .+ ab ≈ cx |> Array
    ax,ab = rand(3,5,4,2),rand(3,1,1,1)
    cx,cb = CuArray.((ax,ab))
    @test ax .+ ab ≈ cudnnAddTensor(cx, cb, format=CUDNN_TENSOR_NHWC) |> Array
end


text:
"""
```
    Canny <: AbstractEdgeDetectionAlgorithm
    Canny(;) #TODO Decide on keywords

    detect_edges([T,] img, f::Canny)
    detect_edges!([out,] img, f::Canny)
```

Returns an image depicting the edges of the input image.

# Details

TODO

# Example

```julia

using TestImages, FileIO, ImageView

img =  testimage("mandril_gray")
img_edges = detect_edges(img, Canny())

imshow(img)
imshow(img_edges)
```

# References
TODO.
"""
@with_kw struct Canny{T₁ <: Union{Real,AbstractGray},
                      T₂ <: Union{Real,AbstractGray}} <: AbstractEdgeDetectionAlgorithm
    # TO BE DECIDED
    example_keyword_1::T₁ = 0
    example_keyword_2::T₂ = 0
end

function (f::Canny)(out::GenericGrayImage, img::GenericGrayImage)
    # TODO
end

function (f::Canny)(out::AbstractArray{<:Color3}, img::AbstractArray{<:Color3})
    # TODO
end

(f::Canny)(out::GenericGrayImage, img::AbstractArray{<:Color3}) =
    f(out, of_eltype(Gray, img))


text:
abstract type Primitive{T} end

struct NormedPrimitive{P,T} <: Primitive{T}
  radius::T
  NormedPrimitive{P,T}(radius::T) where {P,T} = new{P,T}(radius)
  NormedPrimitive{P}(radius::T) where {P,T} = NormedPrimitive{P,T}(radius)
end

norm(p::Point, ::Type{<:NormedPrimitive{P}}) where {P} = norm(p, P)

(np::NormedPrimitive)(p) = norm(p, typeof(np)) - np.radius

(≈)(x::NormedPrimitive, y::NormedPrimitive) = typeof(x) == typeof(y) && x.radius ≈ y.radius

const HyperSphere{T} = NormedPrimitive{2,T}

norm(p::Point, ::Type{<:HyperSphere}) = hypot(p...)

const HyperCube{T} = NormedPrimitive{Inf,T}

"""
    Ellipsoid(semiaxes)

An ellipsoid with semi-axes `semiaxes`.
Is equivalent to a scaled [`HyperSphere`](@ref).
"""
const Ellipsoid{Dim,T} = Scaled{HyperSphere{T},Dim,T}
Ellipsoid(radius::T, transf::Scaling{Dim,T}) where {Dim,T} = Ellipsoid{Dim,T}(HyperSphere(radius), transf)

function Ellipsoid(semiaxes::AbstractVector)
  radius = norm(semiaxes, HyperSphere)
  Ellipsoid(radius, Scaling(semiaxes ./ radius))
end
Ellipsoid(semiaxes::Number...) = Ellipsoid(collect(semiaxes))

(≈)(x::Ellipsoid, y::Ellipsoid) = x.obj.radius .* x.transf.vec ≈ y.obj.radius .* y.transf.vec

Base.show(io::IO, elps::Ellipsoid{Dim,T}) where {Dim,T} = print(io, "Ellipsoid{$Dim, $T}($(elps.transf.vec .* elps.obj.radius))")

# A box around the origin. This is not the same as `Meshes.Box`, it is missing the translation part.
const Box{Dim,T} = Scaled{HyperCube{T},Dim,T}
Box(radius::T, transf::Scaling{Dim,T}) where {Dim,T} = Box{Dim,T}(HyperCube(radius), transf)

apply(t::Scaling, s::NormedPrimitive) = typeof(s)(norm(s, t.vec) * s.radius)
apply(t::Rotation, s::HyperSphere) = s

const Circle{T} = Projection{2,HyperSphere{T}}
const Square{T} = Projection{2,HyperCube{T}}


text:
using Statistics

export testchaos01

"""
    testchaos01(φ::Vector [, cs, N0]) -> chaotic?
Perform the so called "0-1" test for chaos introduced by Gottwald and
Melbourne [1] on the timeseries `φ`.
Return `true` if `φ` is chaotic, `false` otherwise.

## Description
This method tests if the given timeseries is chaotic or not by transforming
it into a two-dimensional diffusive process. If the timeseries is chaotic,
the mean square displacement of the process grows as `sqrt(length(φ))`,
while it stays constant if the timeseries is regular.
The implementation here computes `K`, the correlation coefficient (median
of `Kc for c ∈ cs`), and simply checks if `K > 0.5`.

If you want to access the various `Kc` you should call the method
`testchaos01(φ, c::Real, N0)` which returns `Kc`.

`cs` defaults to `3π/5*rand(10) + π/4` and `N0`, the length
of the two-dimensional process, is `N0 = length(φ)/10`.

Notice that for data sampled from continous dynamical systems, some
care must be taken regarding the values of `cs`, see [1].

## References

[1] : Gottwald & Melbourne, “The 0-1 test for chaos: A review”
[Lect. Notes Phys., vol. 915, pp. 221–247, 2016.](www.doi.org/10.1007/978-3-662-48410-4_7)
"""
function testchaos01(φ::Vector, cs = 3π/5*rand(10) .+ π/4, N0 = Int(length(φ)÷10))
    K = median(testchaos01(φ, c, N0) for c in cs)
    return K > 0.5
end

function testchaos01(φ::Vector, c::Real, N0 = Int(length(φ)÷10))
    N, E = length(φ), mean(φ)
    @assert N0 ≤ N/10
    pc, qc = trigonometric_decomposition(φ, c)
    Dc = mmsd(E, pc, qc, N0, c)
    Kc = cor(Dc, 1:N0)
    return Kc
end

function trigonometric_decomposition(φ, c)
    X = promote_type(eltype(φ), typeof(c))
    N = length(φ)
    pc, qc = zeros(X, N), zeros(X, N)
    @inbounds for n in 1:N-1
        si, co = sincos(c*n)
        pc[n+1] = pc[n] + φ[n]*co
        qc[n+1] = qc[n] + φ[n]*si
    end
    return pc, qc
end

"modified mean square displacement"
function mmsd(mf, pc::Vector{T}, qc, N0, c) where T
    N, Dc = length(pc), zeros(T, N0)
    f = mf^2/(1-cos(c)) # constant factor for Dc
    for n in 1:N0
        Mc = sum((pc[n+j] - pc[j])^2 + (qc[n+j] - qc[j])^2 for j in 1:N-n)/(N-n)
        Dc[n] = Mc - f*(1 - cos(n*c))
    end
    return Dc
end


text:
using Ansillary

println("Press Ctrl+c to exit.")

Screen.raw() do
    for event in Inputs.Events()
        if event == Inputs.CTRL_C
            break
        end
        println(event)
    end
end

println("Bye!")


text:
# ---
# title: 846. Hand of Straights
# id: problem846
# author: Indigo
# date: 2021-07-09
# difficulty: Medium
# categories: Ordered Map
# link: <https://leetcode.com/problems/hand-of-straights/description/>
# hidden: true
# ---
#
# Alice has a `hand` of cards, given as an array of integers.
#
# Now she wants to rearrange the cards into groups so that each group is size
# `W`, and consists of `W` consecutive cards.
#
# Return `true` if and only if she can.
#
# **Note:** This question is the same as 1296:
# <https://leetcode.comhttps://leetcode.com/problems/divide-array-in-sets-of-k-consecutive-numbers/>
#
#
#
# **Example 1:**
#
#
#
#     Input: hand = [1,2,3,6,2,3,4,7,8], W = 3
#     Output: true
#     Explanation: Alice's hand can be rearranged as [1,2,3],[2,3,4],[6,7,8]
#
#
# **Example 2:**
#
#
#
#     Input: hand = [1,2,3,4,5], W = 4
#     Output: false
#     Explanation: Alice's hand can't be rearranged into groups of 4.
#
#
#
#
#
# **Constraints:**
#
#   * `1 <= hand.length <= 10000`
#   * `0 <= hand[i] <= 10^9`
#   * `1 <= W <= hand.length`
#
#
## @lc code=start
using LeetCode

function is_n_straight_hand(hand::Vector{Int}, group_size::Integer)::Bool
    length(hand) % group_size == 0 || return false
    cntr = counter(hand)
    while !isempty(cntr)
        k = minimum(keys(cntr))
        v = cntr[k]
        for i in (0:group_size-1) .+ k
            cntr[i] -= v
            cntr[i] < 0 && return false
            cntr[i] == 0 && pop!(cntr.map, i)
        end
    end
    return true
end
## @lc code=end


text:
# # Minimal Robust Positively Invariant Set

#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/generated/Minimal Robust Positively Invariant Set.ipynb)
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/generated/Minimal Robust Positively Invariant Set.ipynb)

# ### Introduction

# In this example, we implement the method presented in [RKKM05] to compute a robust positively invariant polytope of a linear system under a disturbance bounded by a polytopic set.
#
# We consider the example given in equation (15) of [RKKM05]:
# ```math
# x^+ =
# \begin{bmatrix}
#   1 & 1\\
#   0 & 1
# \end{bmatrix}x +
# \begin{bmatrix}
#   1\\
#   1
# \end{bmatrix} u
# + w
# ```
# with the state feedback control $u(x) = -\begin{bmatrix}1.17 & 1.03\end{bmatrix} x$.
# The controlled system is therefore
# ```math
# x^+ =
# \left(\begin{bmatrix}
#   1 & 1\\
#   0 & 1
# \end{bmatrix} -
# \begin{bmatrix}
#   1\\
#   1
# \end{bmatrix}
# \begin{bmatrix}1.17 & 1.03\end{bmatrix}\right)x
# + w =
# \begin{bmatrix}
#   -0.17 & -0.03\\
#   -1.17 & -0.03
# \end{bmatrix}x
# + w.
# ```
#
# [RKKM05] Sasa V. Rakovic, Eric C. Kerrigan, Konstantinos I. Kouramas, David Q. Mayne *Invariant approximations of the minimal robust positively Invariant set*. IEEE Transactions on Automatic Control 50 (**2005**): 406-410.

A = [1 1; 0 1] - [1; 1] * [1.17 1.03]

# The set of disturbance is the unit ball of the infinity norm.

using Test #jl
using Polyhedra
Wv = vrep([[x, y] for x in [-1.0, 1.0] for y in [-1.0, 1.0]])

# We will use the default library for this example but feel free to pick any other library from [this list of available libraries](https://juliapolyhedra.github.io/) such as [CDDLib](https://github.com/JuliaPolyhedra/CDDLib.jl).
# The LP solver used to detect redundant points in the V-representation is [GLPK](https://github.com/JuliaOpt/GLPK.jl). Again, you can replace it with any other solver listed [here](http://jump.dev/JuMP.jl/dev/installation/#Getting-Solvers-1) that supports LP.

using GLPK
using JuMP
lib = DefaultLibrary{Float64}(GLPK.Optimizer)
W = polyhedron(Wv, lib)

# The $F_s$ function of equation (2) of [RKKM05] can be implemented as follows.

function Fs(s::Integer, verbose=1)
    @assert s ≥ 1
    F = W
    A_W = W
    for i in 1:(s-1)
        A_W = A * A_W
        F += A_W
        if verbose ≥ 1
            println("Number of points after adding A^$i * W: ", npoints(F))
        end
        removevredundancy!(F)
        if verbose ≥ 1
            println("Number of points after removing redundant ones: ", npoints(F))
        end
    end
    return F
end

# We can see below that only the V-representation is computed. In fact, no H-representation was ever computed during `Fs`. Computing $AW$ is done by multiplying all the points by $A$ and doing the Minkowski sum is done by summing each pair of points. The redundancy removal is carried out by CDD's internal LP solver.

@time Fs(4) #!jl
@test npoints(Fs(4)) == 16 #jl

# The Figure 1 of [RKKM05] can be reproduced as follows:

using Plots          #!jl
plot()               #!jl
for i in 10:-1:1     #!jl
    plot!(Fs(i, 0))  #!jl
end                  #!jl
# The cell needs to return the plot for it to be displayed
# but the `for` loop returns `nothing` so we add this dummy `plot!` that returns the plot
plot!()              #!jl

# Now, suppose we want to compute an invariant set by scaling $F_s$ by the appropriate $\alpha$.
# In equation (11) of [RKKM05], we want to check whether $A^s W \subseteq \alpha W$ which is equivalent to $W \subseteq \alpha A^{-s} W$.
# Note that `A^s \ W` triggers the computation of the H-representation of `W` and `A_W` is H-represented.

function αo(s)
    A_W = A^s \ W
    hashyperplanes(A_W) && error("HyperPlanes not supported")
    return maximum([Polyhedra.support_function(h.a, W) / h.β for h in halfspaces(A_W)])
end

# We obtain $\alpha \approx 1.9 \cdot 10^{-5}$ like in [RKKM05].

α = αo(10)
@test α ≈ 1.91907e-5 #jl

# The scaled set is is the following:

using Plots                    #!jl
plot((1 - α)^(-1) * Fs(10, 0)) #!jl

max_stars_repo_path: examples/Minimal Robust Positively Invariant Set.jl

max_stars_repo_name: blegat/Polyhedra

max_issues_repo_name: blegat/Polyhedra

max_forks_repo_name: blegat/Polyhedra


text:
using DelimitedFiles, Printf
using Plots
using CuArrays

using Oceananigans
using Oceananigans.Diagnostics
using Oceananigans.OutputWriters
using Oceananigans.AbstractOperations
using Oceananigans.Utils

using Oceananigans: Face, Cell

# setting up a 2-d model in the style of NG17

#####
##### Some useful constants
#####

const km = 1000
const Ω_Earth = 7.292115e-5  # [s⁻¹]
const φ = -75  # degrees latitude

#####
##### Model grid and domain size
#####

arch = CPU()
FT = Float64

Nx = 1
Ny = 32
Nz = 32

Lx = 5km/Ny
Ly = 5km
Lz = 300

end_time = 7day

# Make up temperature profiles
zC = collect(((-Lz:Lz/Nz:0).+Lz/(2*Nz))[1:end-1])
T₀ = collect(1:2/(Nz-1):3)
S₀ = 34*ones(Nz)

# convert to CuArray if we are running on a GPU
if arch == GPU()
    T₀ = CuArray(T₀)
    S₀ = CuArray(S₀)
end

#####
##### Set up relaxation areas for the meltwater source and for the northern boundary
#####

# Meltwater source location - implemented as a box
source_corners_m = ((1,1,1),(Lx,50,10))
N = (Nx,Ny,Nz)
L = (Lx,Ly,Lz)
source_corners = (Int.(ceil.(source_corners_m[1].*N./L)),Int.(ceil.(source_corners_m[2].*N./L)))

λ = 1/(60)  # Relaxation timescale [s⁻¹].

# Temperature and salinity of the meltwater outflow.
T_source = 2.0
S_source = 34.0

# Specify width of stable relaxation area
stable_relaxation_width_m = 400.0
stable_relaxation_width = Int(ceil(stable_relaxation_width_m.*Ny./Ly))

# Forcing functions
@inline T_relax(i, j, k, grid, time, U, C, p) =
    @inbounds ifelse(p.source_corners[1][1]<=i<=p.source_corners[2][1] && p.source_corners[1][2]<=j<=p.source_corners[2][2] && p.source_corners[1][3]<=k<=p.source_corners[2][3], -p.λ * (C.T[i, j, k] - p.T_source), 0) +
              ifelse(j>grid.Ny-p.stable_relaxation_width, -p.λ * (C.T[i, j, k] - p.T₀[k]), 0)

@inline S_relax(i, j, k, grid, time, U, C, p) =
    @inbounds ifelse(p.source_corners[1][1]<=i<=p.source_corners[2][1] && p.source_corners[1][2]<=j<=p.source_corners[2][2] && p.source_corners[1][3]<=k<=p.source_corners[2][3], -p.λ * (C.S[i, j, k] - p.S_source), 0) +
              ifelse(j>grid.Ny-p.stable_relaxation_width, -p.λ * (C.S[i, j, k] - p.S₀[k]), 0)

params = (source_corners=source_corners, T_source=T_source, S_source=S_source, λ=λ, stable_relaxation_width=stable_relaxation_width, T₀=T₀,S₀=S₀)

forcing = ModelForcing(T = T_relax, S = S_relax)

#####
##### Set up model and simulation
#####

topology = (Periodic, Bounded, Bounded)
grid = RegularCartesianGrid(topology=topology, size=(Nx, Ny, Nz), x=(-Lx/2, Lx/2), y=(0, Ly), z=(-Lz, 0))

# eos = LinearEquationOfState()
eos = RoquetIdealizedNonlinearEquationOfState(:freezing)

model = IncompressibleModel(
           architecture = arch,
             float_type = FT,
                   grid = grid,
                tracers = (:T, :S, :meltwater),
               coriolis = FPlane(rotation_rate=Ω_Earth, latitude=φ),
               buoyancy = SeawaterBuoyancy(equation_of_state=eos),
                closure = AnisotropicMinimumDissipation(),
                forcing = forcing,
             parameters = params
)

#####
##### Setting up initial conditions
#####

T₀_3D = repeat(reshape(T₀, 1, 1, Nz), Nx, Ny, 1)
S₀_3D = repeat(reshape(S₀, 1, 1, Nz), Nx, Ny, 1)

set!(model.tracers.T, T₀_3D)
set!(model.tracers.S, S₀_3D)

# Set meltwater concentration to 1 at the source.
model.tracers.meltwater.data[source_corners[1][1]:source_corners[2][1],source_corners[1][2]:source_corners[2][2],source_corners[1][3]:source_corners[2][3]] .= 1

#####
##### Write out 3D fields and slices to NetCDF files.
#####

# Define vorticity computation
u, v, w = model.velocities
vorticity_operation = ∂x(v) - ∂y(u)
ω = Field(Face, Face, Cell, model.architecture, model.grid, TracerBoundaryConditions(grid))
vorticity_computation = Computation(vorticity_operation, ω)

function get_vorticity(model)
    compute!(vorticity_computation)
    return Array(interior(ω))
end

fields = Dict(
        "u" => model.velocities.u,
        "v" => model.velocities.v,
        "w" => model.velocities.w,
        "T" => model.tracers.T,
        "S" => model.tracers.S,
"meltwater" => model.tracers.meltwater,
       "nu" => model.diffusivities.νₑ,
   "kappaT" => model.diffusivities.κₑ.T,
   "kappaS" => model.diffusivities.κₑ.S,
   "vorticity" => get_vorticity
)

dimensions = Dict(
    "vorticity" => ("xF", "yF", "zC")
)

output_attributes = Dict(
"meltwater" => Dict("longname" => "Meltwater concentration"),
       "nu" => Dict("longname" => "Nonlinear LES viscosity", "units" => "m^2/s"),
   "kappaT" => Dict("longname" => "Nonlinear LES diffusivity for temperature", "units" => "m^2/s"),
   "kappaS" => Dict("longname" => "Nonlinear LES diffusivity for salinity", "units" => "m^2/s"),
   "vorticity" => Dict("longname" => "Vorticity", "units" => "1/s")
)

eos_name(::LinearEquationOfState) = "LinearEOS"
eos_name(::RoquetIdealizedNonlinearEquationOfState) = "RoquetEOS"
prefix = "ice_shelf_meltwater_outflow_2d_$(eos_name(eos))_"

#####
##### Print banner
#####

@printf("""

    Simulating ocean dynamics of meltwater outflow from beneath Antarctic ice shelves
        N : %d, %d, %d
        L : %.3g, %.3g, %.3g [km]
        Δ : %.3g, %.3g, %.3g [m]
        φ : %.3g [latitude]
        f : %.3e [s⁻¹]
     days : %d
 T_source : %.2f [°C]
 S_source : %.2f [g/kg]
  closure : %s
      EoS : %s

""", model.grid.Nx, model.grid.Ny, model.grid.Nz,
     model.grid.Lx / km, model.grid.Ly / km, model.grid.Lz / km,
     model.grid.Δx, model.grid.Δy, model.grid.Δz,
     φ, model.coriolis.f, end_time / day,
     T_source, S_source,
     typeof(model.closure), typeof(model.buoyancy.equation_of_state))

#####
##### Time step!
#####

# Wizard utility that calculates safe adaptive time steps.
wizard = TimeStepWizard(cfl=0.1, Δt=1second, max_change=1.2, max_Δt=10second)

# Number of time steps to perform at a time before printing a progress
# statement and updating the adaptive time step.
Ni = 20

# CFL utilities for reporting stability criterions.
cfl = AdvectiveCFL(wizard)
dcfl = DiffusiveCFL(wizard)

function progress_statement(simulation)
    model = simulation.model
    C_mw = model.tracers.meltwater  # Convenient alias

    # add passive meltwater tracer at source, remove at boundary
    C_mw.data[source_corners[1][1]:source_corners[2][1],source_corners[1][2]:source_corners[2][2],source_corners[1][3]:source_corners[2][3]] .= 1
    C_mw.data[:,Ny-stable_relaxation_width:Ny,:] .= 0

    ## Normalize meltwater concentration to be 0 <= C_mw <= 1.
    #C_mw.data .= max.(0, C_mw.data)
    #C_m w.data .= C_mw.data ./ maximum(C_mw.data)

    # Calculate simulation progress in %.
    progress = 100 * (model.clock.time / end_time)

    # Find maximum velocities.
    umax = maximum(abs, model.velocities.u.data.parent)
    vmax = maximum(abs, model.velocities.v.data.parent)
    wmax = maximum(abs, model.velocities.w.data.parent)

    # Find maximum ν and κ.
    νmax = maximum(model.diffusivities.νₑ.data.parent)
    κmax = maximum(model.diffusivities.κₑ.T.data.parent)

    # Print progress statement.
    i, t = model.clock.iteration, model.clock.time
    @printf("[%06.2f%%] i: %d, t: %5.2f days, umax: (%6.3g, %6.3g, %6.3g) m/s, CFL: %6.4g, νκmax: (%6.3g, %6.3g), νκCFL: %6.4g, next Δt: %8.5g s\n",
            progress, i, t / day, umax, vmax, wmax, cfl(model), νmax, κmax, dcfl(model), wizard.Δt)
end

_ω = get_vorticity(model)
# display(_ω)
# @show typeof(_ω)

# Simulation that manages time stepping.
simulation = Simulation(model, Δt=wizard, stop_time=end_time, progress=progress_statement, progress_frequency=Ni)

simulation.output_writers[:fields] =
    NetCDFOutputWriter(model, fields, filename = prefix * "fields_vort_test.nc",
                       interval = 6hour, output_attributes = output_attributes,
		       dimensions = dimensions)

simulation.output_writers[:along_channel_slice] =
    NetCDFOutputWriter(model, fields, filename = prefix * "along_channel_yz_slice_vort_test.nc",
                       interval = 5minute, output_attributes = output_attributes,
                       dimensions = dimensions, xC = 1, xF = 1)

run!(simulation)

for ow in simulation.output_writers
    ow isa NetCDFOutputWriter && close(ow)
end


text:
module ArnoldiMethod

using LinearAlgebra
using StaticArrays

using Base: RefValue, OneTo

export partialschur, LM, SR, LR, SI, LI, partialeigen

"""
    Arnoldi(n, k) → Arnoldi

Pre-allocated Arnoldi relation of the Vₖ₊₁ and Hₖ matrices that satisfy
A * Vₖ = Vₖ₊₁ * Hₖ, where Vₖ₊₁ is orthonormal of size n × (k+1) and Hₖ upper
Hessenberg of size (k+1) × k. The constructor will just allocate sufficient
space, but will *not* initialize the first vector of `v₁`. For the latter see
`reinitialize!`.
"""
struct Arnoldi{T,TV<:StridedMatrix{T},TH<:StridedMatrix{T}}
    V::TV
    H::TH

    function Arnoldi{T}(matrix_order::Int, krylov_dimension::Int) where {T}
        krylov_dimension <= matrix_order || throw(ArgumentError("Krylov dimension should be less than matrix order."))
        V = Matrix{T}(undef, matrix_order, krylov_dimension + 1)
        H = zeros(T, krylov_dimension + 1, krylov_dimension)
        return new{T,typeof(V),typeof(H)}(V, H)
    end
end

"""
    RitzValues(maxdim) → RitzValues

Convenience wrapper for Ritz values + residual norms and some permutation of
these values. The Ritz values are computed from the active part of the
Hessenberg matrix `H[active:maxdim,active:maxdim]`.

When computing exact shifts in the implicit restart, we need to reorder the Ritz
values in some way. For convenience we simply keep track of a permutation `ord`
of the Ritz values rather than moving the Ritz values themselves around. That
way we don't lose the order of the residual norms.
"""
struct RitzValues{Tv,Tr}
    λs::Vector{Tv}
    rs::Vector{Tr}
    ord::Vector{Int}

    function RitzValues{T}(maxdim::Int) where {T}
        λs = Vector{complex(T)}(undef, maxdim)
        rs = Vector{real(T)}(undef, maxdim)
        ord = Vector{Int}(undef, maxdim)
        return new{complex(T),real(T)}(λs, rs, ord)
    end
end

"""
    PartialSchur(Q, R, eigenvalues)

Holds an orthonormal basis `Q` and a (quasi) upper triangular matrix `R`.

For convenience the eigenvalues that appear on the diagonal of `R` are also
listed as `eigenvalues`, which is in particular useful in the case of real
matrices with complex eigenvalues. Note that the eigenvalues are always a
complex, even when the matrix `R` is real.
"""
struct PartialSchur{TQ,TR,Tλ<:Complex}
    "Orthonormal matrix"
    Q::TQ
    "Quasi upper triangular matrix"
    R::TR
    "Complex-valued vector of eigenvalues"
    eigenvalues::Vector{Tλ}
end

include("targets.jl")
include("partition.jl")
include("schurfact.jl")
include("schursort.jl")
include("restore_hessenberg.jl")
include("expansion.jl")
include("run.jl")
include("eigvals.jl")
include("eigenvector_uppertriangular.jl")
include("show.jl")

end


text:
include("test_params.jl")
println("Testing whether probPickFirst works.")
using SymbolicRegression, Test

n = 10

options = Options(;
    default_params...,
    binary_operators=(+, -, *, /),
    unary_operators=(cos, sin),
    probPickFirst=0.999,
    ns=n,
)

for reverse in [false, true]
    members = PopMember{Float32}[]

    # Generate members with scores from 0 to 1:
    for i in 1:n
        tree = Node("x1") * 3.2f0
        score = Float32(i - 1) / (n - 1)
        if reverse
            score = 1 - score
        end
        test_loss = 1.0f0  # (arbitrary for this test)
        push!(members, PopMember(tree, score, test_loss))
    end

    pop = Population(members, n)

    dummy_frequencies = [0.0f-10 for i in 1:100]
    best_pop_member = [
        SymbolicRegression.best_of_sample(pop, dummy_frequencies, options).score for
        j in 1:100
    ]

    mean_value = sum(best_pop_member) / length(best_pop_member)

    # Make sure average score is small
    @test mean_value < 0.1
end

println("Passed.")


text:
using Statistics

abstract type CorMethod end
struct CorAll <: CorMethod end
struct CorComplete <: CorMethod end
struct CorPairwise <: CorMethod end

function Statistics.cor(df::DataFrame; method::CorMethod=CorAll())
    cor1(i, j) = nrow(df2) == 0 ? missing : cor(df2[i], df2[j])

    function cor2(i, j)
       x = dropmissing(DataFrame([df2[i], df2[j]]))
       nrow(x) == 0 ? missing : cor(x[1], x[2])
    end

    use_cor = method == CorPairwise() ? cor2 : cor1
    df2 = method == CorComplete() ? dropmissing(df) : df

    m = Matrix{Union{Float64, Missing}}(undef, ncol(df), ncol(df))
    for i in 1:ncol(df), j in i:ncol(df)
        m[i, j] = use_cor(i, j)
        m[j, i] = m[i, j]
    end
    m
end

max_stars_repo_path: Chapter07/05. Handling missing data/cor.jl


text:
using CSV, DataFrames, StringEncodings
# Let's encode age ranges with UnitRange, and transform conveniently with CombinedParsers
using CombinedParsers
jahrp = Either("unter 1 Jahr" => 0:1,
               map(a -> a:a,
                   Sequence(1, CombinedParsers.Numeric(Int),"-Jährige")),
               map(a -> a:(MAX-1), Sequence(1, CombinedParsers.Numeric(Int)," Jahre und mehr")),
               "Insgesamt" => 0:MAX-1)

@deprecate popdat() population()
function population()
    pdat = CSV.File(
        open(read, joinpath(datadep"Bevölkerung Deutschland","12411-0006.csv"), enc"ISO-8859-1");
        skipto=9, footerskip=4,
        delim=';', decimal=',',
        dateformat="dd.mm.yyy",
        header=["Stichtag", "Altersgruppe","Männlich","Weiblich","Insgesamt" ]) |> DataFrame
    pdat = pdat[pdat.Altersgruppe .!= "Insgesamt",:]

    pdat[!,"Altersgruppe"] = [ parse(jahrp, v)
                               for v in pdat[!,2] ]
    sort!(pdat, [2, 1])
    pdat
end

@deprecate popdat_bundesländer() population_bundesländer()
function population_bundesländer()
    pdat = vcat([ CSV.File(
        open(read, joinpath(datadep"Bevölkerung Deutschland Bundesländer","12411-0013$suffix.csv"),
             enc"ISO-8859-1");
        skipto=7, footerskip=4,
        delim=';', decimal=',',
        dateformat="dd.mm.yyy",
        header=5:6)  |> DataFrame for suffix in ["", "(1)", "(2)", "(3)", "(4)" ]]...)

    pdat = pdat[pdat[:,2] .!= "Insgesamt",:]

    pdat[!,"Altersgruppe"] = [
        parse(jahrp, v)
        for v in pdat[!,2]
            ];
    sort!(pdat, [1,:Altersgruppe])
    pdat
end

using Dates
"better idea with Intervals?"
iseq(x,y) = [ e.start >= y.start && ( e.stop < y.stop || e.stop == MAX)
              for e in x ]

export population
function population(popdat, deathdat=missing; jahr, geschlecht="Insgesamt", bundesland=missing, alter, monat = missing, kw=missing)
    if monat === missing && kw === missing
        r = if bundesland===missing
            # @assert geschlecht in geschlechter
            ## todo: boundary checks!
            rows = (popdat[:,1] .== Date(jahr-1,12,31)) .&
                iseq(popdat[:,:Altersgruppe], alter)
            sum(popdat[rows, geschlecht])
        else
            rows = (popdat[:,1] .== Date(jahr-1,12,31)) .&
                iseq(popdat[:,"Altersgruppe"], alter)
            ##print(popdat[rows, end])
            sum(popdat[rows, "$(bundesland)_$(lowercase(geschlecht))"])
        end
        convert(Int, round(r))
    elseif kw !== missing
        p1 = population(popdat,deathdat; alter=alter, geschlecht=geschlecht, bundesland = bundesland, jahr=jahr)
        p2 = population(popdat,deathdat; alter=alter, geschlecht=geschlecht, bundesland = bundesland, jahr=jahr+1)
        d = [ deaths(deathdat;bundesland = bundesland, alter=alter, geschlecht=geschlecht, jahr=jahr, kw=w) for w in 1:53 ]
        diedbefore = sum( d[1:(kw)] )
        adjust_kw = p2 == 0 ? 0 : (kw) * (p2-(p1-sum( d ))) / 52.0
        p1 - diedbefore + adjust_kw
    else
        p1 = population(popdat,deathdat; alter=alter, geschlecht=geschlecht, bundesland = bundesland, jahr=jahr)
        p2 = population(popdat,deathdat; alter=alter, geschlecht=geschlecht, bundesland = bundesland, jahr=jahr+1)
        d = [ deaths(deathdat;bundesland = bundesland, alter=alter, geschlecht=geschlecht, jahr=jahr, monat=m) for m in 1:12 ]
        diedbefore = sum( d[1:(monat)] )
        adjust_month = p2 == 0 ? 0 : (monat) * (p2-(p1-sum( d ))) / 12.0
        p1 - diedbefore + adjust_month
    end
end

max_stars_repo_name: gkappler/DeutschlandDaten

max_issues_repo_name: gkappler/DeutschlandDaten

max_forks_repo_name: gkappler/DeutschlandDaten


text:
using Knet,Images;
include(Pkg.dir("Knet","data","mnist.jl"))
global atype = gpu()>=0 ? KnetArray{Float32} : Array{Float32}

#A generic MLP function with customizable activation functions
function mlp(w,x;p=0.0,activation=leakyrelu,outputactivation=sigm)
    for i=1:2:length(w)
        x = w[i]*dropout(mat(x),p) .+ w[i+1]   # mat() used for flatten images to vector.
        i<length(w)-1 && (x = activation.(x))
    end
    return outputactivation.(x) #output layer
end
leakyrelu(x;α=Float32(0.2)) = max(0,x) + α*min(0,x) # LeakyRelu activation

global const 𝜀=Float32(1e-8) #  a small number prevent from getting NaN  in logs
D(w,x;p=0.0) = mlp(w,x;p=p)  #  Discriminator
G(w,z;p=0.0) = mlp(w,z;p=p)  #  Generator
𝑱d(𝗪d,x,Gz) = -mean(log.(D(𝗪d,x)+𝜀)+log.(1-D(𝗪d,Gz)+𝜀))/2 # Discriminator Loss
𝑱g(𝗪g, 𝗪d, z) = -mean(log.(D(𝗪d,G(𝗪g,z))+𝜀))             # Generator Loss
𝒩(input, batch) = atype(randn(Float32, input, batch))      # SampleNoise

∇d  = grad(𝑱d) # Discriminator gradient
∇g  = grad(𝑱g) # Generator gradient

function initweights(hidden,input, output)
    𝗪 = Any[];
    x = input
    for h in [hidden... output]
        push!(𝗪, atype(xavier(h,x)), atype(zeros(h, 1))) #FC Layers weights and bias
        x = h
    end
    return 𝗪  #return model params
end

function generate_and_save(𝗪,number,𝞗;fldr="generations/")
    Gz = G(𝗪[1],𝒩(𝞗[:ginp],number)) .> 0.5
    Gz = permutedims(reshape(Gz,(28,28,number)),(2,1,3))
    [save(fldr*string(i)*".png",Gray.(Gz[:,:,i])) for i=1:number]
end

#(if) train ? it updates model parameters : (else) it print losses
function runmodel(𝗪, data, 𝞗;dtst=nothing,optim=nothing,train=false,saveinterval=10)
    gloss=dloss=counter=0.0;
    B = 𝞗[:batchsize]
    for i=1:(train ? 𝞗[:epochs]:1)
        for (x,_) in data
            counter+=2B
            Gz = G(𝗪[1],𝒩(𝞗[:ginp],B)) #Fake Images
            train ? update!(𝗪[2], ∇d(𝗪[2],x,Gz), optim[2])      : (dloss += 2B*𝑱d(𝗪[2],x,Gz)) #update discriminator
            z=𝒩(𝞗[:ginp],2B) #Sample z from Noise
            train ? update!(𝗪[1], ∇g(𝗪[1], 𝗪[2], z), optim[1]) : (gloss += 2B*𝑱g(𝗪[1],𝗪[2],z)) #update generator
        end
        train ? runmodel(𝗪,dtst,𝞗;train=false) : println((gloss/counter,dloss/counter)) #Print average losses in each epoch
        i % saveinterval == 0 && generate_and_save(𝗪,10,𝞗)  # save 10 images to generations folder
    end
end

function main()
    𝞗 = Dict(:batchsize=>32,:epochs=>75,:ginp=>256,:genh=>[512],:disch=>[512],:optim=>Adam,:lr=>0.002);
    xtrn,ytrn,xtst,ytst = mnist()
    global dtrn = minibatch(xtrn, ytrn, 𝞗[:batchsize]; xtype=atype)
    global dtst = minibatch(xtst, ytst, 𝞗[:batchsize]; xtype=atype)
    𝗪 = (𝗪g,𝗪d)   = initweights(𝞗[:genh], 𝞗[:ginp], 784), initweights(𝞗[:disch], 784, 1)
    𝚶 = (𝚶pg,𝚶pd)  = optimizers(𝗪g,𝞗[:optim];lr=𝞗[:lr]), optimizers(𝗪d,𝞗[:optim];lr=𝞗[:lr])
    runmodel(𝗪, dtst, 𝞗;optim=𝚶, train=false) # initial losses
    runmodel(𝗪, dtrn, 𝞗;optim=𝚶,train=true, dtst=dtst)  # training
    𝗪,𝚶,𝞗,(dtrn,dtst)    # return weights,optimizers,options and dataset
end
main() #enjoy!


text:
struct vcovIID <: vcov
end

"""
    inference(N, K, resid, β, XX_inv, vcov)

Perform inference given model primitives. This function doesn't
use model objects but residuals, βs, and XXs so we can method dispatch
the appropriate "meat" in the vcv matrix depending on the model passed.

These are defined later:
e.g. TSLS will XX_inv = (X' P_Z X)^{-1} whereas OLS will just pass (X'X)^{-1}
"""
function inference(N::Int,
                   K::Int,
                   resid::Vector,
                   β::Vector,
                   XX_inv::Matrix,
                   X, # not used for IID but required for cluster/hetero
                      # Need a better design pattern here.
                   vcov::vcovIID)
    σ_sq = sum(resid.^2) / (N - K)
    vcov_matrix = XX_inv * σ_sq
    se = sqrt.(diag(vcov_matrix))
    pval = 2 .* cdf.(TDist(N - K), -abs.(β ./ se))

    return se, pval, σ_sq, vcov_matrix
end

function Ftest(fit_obj::LinearModelFit; signif_level = 0.05)
    resid_2 = fit_obj.modelfit.resid
    RSS_1 = sum((fit_obj.y .- mean(fit_obj.y)).^2)
    RSS_2 = sum(resid_2.^2)
    K = fit_obj.K
    N = fit_obj.N
    F = ((RSS_1 - RSS_2)/(K - 1)) / ((RSS_2)/(N-K))
    F_CV = quantile(FDist(K - 1, N - K), 1 - signif_level)
    reject = F > F_CV
    return F, F_CV, reject
end

# TODO copy code less here
function Ftest(fit_obj::LinearModelFit,
               restricted_fit_obj::LinearModelFit;
               signif_level = 0.05)
    resid_2 = fit_obj.modelfit.resid
    resid_1 = restricted_fit_obj.modelfit.resid
    RSS_1 = sum(resid_1.^2)
    RSS_2 = sum(resid_2.^2)
    K_2 = fit_obj.K
    K_1 = restricted_fit_obj.K
    N = fit_obj.N
    F = ((RSS_1 - RSS_2)/(K_2 - K_1)) / ((RSS_2)/(N-K_2))
    F_CV = quantile(FDist(K_2 - K_1, N - K_2), 1 - signif_level)
    reject = F > F_CV
    return F, F_CV, reject
end

# Now we just use method dispatch to call the necessary inference function
# whilst also handling different ways of passing arguments to inference.

# N.B. generic vcov::vcov. Hopefully we should only need to write a single
# vcov function  for each type of SE using primitives and everything in
# this file will dispatch away
"""
inference taking β and resid as arguments
"""
function inference(resid::Vector,
                   β::Vector, # todo change type to OLS?
                   fit::Union{LinearModel,LinearModelFit},
                   vcov::vcov)
    R = fit.R
    N = fit.N
    K = fit.K
    X = fit.X
    XX_inv = inv(cholesky(R' * R))
    return inference(N, K, resid, β, XX_inv, X, vcov)
end

"""
OLS inference using fitted model to pass arguments to above function
"""
function inference(fit::FittedOLSModel, vcov::vcov)
    resid = fit.modelfit.resid
    β = fit.modelfit.β
    return return inference(resid, β, fit, vcov)
end

"""
TSLS inference passing X ' P_Z X to meat matrix.
"""
function inference(resid::Vector,
                   β::Vector,
                   P_Z::Matrix,
                   fit::TSLSModel, vcov::vcov)
    X = fit.X
    XX_inv =  inv(X' * P_Z * X)
    return inference(fit.N, fit.K, resid, β, XX_inv, P_Z*X, vcov)
end
"""
TSLS using fitted model
"""
function inference(fit::FittedTSLSModel, vcov::vcov)
    X = fit.X
    P_Z = fit.P_Z
    XX_inv =  inv(X' * P_Z * X)
    return inference(fit.N, fit.K, fit.modelfit.resid, fit.modelfit.β, XX_inv, P_Z*X, vcov)
end

function inference(fit::FittedJIVEModel, vcov::vcov)
    X_jive = fit.X_jive
    X = fit.X
    XX_inv =  inv(X_jive' * X)
    return inference(fit.N, fit.K, fit.modelfit.resid, fit.modelfit.β, XX_inv, X_jive, vcov)
end

max_stars_repo_name: EdJeeOnGitHub/jeeves

max_issues_repo_name: EdJeeOnGitHub/jeeves

max_forks_repo_name: EdJeeOnGitHub/jeeves


text:
## 数据来源 https://tianchi.aliyun.com/competition/entrance/231573/introduction

## 加载包
using DataFrames,CSV
## 设置文件路径
data_raw_path="/home/jerrywang/Desktop/user_balance_table.csv"

## 读取CSV文件
raw=CSV.read(data_raw_path)
data=raw[:,[:report_date,:total_purchase_amt,:total_redeem_amt]]

## 时间区间为 20130701 至 20140831
Date=sort(collect(Set(data[:,:report_date])))

## 初始化申购和赎回
purchase=zeros(Int64,length(Date))
redeem=zeros(Int64,length(Date))
for i in 1:length(Date)
    d=Date[i]
    index=findall(x->isequal(x,d) ,data[:,:report_date])
    purchase[i]=sum(data[index,:total_purchase_amt])
    redeem[i]=sum(data[index,:total_redeem_amt])
end

T=1:427

Out_Data=DataFrame(date_order=1:427,Date=Date,purchase=purchase,redeem=redeem)
CSV.write("Data.csv",Out_Data)


text:
@inline function ode_interpolant(Θ,dt,y₀,y₁,k,cache::DiscreteConstantCache,idxs,T::Type{Val{0}})
  y₀
end

@inline function ode_interpolant!(out,Θ,dt,y₀,y₁,k,cache::DiscreteCache,idxs,T::Type{Val{0}})
  recursivecopy!(out,y₀)
end

"""
Hairer Norsett Wanner Solving Ordinary Differential Euations I - Nonstiff Problems Page 192
"""
@inline function ode_interpolant(Θ,dt,y₀,y₁,k,cache::DP5ConstantCache,idxs,T::Type{Val{0}})
  Θ1 = 1-Θ
  y₀ + dt*Θ*(k[1]+Θ1*(k[2]+Θ*(k[3]+Θ1*k[4])))
end

"""
Hairer Norsett Wanner Solving Ordinary Differential Euations I - Nonstiff Problems Page 192
"""
@inline function ode_interpolant!(out,Θ,dt,y₀,y₁,k,cache::DP5Cache,idxs,T::Type{Val{0}})
  Θ1 = 1-Θ
  @inbounds for (j,i) in enumerate(idxs)
    out[j] = y₀[i] + dt*Θ*(k[1][i]+Θ1*(k[2][i]+Θ*(k[3][i]+Θ1*k[4][i])))
  end
end

"""
Hairer Norsett Wanner Solving Ordinary Differential Euations I - Nonstiff Problems Page 192
"""
@inline function ode_interpolant!(out,Θ,dt,y₀,y₁,k,cache::DP5ThreadedCache,idxs,T::Type{Val{0}})
  Θ1 = 1-Θ
  @inbounds for (j,i) in enumerate(idxs)
    out[j] = y₀[i] + dt*Θ*(k[1][i]+Θ1*(k[2][i]+Θ*(k[3][i]+Θ1*k[4][i])))
  end
end

"""
From MATLAB ODE Suite by Shampine
"""
@inline function ode_interpolant(Θ,dt,y₀,y₁,k,cache::Rosenbrock23ConstantCache,idxs,T::Type{Val{0}})
  d = cache.d
  c1 = Θ*(1-Θ)/(1-2d)
  c2 = Θ*(Θ-2d)/(1-2d)
  y₀ + dt*(c1*k[1] + c2*k[2])
end

"""
From MATLAB ODE Suite by Shampine
"""
@inline function ode_interpolant!(out,Θ,dt,y₀,y₁,k,cache::Rosenbrock23Cache,idxs,T::Type{Val{0}})
  d = cache.tab.d
  c1 = Θ*(1-Θ)/(1-2d)
  c2 = Θ*(Θ-2d)/(1-2d)
  @inbounds for (j,i) in enumerate(idxs)
    out[j] = y₀[i] + dt*(c1*k[1][i] + c2*k[2][i])
  end
end

"""
From MATLAB ODE Suite by Shampine
"""
@inline function ode_interpolant(Θ,dt,y₀,y₁,k,cache::Rosenbrock32ConstantCache,idxs,T::Type{Val{0}})
  d = cache.d
  c1 = Θ*(1-Θ)/(1-2d)
  c2 = Θ*(Θ-2d)/(1-2d)
  y₀ + dt*(c1*k[1] + c2*k[2])
end

"""
From MATLAB ODE Suite by Shampine
"""
@inline function ode_interpolant!(out,Θ,dt,y₀,y₁,k,cache::Rosenbrock32Cache,idxs,T::Type{Val{0}})
  d = cache.tab.d
  c1 = Θ*(1-Θ)/(1-2d)
  c2 = Θ*(Θ-2d)/(1-2d)
  y₀ + dt*(c1*k[1] + c2*k[2])
  @inbounds for (j,i) in enumerate(idxs)
    out[j] = y₀[i] + dt*(c1*k[1][i] + c2*k[2][i])
  end
end

"""
Runge–Kutta pairs of order 5(4) satisfying only the first column
simplifying assumption

Ch. Tsitouras
"""
@inline function ode_interpolant!(out,Θ,dt,y₀,y₁,k,cache::Tsit5Cache,idxs,T::Type{Val{0}})
  b1Θ = -1.0530884977290216Θ * (Θ - 1.3299890189751412)*(Θ^2 - 1.4364028541716351Θ + 0.7139816917074209)
  b2Θ = 0.1017Θ^2 * (Θ^2 - 2.1966568338249754Θ + 1.2949852507374631)
  b3Θ = 2.490627285651252793Θ^2 * (Θ^2 - 2.38535645472061657Θ + 1.57803468208092486)
  b4Θ = -16.54810288924490272*(Θ - 1.21712927295533244)*(Θ - 0.61620406037800089)*Θ^2
  b5Θ = 47.37952196281928122*(Θ - 1.203071208372362603)*(Θ - 0.658047292653547382)*Θ^2
  b6Θ = -34.87065786149660974*(Θ - 1.2)*(Θ - 0.666666666666666667)*Θ^2
  b7Θ = 2.5*(Θ - 1)*(Θ - 0.6)*Θ^2
  @inbounds for (j,i) in enumerate(idxs)
    out[j] = y₀[i] + dt*(k[1][i]*b1Θ + k[2][i]*b2Θ + k[3][i]*b3Θ + k[4][i]*b4Θ + k[5][i]*b5Θ + k[6][i]*b6Θ + k[7][i]*b7Θ)
  end
end

"""
Runge–Kutta pairs of order 5(4) satisfying only the first column
simplifying assumption

Ch. Tsitouras
"""
@inline function ode_interpolant(Θ,dt,y₀,y₁,k,cache::Tsit5ConstantCache,idxs,T::Type{Val{0}})
  b1Θ = -1.0530884977290216Θ * (Θ - 1.3299890189751412)*(Θ^2 - 1.4364028541716351Θ + 0.7139816917074209)
  b2Θ = 0.1017Θ^2 * (Θ^2 - 2.1966568338249754Θ + 1.2949852507374631)
  b3Θ = 2.490627285651252793Θ^2 * (Θ^2 - 2.38535645472061657Θ + 1.57803468208092486)
  b4Θ = -16.54810288924490272*(Θ - 1.21712927295533244)*(Θ - 0.61620406037800089)*Θ^2
  b5Θ = 47.37952196281928122*(Θ - 1.203071208372362603)*(Θ - 0.658047292653547382)*Θ^2
  b6Θ = -34.87065786149660974*(Θ - 1.2)*(Θ - 0.666666666666666667)*Θ^2
  b7Θ = 2.5*(Θ - 1)*(Θ - 0.6)*Θ^2
  y₀ + dt*(k[1]*b1Θ + k[2]*b2Θ + k[3]*b3Θ + k[4]*b4Θ + k[5]*b5Θ + k[6]*b6Θ + k[7]*b7Θ)
end

"""
Coefficients taken from RKSuite
"""
@inline function ode_interpolant(Θ,dt,y₀,y₁,k,cache::BS5ConstantCache,idxs,T::Type{Val{0}})
  @unpack r016,r015,r014,r013,r012,r036,r035,r034,r033,r032,r046,r045,r044,r043,r042,r056,r055,r054,r053,r052,r066,r065,r064,r063,r062,r076,r075,r074,r073,r072,r086,r085,r084,r083,r082,r096,r095,r094,r093,r106,r105,r104,r103,r102,r116,r115,r114,r113,r112 = cache
  Θ2 = Θ^2
  Θ3 = Θ2*Θ
  Θ4 = Θ3*Θ
  Θ5 = Θ4*Θ
  Θ6 = Θ5*Θ
#=
  bΘ6 = r056*k[5] +
  ((r106*k[10]+r086*k[8])+
  (r076*k[7]+r066*k[6])) +
  ((r046*k[4]+r096*k[9])+
  (r036*k[3]+r116*k[11])+
  r016*k[1])

  bΘ5 = (r105*k[10]+r095*k[9]) +
   ((r075*k[7]+r065*k[6])+
   r055*k[5]) + ((r045*k[4]+
   r085*k[8])+(r035*k[3]+r115*k[11])+r015*k[1])

 bΘ4 = ((r044*k[4]+r084*k[8])+
         (r074*k[7]+r064*k[6])+
         r054*k[5]) + ((r104*k[10]+
         r094*k[9])+(r034*k[3]+
         r114*k[11])+r014*k[1])

  bΘ3 = r053*k[5] + r063*k[6] +
         ((r033*k[3]+r093*k[9])+
         (r103*k[10]+r083*k[8])+r013*
         k[1])+((r043*k[4]+r113*
         k[11])+r073*k[7])

  bΘ2 = r052*k[5] + ((r062*k[6]+
         r082*k[8])+r012*k[1]) +
         ((r032*k[3]+r092*k[9])+
         r102*k[10]) + ((r042*k[4]+
         r112*k[11])+r072*k[7])

  return y₀ + dt*Θ*k[1] + dt*(bΘ2*Θ2 + Θ3*bΘ3 + Θ4*bΘ4 + Θ5*bΘ5 + Θ6*bΘ6)
=#

  b1Θ =           r012*Θ2 + r013*Θ3 + r014*Θ4 + r015*Θ5 + r016*Θ6
  b3Θ =           r032*Θ2 + r033*Θ3 + r034*Θ4 + r035*Θ5 + r036*Θ6
  b4Θ =           r042*Θ2 + r043*Θ3 + r044*Θ4 + r045*Θ5 + r046*Θ6
  b5Θ =           r052*Θ2 + r053*Θ3 + r054*Θ4 + r055*Θ5 + r056*Θ6
  b6Θ =           r062*Θ2 + r063*Θ3 + r064*Θ4 + r065*Θ5 + r066*Θ6
  b7Θ =           r072*Θ2 + r073*Θ3 + r074*Θ4 + r075*Θ5 + r076*Θ6
  b8Θ =           r082*Θ2 + r083*Θ3 + r084*Θ4 + r085*Θ5 + r086*Θ6
  b9Θ =                     r093*Θ3 + r094*Θ4 + r095*Θ5 + r096*Θ6
  b10Θ=           r102*Θ2 + r103*Θ3 + r104*Θ4 + r105*Θ5 + r106*Θ6
  b11Θ=           r112*Θ2 + r113*Θ3 + r114*Θ4 + r115*Θ5 + r116*Θ6
  y₀ + dt*Θ*k[1] + dt*(k[1]*b1Θ  + k[3]*b3Θ + k[4]*b4Θ  + k[5]*b5Θ + k[6]*b6Θ + k[7]*b7Θ + k[8]*b8Θ + k[9]*b9Θ + k[10]*b10Θ + k[11]*b11Θ)
end

"""
Coefficients taken from RKSuite
"""
@inline function ode_interpolant!(out,Θ,dt,y₀,y₁,k,cache::BS5Cache,idxs,T::Type{Val{0}})
  @unpack r016,r015,r014,r013,r012,r036,r035,r034,r033,r032,r046,r045,r044,r043,r042,r056,r055,r054,r053,r052,r066,r065,r064,r063,r062,r076,r075,r074,r073,r072,r086,r085,r084,r083,r082,r096,r095,r094,r093,r106,r105,r104,r103,r102,r116,r115,r114,r113,r112 = cache.tab
  Θ2 = Θ^2
  Θ3 = Θ2*Θ
  Θ4 = Θ3*Θ
  Θ5 = Θ4*Θ
  Θ6 = Θ5*Θ
#=
  bΘ6 = r056*k[5] +
  ((r106*k[10]+r086*k[8])+
  (r076*k[7]+r066*k[6])) +
  ((r046*k[4]+r096*k[9])+
  (r036*k[3]+r116*k[11])+
  r016*k[1])

  bΘ5 = (r105*k[10]+r095*k[9]) +
   ((r075*k[7]+r065*k[6])+
   r055*k[5]) + ((r045*k[4]+
   r085*k[8])+(r035*k[3]+r115*k[11])+r015*k[1])

 bΘ4 = ((r044*k[4]+r084*k[8])+
         (r074*k[7]+r064*k[6])+
         r054*k[5]) + ((r104*k[10]+
         r094*k[9])+(r034*k[3]+
         r114*k[11])+r014*k[1])

  bΘ3 = r053*k[5] + r063*k[6] +
         ((r033*k[3]+r093*k[9])+
         (r103*k[10]+r083*k[8])+r013*
         k[1])+((r043*k[4]+r113*
         k[11])+r073*k[7])

  bΘ2 = r052*k[5] + ((r062*k[6]+
         r082*k[8])+r012*k[1]) +
         ((r032*k[3]+r092*k[9])+
         r102*k[10]) + ((r042*k[4]+
         r112*k[11])+r072*k[7])

  return y₀ + dt*Θ*k[1] + dt*(bΘ2*Θ2 + Θ3*bΘ3 + Θ4*bΘ4 + Θ5*bΘ5 + Θ6*bΘ6)
=#

  b1Θ =           r012*Θ2 + r013*Θ3 + r014*Θ4 + r015*Θ5 + r016*Θ6
  b3Θ =           r032*Θ2 + r033*Θ3 + r034*Θ4 + r035*Θ5 + r036*Θ6
  b4Θ =           r042*Θ2 + r043*Θ3 + r044*Θ4 + r045*Θ5 + r046*Θ6
  b5Θ =           r052*Θ2 + r053*Θ3 + r054*Θ4 + r055*Θ5 + r056*Θ6
  b6Θ =           r062*Θ2 + r063*Θ3 + r064*Θ4 + r065*Θ5 + r066*Θ6
  b7Θ =           r072*Θ2 + r073*Θ3 + r074*Θ4 + r075*Θ5 + r076*Θ6
  b8Θ =           r082*Θ2 + r083*Θ3 + r084*Θ4 + r085*Θ5 + r086*Θ6
  b9Θ =                     r093*Θ3 + r094*Θ4 + r095*Θ5 + r096*Θ6
  b10Θ=           r102*Θ2 + r103*Θ3 + r104*Θ4 + r105*Θ5 + r106*Θ6
  b11Θ=           r112*Θ2 + r113*Θ3 + r114*Θ4 + r115*Θ5 + r116*Θ6
  @inbounds for (j,i) in enumerate(idxs)
    out[j] = y₀[i] + dt*Θ*k[1][i] + dt*(k[1][i]*b1Θ  + k[3][i]*b3Θ + k[4][i]*b4Θ  + k[5][i]*b5Θ + k[6][i]*b6Θ + k[7][i]*b7Θ + k[8][i]*b8Θ + k[9][i]*b9Θ + k[10][i]*b10Θ + k[11][i]*b11Θ)
  end
end

"""

"""
@inline function ode_interpolant!(out,Θ,dt,y₀,y₁,k,cache::Vern6Cache,idxs,T::Type{Val{0}})
  @unpack r011,r012,r013,r014,r015,r016,r042,r043,r044,r045,r046,r052,r053,r054,r055,r056,r062,r063,r064,r065,r066,r072,r073,r074,r075,r076,r082,r083,r084,r085,r086,r092,r093,r094,r095,r096,r102,r103,r104,r105,r106,r112,r113,r114,r115,r116,r122,r123,r124,r125,r126 = cache.tab
  Θ2 = Θ^2
  Θ3 = Θ2*Θ
  Θ4 = Θ3*Θ
  Θ5 = Θ4*Θ
  Θ6 = Θ5*Θ
  b1Θ = r011*Θ + r012*Θ2 + r013*Θ3 + r014*Θ4 + r015*Θ5 + r016*Θ6
  b4Θ =          r042*Θ2 + r043*Θ3 + r044*Θ4 + r045*Θ5 + r046*Θ6
  b5Θ =          r052*Θ2 + r053*Θ3 + r054*Θ4 + r055*Θ5 + r056*Θ6
  b6Θ =          r062*Θ2 + r063*Θ3 + r064*Θ4 + r065*Θ5 + r066*Θ6
  b7Θ =          r072*Θ2 + r073*Θ3 + r074*Θ4 + r075*Θ5 + r076*Θ6
  b8Θ =          r082*Θ2 + r083*Θ3 + r084*Θ4 + r085*Θ5 + r086*Θ6
  b9Θ =          r092*Θ2 + r093*Θ3 + r094*Θ4 + r095*Θ5 + r096*Θ6
  b10Θ=          r102*Θ2 + r103*Θ3 + r104*Θ4 + r105*Θ5 + r106*Θ6
  b11Θ=          r112*Θ2 + r113*Θ3 + r114*Θ4 + r115*Θ5 + r116*Θ6
  b12Θ=          r122*Θ2 + r123*Θ3 + r124*Θ4 + r125*Θ5 + r126*Θ6
  @inbounds for (j,i) in enumerate(idxs)
    out[j] = y₀[i] + dt*(k[1][i]*b1Θ + k[4][i]*b4Θ + k[5][i]*b5Θ + k[6][i]*b6Θ + k[7][i]*b7Θ + k[8][i]*b8Θ + k[9][i]*b9Θ + k[10][i]*b10Θ + k[11][i]*b11Θ + k[12][i]*b12Θ)
  end
end

"""

"""
@inline function ode_interpolant(Θ,dt,y₀,y₁,k,cache::Vern6ConstantCache,idxs,T::Type{Val{0}})
  @unpack r011,r012,r013,r014,r015,r016,r042,r043,r044,r045,r046,r052,r053,r054,r055,r056,r062,r063,r064,r065,r066,r072,r073,r074,r075,r076,r082,r083,r084,r085,r086,r092,r093,r094,r095,r096,r102,r103,r104,r105,r106,r112,r113,r114,r115,r116,r122,r123,r124,r125,r126 = cache
  Θ2 = Θ^2
  Θ3 = Θ2*Θ
  Θ4 = Θ3*Θ
  Θ5 = Θ4*Θ
  Θ6 = Θ5*Θ
  b1Θ = r011*Θ + r012*Θ2 + r013*Θ3 + r014*Θ4 + r015*Θ5 + r016*Θ6
  b4Θ =          r042*Θ2 + r043*Θ3 + r044*Θ4 + r045*Θ5 + r046*Θ6
  b5Θ =          r052*Θ2 + r053*Θ3 + r054*Θ4 + r055*Θ5 + r056*Θ6
  b6Θ =          r062*Θ2 + r063*Θ3 + r064*Θ4 + r065*Θ5 + r066*Θ6
  b7Θ =          r072*Θ2 + r073*Θ3 + r074*Θ4 + r075*Θ5 + r076*Θ6
  b8Θ =          r082*Θ2 + r083*Θ3 + r084*Θ4 + r085*Θ5 + r086*Θ6
  b9Θ =          r092*Θ2 + r093*Θ3 + r094*Θ4 + r095*Θ5 + r096*Θ6
  b10Θ=          r102*Θ2 + r103*Θ3 + r104*Θ4 + r105*Θ5 + r106*Θ6
  b11Θ=          r112*Θ2 + r113*Θ3 + r114*Θ4 + r115*Θ5 + r116*Θ6
  b12Θ=          r122*Θ2 + r123*Θ3 + r124*Θ4 + r125*Θ5 + r126*Θ6
  y₀ + dt*(k[1]*b1Θ + k[4]*b4Θ + k[5]*b5Θ + k[6]*b6Θ + k[7]*b7Θ + k[8]*b8Θ + k[9]*b9Θ + k[10]*b10Θ + k[11]*b11Θ + k[12]*b12Θ)
end

"""

"""
@inline function ode_interpolant(Θ,dt,y₀,y₁,k,cache::Vern7ConstantCache,idxs,T::Type{Val{0}})
  @unpack r011,r012,r013,r014,r015,r016,r017,r042,r043,r044,r045,r046,r047,r052,r053,r054,r055,r056,r057,r062,r063,r064,r065,r066,r067,r072,r073,r074,r075,r076,r077,r082,r083,r084,r085,r086,r087,r092,r093,r094,r095,r096,r097,r112,r113,r114,r115,r116,r117,r122,r123,r124,r125,r126,r127,r132,r133,r134,r135,r136,r137,r142,r143,r144,r145,r146,r147,r152,r153,r154,r155,r156,r157,r162,r163,r164,r165,r166,r167 = cache
  Θ2 = Θ^2
  Θ3 = Θ2*Θ
  Θ4 = Θ3*Θ
  Θ5 = Θ4*Θ
  Θ6 = Θ5*Θ
  Θ7 = Θ6*Θ
  b1Θ = r011*Θ + r012*Θ2 + r013*Θ3 + r014*Θ4 + r015*Θ5 + r016*Θ6 + r017*Θ7
  b4Θ =          r042*Θ2 + r043*Θ3 + r044*Θ4 + r045*Θ5 + r046*Θ6 + r047*Θ7
  b5Θ =          r052*Θ2 + r053*Θ3 + r054*Θ4 + r055*Θ5 + r056*Θ6 + r057*Θ7
  b6Θ =          r062*Θ2 + r063*Θ3 + r064*Θ4 + r065*Θ5 + r066*Θ6 + r067*Θ7
  b7Θ =          r072*Θ2 + r073*Θ3 + r074*Θ4 + r075*Θ5 + r076*Θ6 + r077*Θ7
  b8Θ =          r082*Θ2 + r083*Θ3 + r084*Θ4 + r085*Θ5 + r086*Θ6 + r087*Θ7
  b9Θ =          r092*Θ2 + r093*Θ3 + r094*Θ4 + r095*Θ5 + r096*Θ6 + r097*Θ7
  b11Θ=          r112*Θ2 + r113*Θ3 + r114*Θ4 + r115*Θ5 + r116*Θ6 + r117*Θ7
  b12Θ=          r122*Θ2 + r123*Θ3 + r124*Θ4 + r125*Θ5 + r126*Θ6 + r127*Θ7
  b13Θ=          r132*Θ2 + r133*Θ3 + r134*Θ4 + r135*Θ5 + r136*Θ6 + r137*Θ7
  b14Θ=          r142*Θ2 + r143*Θ3 + r144*Θ4 + r145*Θ5 + r146*Θ6 + r147*Θ7
  b15Θ=          r152*Θ2 + r153*Θ3 + r154*Θ4 + r155*Θ5 + r156*Θ6 + r157*Θ7
  b16Θ=          r162*Θ2 + r163*Θ3 + r164*Θ4 + r165*Θ5 + r166*Θ6 + r167*Θ7
  y₀ + dt*(k[1]*b1Θ + k[4]*b4Θ + k[5]*b5Θ + k[6]*b6Θ + k[7]*b7Θ + k[8]*b8Θ + k[9]*b9Θ + k[11]*b11Θ + k[12]*b12Θ + k[13]*b13Θ + k[14]*b14Θ + k[15]*b15Θ + k[16]*b16Θ)
end

"""

"""
@inline function ode_interpolant!(out,Θ,dt,y₀,y₁,k,cache::Vern7Cache,idxs,T::Type{Val{0}})
  @unpack r011,r012,r013,r014,r015,r016,r017,r042,r043,r044,r045,r046,r047,r052,r053,r054,r055,r056,r057,r062,r063,r064,r065,r066,r067,r072,r073,r074,r075,r076,r077,r082,r083,r084,r085,r086,r087,r092,r093,r094,r095,r096,r097,r112,r113,r114,r115,r116,r117,r122,r123,r124,r125,r126,r127,r132,r133,r134,r135,r136,r137,r142,r143,r144,r145,r146,r147,r152,r153,r154,r155,r156,r157,r162,r163,r164,r165,r166,r167 = cache.tab
  Θ2 = Θ^2
  Θ3 = Θ2*Θ
  Θ4 = Θ3*Θ
  Θ5 = Θ4*Θ
  Θ6 = Θ5*Θ
  Θ7 = Θ6*Θ
  b1Θ = r011*Θ + r012*Θ2 + r013*Θ3 + r014*Θ4 + r015*Θ5 + r016*Θ6 + r017*Θ7
  b4Θ =          r042*Θ2 + r043*Θ3 + r044*Θ4 + r045*Θ5 + r046*Θ6 + r047*Θ7
  b5Θ =          r052*Θ2 + r053*Θ3 + r054*Θ4 + r055*Θ5 + r056*Θ6 + r057*Θ7
  b6Θ =          r062*Θ2 + r063*Θ3 + r064*Θ4 + r065*Θ5 + r066*Θ6 + r067*Θ7
  b7Θ =          r072*Θ2 + r073*Θ3 + r074*Θ4 + r075*Θ5 + r076*Θ6 + r077*Θ7
  b8Θ =          r082*Θ2 + r083*Θ3 + r084*Θ4 + r085*Θ5 + r086*Θ6 + r087*Θ7
  b9Θ =          r092*Θ2 + r093*Θ3 + r094*Θ4 + r095*Θ5 + r096*Θ6 + r097*Θ7
  b11Θ=          r112*Θ2 + r113*Θ3 + r114*Θ4 + r115*Θ5 + r116*Θ6 + r117*Θ7
  b12Θ=          r122*Θ2 + r123*Θ3 + r124*Θ4 + r125*Θ5 + r126*Θ6 + r127*Θ7
  b13Θ=          r132*Θ2 + r133*Θ3 + r134*Θ4 + r135*Θ5 + r136*Θ6 + r137*Θ7
  b14Θ=          r142*Θ2 + r143*Θ3 + r144*Θ4 + r145*Θ5 + r146*Θ6 + r147*Θ7
  b15Θ=          r152*Θ2 + r153*Θ3 + r154*Θ4 + r155*Θ5 + r156*Θ6 + r157*Θ7
  b16Θ=          r162*Θ2 + r163*Θ3 + r164*Θ4 + r165*Θ5 + r166*Θ6 + r167*Θ7
  @inbounds for (j,i) in enumerate(idxs)
    out[j] = y₀[i] + dt*(k[1][i]*b1Θ + k[4][i]*b4Θ + k[5][i]*b5Θ + k[6][i]*b6Θ + k[7][i]*b7Θ + k[8][i]*b8Θ + k[9][i]*b9Θ + k[11][i]*b11Θ + k[12][i]*b12Θ + k[13][i]*b13Θ + k[14][i]*b14Θ + k[15][i]*b15Θ + k[16][i]*b16Θ)
  end
end

"""

"""
@inline function ode_interpolant(Θ,dt,y₀,y₁,k,cache::Vern8ConstantCache,idxs,T::Type{Val{0}})
  @unpack r011,r012,r013,r014,r015,r016,r017,r018,r062,r063,r064,r065,r066,r067,r068,r072,r073,r074,r075,r076,r077,r078,r082,r083,r084,r085,r086,r087,r088,r092,r093,r094,r095,r096,r097,r098,r102,r103,r104,r105,r106,r107,r108,r112,r113,r114,r115,r116,r117,r118,r122,r123,r124,r125,r126,r127,r128,r142,r143,r144,r145,r146,r147,r148,r152,r153,r154,r155,r156,r157,r158,r162,r163,r164,r165,r166,r167,r168,r172,r173,r174,r175,r176,r177,r178,r182,r183,r184,r185,r186,r187,r188,r192,r193,r194,r195,r196,r197,r198,r202,r203,r204,r205,r206,r207,r208,r212,r213,r214,r215,r216,r217,r218 = cache
  Θ2 = Θ^2
  Θ3 = Θ2*Θ
  Θ4 = Θ3*Θ
  Θ5 = Θ4*Θ
  Θ6 = Θ5*Θ
  Θ7 = Θ6*Θ
  Θ8 = Θ7*Θ
  b1Θ = r011*Θ + r012*Θ2 + r013*Θ3 + r014*Θ4 + r015*Θ5 + r016*Θ6 + r017*Θ7 + r018*Θ8
  b6Θ =          r062*Θ2 + r063*Θ3 + r064*Θ4 + r065*Θ5 + r066*Θ6 + r067*Θ7 + r068*Θ8
  b7Θ =          r072*Θ2 + r073*Θ3 + r074*Θ4 + r075*Θ5 + r076*Θ6 + r077*Θ7 + r078*Θ8
  b8Θ =          r082*Θ2 + r083*Θ3 + r084*Θ4 + r085*Θ5 + r086*Θ6 + r087*Θ7 + r088*Θ8
  b9Θ =          r092*Θ2 + r093*Θ3 + r094*Θ4 + r095*Θ5 + r096*Θ6 + r097*Θ7 + r098*Θ8
  b10Θ=          r102*Θ2 + r103*Θ3 + r104*Θ4 + r105*Θ5 + r106*Θ6 + r107*Θ7 + r108*Θ8
  b11Θ=          r112*Θ2 + r113*Θ3 + r114*Θ4 + r115*Θ5 + r116*Θ6 + r117*Θ7 + r118*Θ8
  b12Θ=          r122*Θ2 + r123*Θ3 + r124*Θ4 + r125*Θ5 + r126*Θ6 + r127*Θ7 + r128*Θ8
  b14Θ=          r142*Θ2 + r143*Θ3 + r144*Θ4 + r145*Θ5 + r146*Θ6 + r147*Θ7 + r148*Θ8
  b15Θ=          r152*Θ2 + r153*Θ3 + r154*Θ4 + r155*Θ5 + r156*Θ6 + r157*Θ7 + r158*Θ8
  b16Θ=          r162*Θ2 + r163*Θ3 + r164*Θ4 + r165*Θ5 + r166*Θ6 + r167*Θ7 + r168*Θ8
  b17Θ=          r172*Θ2 + r173*Θ3 + r174*Θ4 + r175*Θ5 + r176*Θ6 + r177*Θ7 + r178*Θ8
  b18Θ=          r182*Θ2 + r183*Θ3 + r184*Θ4 + r185*Θ5 + r186*Θ6 + r187*Θ7 + r188*Θ8
  b19Θ=          r192*Θ2 + r193*Θ3 + r194*Θ4 + r195*Θ5 + r196*Θ6 + r197*Θ7 + r198*Θ8
  b20Θ=          r202*Θ2 + r203*Θ3 + r204*Θ4 + r205*Θ5 + r206*Θ6 + r207*Θ7 + r208*Θ8
  b21Θ=          r212*Θ2 + r213*Θ3 + r214*Θ4 + r215*Θ5 + r216*Θ6 + r217*Θ7 + r218*Θ8
  y₀ + dt*(k[1]*b1Θ + k[6]*b6Θ + k[7]*b7Θ + k[8]*b8Θ + k[9]*b9Θ + k[10]*b10Θ + k[11]*b11Θ + k[12]*b12Θ + k[14]*b14Θ + k[15]*b15Θ + k[16]*b16Θ + k[17]*b17Θ + k[18]*b18Θ + k[19]*b19Θ + k[20]*b20Θ + k[21]*b21Θ)
end

"""

"""
@inline function ode_interpolant!(out,Θ,dt,y₀,y₁,k,cache::Vern8Cache,idxs,T::Type{Val{0}})
  @unpack r011,r012,r013,r014,r015,r016,r017,r018,r062,r063,r064,r065,r066,r067,r068,r072,r073,r074,r075,r076,r077,r078,r082,r083,r084,r085,r086,r087,r088,r092,r093,r094,r095,r096,r097,r098,r102,r103,r104,r105,r106,r107,r108,r112,r113,r114,r115,r116,r117,r118,r122,r123,r124,r125,r126,r127,r128,r142,r143,r144,r145,r146,r147,r148,r152,r153,r154,r155,r156,r157,r158,r162,r163,r164,r165,r166,r167,r168,r172,r173,r174,r175,r176,r177,r178,r182,r183,r184,r185,r186,r187,r188,r192,r193,r194,r195,r196,r197,r198,r202,r203,r204,r205,r206,r207,r208,r212,r213,r214,r215,r216,r217,r218 = cache.tab
  Θ2 = Θ^2
  Θ3 = Θ2*Θ
  Θ4 = Θ3*Θ
  Θ5 = Θ4*Θ
  Θ6 = Θ5*Θ
  Θ7 = Θ6*Θ
  Θ8 = Θ7*Θ
  b1Θ = r011*Θ + r012*Θ2 + r013*Θ3 + r014*Θ4 + r015*Θ5 + r016*Θ6 + r017*Θ7 + r018*Θ8
  b6Θ =          r062*Θ2 + r063*Θ3 + r064*Θ4 + r065*Θ5 + r066*Θ6 + r067*Θ7 + r068*Θ8
  b7Θ =          r072*Θ2 + r073*Θ3 + r074*Θ4 + r075*Θ5 + r076*Θ6 + r077*Θ7 + r078*Θ8
  b8Θ =          r082*Θ2 + r083*Θ3 + r084*Θ4 + r085*Θ5 + r086*Θ6 + r087*Θ7 + r088*Θ8
  b9Θ =          r092*Θ2 + r093*Θ3 + r094*Θ4 + r095*Θ5 + r096*Θ6 + r097*Θ7 + r098*Θ8
  b10Θ=          r102*Θ2 + r103*Θ3 + r104*Θ4 + r105*Θ5 + r106*Θ6 + r107*Θ7 + r108*Θ8
  b11Θ=          r112*Θ2 + r113*Θ3 + r114*Θ4 + r115*Θ5 + r116*Θ6 + r117*Θ7 + r118*Θ8
  b12Θ=          r122*Θ2 + r123*Θ3 + r124*Θ4 + r125*Θ5 + r126*Θ6 + r127*Θ7 + r128*Θ8
  b14Θ=          r142*Θ2 + r143*Θ3 + r144*Θ4 + r145*Θ5 + r146*Θ6 + r147*Θ7 + r148*Θ8
  b15Θ=          r152*Θ2 + r153*Θ3 + r154*Θ4 + r155*Θ5 + r156*Θ6 + r157*Θ7 + r158*Θ8
  b16Θ=          r162*Θ2 + r163*Θ3 + r164*Θ4 + r165*Θ5 + r166*Θ6 + r167*Θ7 + r168*Θ8
  b17Θ=          r172*Θ2 + r173*Θ3 + r174*Θ4 + r175*Θ5 + r176*Θ6 + r177*Θ7 + r178*Θ8
  b18Θ=          r182*Θ2 + r183*Θ3 + r184*Θ4 + r185*Θ5 + r186*Θ6 + r187*Θ7 + r188*Θ8
  b19Θ=          r192*Θ2 + r193*Θ3 + r194*Θ4 + r195*Θ5 + r196*Θ6 + r197*Θ7 + r198*Θ8
  b20Θ=          r202*Θ2 + r203*Θ3 + r204*Θ4 + r205*Θ5 + r206*Θ6 + r207*Θ7 + r208*Θ8
  b21Θ=          r212*Θ2 + r213*Θ3 + r214*Θ4 + r215*Θ5 + r216*Θ6 + r217*Θ7 + r218*Θ8
  @inbounds for (j,i) in enumerate(idxs)
    out[j] = y₀[i] + dt*(k[1][i]*b1Θ + k[6][i]*b6Θ + k[7][i]*b7Θ + k[8][i]*b8Θ + k[9][i]*b9Θ + k[10][i]*b10Θ + k[11][i]*b11Θ + k[12][i]*b12Θ + k[14][i]*b14Θ + k[15][i]*b15Θ + k[16][i]*b16Θ + k[17][i]*b17Θ + k[18][i]*b18Θ + k[19][i]*b19Θ + k[20][i]*b20Θ + k[21][i]*b21Θ)
  end
end

"""

"""
@inline function ode_interpolant(Θ,dt,y₀,y₁,k,cache::Vern9ConstantCache,idxs,T::Type{Val{0}})
  @unpack r011,r012,r013,r014,r015,r016,r017,r018,r019,r082,r083,r084,r085,r086,r087,r088,r089,r092,r093,r094,r095,r096,r097,r098,r099,r102,r103,r104,r105,r106,r107,r108,r109,r112,r113,r114,r115,r116,r117,r118,r119,r122,r123,r124,r125,r126,r127,r128,r129,r132,r133,r134,r135,r136,r137,r138,r139,r142,r143,r144,r145,r146,r147,r148,r149,r152,r153,r154,r155,r156,r157,r158,r159,r172,r173,r174,r175,r176,r177,r178,r179,r182,r183,r184,r185,r186,r187,r188,r189,r192,r193,r194,r195,r196,r197,r198,r199,r202,r203,r204,r205,r206,r207,r208,r209,r212,r213,r214,r215,r216,r217,r218,r219,r222,r223,r224,r225,r226,r227,r228,r229,r232,r233,r234,r235,r236,r237,r238,r239,r242,r243,r244,r245,r246,r247,r248,r249,r252,r253,r254,r255,r256,r257,r258,r259,r262,r263,r264,r265,r266,r267,r268,r269 = cache
  Θ2 = Θ^2
  Θ3 = Θ2*Θ
  Θ4 = Θ3*Θ
  Θ5 = Θ4*Θ
  Θ6 = Θ5*Θ
  Θ7 = Θ6*Θ
  Θ8 = Θ7*Θ
  Θ9 = Θ8*Θ
  b1Θ = r011*Θ + r012*Θ2 + r013*Θ3 + r014*Θ4 + r015*Θ5 + r016*Θ6 + r017*Θ7 + r018*Θ8 + r019*Θ9
  b8Θ =          r082*Θ2 + r083*Θ3 + r084*Θ4 + r085*Θ5 + r086*Θ6 + r087*Θ7 + r088*Θ8 + r089*Θ9
  b9Θ =          r092*Θ2 + r093*Θ3 + r094*Θ4 + r095*Θ5 + r096*Θ6 + r097*Θ7 + r098*Θ8 + r099*Θ9
  b10Θ=          r102*Θ2 + r103*Θ3 + r104*Θ4 + r105*Θ5 + r106*Θ6 + r107*Θ7 + r108*Θ8 + r109*Θ9
  b11Θ=          r112*Θ2 + r113*Θ3 + r114*Θ4 + r115*Θ5 + r116*Θ6 + r117*Θ7 + r118*Θ8 + r119*Θ9
  b12Θ=          r122*Θ2 + r123*Θ3 + r124*Θ4 + r125*Θ5 + r126*Θ6 + r127*Θ7 + r128*Θ8 + r129*Θ9
  b13Θ=          r132*Θ2 + r133*Θ3 + r134*Θ4 + r135*Θ5 + r136*Θ6 + r137*Θ7 + r138*Θ8 + r139*Θ9
  b14Θ=          r142*Θ2 + r143*Θ3 + r144*Θ4 + r145*Θ5 + r146*Θ6 + r147*Θ7 + r148*Θ8 + r149*Θ9
  b15Θ=          r152*Θ2 + r153*Θ3 + r154*Θ4 + r155*Θ5 + r156*Θ6 + r157*Θ7 + r158*Θ8 + r159*Θ9
  b17Θ=          r172*Θ2 + r173*Θ3 + r174*Θ4 + r175*Θ5 + r176*Θ6 + r177*Θ7 + r178*Θ8 + r179*Θ9
  b18Θ=          r182*Θ2 + r183*Θ3 + r184*Θ4 + r185*Θ5 + r186*Θ6 + r187*Θ7 + r188*Θ8 + r189*Θ9
  b19Θ=          r192*Θ2 + r193*Θ3 + r194*Θ4 + r195*Θ5 + r196*Θ6 + r197*Θ7 + r198*Θ8 + r199*Θ9
  b20Θ=          r202*Θ2 + r203*Θ3 + r204*Θ4 + r205*Θ5 + r206*Θ6 + r207*Θ7 + r208*Θ8 + r209*Θ9
  b21Θ=          r212*Θ2 + r213*Θ3 + r214*Θ4 + r215*Θ5 + r216*Θ6 + r217*Θ7 + r218*Θ8 + r219*Θ9
  b22Θ=          r222*Θ2 + r223*Θ3 + r224*Θ4 + r225*Θ5 + r226*Θ6 + r227*Θ7 + r228*Θ8 + r229*Θ9
  b23Θ=          r232*Θ2 + r233*Θ3 + r234*Θ4 + r235*Θ5 + r236*Θ6 + r237*Θ7 + r238*Θ8 + r239*Θ9
  b24Θ=          r242*Θ2 + r243*Θ3 + r244*Θ4 + r245*Θ5 + r246*Θ6 + r247*Θ7 + r248*Θ8 + r249*Θ9
  b25Θ=          r252*Θ2 + r253*Θ3 + r254*Θ4 + r255*Θ5 + r256*Θ6 + r257*Θ7 + r258*Θ8 + r259*Θ9
  b26Θ=          r262*Θ2 + r263*Θ3 + r264*Θ4 + r265*Θ5 + r266*Θ6 + r267*Θ7 + r268*Θ8 + r269*Θ9
  y₀ + dt*(k[1]*b1Θ + k[8]*b8Θ + k[9]*b9Θ + k[10]*b10Θ + k[11]*b11Θ + k[12]*b12Θ + k[13]*b13Θ + k[14]*b14Θ + k[15]*b15Θ + k[17]*b17Θ + k[18]*b18Θ + k[19]*b19Θ + k[20]*b20Θ + k[21]*b21Θ + k[22]*b22Θ + k[23]*b23Θ + k[24]*b24Θ + k[25]*b25Θ + k[26]*b26Θ)
end

"""

"""
@inline function ode_interpolant!(out,Θ,dt,y₀,y₁,k,cache::Vern9Cache,idxs,T::Type{Val{0}})
  @unpack r011,r012,r013,r014,r015,r016,r017,r018,r019,r082,r083,r084,r085,r086,r087,r088,r089,r092,r093,r094,r095,r096,r097,r098,r099,r102,r103,r104,r105,r106,r107,r108,r109,r112,r113,r114,r115,r116,r117,r118,r119,r122,r123,r124,r125,r126,r127,r128,r129,r132,r133,r134,r135,r136,r137,r138,r139,r142,r143,r144,r145,r146,r147,r148,r149,r152,r153,r154,r155,r156,r157,r158,r159,r172,r173,r174,r175,r176,r177,r178,r179,r182,r183,r184,r185,r186,r187,r188,r189,r192,r193,r194,r195,r196,r197,r198,r199,r202,r203,r204,r205,r206,r207,r208,r209,r212,r213,r214,r215,r216,r217,r218,r219,r222,r223,r224,r225,r226,r227,r228,r229,r232,r233,r234,r235,r236,r237,r238,r239,r242,r243,r244,r245,r246,r247,r248,r249,r252,r253,r254,r255,r256,r257,r258,r259,r262,r263,r264,r265,r266,r267,r268,r269 = cache.tab
  Θ2 = Θ^2
  Θ3 = Θ2*Θ
  Θ4 = Θ3*Θ
  Θ5 = Θ4*Θ
  Θ6 = Θ5*Θ
  Θ7 = Θ6*Θ
  Θ8 = Θ7*Θ
  Θ9 = Θ8*Θ
  b1Θ = r011*Θ + r012*Θ2 + r013*Θ3 + r014*Θ4 + r015*Θ5 + r016*Θ6 + r017*Θ7 + r018*Θ8 + r019*Θ9
  b8Θ =          r082*Θ2 + r083*Θ3 + r084*Θ4 + r085*Θ5 + r086*Θ6 + r087*Θ7 + r088*Θ8 + r089*Θ9
  b9Θ =          r092*Θ2 + r093*Θ3 + r094*Θ4 + r095*Θ5 + r096*Θ6 + r097*Θ7 + r098*Θ8 + r099*Θ9
  b10Θ=          r102*Θ2 + r103*Θ3 + r104*Θ4 + r105*Θ5 + r106*Θ6 + r107*Θ7 + r108*Θ8 + r109*Θ9
  b11Θ=          r112*Θ2 + r113*Θ3 + r114*Θ4 + r115*Θ5 + r116*Θ6 + r117*Θ7 + r118*Θ8 + r119*Θ9
  b12Θ=          r122*Θ2 + r123*Θ3 + r124*Θ4 + r125*Θ5 + r126*Θ6 + r127*Θ7 + r128*Θ8 + r129*Θ9
  b13Θ=          r132*Θ2 + r133*Θ3 + r134*Θ4 + r135*Θ5 + r136*Θ6 + r137*Θ7 + r138*Θ8 + r139*Θ9
  b14Θ=          r142*Θ2 + r143*Θ3 + r144*Θ4 + r145*Θ5 + r146*Θ6 + r147*Θ7 + r148*Θ8 + r149*Θ9
  b15Θ=          r152*Θ2 + r153*Θ3 + r154*Θ4 + r155*Θ5 + r156*Θ6 + r157*Θ7 + r158*Θ8 + r159*Θ9
  b17Θ=          r172*Θ2 + r173*Θ3 + r174*Θ4 + r175*Θ5 + r176*Θ6 + r177*Θ7 + r178*Θ8 + r179*Θ9
  b18Θ=          r182*Θ2 + r183*Θ3 + r184*Θ4 + r185*Θ5 + r186*Θ6 + r187*Θ7 + r188*Θ8 + r189*Θ9
  b19Θ=          r192*Θ2 + r193*Θ3 + r194*Θ4 + r195*Θ5 + r196*Θ6 + r197*Θ7 + r198*Θ8 + r199*Θ9
  b20Θ=          r202*Θ2 + r203*Θ3 + r204*Θ4 + r205*Θ5 + r206*Θ6 + r207*Θ7 + r208*Θ8 + r209*Θ9
  b21Θ=          r212*Θ2 + r213*Θ3 + r214*Θ4 + r215*Θ5 + r216*Θ6 + r217*Θ7 + r218*Θ8 + r219*Θ9
  b22Θ=          r222*Θ2 + r223*Θ3 + r224*Θ4 + r225*Θ5 + r226*Θ6 + r227*Θ7 + r228*Θ8 + r229*Θ9
  b23Θ=          r232*Θ2 + r233*Θ3 + r234*Θ4 + r235*Θ5 + r236*Θ6 + r237*Θ7 + r238*Θ8 + r239*Θ9
  b24Θ=          r242*Θ2 + r243*Θ3 + r244*Θ4 + r245*Θ5 + r246*Θ6 + r247*Θ7 + r248*Θ8 + r249*Θ9
  b25Θ=          r252*Θ2 + r253*Θ3 + r254*Θ4 + r255*Θ5 + r256*Θ6 + r257*Θ7 + r258*Θ8 + r259*Θ9
  b26Θ=          r262*Θ2 + r263*Θ3 + r264*Θ4 + r265*Θ5 + r266*Θ6 + r267*Θ7 + r268*Θ8 + r269*Θ9
  @inbounds for (j,i) in enumerate(idxs)
    out[j] = y₀[i] + dt*(k[1][i]*b1Θ + k[8][i]*b8Θ + k[9][i]*b9Θ + k[10][i]*b10Θ + k[11][i]*b11Θ + k[12][i]*b12Θ + k[13][i]*b13Θ + k[14][i]*b14Θ + k[15][i]*b15Θ + k[17][i]*b17Θ + k[18][i]*b18Θ + k[19][i]*b19Θ + k[20][i]*b20Θ + k[21][i]*b21Θ + k[22][i]*b22Θ + k[23][i]*b23Θ + k[24][i]*b24Θ + k[25][i]*b25Θ + k[26][i]*b26Θ)
  end
end

"""

"""
@inline function ode_interpolant(Θ,dt,y₀,y₁,k,cache::DP8ConstantCache,idxs,T::Type{Val{0}})
  Θ1 = 1-Θ
  conpar = k[4] + Θ*(k[5] + Θ1*(k[6]+Θ*k[7]))
  y₀ + dt*Θ*(k[1] + Θ1*(k[2] + Θ*(k[3]+Θ1*conpar)))
end

"""

"""
@inline function ode_interpolant!(out,Θ,dt,y₀,y₁,k,cache::DP8Cache,idxs,T::Type{Val{0}})
  Θ1 = 1-Θ
  @inbounds for (j,i) in enumerate(idxs)
    out[j] = y₀[i] + dt*Θ*(k[1][i] + Θ1*(k[2][i] + Θ*(k[3][i]+Θ1*(k[4][i] + Θ*(k[5][i] + Θ1*(k[6][i]+Θ*k[7][i]))))))
  end
end


text:
function nlp_cvx_104_010(
    optimizer,
    objective_tol,
    primal_tol,
    dual_tol,
    termination_target = TERMINATION_TARGET_LOCAL,
    primal_target = PRIMAL_TARGET_LOCAL,
)
    # Test Goals:
    # - redundant nonlinear constraint (feasible set is inside the circle constraint)
    # Variants
    #   010 - redundant constraint

    model = Model(optimizer)

    @variable(model, x)
    @variable(model, y)

    @objective(model, Min, -x)
    @NLconstraint(model, x^2 <= y)
    @NLconstraint(model, -x^2 + 1 >= y)
    @NLconstraint(model, x^2 + (y - 0.5)^2 <= 1.0)

    optimize!(model)

    check_status(model, FEASIBLE_PROBLEM, termination_target, primal_target)
    check_objective(model, -1 / sqrt(2), tol = objective_tol)
    return check_solution([x, y], [1 / sqrt(2), 1 / 2], tol = primal_tol)
end


text:
#=
# This experiment file establishes the initial conditions, boundary conditions,
# source terms and simulation parameters (domain size + resolution) for the
# BOMEX LES case. The set of parameters presented in the `master` branch copy
# include those that have passed offline tests at the full simulation time of
# 6 hours. Suggested offline tests included plotting horizontal-domain averages
# of key properties (see AtmosDiagnostics). The timestepper configuration is in
# `src/Driver/solver_configs.jl` while the `AtmosModel` defaults can be found in
# `src/Atmos/Model/AtmosModel.jl` and `src/Driver/driver_configs.jl`
#
# This setup works in both Float32 and Float64 precision. `FT`
#
# To simulate the full 6 hour experiment, change `timeend` to (3600*6) and type in
#
# julia --project experiments/AtmosLES/bomex.jl
#
# See `src/Driver/driver_configs.jl` for additional flags (e.g. VTK, diagnostics,
# update-interval, output directory settings)
#
# Upcoming changes:
# 1) Atomic sources
# 2) Improved boundary conditions
# 3) Collapsed experiment design
# 4) Updates to generally keep this in sync with master

@article{doi:10.1175/1520-0469(2003)60<1201:ALESIS>2.0.CO;2,
author = {Siebesma, A. Pier and Bretherton,
          Christopher S. and Brown,
          Andrew and Chlond,
          Andreas and Cuxart,
          Joan and Duynkerke,
          Peter G. and Jiang,
          Hongli and Khairoutdinov,
          Marat and Lewellen,
          David and Moeng,
          Chin-Hoh and Sanchez,
          Enrique and Stevens,
          Bjorn and Stevens,
          David E.},
title = {A Large Eddy Simulation Intercomparison Study of Shallow Cumulus Convection},
journal = {Journal of the Atmospheric Sciences},
volume = {60},
number = {10},
pages = {1201-1219},
year = {2003},
doi = {10.1175/1520-0469(2003)60<1201:ALESIS>2.0.CO;2},
URL = {https://journals.ametsoc.org/doi/abs/10.1175/1520-0469%282003%2960%3C1201%3AALESIS%3E2.0.CO%3B2},
eprint = {https://journals.ametsoc.org/doi/pdf/10.1175/1520-0469%282003%2960%3C1201%3AALESIS%3E2.0.CO%3B2}
=#

using Distributions
using Random
using StaticArrays
using Test
using DocStringExtensions
using LinearAlgebra

using CLIMA
using CLIMA.Atmos
using CLIMA.ConfigTypes
using CLIMA.DGmethods.NumericalFluxes
using CLIMA.Diagnostics
using CLIMA.GenericCallbacks
using CLIMA.Mesh.Filters
using CLIMA.ODESolvers
using CLIMA.MoistThermodynamics
using CLIMA.VariableTemplates

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

import CLIMA.DGmethods: vars_state, vars_aux
import CLIMA.Atmos: source!, atmos_source!, altitude
import CLIMA.Atmos: flux_diffusive!, thermo_state

"""
  Bomex Geostrophic Forcing (Source)
"""
struct BomexGeostrophic{FT} <: Source
    "Coriolis parameter [s⁻¹]"
    f_coriolis::FT
    "Eastward geostrophic velocity `[m/s]` (Base)"
    u_geostrophic::FT
    "Eastward geostrophic velocity `[m/s]` (Slope)"
    u_slope::FT
    "Northward geostrophic velocity `[m/s]`"
    v_geostrophic::FT
end
function atmos_source!(
    s::BomexGeostrophic,
    atmos::AtmosModel,
    source::Vars,
    state::Vars,
    diffusive::Vars,
    aux::Vars,
    t::Real,
    direction,
)

    f_coriolis = s.f_coriolis
    u_geostrophic = s.u_geostrophic
    u_slope = s.u_slope
    v_geostrophic = s.v_geostrophic

    z = altitude(atmos, aux)
    # Note z dependence of eastward geostrophic velocity
    u_geo = SVector(u_geostrophic + u_slope * z, v_geostrophic, 0)
    ẑ = vertical_unit_vector(atmos, aux)
    fkvector = f_coriolis * ẑ
    # Accumulate sources
    source.ρu -= fkvector × (state.ρu .- state.ρ * u_geo)
    return nothing
end

"""
  Bomex Sponge (Source)
"""
struct BomexSponge{FT} <: Source
    "Maximum domain altitude (m)"
    z_max::FT
    "Altitude at with sponge starts (m)"
    z_sponge::FT
    "Sponge Strength 0 ⩽ α_max ⩽ 1"
    α_max::FT
    "Sponge exponent"
    γ::FT
    "Eastward geostrophic velocity `[m/s]` (Base)"
    u_geostrophic::FT
    "Eastward geostrophic velocity `[m/s]` (Slope)"
    u_slope::FT
    "Northward geostrophic velocity `[m/s]`"
    v_geostrophic::FT
end
function atmos_source!(
    s::BomexSponge,
    atmos::AtmosModel,
    source::Vars,
    state::Vars,
    diffusive::Vars,
    aux::Vars,
    t::Real,
    direction,
)

    z_max = s.z_max
    z_sponge = s.z_sponge
    α_max = s.α_max
    γ = s.γ
    u_geostrophic = s.u_geostrophic
    u_slope = s.u_slope
    v_geostrophic = s.v_geostrophic

    z = altitude(atmos, aux)
    u_geo = SVector(u_geostrophic + u_slope * z, v_geostrophic, 0)
    ẑ = vertical_unit_vector(atmos, aux)
    # Accumulate sources
    if z_sponge <= z
        r = (z - z_sponge) / (z_max - z_sponge)
        β_sponge = α_max * sinpi(r / 2)^s.γ
        source.ρu -= β_sponge * (state.ρu .- state.ρ * u_geo)
    end
    return nothing
end

"""
  BomexTendencies (Source)
Moisture, Temperature and Subsidence tendencies
"""
struct BomexTendencies{FT} <: Source
    "Advection tendency in total moisture `[s⁻¹]`"
    ∂qt∂t_peak::FT
    "Lower extent of piecewise profile (moisture term) `[m]`"
    zl_moisture::FT
    "Upper extent of piecewise profile (moisture term) `[m]`"
    zh_moisture::FT
    "Cooling rate `[K/s]`"
    ∂θ∂t_peak::FT
    "Lower extent of piecewise profile (subsidence term) `[m]`"
    zl_sub::FT
    "Upper extent of piecewise profile (subsidence term) `[m]`"
    zh_sub::FT
    "Subsidence peak velocity"
    w_sub::FT
    "Max height in domain"
    z_max::FT
end
function atmos_source!(
    s::BomexTendencies,
    atmos::AtmosModel,
    source::Vars,
    state::Vars,
    diffusive::Vars,
    aux::Vars,
    t::Real,
    direction,
)
    FT = eltype(state)
    ρ = state.ρ
    z = altitude(atmos, aux)
    _e_int_v0 = FT(e_int_v0(atmos.param_set))

    # Establish thermodynamic state
    TS = thermo_state(atmos, state, aux)

    # Moisture tendencey (sink term)
    # Temperature tendency (Radiative cooling)
    # Large scale subsidence
    # Unpack struct
    zl_moisture = s.zl_moisture
    zh_moisture = s.zh_moisture
    z_max = s.z_max
    zl_sub = s.zl_sub
    zh_sub = s.zh_sub
    zl_temperature = zl_sub
    w_sub = s.w_sub
    ∂qt∂t_peak = s.∂qt∂t_peak
    ∂θ∂t_peak = s.∂θ∂t_peak
    k̂ = vertical_unit_vector(atmos, aux)

    # Thermodynamic state identification
    q_pt = PhasePartition(TS)
    cvm = cv_m(TS)

    # Piecewise term for moisture tendency
    linscale_moisture = (z - zl_moisture) / (zh_moisture - zl_moisture)
    if z <= zl_moisture
        ρ∂qt∂t = ρ * ∂qt∂t_peak
    elseif zl_moisture < z <= zh_moisture
        ρ∂qt∂t = ρ * (∂qt∂t_peak - ∂qt∂t_peak * linscale_moisture)
    else
        ρ∂qt∂t = -zero(FT)
    end

    # Piecewise term for internal energy tendency
    linscale_temp = (z - zl_sub) / (z_max - zl_sub)
    if z <= zl_sub
        ρ∂θ∂t = ρ * ∂θ∂t_peak
    elseif zl_temperature < z <= z_max
        ρ∂θ∂t = ρ * (∂θ∂t_peak - ∂θ∂t_peak * linscale_temp)
    else
        ρ∂θ∂t = -zero(FT)
    end

    # Piecewise terms for subsidence
    linscale_sub = (z - zl_sub) / (zh_sub - zl_sub)
    w_s = -zero(FT)
    if z <= zl_sub
        w_s = -zero(FT) + z * (w_sub) / (zl_sub)
    elseif zl_sub < z <= zh_sub
        w_s = w_sub - (w_sub) * linscale_sub
    else
        w_s = -zero(FT)
    end

    # Collect Sources
    source.moisture.ρq_tot += ρ∂qt∂t
    source.ρe += cvm * ρ∂θ∂t * exner(TS) + _e_int_v0 * ρ∂qt∂t
    source.ρe -= ρ * w_s * dot(k̂, diffusive.∇h_tot)
    source.moisture.ρq_tot -= ρ * w_s * dot(k̂, diffusive.moisture.∇q_tot)
    return nothing
end

"""
  Initial Condition for BOMEX LES
"""
seed = MersenneTwister(0)
function init_bomex!(bl, state, aux, (x, y, z), t)
    # This experiment runs the BOMEX LES Configuration
    # (Shallow cumulus cloud regime)
    # x,y,z imply eastward, northward and altitude coordinates in `[m]`

    # Problem floating point precision
    FT = eltype(state)

    P_sfc::FT = 1.015e5 # Surface air pressure
    qg::FT = 22.45e-3 # Total moisture at surface
    q_pt_sfc = PhasePartition(qg) # Surface moisture partitioning
    Rm_sfc = gas_constant_air(q_pt_sfc, bl.param_set) # Moist gas constant
    θ_liq_sfc = FT(299.1) # Prescribed θ_liq at surface
    T_sfc = FT(300.4) # Surface temperature
    _grav = FT(grav(bl.param_set))

    # Initialise speeds [u = Eastward, v = Northward, w = Vertical]
    u::FT = 0
    v::FT = 0
    w::FT = 0

    # Prescribed altitudes for piece-wise profile construction
    zl1::FT = 520
    zl2::FT = 1480
    zl3::FT = 2000
    zl4::FT = 3000

    # Assign piecewise quantities to θ_liq and q_tot
    θ_liq::FT = 0
    q_tot::FT = 0

    # Piecewise functions for potential temperature and total moisture
    if FT(0) <= z <= zl1
        # Well mixed layer
        θ_liq = 298.7
        q_tot = 17.0 + (z / zl1) * (16.3 - 17.0)
    elseif z > zl1 && z <= zl2
        # Conditionally unstable layer
        θ_liq = 298.7 + (z - zl1) * (302.4 - 298.7) / (zl2 - zl1)
        q_tot = 16.3 + (z - zl1) * (10.7 - 16.3) / (zl2 - zl1)
    elseif z > zl2 && z <= zl3
        # Absolutely stable inversion
        θ_liq = 302.4 + (z - zl2) * (308.2 - 302.4) / (zl3 - zl2)
        q_tot = 10.7 + (z - zl2) * (4.2 - 10.7) / (zl3 - zl2)
    else
        θ_liq = 308.2 + (z - zl3) * (311.85 - 308.2) / (zl4 - zl3)
        q_tot = 4.2 + (z - zl3) * (3.0 - 4.2) / (zl4 - zl3)
    end

    # Set velocity profiles - piecewise profile for u
    zlv::FT = 700
    if z <= zlv
        u = -8.75
    else
        u = -8.75 + (z - zlv) * (-4.61 + 8.75) / (zl4 - zlv)
    end

    # Convert total specific humidity to kg/kg
    q_tot /= 1000
    # Scale height based on surface parameters
    H = Rm_sfc * T_sfc / _grav
    # Pressure based on scale height
    P = P_sfc * exp(-z / H)

    # Establish thermodynamic state and moist phase partitioning
    TS = LiquidIcePotTempSHumEquil_given_pressure(θ_liq, P, q_tot, bl.param_set)
    T = air_temperature(TS)
    ρ = air_density(TS)
    q_pt = PhasePartition(TS)

    # Compute momentum contributions
    ρu = ρ * u
    ρv = ρ * v
    ρw = ρ * w

    # Compute energy contributions
    e_kin = FT(1 // 2) * (u^2 + v^2 + w^2)
    e_pot = _grav * z
    ρe_tot = ρ * total_energy(e_kin, e_pot, TS)

    # Assign initial conditions for prognostic state variables
    state.ρ = ρ
    state.ρu = SVector(ρu, ρv, ρw)
    state.ρe = ρe_tot
    state.moisture.ρq_tot = ρ * q_tot

    if z <= FT(400) # Add random perturbations to bottom 400m of model
        state.ρe += rand(seed) * ρe_tot / 100
        state.moisture.ρq_tot += rand(seed) * ρ * q_tot / 100
    end
end

function config_bomex(FT, N, resolution, xmax, ymax, zmax)

    ics = init_bomex!     # Initial conditions

    C_smag = FT(0.18)     # Smagorinsky coefficient

    u_star = FT(0.28)     # Friction velocity

    T_sfc = FT(300.4)     # Surface temperature `[K]`
    LHF = FT(147.2)       # Latent heat flux `[W/m²]`
    SHF = FT(9.5)         # Sensible heat flux `[W/m²]`

    ∂qt∂t_peak = FT(-1.2e-8)  # Moisture tendency (energy source)
    zl_moisture = FT(300)     # Low altitude limit for piecewise function (moisture source)
    zh_moisture = FT(500)     # High altitude limit for piecewise function (moisture source)
    ∂θ∂t_peak = FT(-2 / FT(day(param_set)))  # Potential temperature tendency (energy source)

    z_sponge = FT(2400)     # Start of sponge layer
    α_max = FT(0.75)        # Strength of sponge layer (timescale)
    γ = 2              # Strength of sponge layer (exponent)

    u_geostrophic = FT(-10)        # Eastward relaxation speed
    u_slope = FT(1.8e-3)     # Slope of altitude-dependent relaxation speed
    v_geostrophic = FT(0)          # Northward relaxation speed

    zl_sub = FT(1500)         # Low altitude for piecewise function (subsidence source)
    zh_sub = FT(2100)         # High altitude for piecewise function (subsidence source)
    w_sub = FT(-0.65e-2)     # Subsidence velocity peak value

    f_coriolis = FT(0.376e-4) # Coriolis parameter

    # Assemble source components
    source = (
        Gravity(),
        BomexTendencies{FT}(
            ∂qt∂t_peak,
            zl_moisture,
            zh_moisture,
            ∂θ∂t_peak,
            zl_sub,
            zh_sub,
            w_sub,
            zmax,
        ),
        BomexSponge{FT}(
            zmax,
            z_sponge,
            α_max,
            γ,
            u_geostrophic,
            u_slope,
            v_geostrophic,
        ),
        BomexGeostrophic{FT}(f_coriolis, u_geostrophic, u_slope, v_geostrophic),
    )

    # Assemble timestepper components
    ode_solver_type = CLIMA.DefaultSolverType()

    # Assemble model components
    model = AtmosModel{FT}(
        AtmosLESConfigType,
        param_set;
        turbulence = SmagorinskyLilly{FT}(C_smag),
        moisture = EquilMoist{FT}(; maxiter = 5, tolerance = FT(0.1)),
        source = source,
        boundarycondition = (
            AtmosBC(
                momentum = Impenetrable(DragLaw(
                    # normPu_int is the internal horizontal speed
                    # P represents the projection onto the horizontal
                    (state, aux, t, normPu_int) -> (u_star / normPu_int)^2,
                )),
                energy = PrescribedEnergyFlux((state, aux, t) -> LHF + SHF),
                moisture = PrescribedMoistureFlux(
                    (state, aux, t) -> LHF / latent_heat_vapor(T_sfc),
                ),
            ),
            AtmosBC(),
        ),
        init_state = ics,
    )

    # Assemble configuration
    config = CLIMA.AtmosLESConfiguration(
        "BOMEX",
        N,
        resolution,
        xmax,
        ymax,
        zmax,
        param_set,
        init_bomex!,
        solver_type = ode_solver_type,
        model = model,
    )
    return config
end

function config_diagnostics(driver_config)
    interval = 10000 # in time steps
    dgngrp = setup_atmos_default_diagnostics(interval, driver_config.name)
    return CLIMA.DiagnosticsConfiguration([dgngrp])
end

function main()
    CLIMA.init()

    FT = Float32

    # DG polynomial order
    N = 4
    # Domain resolution and size
    Δh = FT(100)
    Δv = FT(40)

    resolution = (Δh, Δh, Δv)

    # Prescribe domain parameters
    xmax = FT(6400)
    ymax = FT(6400)
    zmax = FT(3000)

    t0 = FT(0)

    # For a full-run, please set the timeend to 3600*6 seconds
    # For the test we set this to == 30 minutes
    timeend = FT(1800)
    #timeend = FT(3600 * 6)
    CFLmax = FT(1.0)

    driver_config = config_bomex(FT, N, resolution, xmax, ymax, zmax)
    solver_config = CLIMA.SolverConfiguration(
        t0,
        timeend,
        driver_config,
        init_on_cpu = true,
        Courant_number = CFLmax,
    )
    dgn_config = config_diagnostics(driver_config)

    cbtmarfilter = GenericCallbacks.EveryXSimulationSteps(1) do (init = false)
        Filters.apply!(solver_config.Q, 6, solver_config.dg.grid, TMARFilter())
        nothing
    end

    result = CLIMA.invoke!(
        solver_config;
        diagnostics_config = dgn_config,
        user_callbacks = (cbtmarfilter,),
        check_euclidean_distance = true,
    )

    @test isapprox(result, FT(1); atol = 2e-3)
end

main()

max_stars_repo_name: yairchn/CLIMA

max_issues_repo_name: yairchn/CLIMA

max_forks_repo_name: yairchn/CLIMA


text:
using Ising2D
using Random: seed!
using Test

@testset "Ising2D.jl" begin
    N = 100

    seed!(4649373)
    S_ifelse = [ising2d!(IfElse(), rand_ising2d(), β_ising2d, 100) for _ in 1:N]
    @test abs(sum(energy_density_ising2d.(S_ifelse))/N + 1.37) < 0.015
    @test abs(sum(abs.(magnetization_ising2d.(S_ifelse)))/N - 0.184) < 0.07

    seed!(4649373)
    S_multifor = [ising2d!(MultiFor(), rand_ising2d(), β_ising2d, 100) for _ in 1:N]
    @test abs(sum(energy_density_ising2d.(S_multifor))/N + 1.37) < 0.015
    @test abs(sum(abs.(magnetization_ising2d.(S_multifor)))/N - 0.184) < 0.07

    @test S_ifelse == S_multifor
end


text:
# This file is a part of TextSearch.jl

import SparseArrays: sparsevec, sparse

export dvec

"""
    dvec(x::AbstractSparseVector)

Converts an sparse vector into a DVEC sparse vector
"""
function dvec(x::AbstractSparseVector)
    DVEC{eltype(x.nzind), eltype(x.nzval)}(t => v for (t, v) in zip(x.nzind, x.nzval))
end

"""
    sparsevec(vec::DVEC{Ti,Tv}, m=0) where {Ti<:Integer,Tv<:Number}

Creates a sparse vector from a DVEC sparse vector
"""
function sparsevec(vec::DVEC{Ti,Tv}, m=0) where {Ti<:Integer,Tv<:Number}
    I = Ti[]
    F = Tv[]

    for (t, weight) in vec
        if t > 0
            push!(I, t)
            push!(F, weight)
        end
    end

    if m == 0
        sparsevec(I, F)
    else
        sparsevec(I, F, m)
    end
end

"""
    sparse(cols::AbstractVector{S}, m=0; minweight=1e-9) where S<:DVEC{Ti,Tv} where {Ti<:Integer,Tv<:Number}

Creates a sparse matrix from an array of DVEC sparse vectors.
"""
function sparse(cols::AbstractVector{DVEC{Ti,Tv}}, m=0; minweight=1e-9) where {Ti<:Integer,Tv<:Number}
    I = Ti[]
    J = Ti[]
    F = Tv[]

    for j in eachindex(cols)
        for (term, weight) in cols[j]
            if term > 0 && weight >= minweight
                push!(I, term)
                push!(J, j)
                push!(F, weight)
            end
        end
    end

    if m == 0
        sparse(I, J, F)
    else
        sparse(I, J, F, m, length(cols))
    end
end

max_stars_repo_name: sadit/TextModel

max_issues_repo_name: sadit/TextModel

max_forks_repo_name: sadit/TextModel