text:
## --------------------------------------------------------------------------
function extract_input_file(line)

    # extract uncomment
    uncomment, _ = split_comment(line)

    reg = REGEXS[:TEX_INPUT_CMD_REGEX]
    inputm = match(reg, uncomment)
    return isnothing(inputm) ? "" : inputm[:label]
end

## --------------------------------------------------------------------------
# tex section
function _extract_tex_label(line, secregex)

    # extract uncomment
    uncomment, _ = split_comment(line)

    match_ = match(secregex, uncomment)
    return isnothing(match_) ? "" : match_[:label]
end

extract_tex_section_label(line) = _extract_tex_label(line, REGEXS[:TEX_SECTION_CMD_REGEX])
extract_tex_subsection_label(line) = _extract_tex_label(line, REGEXS[:TEX_SUBSECTION_CMD_REGEX])
extract_tex_subsubsection_label(line) = _extract_tex_label(line, REGEXS[:TEX_SUBSUBSECTION_CMD_REGEX])

extract_tex_begin_label(line) = _extract_tex_label(line, REGEXS[:TEX_BEGIN_REGEX])
extract_tex_end_label(line) = _extract_tex_label(line, REGEXS[:TEX_END_REGEX])

## --------------------------------------------------------------------------
# texass
function _extract_texass_label(line, reg)
    match_ = match(reg, line)
    return isnothing(match_) ? "" : match_[:label]
end
extract_texass_sec_name(line) = _extract_texass_label(line, REGEXS[:TEXASS_SECTION_REGEX])
extract_texass_end_label(line) = _extract_texass_label(line, REGEXS[:TEXASS_END_REGEX])

## --------------------------------------------------------------------------
# utils
function extract_indent(line)
    match_ = match(REGEXS[:INDENT], line)
    return isnothing(match_) ? "" : match_[:indent]
end


text:
### A Pluto.jl notebook ###
# v0.16.1

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ ab02837b-79ec-40d7-bff1-c1d2dd7362ef
md"""
> # PlutoTest.jl
>
> _Visual, reactive testing library for Julia_

A macro `@test` that you can use to verify your code's correctness.

**But instead of just saying _"Pass"_ or _"Fail"_, let's try to show you _why_ a test failed.**
- ✨ _time travel_ ✨ to replay the execution step-by-step
- ⭐ _multimedia display_ ⭐ to make results easy to read
"""

# ╔═╡ 56347b7e-5007-45f8-8f6d-8ac8cc719637
md"""
Tests have _time-travel_ functionality built in! **Click on the test results above.**
"""

# ╔═╡ fd8428a3-9fa3-471a-8b2d-5bbb8fdb3137
is_good_boy(x) = true;

# ╔═╡ 191f1f04-18d4-485b-af8b-a2f073b7043b
md"""
## Installation and more info

> [github.com/JuliaPluto/PlutoTest.jl](https://github.com/JuliaPluto/PlutoTest.jl)
"""

# ╔═╡ ec1fd70a-d92a-4688-98b2-135879f07141
md"""
### (You need `Pluto ≥ 0.14.5` to run this notebook)
"""

# ╔═╡ 80b2fb3f-94b2-4024-94ff-d111a249c8b0

# ╔═╡ 9d49ea50-8158-4d8b-97af-edba1f7dc38b
x = [1,3]

# ╔═╡ 1aa24b1c-e8ca-4de7-b614-7a3f02b4833d
always_false(args...; kwargs...) = true

# ╔═╡ 8a2e8348-49cf-4855-b5b3-cdee33e5ed67
const pluto_test_css = """
pt-dot {
	flex: 0 0 auto;
	background: grey;
	width: 1em;
	height: 1em;
	bottom: -.1em;
	border-radius: 100%;
	margin-right: .7em;
	display: block;
	position: relative;
	cursor: pointer;
}

pt-dot.floating {
	position: fixed;
	z-index: 60;
	visibility: hidden;
	transition: transform linear 240ms;
	opacity: .8;
}
.show-top-float > pt-dot.floating.top,
.show-bottom-float > pt-dot.floating.bottom {
	visibility: visible;
}

pt-dot.floating.top {
	top: 5px;
}
pt-dot.floating.bottom {
	bottom: 5px;
}

.fail > pt-dot {
	background: #f75d5d;

}
.pass > pt-dot {
	background: #56a038;
}

@keyframes fadeout {
    0% { opacity: 1;}
    100% { opacity: 0; pointer-events: none;}
}

.pass > pt-dot.floating {

    animation: fadeout 2s;

	animation-fill-mode: both;
	animation-delay: 2s;

	/*opacity: 0.4;*/

}

.pluto-test {
	font-family: "JuliaMono", monospace;
	font-size: 0.75rem;
padding: 4px;

min-height: 25px;
}

.pluto-test.pass {
	color: rgba(0,0,0,.5);
}

.pluto-test.fail {
background: linear-gradient(90deg, #ff2e2e14, transparent);
border-radius: 7px;
}

.pluto-test>.arg_result {
	flex: 0 0 auto;
}

.pluto-test>.arg_result>div,
.pluto-test>.arg_result>div>pluto-display>div {
	display: inline-flex;
}

.pluto-test>.comma {
	margin-right: .5em;
}

.pluto-test.call>code {
	padding: 0px;
}

.pluto-test.call.infix-operator>div {
	overflow-x: auto;
}

.pluto-test {
	display: flex;
	align-items: baseline;
}

.pluto-test.call.infix-operator>.fname {
	margin: 0px .6em;
	/*color: darkred;*/
}

/* expanding */

.pluto-test:not(.expanded) {
	cursor: pointer;
}

.pluto-test:not(.expanded) > p-frame-viewer > p-frame-controls {
	display: none;

}

.pluto-test.expanded > p-frame-viewer {
    max-width: 100%;
}
.pluto-test.expanded > p-frame-viewer > p-frames > slotted-code > line-like {
	flex-wrap: wrap;
}
.pluto-test.expanded > p-frame-viewer > p-frames > slotted-code > line-like > pluto-display[mime="application/vnd.pluto.tree+object"] {
	/*flex-basis: 100%;*/
}
""";

# ╔═╡ 42671258-07a0-4015-8f47-4b3032595f08
const frames_css = """
p-frame-viewer {
	display: inline-flex;
	flex-direction: column;
}
p-frames,
p-frame-controls {
	display: inline-flex;
}
""";

# ╔═╡ 0d70962a-3880-4dee-a439-35068d019f5a
md"""
# Type definitions
"""

# ╔═╡ 113cc425-e224-4f77-bfbd-ef4eb1d1ed70
abstract type TestResult end

# ╔═╡ 6188f559-bcab-4da6-84b2-a3fe522a5c3c
abstract type Fail <: TestResult end

# ╔═╡ c24b46ce-bcbb-4dc9-8a59-b5b1bd2cd617
abstract type Pass <: TestResult end

# ╔═╡ 5041085e-a406-4ed4-ab82-84d8f126cf0f
const Code = Any

# ╔═╡ 8c92bad9-234e-47dd-a599-b75dc6d5db89
# struct Correct <: Pass
# 	expr::Code
# end

# ╔═╡ 03ccd498-83c3-41bb-84d7-625adabd7aee
struct CorrectCall <: Pass
	expr::Code
	steps::Vector
end

# ╔═╡ 1bcf8bd1-c8a3-49a1-9791-d813aa856399
# struct Error <: Fail
# 	expr::Code
# 	error
# end

# ╔═╡ 656c4190-b49e-4225-869d-eeb7e8e41e72
# struct Wrong <: Fail
# 	expr::Code
# 	step
# end

# ╔═╡ 14c525a1-eca1-466b-8e63-3a90d7d7111c
struct WrongCall <: Fail
	expr::Code
	steps::Vector
end

# ╔═╡ a2efc968-246c-40c2-b285-2ec94b185a44
md"""
# Test macro
"""

# ╔═╡ b6b1e037-754e-4c9c-a79d-2a1517cb409d
remove_quotes(ex::Expr) = Expr(if ex.head === :$
		:block
	else
		ex.head
	end, remove_quotes.(ex.args)...)

# ╔═╡ 8020701d-e23a-4d95-b4ad-ebeddd374fb5
remove_quotes(x) = x

# ╔═╡ fc1ec730-4e11-452c-94a7-26605a87a89e
let
	e = """
		map(1:15) do i
			@test 2 * \$i > 0.19
		end
	""" |> Meta.parse

	e, remove_quotes(e)
end

# ╔═╡ bfe4dc61-9160-4c7e-8897-9c723b309adc
# function test(expr)
# 	if Meta.isexpr(expr, :call, 3) && expr.args[1] === :(==)
# 	quote
# 		expr_raw = $(QuoteNode(expr))
# 		try
# 			left = $(expr.args[2] |> esc)
# 			right = $(expr.args[3] |> esc)

# 			result = left == right

# 			if result === true
# 				Correct(expr_raw)
# 			elseif result === false
# 				WrongEquality(expr_raw, left, right)
# 			else
# 				Wrong(expr_raw, result)
# 			end
# 		catch e
# 			rethrow(e)
# 			# Error(expr_raw, e)
# 		end
# 	end
# 	end
# end

# ╔═╡ 4bc1b7a4-0a36-4a07-b7ee-3d5be50350e1
("asd"*"asd","asd","asd"*" asd","asd","asd"*"asd","asd","asd"*"asd","asd","asd"*"asd","asd","asd"*"asd","asd","asd"*"asd","asd")

# ╔═╡ 3b2e8f55-1d4b-4a36-83f6-26becbd79e4b
# @test (@test t isa Pass) isa Pass

# ╔═╡ ec2ed42c-1227-4e0d-b642-20e6f3503d2a
embed_display(x) = if isdefined(Main, :PlutoRunner) && isdefined(Main.PlutoRunner, :embed_display)
	# if this package is used inside Pluto, and Pluto is new enough
	Main.PlutoRunner.embed_display(x)
else
	identity(x)
end

# ╔═╡ 1ac164c8-88fc-4a87-a194-60ef616fb399
flatmap(args...) = vcat(map(args...)...)

# ╔═╡ 98ac4c36-49c7-4f65-982d-0b8bf6c372c0
emb = embed_display

# ╔═╡ 4d5f44e4-85e9-4985-9b76-73be5e097186
remove_linenums(e::Expr) = if e.head === :macrocall
	Expr(e.head, (remove_linenums(x) for x in e.args)...)
else
	Expr(e.head, (remove_linenums(x) for x in e.args if !(x isa LineNumberNode))...)
end

# ╔═╡ dd495e00-d74d-47d4-a5d5-422fb147ec3b
remove_linenums(x) = x

# ╔═╡ a83a6a4c-664c-46fa-a07f-81088493dc35
const ObjectID = typeof(objectid("hello computer"))

# ╔═╡ 5cb03161-2cbc-4080-ba59-f94efd3b620c
expr_hash(e::Expr) = objectid(e.head) + mapreduce(p -> objectid((p[1], expr_hash(p[2]))), +, enumerate(e.args); init=zero(ObjectID))

# ╔═╡ 0611a36b-b4be-4b17-a485-7c4a8fa04927
expr_hash(x) = objectid(x)

# ╔═╡ dbfbcc16-c740-436c-bbf0-fee16b0a20c5
md"""
# $(html"") _Time travel_ evaluation

In Julia, expressions are objects! This means that, before evaluation, code is expressed as a Julia object:
"""

# ╔═╡ d97987a0-bdc0-46ed-a6a5-f35c1ce961dc
ex1 = :(first([56,sqrt(9)]))

# ╔═╡ 7c2bab29-8609-4575-b2ca-7feb34915645
md"""
You can use `Core.eval` to evaluate expressions at runtime:
"""

# ╔═╡ 838b5904-1de2-4d9f-8f3c-a93ec224d40e
md"""
But _did you know_ that you can also **partially evaluate** expressions?
"""

# ╔═╡ a3c41025-2f4a-4f9c-8577-72e4b7abbb98
ex2_inner = ex1.args[2]

# ╔═╡ b056a99d-5b13-47ba-a199-d788410e3c99
md"""
Here, `ex2` is not a raw `Expr` — it _contains_ an evaluated array!
"""

# ╔═╡ 5b093e83-78c1-4187-b406-56e79800e1be
md"""
#### `Computed` struct

Our time travel mechanism will be based on the partial evaluation principle introduced above. To differentiate between computed results and the original expression, we will wrap all computed results in a `struct`.
"""

# ╔═╡ a461f1fd-b5a5-4ae3-a47c-067a6081fb24
struct Computed
	x
end

# ╔═╡ b765dbfe-4e58-4bb9-b1d6-aa4378d4e9c9
expr_to_str(e, mod=@__MODULE__()) = let
	Computed
	sprint() do io
		Base.print(IOContext(io, :module => @__MODULE__), remove_linenums(e))
	end
end

# ╔═╡ ea73a35e-a34f-4708-acc1-858f2466e9ba
expr_to_str(:(x+1))

# ╔═╡ ef6fc423-f1b1-4dcb-a059-276121391bc6
prettycolors(e) = Markdown.MD([Markdown.Code("julia", expr_to_str(e))])

# ╔═╡ f9c81ab1-556c-4d81-bee8-2897c20e324d
md"""
We also add a function to recursively _unwrap_ an expression with `Computed` entries:
"""

# ╔═╡ a392d2d6-5a16-4383-b0ef-5003aa2de9fa
unwrap_computed(x) = x

# ╔═╡ ae95b691-f54b-4bf5-b17b-3e5bd1edf75e
unwrap_computed(c::Computed) = c.x

# ╔═╡ 12119016-fa61-4d38-8c58-821ea435df7d
unwrap_computed(e::Expr) = Expr(e.head, unwrap_computed.(e.args)...)

# ╔═╡ 2c1b906d-71b9-430e-83ed-d4c8c0018632
md"""
## Stepping function

"""

# ╔═╡ 450a36ea-2c43-4f01-a775-0b8c59bf6dca
function onestep_light(e::Expr; m=Module())::Vector
	results = Any[]

	seval(ex) = Computed(Core.eval(m,ex))

	# will be modified
	arg_results = Any[a for a in e.args]

	push_intermediate() = push!(results, Expr(e.head, arg_results...))

	# we only step for expressions where this is possible/easy
	if e.head === :call || e.head === :begin || e.head === :block# || e.head === :vect

		for (i,a) in enumerate(e.args)

			if a isa QuoteNode
				a

			elseif (Meta.isexpr(e, :call) || Meta.isexpr(e, :let)) && i == 1
				a

			elseif a isa Symbol
				arg_results[i] = seval(a)
				push_intermediate()

			elseif a isa Expr
				inner_results = onestep_light(a; m=m)
				for ir in inner_results
					arg_results[i] = ir
					push_intermediate()
				end

				arg_results[i] = inner_results[end]
			else
				a
			end
		end
	end

	push!(results, seval(unwrap_computed(Expr(e.head, arg_results...))))

	results
end

# ╔═╡ ef45af85-11f8-4505-9f1e-3ffb15a47142
onestep_light(x::Any; m=Module()) = [Computed(Core.eval(m,x))]

# ╔═╡ 97a05d7c-2e08-458d-a628-7992f204b4ea
md"""

"""

# ╔═╡ 88f6a040-07cf-47e0-a8be-2478ea350aa7
can_interpret(x) = true

# ╔═╡ d36a8a72-eced-4e63-9130-7fcb6c86df76
can_interpret_call_arg(e::Expr) = !(e.head === :(...) || e.head === :kw || e.head === :parameters)

# ╔═╡ e9659020-d433-4357-9099-71a65b66a091
can_interpret_call_arg(x) = true

# ╔═╡ c3fc3292-b7eb-4b01-8fba-159c86228de9
Meta.isbinaryoperator(:(==))

# ╔═╡ 89578bff-16b9-4eb2-b8ee-b2839ff2d74c
can_interpret(e::Expr) = if false
	false
elseif e.head === :call && !all(can_interpret_call_arg, e.args)
	false
# elseif e.head === :(=) || e.head === :macrocall
# 	false
else
	all(can_interpret, e.args)
end

# ╔═╡ e1c306e3-0a47-4149-a9fb-ec7ab380fa11
function step_by_step(expr; __module__)
	Computed
	onestep_light
	if can_interpret(expr)
		quote
			Any[$(QuoteNode(expr)), $(onestep_light)($(esc(Expr(:quote,expr))); m=$(__module__))...]
		end
	else
		quote
			[$(QuoteNode(expr)), Computed($(esc(expr)))]
		end
	end
end

# ╔═╡ a661e172-6afb-42ff-bd43-bb5b787ee5ed
macro eval_step_by_step(e)
	step_by_step
	quote
		if false
			$(esc(remove_quotes(e)))
		end
		$(step_by_step(e; __module__=__module__))
	end
end

# ╔═╡ 930f8244-cf25-4c1a-95f6-5c8963559c62
@macroexpand @eval_step_by_step x == [1,2]

# ╔═╡ 68ba60db-44ad-43e4-b33e-d27696babc99
@eval_step_by_step sqrt(sqrt(length([1,2])))

# ╔═╡ 807bcd72-26c3-44d3-a295-56874cb51a89
@eval_step_by_step xasdf = 123

# ╔═╡ 8a5a4c26-e36c-4061-b32f-4448625ce4a6
try
	xasdf
catch e
	e
end

# ╔═╡ 21d4560e-721f-4ed4-9db7-86a8151ab22c
md"""
## Displaying objects inside code
"""

# ╔═╡ 99afc7f4-727c-4277-8311-f2ffa94830ae
md"""
#### Slotting

We walk through the expression tree. Whenever we find a `Computed` object, we generate a random key (e.g. `iosjddfo`), we add it to our dictionary (`found`). In the expression, we replace the `Computed` object with a placeholder symbol `__slotiosjddfo__`. We will later be able to match the object to this slot.
"""

# ╔═╡ 4956526a-daf9-43c9-bff3-ff2446016e2e
slot!(found, c::Computed) = let
	k = Symbol("__slot", join(rand('a':'z', 16)), "__")
	found[k] = c
	k
end

# ╔═╡ 84ff6a23-c134-4910-b630-a7ad45f3bf29
slot!(found, x) = x

# ╔═╡ 318363d0-6d9e-4144-b478-b775f437edaf
slot!(found, e::Expr) = Expr(e.head, slot!.([found], e.args)...)

# ╔═╡ 67fd07b7-340b-4e24-bc06-e4c85b186872
slot(e) = let
	d = Dict{Symbol,Any}()
	new_e = slot!(d, e)
	d, new_e
end

# ╔═╡ c6d5597c-d505-4125-88c4-10415934d2a4
md"""
#### SlottedDisplay

We use `print` to turn the expression into source code.

For each line, we regex-search for slot variables, and we split the line around those. The code segments around slots are rendered inside `<pre-ish>` tags (like `<pre>` but inline), and the slots are replaced by [embedded displays](https://github.com/fonsp/Pluto.jl/pull/1126) of the objects.
"""

# ╔═╡ c877c109-db16-468c-8f3c-8294db859d6d
begin
	struct SlottedDisplay
		d
		e
	end
	SlottedDisplay(expr) = SlottedDisplay(slot(expr)...)
end

# ╔═╡ 8480d0d7-bdf7-468d-9344-5b789e33921c
const slotted_code_css = """
slotted-code {
	font-family: "JuliaMono", monospace;
	font-size: .75rem;
	display: flex;
	flex-direction: column;
}
pre-ish {
	white-space: pre;
}

line-like {
	display: flex;
	align-items: baseline;
}
"""

# ╔═╡ b5763c10-e11c-4389-b6fc-421d2c9682f1
md"""
#### Frame viewer

A widget that takes a series of elements and displays them as 'video frames' with a timeline scrubber.
"""

# ╔═╡ 3d5abd58-02ab-4b91-a7a3-d9068d4df017
md"""
#### Macro to test frames
"""

# ╔═╡ 34f613a3-85fb-45a8-be3b-cd8e6b3cb5a2

# ╔═╡ f9ed2487-a7f6-4ce9-b673-f8a298cd5fc3
md"""
# Appendix
"""

# ╔═╡ 35f63c4e-3583-4ea8-a057-31f18f8a09d6
md"""
## `@skip_as_script`
"""

# ╔═╡ 35b2770e-1db6-4327-bf86-c27a4b61dbd3
function is_inside_pluto(m::Module)::Bool
	if isdefined(m, :PlutoForceDisplay)
		return m.PlutoForceDisplay
	else
		isdefined(m, :PlutoRunner) && parentmodule(m) === Main
	end
end

# ╔═╡ 22640a2f-ea38-4517-a4f3-7a65e60ffebe
"""
	@displayonly expression

Marks a expression as Pluto-only, which means that it won't be executed when running outside Pluto. Do not use this for your own projects.
"""
macro skip_as_script(ex) is_inside_pluto(__module__) ? esc(ex) : nothing end

# ╔═╡ cf314b21-3f4f-4637-b1ce-ec1d5d5af966
begin
	if is_inside_pluto(@__MODULE__)
		import Pkg
		Pkg.activate("..")
		Pkg.instantiate()
	end
	import HypertextLiteral: @htl
	import Test
end

# ╔═╡ 05f42764-acfe-4370-b85b-ce0e7c4270d0
begin
	export @test_nowarn, @test_warn, @test_logs, @test_skip, @test_broken, @test_throws, @test_deprecated

	var"@test_warn" = Test.var"@test_warn"
	var"@test_nowarn" = Test.var"@test_nowarn"
	var"@test_logs" = Test.var"@test_logs"
	var"@test_skip" = Test.var"@test_skip"
	var"@test_broken" = Test.var"@test_broken"
	var"@test_throws" = Test.var"@test_throws"
	var"@test_deprecated" = Test.var"@test_deprecated"
end

# ╔═╡ b6e8a170-12cc-4d97-905d-274e2609bfd8
function test(expr, extra_args...; __module__)
	step_by_step
	Test.test_expr!("", expr, extra_args...)

	quote
		expr_raw = $(QuoteNode(expr))
		try
			# steps = @eval_step_by_step($(expr))

			steps = $(step_by_step(expr; __module__=__module__))

# 			arg_results = [$((expr.args[2:end] .|> esc)...)]

# 			result = $(esc(:eval))(Expr(:call, $(expr.args[1] |> QuoteNode), arg_results...))

			result = unwrap_computed(last(steps))

			if result === true
				CorrectCall(expr_raw, steps)
			# elseif result === false
			# 	WrongCall(expr_raw, steps)
			else
				WrongCall(expr_raw, steps)
			end
		catch e
			rethrow(e)
			# Error(expr_raw, e)
		end
	end
end

# ╔═╡ 0fcc6cb0-2711-4609-9bf3-634cf9407840
div(x; class="", style="") = @htl("
$(x)
")

# ╔═╡ 69200d7c-b7bc-4c7e-a9a1-5e26979179a3
div(; class="", style="") = x -> @htl("
$(x)
")

# ╔═╡ 872b4877-30dd-4a92-a3c8-69eb50675dcb
preish(x) = @htl("
$(x)
")

# ╔═╡ ab0a19b8-cf7c-4c4f-802a-f85eef81fc02
function Base.show(io::IO, m::MIME"text/html", sd::SlottedDisplay)

	d, e = sd.d, sd.e

	s = sprint() do iobuffer
		print(IOContext(iobuffer, io), e |> remove_linenums)
	end

	lines = split(s, "\n")

	r = r"\_\_slot[a-z]{16}\_\_"
	embed_display
	h = @htl("""<slotted-code>
		$(
	map(lines) do l
		keys = [Symbol(m.match) for m in eachmatch(r, l)]
		rest = split(l, r; keepempty=true)

		result = vcat((
			[(isempty(r) ? @htl("") : preish(r)), embed_display(d[k].x)]
			for (r,k) in zip(rest, keys)
			)...)

		push!(result, preish(last(rest)))

		@htl("<line-like>$(result)</line-like>")
	end
	)
		</slotted-code>""")
	show(io, m, h)
end

# ╔═╡ e968fc57-d850-4e2d-9410-8777d03b7b3c
function frames(fs::Vector)
	l = length(fs)

	startframe = l > 2 ? l - 1 : l

	@htl("""

		$(fs)

		<input class="timescrub" style="filter: hue-rotate(149deg) grayscale(.9);" type=range min=1 max=$(l) value=$(startframe)>

		""")

end

# ╔═╡ b273d3d3-648f-4d34-94e7-e49277d4ba29
with_slotted_css(x) = @htl("""
	$(x)

	""")

# ╔═╡ 326f7661-3482-4bf2-a97b-57cc7ac60ee2
macro visual_debug(expr)
	frames
	SlottedDisplay
	var"@eval_step_by_step"
	with_slotted_css
	quote
		@eval_step_by_step($(expr)) .|> SlottedDisplay |> frames |> with_slotted_css
	end
end

# ╔═╡ 69bfb438-7ecf-4f9b-8bc4-51e07aa46ef1
@skip_as_script Core.eval(Module(), ex1)

# ╔═╡ 3e79ff61-6532-4879-9402-86473aa7d960
ex2_inner_result = @skip_as_script Core.eval(Module(), ex2_inner)

# ╔═╡ 275c5f57-623d-439f-b09d-f7c745e0bed6
ex2 = Expr(:call, :first, ex2_inner_result)

# ╔═╡ 38e54516-cdf4-4c1d-815b-68e1e7a7f6f7
ex3 = Expr(:call, :first, Computed(ex2_inner_result))

# ╔═╡ 9bed78b6-5a8f-44ce-ab66-cab685daf264
unwrap_computed(ex3)

# ╔═╡ b5a9b7c1-43fb-4ec6-aaeb-0ec55580ef64
@skip_as_script quote
	sqrt(sqrt(4)) + 2
end |> onestep_light .|> prettycolors

# ╔═╡ 8ef356ea-7d54-43e6-a936-7c8be04c595f
@skip_as_script onestep_light(quote
		1+2
		2+3
		4+5
		sqrt(sqrt(sqrt(5)))
	end) .|> prettycolors

# ╔═╡ 4edf747b-3838-4315-a397-e452ac9b5465
@skip_as_script onestep_light(quote
		(1+2) + (7-6)
		2+3
		4+5
		sqrt(sqrt(sqrt(5)))
	end |> remove_linenums) .|> slot

# ╔═╡ d414f840-4952-4de5-a565-7fdc81a94817
"The opposite of `@skip_as_script`"
macro only_as_script(ex) is_inside_pluto(__module__) ? nothing : esc(ex) end

# ╔═╡ 326825b0-a17f-427a-9056-8e8156098418
@skip_as_script "hello"

# ╔═╡ 64bf02a4-4fe3-424d-ae6e-5906c3395278
md"""
## PlutoUI favourites
"""

# ╔═╡ f3916810-1911-48bd-936b-776206fcad54
const toc_css = """
@media screen and (min-width: 1081px) {
	.plutoui-toc.aside {
		position:fixed;
		right: 1rem;
		top: 5rem;
		width:25%;
		padding: 10px;
		border: 3px solid rgba(0, 0, 0, 0.15);
		border-radius: 10px;
		box-shadow: 0 0 11px 0px #00000010;
		/* That is, viewport minus top minus Live Docs */
		max-height: calc(100vh - 5rem - 56px);
		overflow: auto;
		z-index: 50;
		background: white;
	}
}

.plutoui-toc header {
	display: block;
	font-size: 1.5em;
	margin-top: 0.67em;
	margin-bottom: 0.67em;
	margin-left: 0;
	margin-right: 0;
	font-weight: bold;
	border-bottom: 2px solid rgba(0, 0, 0, 0.15);
}

.plutoui-toc section .toc-row {
	white-space: nowrap;
	overflow: hidden;
	text-overflow: ellipsis;
	padding-bottom: 2px;
}

.highlight-pluto-cell-shoulder {
	background: rgba(0, 0, 0, 0.05);
	background-clip: padding-box;
}

.plutoui-toc section a {
	text-decoration: none;
	font-weight: normal;
	color: gray;
}
.plutoui-toc section a:hover {
	color: black;
}

.plutoui-toc.indent section a.H1 {
	font-weight: 700;
	line-height: 1em;
}

.plutoui-toc.indent section a.H1 {
	padding-left: 0px;
}
.plutoui-toc.indent section a.H2 {
	padding-left: 10px;
}
.plutoui-toc.indent section a.H3 {
	padding-left: 20px;
}
.plutoui-toc.indent section a.H4 {
	padding-left: 30px;
}
.plutoui-toc.indent section a.H5 {
	padding-left: 40px;
}
.plutoui-toc.indent section a.H6 {
	padding-left: 50px;
}
"""

# ╔═╡ 122c27a5-a6e8-45ef-a968-b9b4b3f9ad09
const toc_js = toc -> """
const getParentCell = el => el.closest("pluto-cell")

const getHeaders = () => {
	const depth = Math.max(1, Math.min(6, $(toc.depth))) // should be in range 1:6
	const range = Array.from({length: depth}, (x, i) => i+1) // [1, ..., depth]

	const selector = range.map(i => `pluto-notebook pluto-cell h\${i}`).join(",")
	return Array.from(document.querySelectorAll(selector))
}

const indent = $(repr(toc.indent))
const aside = $(repr(toc.aside))

const render = (el) => html`\${el.map(h => {
	const parent_cell = getParentCell(h)

	const a = html`\${h.innerText}`
	/* a.onmouseover=()=>{
		parent_cell.firstElementChild.classList.add(
			'highlight-pluto-cell-shoulder'
		)
	}
	a.onmouseout=() => {
		parent_cell.firstElementChild.classList.remove(
			'highlight-pluto-cell-shoulder'
		)
	} */
	a.onclick=(e) => {
		e.preventDefault();
		h.scrollIntoView({
			behavior: 'smooth',
			block: 'center'
		})
	}

	return html`<div class="toc-row">\${a}</div>`
})}`

const tocNode = html`
	$(toc.title)

`
tocNode.classList.toggle("aside", aside)
tocNode.classList.toggle("indent", aside)

const updateCallback = () => {
	tocNode.querySelector("section").replaceWith(
		html`<section>\${render(getHeaders())}</section>`
	)
}
updateCallback()

const notebook = document.querySelector("pluto-notebook")

// We have a mutationobserver for each cell:
const observers = {
	current: [],
}

const createCellObservers = () => {
	observers.current.forEach((o) => o.disconnect())
	observers.current = Array.from(notebook.querySelectorAll("pluto-cell")).map(el => {
		const o = new MutationObserver(updateCallback)
		o.observe(el, {attributeFilter: ["class"]})
		return o
	})
}
createCellObservers()

// And one for the notebook's child list, which updates our cell observers:
const notebookObserver = new MutationObserver(() => {
	updateCallback()
	createCellObservers()
})
notebookObserver.observe(notebook, {childList: true})

// And finally, an observer for the document.body classList, to make sure that the toc also works when if is loaded during notebook initialization
const bodyClassObserver = new MutationObserver(updateCallback)
bodyClassObserver.observe(document.body, {attributeFilter: ["class"]})

invalidation.then(() => {
	notebookObserver.disconnect()
	bodyClassObserver.disconnect()
	observers.current.forEach((o) => o.disconnect())
})

return tocNode
"""

# ╔═╡ 9126f47d-cbc7-411f-93bd-8684ba06c9e9
toc() = HTML("""

	""")

# ╔═╡ c763ed72-82c9-445c-a8f7-a0c40982e4d9
toc()

# ╔═╡ 955705f9-c90d-495d-86b4-5f3b5bc9fc8e
begin
	struct Slider
		xs
	end

	Base.get(s::Slider) = first(s.xs)

	Base.show(io::IO, m::MIME"text/html", s::Slider) = show(io, m, @htl("<input type=range min=$(minimum(s.xs)) max=$(maximum(s.xs)) value=$(first(s.xs))>"))

	Slider
end

# ╔═╡ a6f82260-3519-4254-a21e-abc7bb19ec4e
@bind howmuch Slider(0:100)

# ╔═╡ c369b4b5-2fcf-4029-a1f6-352120b2fc4b
@bind n Slider(1:10)

# ╔═╡ 98992db9-4f14-4aa6-a7c5-477622266112
@bind k Slider(0:15)

# ╔═╡ 187c3005-cd43-45a0-8cbd-bc96b9cb39da
Dump(x; maxdepth=8) = sprint(io -> dump(io, x; maxdepth=maxdepth)) |> Text

# ╔═╡ a6709e08-964d-46ea-9813-2c70a834824b
Dump(ex1)

# ╔═╡ 10803c0d-d0a5-45c5-b7ef-9659e441df69
Dump(ex2)

# ╔═╡ 411271a6-4236-45e2-ab34-f26410108821
Dump(ex3)

# ╔═╡ daee414b-3e3c-4e2a-a25a-429a1e7275d5
Dump(Ref(1))

# ╔═╡ 6c0156a9-7281-4326-9e1f-989efa73bb7b
begin
	struct Show{M <: MIME}
		mime::M
		data
	end

	Base.show(io::IO, ::M, x::Show{M}) where M <: MIME = write(io, x.data)

	Show
end

# ╔═╡ e46cf3e0-aa15-4c17-a925-3e9fc5109d54
Hannes = let
	url = "https://user-images.githubusercontent.com/6933510/116753174-fa40ab80-aa06-11eb-94d7-88f4171970b2.jpeg"
	data = read(download(url))
	Show(MIME"image/jpg"(), data)
end;

# ╔═╡ 6f5ba692-4b6a-405a-8cd3-1a8f9cc06611
plot(args...; kwargs...) = Hannes

# ╔═╡ b4b317d7-bed1-489c-9650-8d336e330689
rs = @eval_step_by_step(begin
		(1+2) + (7-6)
		plot(2000 .+ 30 .* rand(2+2))
		4+5
		sqrt(sqrt(sqrt(5)))
	end) .|> SlottedDisplay

# ╔═╡ 93ed973f-daf6-408b-9d4b-d53495418610
@bind rindex Slider(eachindex(rs))

# ╔═╡ dea898a0-1904-4d09-ad0b-6915008fe946
rs[rindex]

# ╔═╡ 74c19786-1ba7-4865-a993-590a779ae564
frames(rs)

# ╔═╡ a2cbb0c3-23b9-4091-9ca7-5ba96e85e3a3
@visual_debug begin
	(1+2) + (7-6)
	plot(2000 .+ 30 .* rand(2+2))
	4+5
	sqrt(sqrt(sqrt(5)))
	md"# Wow"
end

# ╔═╡ 5b70aaf1-9623-4f55-b055-4263ed8be31d
Floep = let
	url = "https://user-images.githubusercontent.com/6933510/116753861-142ebe00-aa08-11eb-8ce8-684af1098935.jpeg"
	data = read(download(url))
	Show(MIME"image/jpg"(), data)
end;

# ╔═╡ bf2abe01-6ae0-4066-8704-12f64e04511b
friends = Any[Hannes, Floep];

# ╔═╡ 9c3f6eab-b1c3-4607-add8-d6d7e468c11a
begin
	export @test

	macro test(expr...)
		quote
			if false
				$(esc(remove_quotes(expr[1])))
			end
			$(test(expr...; __module__=__module__))
		end
	end

	function Base.show(io::IO, m::MIME"text/html", call::Union{WrongCall,CorrectCall})

		infix = if Meta.isexpr(call.expr, :call)
			fname = call.expr.args[1]
			Meta.isbinaryoperator(fname)
		else
			false
		end

		classes = [
			"pluto-test",
			"call",
			(isa(call,CorrectCall) ? "correct" : "wrong"),
			(isa(call,Pass) ? "pass" : "fail"),
			infix ? "infix-operator" : "prefix-operator",
			]

		result = @htl("""

			<pt-dot></pt-dot>
			<pt-dot class="floating top"></pt-dot>
			<pt-dot class="floating bottom"></pt-dot>

			$(frames(SlottedDisplay.( call.steps)))

		""")
		Base.show(io, m, result)
	end
end

# ╔═╡ 73d74146-8f60-4388-aaba-0dfe4215cb5d
@test sqrt(20-11) == 3

# ╔═╡ 71b22e76-2b50-4d16-85f6-9dad0415630e
@test iseven(123 + 7^3)

# ╔═╡ 6762ed72-f422-43a9-a782-de78f739c0ae
@test 4+4 ∈ [1:7...]

# ╔═╡ f77275b9-90aa-4e07-a608-981b5df727af
@test is_good_boy(first(friends))

# ╔═╡ eab4ba31-c787-46dd-8024-693eca7fd1a0
@test x == [1,2+2]

# ╔═╡ 26b0faf0-9016-48d7-8667-c1c1cfce655e
@test missing == 2

# ╔═╡ 26b4fb86-892f-415c-8046-6a5449052fd7
@test 2+2 == 2+2

# ╔═╡ 96dc7b01-3766-4206-88ba-eca1665bc5cb
@test rand(50) == [rand(50),2]

# ╔═╡ 7c6ce205-053d-434c-b5b1-500babb8ec02
@test always_false(rand(howmuch), rand(howmuch),123)

# ╔═╡ fe7f8cce-a706-476d-8680-a2fe793b474f
@test always_false(rand(2), rand(2),123)

# ╔═╡ 1872f785-1ae5-43ea-bce1-6c5cd893f3a8
@test !!always_false(rand(2), rand(2),123; r=123)

# ╔═╡ 5570972e-9309-4458-99a6-ea718ec2c3ab
@test always_false([1,2,3]...)

# ╔═╡ 8d340983-ea07-4038-872f-22a165003ed2
@test isless(2+2,1)

# ╔═╡ ea5a4fc0-db62-41dd-9600-a21d4eabf822
@test isless(1,2+2)

# ╔═╡ 4509cdbf-8b8b-4f70-9e63-bb972eb88c93
@test iseven(n^2)

# ╔═╡ 8360d1bc-b1f4-4263-a042-724cbd120227
@test 4+4 ∈ [1:k...]

# ╔═╡ 064e28de-0c22-48b5-b427-6eb343880287
@test isempty((1:k) .^ 2)

# ╔═╡ 3a6a6ee1-c619-4044-b9b1-68e5ae9d2463
map(1:10) do i
	@test sqrt($i) < 3 && always_false()
end

# ╔═╡ be93a6f4-b626-43db-a2fe-4e754e79c030
@test isempty([1,sqrt(2)])

# ╔═╡ 17bd5cd9-212f-4656-ab79-590dd6c64ff8
@test 1 ∈ [sqrt(20), 5:9...]

# ╔═╡ 539e2c38-993b-4b3b-8aa0-f02d46d79839
@test 1 ∈ rand(60)

# ╔═╡ 3d3f3592-e056-4e7b-8896-a75e5b5dcad6
@test rand(60) ∋ 1

# ╔═╡ c39021dc-157c-4bcb-a3a9-fec8d9286b48
map(1:15) do i
	@test 2 * $i > 0.19
end

# ╔═╡ ac02b12a-3982-4526-a51c-0bf85198b81b
var"@test"; macroexpand(@__MODULE__, :(@test x == [1,2+i]); recursive=false) |> prettycolors

# ╔═╡ bb770f3f-72dd-4a71-8d71-9e773224df05
t = @test always_false(rand(20), rand(20),123)

# ╔═╡ 22a33c8c-e07f-445e-9d8d-a676f704ec45
@test begin
	(1+2) - (3+4)
	false
end

# ╔═╡ 176f39f1-fa36-4ce1-86ba-76248848a834
@test (@test always_false(rand(30),123)) isa Fail

# ╔═╡ 8150cb7e-b2e2-4ee8-a475-db4454c954f0
embed_display(@test false)

# ╔═╡ 716bba60-bfbb-48a4-8924-8bf4e8958cb1
@test always_false("asd"*"asd","asd","asd"*" asd","asd","asd"*"asd","asd","asd"*"asd","asd","asd"*"asd","asd","asd"*"asd","asd","asd"*"asd","asd")

# ╔═╡ 7c1aa057-dff2-48cd-aad5-1bbc1c0a729b
@test π ≈ 3.14 atol=0.01 rtol=1

# ╔═╡ de62537b-a428-48d3-a866-151127b3255b

# ╔═╡ Cell order:
# ╟─ab02837b-79ec-40d7-bff1-c1d2dd7362ef
# ╠═73d74146-8f60-4388-aaba-0dfe4215cb5d
# ╠═71b22e76-2b50-4d16-85f6-9dad0415630e
# ╠═6762ed72-f422-43a9-a782-de78f739c0ae
# ╠═f77275b9-90aa-4e07-a608-981b5df727af
# ╟─56347b7e-5007-45f8-8f6d-8ac8cc719637
# ╟─bf2abe01-6ae0-4066-8704-12f64e04511b
# ╟─e46cf3e0-aa15-4c17-a925-3e9fc5109d54
# ╟─5b70aaf1-9623-4f55-b055-4263ed8be31d
# ╟─fd8428a3-9fa3-471a-8b2d-5bbb8fdb3137
# ╟─191f1f04-18d4-485b-af8b-a2f073b7043b
# ╟─ec1fd70a-d92a-4688-98b2-135879f07141
# ╠═cf314b21-3f4f-4637-b1ce-ec1d5d5af966
# ╠═80b2fb3f-94b2-4024-94ff-d111a249c8b0
# ╠═c763ed72-82c9-445c-a8f7-a0c40982e4d9
# ╠═9d49ea50-8158-4d8b-97af-edba1f7dc38b
# ╠═eab4ba31-c787-46dd-8024-693eca7fd1a0
# ╠═26b0faf0-9016-48d7-8667-c1c1cfce655e
# ╠═26b4fb86-892f-415c-8046-6a5449052fd7
# ╠═96dc7b01-3766-4206-88ba-eca1665bc5cb
# ╠═7c6ce205-053d-434c-b5b1-500babb8ec02
# ╠═a6f82260-3519-4254-a21e-abc7bb19ec4e
# ╠═fe7f8cce-a706-476d-8680-a2fe793b474f
# ╠═1872f785-1ae5-43ea-bce1-6c5cd893f3a8
# ╠═5570972e-9309-4458-99a6-ea718ec2c3ab
# ╠═8d340983-ea07-4038-872f-22a165003ed2
# ╠═ea5a4fc0-db62-41dd-9600-a21d4eabf822
# ╠═c369b4b5-2fcf-4029-a1f6-352120b2fc4b
# ╠═4509cdbf-8b8b-4f70-9e63-bb972eb88c93
# ╠═98992db9-4f14-4aa6-a7c5-477622266112
# ╠═8360d1bc-b1f4-4263-a042-724cbd120227
# ╠═064e28de-0c22-48b5-b427-6eb343880287
# ╠═3a6a6ee1-c619-4044-b9b1-68e5ae9d2463
# ╠═be93a6f4-b626-43db-a2fe-4e754e79c030
# ╟─17bd5cd9-212f-4656-ab79-590dd6c64ff8
# ╟─539e2c38-993b-4b3b-8aa0-f02d46d79839
# ╟─3d3f3592-e056-4e7b-8896-a75e5b5dcad6
# ╠═1aa24b1c-e8ca-4de7-b614-7a3f02b4833d
# ╠═8a2e8348-49cf-4855-b5b3-cdee33e5ed67
# ╠═42671258-07a0-4015-8f47-4b3032595f08
# ╠═05f42764-acfe-4370-b85b-ce0e7c4270d0
# ╟─0d70962a-3880-4dee-a439-35068d019f5a
# ╠═113cc425-e224-4f77-bfbd-ef4eb1d1ed70
# ╠═6188f559-bcab-4da6-84b2-a3fe522a5c3c
# ╠═c24b46ce-bcbb-4dc9-8a59-b5b1bd2cd617
# ╠═5041085e-a406-4ed4-ab82-84d8f126cf0f
# ╠═8c92bad9-234e-47dd-a599-b75dc6d5db89
# ╠═03ccd498-83c3-41bb-84d7-625adabd7aee
# ╠═1bcf8bd1-c8a3-49a1-9791-d813aa856399
# ╠═656c4190-b49e-4225-869d-eeb7e8e41e72
# ╠═14c525a1-eca1-466b-8e63-3a90d7d7111c
# ╟─a2efc968-246c-40c2-b285-2ec94b185a44
# ╠═c39021dc-157c-4bcb-a3a9-fec8d9286b48
# ╠═b6b1e037-754e-4c9c-a79d-2a1517cb409d
# ╠═8020701d-e23a-4d95-b4ad-ebeddd374fb5
# ╠═fc1ec730-4e11-452c-94a7-26605a87a89e
# ╠═e1c306e3-0a47-4149-a9fb-ec7ab380fa11
# ╠═b6e8a170-12cc-4d97-905d-274e2609bfd8
# ╟─bfe4dc61-9160-4c7e-8897-9c723b309adc
# ╠═bb770f3f-72dd-4a71-8d71-9e773224df05
# ╠═22a33c8c-e07f-445e-9d8d-a676f704ec45
# ╠═176f39f1-fa36-4ce1-86ba-76248848a834
# ╠═8150cb7e-b2e2-4ee8-a475-db4454c954f0
# ╠═716bba60-bfbb-48a4-8924-8bf4e8958cb1
# ╠═4bc1b7a4-0a36-4a07-b7ee-3d5be50350e1
# ╠═3b2e8f55-1d4b-4a36-83f6-26becbd79e4b
# ╠═7c1aa057-dff2-48cd-aad5-1bbc1c0a729b
# ╠═ec2ed42c-1227-4e0d-b642-20e6f3503d2a
# ╠═ac02b12a-3982-4526-a51c-0bf85198b81b
# ╠═9c3f6eab-b1c3-4607-add8-d6d7e468c11a
# ╠═1ac164c8-88fc-4a87-a194-60ef616fb399
# ╠═98ac4c36-49c7-4f65-982d-0b8bf6c372c0
# ╠═0fcc6cb0-2711-4609-9bf3-634cf9407840
# ╠═69200d7c-b7bc-4c7e-a9a1-5e26979179a3
# ╠═ea73a35e-a34f-4708-acc1-858f2466e9ba
# ╠═b765dbfe-4e58-4bb9-b1d6-aa4378d4e9c9
# ╠═ef6fc423-f1b1-4dcb-a059-276121391bc6
# ╠═4d5f44e4-85e9-4985-9b76-73be5e097186
# ╠═dd495e00-d74d-47d4-a5d5-422fb147ec3b
# ╟─a83a6a4c-664c-46fa-a07f-81088493dc35
# ╟─5cb03161-2cbc-4080-ba59-f94efd3b620c
# ╟─0611a36b-b4be-4b17-a485-7c4a8fa04927
# ╟─dbfbcc16-c740-436c-bbf0-fee16b0a20c5
# ╠═d97987a0-bdc0-46ed-a6a5-f35c1ce961dc
# ╠═a6709e08-964d-46ea-9813-2c70a834824b
# ╟─7c2bab29-8609-4575-b2ca-7feb34915645
# ╠═69bfb438-7ecf-4f9b-8bc4-51e07aa46ef1
# ╟─838b5904-1de2-4d9f-8f3c-a93ec224d40e
# ╠═a3c41025-2f4a-4f9c-8577-72e4b7abbb98
# ╠═3e79ff61-6532-4879-9402-86473aa7d960
# ╠═275c5f57-623d-439f-b09d-f7c745e0bed6
# ╠═10803c0d-d0a5-45c5-b7ef-9659e441df69
# ╟─b056a99d-5b13-47ba-a199-d788410e3c99
# ╟─5b093e83-78c1-4187-b406-56e79800e1be
# ╠═a461f1fd-b5a5-4ae3-a47c-067a6081fb24
# ╠═38e54516-cdf4-4c1d-815b-68e1e7a7f6f7
# ╠═411271a6-4236-45e2-ab34-f26410108821
# ╟─f9c81ab1-556c-4d81-bee8-2897c20e324d
# ╠═a392d2d6-5a16-4383-b0ef-5003aa2de9fa
# ╠═ae95b691-f54b-4bf5-b17b-3e5bd1edf75e
# ╠═12119016-fa61-4d38-8c58-821ea435df7d
# ╠═9bed78b6-5a8f-44ce-ab66-cab685daf264
# ╟─2c1b906d-71b9-430e-83ed-d4c8c0018632
# ╠═b5a9b7c1-43fb-4ec6-aaeb-0ec55580ef64
# ╠═450a36ea-2c43-4f01-a775-0b8c59bf6dca
# ╠═ef45af85-11f8-4505-9f1e-3ffb15a47142
# ╠═97a05d7c-2e08-458d-a628-7992f204b4ea
# ╠═a661e172-6afb-42ff-bd43-bb5b787ee5ed
# ╠═930f8244-cf25-4c1a-95f6-5c8963559c62
# ╠═68ba60db-44ad-43e4-b33e-d27696babc99
# ╠═807bcd72-26c3-44d3-a295-56874cb51a89
# ╠═8a5a4c26-e36c-4061-b32f-4448625ce4a6
# ╠═8ef356ea-7d54-43e6-a936-7c8be04c595f
# ╠═88f6a040-07cf-47e0-a8be-2478ea350aa7
# ╠═d36a8a72-eced-4e63-9130-7fcb6c86df76
# ╠═e9659020-d433-4357-9099-71a65b66a091
# ╠═c3fc3292-b7eb-4b01-8fba-159c86228de9
# ╠═89578bff-16b9-4eb2-b8ee-b2839ff2d74c
# ╟─21d4560e-721f-4ed4-9db7-86a8151ab22c
# ╟─99afc7f4-727c-4277-8311-f2ffa94830ae
# ╠═4956526a-daf9-43c9-bff3-ff2446016e2e
# ╠═84ff6a23-c134-4910-b630-a7ad45f3bf29
# ╠═318363d0-6d9e-4144-b478-b775f437edaf
# ╠═67fd07b7-340b-4e24-bc06-e4c85b186872
# ╠═4edf747b-3838-4315-a397-e452ac9b5465
# ╟─c6d5597c-d505-4125-88c4-10415934d2a4
# ╠═872b4877-30dd-4a92-a3c8-69eb50675dcb
# ╠═c877c109-db16-468c-8f3c-8294db859d6d
# ╠═ab0a19b8-cf7c-4c4f-802a-f85eef81fc02
# ╠═8480d0d7-bdf7-468d-9344-5b789e33921c
# ╠═6f5ba692-4b6a-405a-8cd3-1a8f9cc06611
# ╠═b4b317d7-bed1-489c-9650-8d336e330689
# ╠═93ed973f-daf6-408b-9d4b-d53495418610
# ╠═dea898a0-1904-4d09-ad0b-6915008fe946
# ╟─b5763c10-e11c-4389-b6fc-421d2c9682f1
# ╠═74c19786-1ba7-4865-a993-590a779ae564
# ╠═e968fc57-d850-4e2d-9410-8777d03b7b3c
# ╟─3d5abd58-02ab-4b91-a7a3-d9068d4df017
# ╠═326f7661-3482-4bf2-a97b-57cc7ac60ee2
# ╟─b273d3d3-648f-4d34-94e7-e49277d4ba29
# ╠═a2cbb0c3-23b9-4091-9ca7-5ba96e85e3a3
# ╟─34f613a3-85fb-45a8-be3b-cd8e6b3cb5a2
# ╟─f9ed2487-a7f6-4ce9-b673-f8a298cd5fc3
# ╟─35f63c4e-3583-4ea8-a057-31f18f8a09d6
# ╟─35b2770e-1db6-4327-bf86-c27a4b61dbd3
# ╟─22640a2f-ea38-4517-a4f3-7a65e60ffebe
# ╟─d414f840-4952-4de5-a565-7fdc81a94817
# ╠═326825b0-a17f-427a-9056-8e8156098418
# ╟─64bf02a4-4fe3-424d-ae6e-5906c3395278
# ╟─f3916810-1911-48bd-936b-776206fcad54
# ╟─122c27a5-a6e8-45ef-a968-b9b4b3f9ad09
# ╟─9126f47d-cbc7-411f-93bd-8684ba06c9e9
# ╟─955705f9-c90d-495d-86b4-5f3b5bc9fc8e
# ╟─187c3005-cd43-45a0-8cbd-bc96b9cb39da
# ╠═daee414b-3e3c-4e2a-a25a-429a1e7275d5
# ╟─6c0156a9-7281-4326-9e1f-989efa73bb7b
# ╠═de62537b-a428-48d3-a866-151127b3255b

ext: jl


text:
"""
    Vasicek()

[Vasicek model](https://en.wikipedia.org/wiki/Vasicek_model).
"""
struct Vasicek <: Model end

"""
    evaluate(IRD, Vasicek(), n = 12)

Evaluate interest rate derivative `IRD` using `Vasicek` model.

# Arguments
- `n`: number of paths to simulate

"""
function evaluate(IRD::InterestRateDerivative, m::Vasicek, n::Int64 = 12)
    Δt = IRD.t / n
    rates = [IRD.r]
    for i = 0:n
        dr = IRD.k * (IRD.θ - rates[end]) * Δt + IRD.σ * rates[end] * randn()
        append!(rates, rates[end] + dr)
    end
    return rates
end


text:
module InteractiveChaos

using Reexport
@reexport using DynamicalSystems
using AbstractPlotting, Observables

"""
    subscript(i::Int)
Transform `i` to a string that has `i` as a subscript.
"""
function subscript(i::Int)
    if i < 0
        "₋"*subscript(-i)
    elseif i == 1
        "₁"
    elseif i == 2
        "₂"
    elseif i == 3
        "₃"
    elseif i == 4
        "₄"
    elseif i == 5
        "₅"
    elseif i == 6
        "₆"
    elseif i == 7
        "₇"
    elseif i == 8
        "₈"
    elseif i == 9
        "₉"
    elseif i == 0
        "₀"
    elseif i ≥ 10
        join(subscript.(digits(i)))
    else
        string(i)
    end
end

include("plot_dataset.jl")
include("orbitdiagram.jl")
include("poincareclick.jl")
include("trajectory_highlighter.jl")

export interactive_orbitdiagram

end


text:
#TODO: Renderable type
abstract type Renderable end

"""
Return true if an object or type is *directly renderable*, false otherwise.

New types should implement the `isrenderable(t::Type{NewType})` method.
"""
function isrenderable end

isrenderable(object) = isrenderable(typeof(object))
isrenderable(::Type{R}) where R <: Renderable = true
isrenderable(t::Type) = hasmethod(render!, Tuple{RenderModel, t})
isrenderable(t::Type{Roadway}) = true

"""
    render(scene)
    render(scene; kwargs...)

Render all the items in `scene` to a Cairo surface and return it.

`scene` is simply an iterable object (e.g. a vector) of items that are either directly renderable or renderable by conversion. See the AutoViz README for more details.
"""
function render(scene; # iterable of renderable objects
                overlays=[],
                rendermodel::RenderModel=RenderModel(),
                cam::Camera=FitToContentCamera(),
                canvas_height::Int=DEFAULT_CANVAS_HEIGHT,
                canvas_width::Int=DEFAULT_CANVAS_WIDTH,
                surface=CairoRGBSurface(canvas_width, canvas_height)
               )

    ctx = creategc(surface)
    clear_setup!(rendermodel)

    for x in scene
        if isrenderable(x)
            render!(rendermodel, x)
        else
            render!(rendermodel, convert(Renderable, x))
        end
    end

    for o in overlays
        render!(rendermodel, o, scene)
    end

    camera_set!(rendermodel, cam, scene, canvas_width, canvas_height)

    render(rendermodel, ctx, canvas_width, canvas_height)
    return surface
end

function Base.convert(::Type{Renderable}, x::AbstractVector{R}) where R <: Real
    @assert length(x) >= 2
    if length(x) == 2
        ac = ArrowCar(SVector(x[1],x[2]))
    else
        ac = ArrowCar(SVector(x[1],x[2]), x[3])
    end
    ac
end


text:
const abspathpattern =
  @windows? r"([a-zA-Z]+:[\\/][a-zA-Z_\./\\ 0-9]+\.jl)(?::([0-9]*))?" : r"(/[a-zA-Z_\./ 0-9]+\.jl)(?::([0-9]*))?"

# Make the prefix optional, but disallow spaces
const relpathpattern =
  @windows? r"([a-zA-Z_\./\\0-9]+\.jl)(?::([0-9]*))?$" : r"([a-zA-Z_\./0-9]+\.jl)(?::([0-9]*))?$"

function basepath(file)
  path = joinpath(JULIA_HOME,"..","share","julia","base",file) |> normpath
  isfile(path) || (path = nothing)
  return path
end

function githublink(path)
  file, line = match(relpathpattern, path).captures
  """$path"""
end

function baselink(path)
  m = match(relpathpattern, path)
  m == nothing && return path
  file, line = m.captures
  link = basepath(file)
  link == nothing && return githublink(path)
  line == nothing || (link *= ":$line")
  """$path"""
end

function filelink(path)
  """$(basename(path))"""
end

function showerror_html(io, e)
  print(io, """""")
  showerror(io, e)
  print(io, """""")
end

function showerror_html(io, e, bt, top_function = :eval_user_input)
  println(io, """
""")
  showerror_html(io, e)
  showbacktrace_html(io, top_function, bt)
  println(io, """
""")
end

function showerror_html(io::IO, e::LoadError, bt, top_function = :eval_user_input)
  println(io, """
""")
  showerror_html(io, e.error)
  showbacktrace_html(io, top_function, bt)
  println(io, """while loading $(e.file), in expression starting on line $(e.line)""")
  println(io, """
""")
end

function showbacktrace_html(io::IO, top_function::Symbol, t, set = 1:typemax(Int))
  ls = map(strip, lines(sprint(Base.show_backtrace, top_function, t, set)))
  print(io, """
""")

  for i = 2:length(ls)
    print(io, """
""")
    if ismatch(abspathpattern, ls[i])
      print(io, replace(ls[i], abspathpattern, filelink))
    elseif ismatch(relpathpattern, ls[i])
      print(io, replace(ls[i], relpathpattern, baselink))
    else
      println(io, ls[i])
    end
    print(io, """
""")
  end

  println(io, """
""")
end

# Methods

Requires.@init Jewel.@inmodule Base begin

stripparams(t) = replace(t, r"\{([A-Za-z, ]*?)\}", "")

function writemime(io::IO, ::MIME"text/html", m::Method)
  print(io, "
")
  print(io, m.func.code.name)
  tv, decls, file, line = Base.arg_decl_parts(m)
  if !isempty(tv)
    print(io,"")
    # Redundant, since type parameters appear in the parameter list
    #Base.show_delim_array(io, tv, '{', ',', '}', false)
    print(io,"")
  end
  print(io, "(")
  print_joined(io, [isempty(d[2]) ? d[1] : d[1]*"::"*stripparams(d[2])*""
                    for d in decls], ",", ",")
  print(io, ")")
  if line > 0
    file = "$file:$line"
    print(io, "	")
    print(io, isabspath(file) ? Main.LightTable.filelink(file) : Main.LightTable.baselink(file))
    print(io, "
")
  end
end

function writemime(io::IO, m::MIME"text/html", mt::MethodTable)
  defs = Method[]
  d = mt.defs
  while !is(d,())
    push!(defs, d)
    d = d.next
  end
  writemime(io, m, defs)
end

function writemime(io::IO, m::MIME"text/html", mt::Vector{Method})
  print(io, """

""")
  file(m) = m.func.code.file |> string |> basename
  line(m) = m.func.code.line
  map(d->writemime(io, m, d), sort!(mt, lt = (a, b) -> file(a) == file(b) ?
                                      line(a) < line(b) :
                                      file(a) < file(b)))
  print(io, """
""")
  print(io, """
""")
end

end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule SQLCipher_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("SQLCipher")
JLLWrappers.@generate_main_file("SQLCipher", UUID("dc592094-f866-5418-b15b-f7542ad3543c"))
end  # module SQLCipher_jll


text:
# by Lazaro Alonso
# using GLMakie # HIDE
using GLMakie
GLMakie.activate!() # HIDE
let
    x = 1:10
    y = 1:10
    z = 1:10
    f(x,y,z) = x^2 + y^2 + z^2
    positions = vec([(i, j, k) for i in x,j in y, k in z])
    vals = [f(ix,iy,iz) for ix in x, iy in y, iz in z]
    fig, ax, pltobj = meshscatter(positions, color = vec(vals),
        marker = FRect3D(Vec3f0(0), Vec3f0(10)), # here, if you use less than 10, you will see smaller squares.
        colormap = :Spectral_11, colorrange = (minimum(vals), maximum(vals)),
        transparency = true, # set to false, if you don't want the transparency.
        shading= false,
        figure = (; resolution = (800,800)),
        axis=(; type=Axis3, perspectiveness = 0.5,  azimuth = 7.19, elevation = 0.57,
            xlabel = "x label", ylabel = "y label", zlabel = "z label",
            aspect = (1,1,1)))
    cbar = Colorbar(fig, pltobj, label = "f values", height = Relative(0.5))
    xlims!(ax,-1,11)
    ylims!(ax,-1,11)
    zlims!(ax,-1,11)
    fig[1,2] = cbar
    fig
    #save("fileName.png", fig) # here, you save your figure.
    save(joinpath(@__DIR__, "output", "volumeScatters.png"), fig, px_per_unit = 2.0) # HIDE
end
using Pkg # HIDE
Pkg.status(["GLMakie"]) # HIDE

max_stars_repo_name: MisterBiggs/BeautifulMakie

max_issues_repo_name: MisterBiggs/BeautifulMakie

max_forks_repo_name: MisterBiggs/BeautifulMakie


text:
#@testset "CronTrigger" begin
#    @test 1 == 1
#end


text:
function UA_ServerConfig_setMinimal(
    config::Ref{UA_ServerConfig},
    portNumber::Integer,
    certificate::Ref{<:Union{Nothing, UA_ByteString}},
)
    UA_ServerConfig_setMinimalCustomBuffer(config, portNumber, certificate, 0, 0)
end

UA_ServerConfig_setDefault(config::Ref{UA_ServerConfig}) = UA_ServerConfig_setMinimal(config, 4840, C_NULL)

## Read functions
for att in attributes_UA_Server_read
    fun_name = Symbol(att[1])
    attr_name = Symbol(att[2])
    ret_type = Symbol(att[3])
    ret_type_ptr = Symbol("UA_TYPES_", uppercase(String(ret_type)[4:end]))
    ua_attr_name = Symbol("UA_ATTRIBUTEID_", uppercase(att[2]))

    @eval begin
        function $(fun_name)(server::Ptr{UA_Server}, nodeId::Ref{UA_NodeId})
            out = Ref{$(ret_type)}()
            retval = __UA_Server_read(server, nodeId, $(ua_attr_name), out)
            if retval == UA_STATUSCODE_GOOD
                return out[]
            else
                error("Reading ´$(attr_name)´ from UA_Server failed with statuscode \"$(UA_StatusCode_name_print(retval))\".")
            end
        end

        function $(fun_name)(server::Ptr{UA_Server}, nodeId::UA_NodeId)
            return $(fun_name)(server, Ref(nodeId))
        end
    end
end

## Write functions
for att in attributes_UA_Server_write
    fun_name = Symbol(att[1])
    attr_name = Symbol(att[2])
    attr_type = Symbol(att[3])
    attr_type_ptr = Symbol("UA_TYPES_", uppercase(String(attr_type)[4:end]))
    ua_attr_name = Symbol("UA_ATTRIBUTEID_", uppercase(att[2]))

    @eval begin
        function $(fun_name)(server::Ptr{UA_Server}, nodeId::Ref{UA_NodeId}, new_val::$(attr_type))
            data_type_ptr = UA_TYPES_PTRS[$(attr_type_ptr)]
            retval = __UA_Server_write(server, nodeId, $(ua_attr_name), data_type_ptr, new_val)
            if retval == UA_STATUSCODE_GOOD
                return retval
            else
                error("Writing ´$(attr_name)´ on UA_Server failed with statuscode \"$(UA_StatusCode_name_print(retval))\".")
            end
        end

        function $(fun_name)(server::Ptr{UA_Server}, nodeId::UA_NodeId, new_val::$(attr_type))
            return $(fun_name)(server, Ref(nodeId), new_val)
        end
    end
end


text:
using SpikingNN
using Distributions
using Distributions: sample

"""
Delay Synapses Have been replaced by
Queued Alpha Synpases, to satisfy an annoying syntactic problem, regarding syntactic statements, and assigning functions to functions
In the Delayed Synapse function definition.
"""

# simulation parameters
T = 1.0
dt = 0.1f-3
n = Int(ceil(T / dt))

# populations parameters
Ne = 800
Ni = 200

# connectivity parameters
w = 0.1
wexcite = 0.1
winhibit = -0.5
sparsity = 0.1
λ = 1f-2

# neuron parameters
τm = 20f-3
τr = 2f-3
vth = 20f-3
delay = Int(ceil(1.5f-3 / dt))

# input parameters
ρ0 = 20f0
Ninput = Ne + Ni

# create system
D = Bernoulli(1 - sparsity)
inputs = InputPopulation([ConstantRate(ρ0, dt) for _ in 1:Ninput])
W_EE = Float32.(rand(D, Ne, Ne) .* wexcite)
W_II = Float32.(rand(D, Ni, Ni) .* winhibit)
W_EI = Float32.(rand(D, Ne, Ni) .* wexcite)
W_IE = Float32.(rand(D, Ni, Ne) .* winhibit)
W_input_E = Float32.(rand(D, Ninput, Ne) .* wexcite)
W_input_I = Float32.(rand(D, Ninput, Ni) .* wexcite)

# neuron parameters
vᵣ = 0
τᵣ = 1.0
vth = 1.0

E = Population(W_EE; cell = () -> LIF(τᵣ, vᵣ),
                     threshold = () -> Threshold.Ideal(vth),
                     synapse = () -> QueuedSynapse(Synapse.Alpha()))
I = Population(W_II; cell = () -> LIF(τᵣ, vᵣ),
                     threshold = () -> Threshold.Ideal(vth),
                     synapse = () -> QueuedSynapse(Synapse.Alpha()))
net = Network(Dict(:input => inputs, :E => E, :I => I))

connect!(net, :E, :I; weights = W_EI, synapse = () -> QueuedSynapse(Synapse.Alpha()))
connect!(net, :I, :E; weights = W_IE, synapse = () -> QueuedSynapse(Synapse.Alpha()))
connect!(net, :input, :E; weights = W_input_E, synapse = () -> QueuedSynapse(Synapse.Alpha()))
connect!(net, :input, :I; weights = W_input_I, synapse = () -> QueuedSynapse(Synapse.Alpha()))

# recording callback
nsample = min(Ne, Ni, 25)
const record_excite = sample(1:Ne, nsample; replace = false)
const record_inhibit = sample(1:Ni, nsample; replace = false)
Ve = Dict{Int, Vector{Float32}}()
Vi = Dict{Int, Vector{Float32}}()
function record()
    global Ve, Vi, record_excite, record_inhibit
    for idx in record_excite
        push!(get!(Ve, idx, Float32[]), getvoltage(net[:E][idx]))
    end
    for idx in record_inhibit
        push!(get!(Vi, idx, Float32[]), getvoltage(net[:I][idx]))
    end
end

# simulate

# n is the simulation time step vector

@time simulate!(net, n; cb = record, dt = dt)


text:
export TMJets

using IntervalArithmetic: IntervalBox

struct TMJets <: ContinuousPost
    options::TwoLayerOptions

    function TMJets(𝑂::Options)
        𝑂new = validate_and_wrap_options(𝑂, options_TMJets())
        return new(𝑂new)
    end
end

# convenience constructor from pairs of symbols
TMJets(𝑂::Pair{Symbol, <:Any}...) = TMJets(Options(Dict{Symbol, Any}(𝑂)))

# default options (they are added in the function validate_and_wrap_options)
TMJets() = TMJets(Options())

include("init.jl")
include("post.jl")
include("project.jl")
include("reach.jl")


text:
using Glob;
using DelimitedFiles;
using GLM;
using DataFrames;
using Logging;

if length(ARGS) != 2
  println("usage: julia statistical_analysis_mcmc.jl <indir> <subpattern>");
  exit(1);
end

global_logger(ConsoleLogger(stderr, Logging.Info));

indir = ARGS[1];
subpattern = ARGS[2];
pattern = Glob.GlobMatch(subpattern * ".out");

chunks = map(piece -> split(piece, "-"), split(subpattern, "_"));
param_map = Dict();
for chunk in chunks
  if chunk[1] == "run"; continue; end
  param_map[chunk[1]] = Meta.parse(chunk[end]) / 1000.0;
end
@show param_map;

data_serieses = Dict();
for subdir in ["standard"; "clustering"; "umbrella"]
  data_series = Dict();
  for infile in readdir(pattern, joinpath(indir, subdir))
    for line in readlines(infile)
      var, val = map(strip, split(line, "="));
      val = eval(Meta.parse(val));
      if haskey(data_series, var)
        push!(data_series[var], val);
      else
        data_series[var] = [val];
      end
    end
  end
  data_serieses[subdir] = data_series;
end

Ustar = (param_map["E0"])^2 * abs(param_map["K2"] - param_map["K1"]);
rstar = param_map["b"];
μstar = (param_map["E0"])*(2*param_map["K2"] + param_map["K1"]);
for (name, data_series) in data_serieses
  @show name;
  avgs = Dict();
  vars = Dict();
  std_norms = Dict();
  for (k, v) in data_series
    avgs[k] = sum(v) / length(v);
    vars[k] = sum(map(x -> (x - avgs[k]).^2, v)) / length(v);
    std_norms[k] = map(sqrt, vars[k]) ./ avgs[k];
  end
  for k in keys(avgs)
    @info "$k: μ = $(avgs[k]); σ^2 = $(vars[k]); σ / μ = $(std_norms[k])";
  end
end

max_stars_repo_name: grasingerm/polymer-stats

max_issues_repo_name: grasingerm/polymer-stats

max_forks_repo_name: grasingerm/polymer-stats


text:
if Base.HOME_PROJECT[] !== nothing
    # JuliaLang/julia/pull/28625
    Base.HOME_PROJECT[] = abspath(Base.HOME_PROJECT[])
end

using Documenter, TSVD

makedocs()

deploydocs(
    deps   = Deps.pip("mkdocs==0.17.5", "python-markdown-math", "mkdocs-material==2.9.4"),
    repo = "github.com/andreasnoack/TSVD.jl.git",
    julia  = "1.0"
)


text:
# run hopper_in_place_hop.jl

include(joinpath(@__DIR__, "..", "src/dynamics", "hopper_2D", "visuals.jl"))
T = Float64
vis = Visualizer()
open(vis)

# get hopper model
# model_sim = get_model("hopper_2D", surf="sinusoidal")
model_sim = get_model("hopper_2D", surf="flat")
nq = model_sim.dim.q
nu = model_sim.dim.u
nc = model_sim.dim.c
nb = model_sim.dim.b

# time
H = traj.H
h = traj.h
N_sample = 5
H_mpc = 10
h_sim = h / N_sample
H_sim = 3000

# Select the initial speed
v0 = 0.0
Tstance = 0.13 # measure using hop-in-place gait
Tflight = 0.62 # measure using hop-in-place gait
include(joinpath(module_dir(), "src/controller/raibert_policy.jl"))
p = raibert_policy(model, v0=v0, Tstance=Tstance, Tflight=Tflight, h=h)

off0 = SVector{nq,T}([0.0, 0.5, 0.0, 0.0])
off1 = SVector{nq,T}([v0*h_sim, 0.5, 0.0, 0.0])
q_ref = SVector{nq,T}([0.0, 0.5, 0.0, 0.5])
q0_sim = copy(q_ref) + off0
q1_sim = copy(q_ref) + off1

# Disturbances
w = [zeros(nw) for t=1:Int(ceil(H_sim/N_sample))]
w[75] += [5.0, -0.0]
d = open_loop_disturbances(w)

sim = simulator(model_sim, q0_sim, q1_sim, h_sim, H_sim,
    p = p,
    d = d,
    ip_opts = InteriorPointOptions(
        r_tol = 1.0e-8,
        κ_init = 1.0e-8,
        κ_tol = 2.0e-8),
    sim_opts = SimulatorOptions(warmstart = true)
    )

@time status = simulate!(sim)

plt = plot(layout=(3,1), legend=false)
plot!(plt[1,1], hcat(Vector.(vcat([fill(ref_traj.q[i], N_sample) for i=1:H]...))...)',
    color=:red, linewidth=3.0)
plot!(plt[1,1], hcat(Vector.(sim.traj.q)...)', color=:blue, linewidth=1.0)
plot!(plt[2,1], hcat(Vector.(vcat([fill(ref_traj.u[i][1:nu], N_sample) for i=1:H]...))...)',
    color=:red, linewidth=3.0)
plot!(plt[2,1], hcat(Vector.([u[1:nu] for u in sim.traj.u]*N_sample)...)', color=:blue, linewidth=1.0)
plot!(plt[3,1], hcat(Vector.([γ[1:nc] for γ in sim.traj.γ]*N_sample)...)', color=:blue, linewidth=1.0)

plot_lines!(vis, model_sim, sim.traj.q[1:10:end])
plot_surface!(vis, model_sim.env, n=200)
anim = visualize_robot!(vis, model_sim, sim.traj, sample=20)
anim = visualize_force!(vis, model_sim, sim.traj, anim=anim, h=h_sim, sample=20)

plot(hcat([q[[1,2,3,4]] for q in Vector.(sim.traj.q)[1:1200]]...)')
plot(hcat([q[[3]] for q in Vector.(sim.traj.q)[1:end]]...)')
plot(hcat([γ[[1]] for γ in Vector.(sim.traj.γ)[1:end]]...)')
plot(hcat([u[[1,2]] for u in Vector.(sim.traj.u)[1:end]]...)')

# filename = "raibert_hopper_push"
# MeshCat.convert_frames_to_video(
#     "/home/simon/Downloads/$filename.tar",
#     "/home/simon/Documents/$filename.mp4", overwrite=true)
#
# convert_video_to_gif(
#     "/home/simon/Documents/$filename.mp4",
#     "/home/simon/Documents/$filename.gif", overwrite=true)


text:
"""
Description

       The @tux macro scaffolds both method source and test files. The source file includes a documentation template.

Usage

       @tux MyPackage foo

Method Signature(s)

      @tux(pkg, method)

Details

      The tux macro is used for its side=effects, which is to create a documentation=stubbed file for the method to be written, and a pending
      test file. It takes two arguments: the package and method names. The auto-generated method file is placed in the src/ directory of the
      package that is listed, while the auto-generated test file is placed in the test/ directory of the package. The macro determines what
      version of Julia is being used and places the two files in the appropriate .julia/ directory. After creating these files from REPL, an
      editor (currently only vi is supported) is opened with both method files and test files in the buffer.

References

      https://github.com/milktrader/Tuxedo

"""

macro tux(pkg, method)

    v    = match(r"[0-9]", string(VERSION), 2)
    ver  = v.match
    home = homedir()
    cd("$home/.julia/v0.$ver/$pkg/")

    src  = Pkg.dir("Tuxedo/src/src.jl")
    test = Pkg.dir("Tuxedo/src/test.jl")

    run(pipeline(`cp $src src/$method.jl`))
    run(pipeline(`cp $test test/$method.jl`))
    run(pipeline(`ls`, `cat src/$method.jl`, pipeline(`sed s/foo/$method/g`, stdout="src/$method.jl")))
    # TODO: make vi an option
    run(`vi src/$method.jl test/$method.jl`)
end

max_stars_repo_name: milktrader/Tuxedo

max_issues_repo_name: milktrader/Tuxedo

max_forks_repo_name: milktrader/Tuxedo


text:
module day05

include("inputs.jl")
input = inputs.input_d5()

lookup = Dict('B'=>'1', 'F'=>'0', 'R'=>'1', 'L'=>'0')

function parse_bp(bp)
    row = parse(Int, join(lookup[c] for c in bp[1:7]), base=2)
    col = parse(Int, join(lookup[c] for c in bp[8:end]), base=2)
    seat_id = 8*row + col
    return seat_id
end

@assert parse_bp("BFFFBBFRRR") == 567
@assert parse_bp("FFFBBBFRRR") == 119
@assert parse_bp("BBFFBBFRLL") == 820

# Part 1
function highest_id(input)
    seat_ids = parse_bp.(input)
    return maximum(seat_ids)
end
@info "Day 5, Part 1 answer: $(highest_id(input))"

#Part 2
# function find_seat(;input=input())
function find_seat(input)
    seat_ids = parse_bp.(input)
    lo = minimum(seat_ids)
    hi = maximum(seat_ids)
    return [x for x in lo:hi if x ∉ seat_ids]
end
@info "Day 5, Part 2 answer: $(find_seat(input)[1])"

end


text:
using Parameters

abstract type AbstractThermo end
export AbstractThermo

@inline function getGibbs(th::P,T::N) where {N<:Number,P<:AbstractThermo}
    return @fastmath getEnthalpy(th,T)-T*getEntropy(th,T)
end

@with_kw struct NASApolynomial <: AbstractThermo
    coefs::Array{Float64,1}
    Tmin::Float64
    Tmax::Float64
end
export NASApolynomial

@inline function calcHSCpdless(poly::NASApolynomial,T::X) where {X<:Real}
    if length(poly.coefs) != 7
        Tpoly0 = T
        Tpoly1 = T*T
        Tpoly2 = Tpoly1*T
        Tpoly3 = Tpoly2*T
        Tpoly4 = 1.0/T
        Tpoly5 = Tpoly4*Tpoly4
        Tpoly6 = log(T)

        ct0 = poly.coefs[1]*Tpoly5
        ct1 = poly.coefs[2]*Tpoly4
        ct2 = poly.coefs[3]
        ct3 = poly.coefs[4]*Tpoly0
        ct4 = poly.coefs[5]*Tpoly1
        ct5 = poly.coefs[6]*Tpoly2
        ct6 = poly.coefs[7]*Tpoly3

        cpdivR = ct0 + ct1 + ct2 + ct3 + ct4 + ct5 + ct6

        hdivRT = -ct0+Tpoly6*ct1+ct2+0.5*ct3+0.33333333333*ct4+0.25*ct5+0.2*ct6+poly.coefs[8]*Tpoly4
        sdivR = -0.5*ct0 - ct1 + Tpoly6*ct2 + ct3 + 0.5*ct4 + 0.33333333333*ct5 + 0.25*ct6 + poly.coefs[9]
    else
        Tpoly0 = T
        Tpoly1 = T*T
        Tpoly2 = Tpoly1*T
        Tpoly3 = Tpoly2*T
        Tpoly4 = 1.0/T
        #Tpoly5 = Tpoly4*Tpoly4
        Tpoly6 = log(T)

        #ct0 = 0.0
        #ct1 = 0.0
        ct2 = poly.coefs[1]
        ct3 = poly.coefs[2]*Tpoly0
        ct4 = poly.coefs[3]*Tpoly1
        ct5 = poly.coefs[4]*Tpoly2
        ct6 = poly.coefs[5]*Tpoly3

        cpdivR = ct2 + ct3 + ct4 + ct5 + ct6

        hdivRT =ct2+0.5*ct3+0.33333333333*ct4+0.25*ct5+0.2*ct6+poly.coefs[6]*Tpoly4
        sdivR = Tpoly6*ct2 + ct3 + 0.5*ct4 + 0.33333333333*ct5 + 0.25*ct6 + poly.coefs[7]
    end
    return (cpdivR,hdivRT,sdivR)::Tuple{Float64,Float64,Float64}
end
export calcHSCpdless

@inline function getHeatCapacity(poly::NASApolynomial,T::N) where {N<:Number}
    if length(poly.coefs) == 9
        return @views @inbounds @fastmath evalpoly(T,poly.coefs[1:7])/T^2*R
    elseif length(poly.coefs) == 7
        return @views @inbounds @fastmath evalpoly(T,poly.coefs[1:5])*R
    else
        throw(error("NASA polynomial has a number of coefficients not equal to 9 or 7"))
    end
end

@inline function getEntropy(poly::NASApolynomial,T::N) where {N<:Number}
    if length(poly.coefs) == 9
        return @views @inbounds @fastmath ((-poly.coefs[1]/(2*T)-poly.coefs[2])/T+poly.coefs[3]*log(T)+T*evalpoly(T,poly.coefs[4:end-2]./(1:4))+poly.coefs[end])*R
    elseif length(poly.coefs) == 7
        return @views @inbounds @fastmath (poly.coefs[1]*log(T)+T*evalpoly(T,poly.coefs[2:end-2]./(1:4))+poly.coefs[end])*R
    else
        throw(error("NASA polynomial has a number of coefficients not equal to 9 or 7"))
    end
end

@inline function getEnthalpy(poly::NASApolynomial,T::N) where {N<:Number}
    if length(poly.coefs) == 9
        return @views @inbounds @fastmath ((-poly.coefs[1]/T+poly.coefs[2]*log(T))/T+evalpoly(T,poly.coefs[3:end-2]./(1:5)))*R*T+poly.coefs[end-1]*R
    elseif length(poly.coefs) == 7
        return @views @inbounds @fastmath evalpoly(T,poly.coefs[1:end-2]./(1:5))*R*T+poly.coefs[end-1]*R
    else
        throw(error("NASA polynomial has a number of coefficients not equal to 9 or 7"))
    end
end

@with_kw struct NASA{T<:AbstractThermoUncertainty} <: AbstractThermo
    polys::Array{NASApolynomial,1}
    unc::T = EmptyThermoUncertainty()
end
export NASA

@inline function selectPoly(nasa::NASA,T::N) where {N<:Real}
    """
    retrieve the nasa polynomial corresponding to the T range
    """
    for p in nasa.polys
        if T<=p.Tmax
            return p
        end
    end
    return nasa.polys[end]
end
export selectPoly

@inline getHeatCapacity(nasa::NASA,T::N) where {N<:Number} = getHeatCapacity(selectPoly(nasa,T),T)
@inline getEntropy(nasa::NASA,T::N) where {N<:Number} = getEntropy(selectPoly(nasa,T),T)
@inline getEnthalpy(nasa::NASA,T::N) where {N<:Number} = getEnthalpy(selectPoly(nasa,T),T)
@inline getGibbs(nasa::NASA,T::N) where {N<:Number} = getGibbs(selectPoly(nasa,T),T)
@inline calcHSCpdless(nasa::NASA,T::N) where {N<:Real}= calcHSCpdless(selectPoly(nasa,T),T)

@with_kw struct Wilhoit{N,Q,T,P,U,R<:Number,M<:AbstractThermoUncertainty} <: AbstractThermo
    Cp0::N
    Cpinf::T
    coefs::Array{Q,1}
    H0::P
    S0::U
    B::R
    unc::M = EmptyThermoUncertainty()
end
export Wilhoit

@inline function getHeatCapacity(w::Wilhoit,T::N) where {N<:Number}
    @fastmath y = T/(T+w.B)
    return @fastmath w.Cp0 + (w.Cpinf-w.Cp0)*y^2*(1+(y-1)*evalpoly(y,w.coefs))
end

@inline function getEnthalpy(w::Wilhoit,T::N) where {N<:Number}
    @fastmath y = T/(T+w.B)
    return @views @fastmath @inbounds w.H0 + w.Cp0 * T - (w.Cpinf - w.Cp0) * T * (
            y * y * ((3 * w.coefs[1] + sum(w.coefs[2:end])) / 6. +
                     (4 * w.coefs[2] + sum(w.coefs[3:end])) * y / 12. +
                     (5 * w.coefs[3] + w.coefs[4]) * y^2 / 20. +
                     w.coefs[4] * y^3 / 5.) +
            (2 + sum(w.coefs)) * (y / 2. - 1 + (1.0 / y - 1.) * log(w.B + T))
        )
end

@inline function getEntropy(w::Wilhoit,T::N) where {N<:Number}
    @fastmath y = T/(T+w.B)
    return @fastmath w.S0 + w.Cpinf*log(T)-(w.Cpinf-w.Cp0)*(log(y)+y*(1+y*evalpoly(y,w.coefs./(2:5))))
end

export getGibbs, getEntropy, getEnthalpy, getHeatCapacity


text:
abstract type Object end

# Avoid accidental infinite recursion:
Base.show(io::IO, ::MIME"text/plain", x::Object) = invoke(show, Tuple{IO,Any}, io, x)

function Base.show(io::IO, x::Object)
    if get(io, :limit, false)
        # Support `Vector{typeof(x)}` etc. with saner default:
        show(IOContext(io, :compact => true), MIME"text/plain"(), x)
        return
    end
    invoke(show, Tuple{IO,Any}, io, x)
    return
end

# TODO: lock-free
mutable struct Swappable{T,Lock}
    value::T
    lock::Lock
end

Swappable{T,L}(value) where {T,L} = Swappable{T,L}(value, L())

swap!(x::Swappable{T}, value::T) where {T} = something(swapif!(_ -> true, x, value))
function swapif!(f, x::Swappable{T}, value::T) where {T}
    lock(x.lock) do
        old = x.value
        if f(old)
            x.value = value
            return Some(old)
        else
            return nothing
        end
    end
end

function locked(f::F, x::Swappable) where {F}
    lock(x.lock) do
        f(x.value)
    end
end

const DEFAULT_SIZE = 256

const SwappableVector{T,L} = Swappable{Vector{T},L}
SwappableVector{T,L}() where {T,L} =
    SwappableVector{T,L}(empty!(Vector{T}(undef, DEFAULT_SIZE)))

const DefaultSwappableVector{T} = SwappableVector{T,Threads.SpinLock}

struct GlobalRecord{T} <: AbstractVector{T}
    sinks::Vector{DefaultSwappableVector{T}}
    reserves::Vector{Vector{T}}
    tmp::Vector{Vector{T}}
    buffer::Vector{T}
end

function GlobalRecord{T}() where {T}
    sinks = [DefaultSwappableVector{T}() for _ in 1:Threads.nthreads()]
    reserves = [empty!(Vector{T}(undef, DEFAULT_SIZE)) for _ in 1:Threads.nthreads()]
    tmp = empty!(Vector{Vector{T}}(undef, Threads.nthreads()))
    buffer = T[]
    return GlobalRecord{T}(sinks, reserves, tmp, buffer)
end

function init!(gr::GlobalRecord{T}) where {T}
    new = GlobalRecord{T}()
    append!(empty!(gr.sinks), new.sinks)
    append!(empty!(gr.reserves), new.reserves)
    append!(empty!(gr.tmp), new.tmp)
    empty!(gr.buffer)
    return gr
end

function maintain!(gr::GlobalRecord)
    need_merge = false
    for i in eachindex(gr.sinks, gr.reserves)
        result = swapif!(!isempty, gr.sinks[i], gr.reserves[i])
        if result isa Some
            old = gr.reserves[i] = something(result)
            push!(gr.tmp, old)
            need_merge = true
        end
    end
    if need_merge
        foldl(append!, gr.tmp; init = gr.buffer)
        foreach(empty!, gr.tmp)
        empty!(gr.tmp)
    end
end

function Base.empty!(gr::GlobalRecord)
    empty!(gr.buffer)
    for i in eachindex(gr.sinks, gr.reserves)
        gr.reserves[i] = empty!(swap!(gr.sinks[i], gr.reserves[i]))
    end
    return gr
end

function Base.size(gr::GlobalRecord)
    maintain!(gr)
    return size(gr.buffer)
end

function Base.getindex(gr::GlobalRecord, i::Int)
    @boundscheck checkbounds(gr, i)
    return gr.buffer[i]
end

function Base.put!(gr::GlobalRecord{T}, v::T) where {T}
    locked(gr.sinks[Threads.nthreads()]) do sink
        push!(sink, v)
    end
end


text:
using ManifoldsBase
import ManifoldsBase: representation_size, manifold_dimension, inner
using LinearAlgebra
using Test

struct ComplexEuclidean{N} <: Manifold{ℂ} where {N} end
ComplexEuclidean(n::Int) = ComplexEuclidean{n}()
representation_size(::ComplexEuclidean{N}) where {N} = (N,)
manifold_dimension(::ComplexEuclidean{N}) where {N} = 2N
@inline inner(::ComplexEuclidean, x, v, w) = dot(v, w)

@testset "Complex Euclidean" begin
    M = ComplexEuclidean(3)

    x = complex.([1.0, 2.0, 3.0], [4.0, 5.0, 6.0])
    v1 = complex.([0.1, 0.2, 0.3], [0.4, 0.5, 0.6])
    v2 = complex.([0.3, 0.2, 0.1], [0.6, 0.5, 0.4])

    @test norm(M, x, v1) isa Real
    @test norm(M, x, v1) ≈ norm(v1)
    @test angle(M, x, v1, v2) isa Real
    @test angle(M, x, v1, v2) ≈ acos(real(dot(v1, v2)) / norm(v1) / norm(v2))

    vv1 = vcat(reim(v1)...)
    vv2 = vcat(reim(v2)...)
    @test angle(M, x, v1, v2) ≈ acos(dot(normalize(vv1), normalize(vv2)))
end


text:
################################################################################
#  Copyright 2021, Arpan Koirala, Tom Van Acker                                #
################################################################################
# StochasticPowerModels.jl                                                     #
# An extention package of PowerModels.jl for Stochastic (Optimal) Power Flow   #
# See http://github.com/timmyfaraday/StochasticPowerModels.jl                  #
################################################################################
# NOTE: dc lines are omitted from the current formulation                      #
################################################################################

""
function run_sopf_acr(data::Dict, model_constructor::Type, optimizer; aux::Bool=true, deg::Int=1, solution_processors=[sol_data_model!], kwargs...)
    @assert _IM.ismultinetwork(data) == false "The data supplied is multinetwork, it should be single-network"
    @assert model_constructor <: _PM.AbstractACRModel "This problem type only supports the ACRModel"

    sdata = build_stochastic_data(data, deg)
    if aux
        result = _PM.run_model(sdata, model_constructor, optimizer, build_sopf_acr_with_aux; multinetwork=true, solution_processors=solution_processors, kwargs...)
    else
        result = _PM.run_model(sdata, model_constructor, optimizer, build_sopf_acr_without_aux; multinetwork=true, solution_processors=solution_processors, kwargs...)
    end
    result["mop"] = sdata["mop"]
    return result
end

""
function run_sopf_acr(file::String, model_constructor, optimizer; aux::Bool=true, deg::Int=1, solution_processors=[sol_data_model!], kwargs...)
    data = _PM.parse_file(file)
    return run_sopf_acr(data, model_constructor, optimizer; aux=aux, deg=deg, solution_processors=solution_processors, kwargs...)
end

""
function build_sopf_acr_with_aux(pm::AbstractPowerModel)
    for (n, network) in _PM.nws(pm)
        variable_bus_voltage(pm, nw=n, aux=true)

        variable_gen_power(pm, nw=n, bounded=false)

        variable_branch_power(pm, nw=n, bounded=false)
        variable_branch_current(pm, nw=n, bounded=false, aux=true)
        for b in _PM.ids(pm, :branch, nw=n)
            expression_branch_voltage_drop_real(pm, b, nw=n)
            expression_branch_voltage_drop_img(pm, b, nw=n)
        end
    end

    for i in _PM.ids(pm, :bus,nw=1)
        constraint_bus_voltage_squared_cc_limit(pm, i, nw=1)
    end

    for g in _PM.ids(pm, :gen, nw=1)
        constraint_gen_power_cc_limit(pm, g, nw=1)
    end

    for b in _PM.ids(pm, :branch, nw=1)
       constraint_branch_series_current_squared_cc_limit(pm, b, nw=1)
    end

    for (n, network) in _PM.nws(pm)
        for i in _PM.ids(pm, :ref_buses, nw=n)
            constraint_bus_voltage_ref(pm, i, nw=n)
        end

        for i in _PM.ids(pm, :bus, nw=n)
            constraint_power_balance(pm, i, nw=n)
            constraint_gp_bus_voltage_squared(pm, i, nw=n)
        end

        for b in _PM.ids(pm, :branch, nw=n)
            constraint_gp_power_branch_to(pm, b, nw=n)
            constraint_gp_power_branch_from(pm, b, nw=n)

            #constraint_branch_voltage(pm, b, nw=n) #replaced by expression

            constraint_gp_current_squared(pm, b, nw=n)
        end
    end

    objective_min_expected_generation_cost(pm)
end

""
function build_sopf_acr_without_aux(pm::AbstractPowerModel)
    for (n, network) in _PM.nws(pm)
        variable_bus_voltage(pm, nw=n, aux=false)

        variable_gen_power(pm, nw=n, bounded=false)

        variable_branch_power(pm, nw=n, bounded=false)
        variable_branch_current(pm, nw=n, bounded=true, aux=false)
        for b in _PMs.ids(pm, :branch, nw=n)
            expression_branch_voltage_drop_real(pm, b, nw=n)
            expression_branch_voltage_drop_img(pm, b, nw=n)
        end
    end

    for i in _PM.ids(pm, :bus,nw=1)
        constraint_bus_voltage_cc_limit(pm, i, nw=1)
    end

    for g in _PM.ids(pm, :gen, nw=1)
        constraint_gen_power_cc_limit(pm, g, nw=1)
    end

    for b in _PM.ids(pm, :branch, nw=1)
       constraint_branch_series_current_cc_limit(pm, b, nw=1)
    end

    for (n, network) in _PM.nws(pm)

        for i in _PM.ids(pm, :ref_buses, nw=n)
            constraint_bus_voltage_ref(pm, i, nw=n)
        end

        for i in _PM.ids(pm, :bus, nw=n)
            constraint_power_balance(pm, i, nw=n)
        end

        for b in _PM.ids(pm, :branch, nw=n)
            constraint_gp_power_branch_to(pm, b, nw=n)
            constraint_gp_power_branch_from(pm, b, nw=n)

            #constraint_branch_voltage(pm, b, nw=n) #replaced by expression
        end
    end

    objective_min_expected_generation_cost(pm)
end


text:
function water_flood(core_props, fluids, rel_perms, core_flood)
    fw, dfw = fractional_flow_function(rel_perms, fluids)
    sw_init = core_flood.initial_water_saturation
    sw_shock = tangent_line_saturation(rel_perms, fluids, (sw_init, fw(sw_init)))
    # println("sw_shock = $sw_shock")
    t_D_BT = 1/dfw(sw_shock) # breakthrough (BT) time [#PV]
    # println("breakthrough time = $t_D_BT")

    # tracer
    sw_tracer = tangent_line_saturation(rel_perms, fluids, (0.0, 0.0))
    t_D_tracer = 1/dfw(sw_tracer)
    xt_tracer_shock = dfw(sw_tracer)
    # construct the recovery factor curve versus the # of PV
    # R = zeros(1)
    # pv_R = zeros(1)
    # at breakthrough
    # push!(R, (1-fw(sw_init))*t_D_BT/(1-sw_init)) # recovery at BT
    # push!(pv_R, t_D_BT) # BT time
    # after breakthrough
    pv_inj = max(core_flood.injected_pore_volume, 2.0) # at least inject 2 pv
    # sw_max = outlet_saturation(pv_inj, sw_shock, dfw, rel_perms)
    # println(sw_max)
    # sw_tmp = linspace(sw_shock, sw_max, 100)
    # t_D_tmp = 1./dfw.(sw_tmp)

    # s_av_tmp = sw_tmp-(fw.(sw_tmp)-1).*t_D_tmp
    # R_tmp = (s_av_tmp-sw_init)/(1-sw_init)
    # append!(R, R_tmp[2:end])
    # append!(pv_R, t_D_tmp[2:end])

    # saturation profile
    sor = rel_perms.sor
    sw_inj = core_flood.injected_water_saturation
    phi = core_props.porosity
    ut = core_flood.injection_velocity
    L = core_props.length
    pv_to_t = phi*L/ut
    s1 = collect(linspace(min(sw_inj, 1-sor-eps()), sw_shock, 500))
    xt_s1 = dfw.(s1)
    xt_shock = dfw(sw_shock)
    xt_prf=[xt_s1; xt_shock+eps(); 2*xt_shock]
    xt_prf[1]=0.0
    sw_prf=[s1; sw_init; sw_init]

    xt_tracer = [0.0, xt_tracer_shock, xt_tracer_shock+eps(), xt_prf[end]]
    c_tracer = [1.0, 1.0, 0.0, 0.0]

    # calculate the pressure drop
    k = core_props.permeability
    krw, kro, dkrw, dkro = rel_perm_functions(rel_perms)
    muo, muw = fluids.oil_viscosity, fluids.water_viscosity

    xt = copy(xt_prf)
    sw = copy(sw_prf)
    i=1
    while(true)
        if i==length(xt)
            break
        elseif xt[i]>=xt[i+1]
            deleteat!(xt, i)
            deleteat!(sw, i)
        else
            i+=1
        end
    end

    x = collect(linspace(0,L,200))
    sw_int = Spline1D(xt, sw, k=1, bc="nearest")
    # print(sw_int)
    # println(ut/phi.*xt)
    # println(sw)
    t_inj=pv_inj*pv_to_t
    t = collect(linspace(0.0,t_inj, 200)) # [s] time
    # println(x/t[2])

    # println(sw_int.(x/t[2]))
    p_inj = zeros(length(t))
    R_int = zeros(length(t))
    WC_int = zeros(length(t))
    p_inj[1]=L*ut./(k*(kro.(sw_init)/muo+krw.(sw_init)/muw))
    WC_int[1] = (krw(sw_init)/muw)/(kro(sw_init)/muo+krw(sw_init)/muw)
    for i in 2:length(t)
        # try
        xt_real = x/t[i]/(ut/phi)
        p_inj[i] = trapz(x, ut./(k*(kro.(sw_int(xt_real))/muo+krw.(sw_int(xt_real))/muw)))
        R_int[i] = (trapz(x, sw_int(xt_real))/L-sw_init)/(1-sw_init)
        WC_int[i] = (krw(sw_int(xt_real[end]))/muw)/(kro(sw_int(xt_real[end]))/muo+krw(sw_int(xt_real[end]))/muw)
        # catch
            # println(x/t[i])
            # error("interpolator at $i !!!")
        # end
    end

    return FracFlowResults([fw], [Line([sw_init, fw(sw_init)], [sw_shock, fw(sw_shock)])],
                            [t/pv_to_t R_int], [t R_int], [xt sw], [xt_tracer c_tracer],
                            [t/pv_to_t p_inj], [t p_inj], [t/pv_to_t WC_int], [t WC_int])
    # return pv_R, R, xt_prf, sw_prf # for the time being to test the code

end

max_stars_repo_name: simulkade/peteng

max_issues_repo_name: simulkade/peteng

max_forks_repo_name: simulkade/peteng


text:
# ## ghost
ghost(env, x_sol, timesteps=[1, 5, 7, 10, T])

set_floor!(env.vis, z=0.01)
set_camera!(env.vis, cam_pos=[0,-7,0], zoom=2)


text:
include("HelperFunctions.jl")
include("InitialConditions.jl")
include("EM.jl")

using DataFrames
using Dates
using Statistics
using CSV

export estimate_dfm
export export_dfm
export import_dfm
export predict_dfm

"""
    estimate a dynamic factor model.
    parameters:
        Y : DataFrame | Array{Number, 2}
            matrix of variables, size (n_obs, n_variables). Must include a column of type Dates.Date.
        blocks : DataFrame, size (n_variables, n_blocks). Note don't include date column in n_variables
            matrix of 1s or 0s for block loadings, i.e. included in block.
            ex for Y 2 x 3, 2 blocks:
                [1, 0]
                [1, 1]
                [0, 1]
        p : Int
            number of lags in transition equation
        max_iter : Int
            number of iterations for EM
        threshold : Float
            EM loop threshold
    returns: Dict
      :X_smooth => Kalman-smoothed data where missing values are replaced by their expectation
      :Z => smoothed states, rows give time, and columns are organized according to matrix C
      :C => measurement matrix, rows correspond to each series, and the columns are organized as
        - columns 1-20 give the factor loadings. For example, 1-5 give loadings for the first
        and are organized in reverse-chronological order (f^G_t, f^G_t-1, f^G_t-2, f^G_t-3, f^G_t-4)
        - Columns 6-10, 11-15, and 16-20 give loadings for the second, third, and fourth blocks respectively.
      :R => covariance for measurement matrix residuals
      :A => transition matrix, a square matrix that follows the same organization scheme as matrix C's columns.
        Identity matrices are used to account for matching terms on the left and righthand side.
        For example, we place an I4 matrix to account for matching (f_t-1; f_t-2; f_t-3; f_t-4) terms
      :Q => covariance for transition equation residuals
      :means =>  series mean
      :stds => series standard deviation
      :Z0 => initial value of state
      :V0 => initial value of covariance matrix
      :p => number of lags in transition equation
      :num_iter => number of iterations required for convergence
      :convergence => bool if algorithm converged successfully (given max_iter)
      :loglik => log likelihood of last iteration
      :LL => sequence of log likelihoods per iteration
"""
function estimate_dfm(Y; blocks, p, max_iter=5000, threshold=1e-5)
    R_mat = [2 -1 0 0 0; 3 0 -1 0 0; 2 0 0 -1 0; 1 0 0 0 -1] # R*λ = q; constraints on loadings of quarterly variables
    q = zeros(4)
    dates = Y[!, date_col_name(Y)]
    y_tmp = Y[!, Not(date_col_name(Y))]
    monthly_quarterly_array = gen_monthly_quarterly(dates, y_tmp)
    # put all quarterly variables to the end
    y_tmp = [y_tmp[!, .!monthly_quarterly_array]  y_tmp[!, monthly_quarterly_array]]
    monthly_quarterly_array = gen_monthly_quarterly(dates, y_tmp)
    nM = sum(.!monthly_quarterly_array)

    # calculate initial variables
    init_conds = initialize_conditions(y_tmp; dates=dates, p=p, blocks=blocks, R_mat=R_mat)
    A = init_conds[:A]; C = init_conds[:C]; Q = init_conds[:Q]; R = init_conds[:R]; Z0 = init_conds[:Z0]; V0 = init_conds[:V0]

    # initialize EM loop values
    #= EM model
    y = C*Z + e
    z = A*Z(-1) + v
    where y is n_obs x n_variables, Z is (pr) x n_variables, etc.
    =#

    # remove leading and ending rows where all missing, transpose
    y_est = y_tmp[[sum(.!ismissing.(Array(x))) > 0 for x in eachrow(y_tmp)], :] |> Array |> transpose
    n_obs = size(y_est)[2]
    prev_loglik = -1e6
    num_iter = 0
    LL = [-1e6]
    converged = false
    y = transpose(y_tmp |> Array) # y for the estimation is WITH missing data

    # EM loop
    gain = []
    while !converged & (num_iter <= max_iter)
        em_output = EM_step(y_est; A=A, C=C, Q=Q, R=R, Z0=Z0, V0=V0, p=p, blocks=blocks, R_mat=R_mat, q=q, nM=nM, monthly_quarterly_array=monthly_quarterly_array)
        A = em_output[:A_new]; C = em_output[:C_new]; Q = em_output[:Q_new]; R = em_output[:R_new]; Z0 = em_output[:Z0]; V0 = em_output[:V0]; global loglik = em_output[:loglik];

        em_conv = EM_convergence(loglik, prev_loglik, threshold)
        converged = em_conv[:converged] == 1

        if ((mod(num_iter, 10) == 0) & (num_iter > 0))
            println("Now running iteration number $num_iter out of a maximum of $max_iter")
            println("Loglik: $(round(loglik, digits=4)); % change: $(round(100 * (loglik - prev_loglik)/prev_loglik, digits=10))%")
        end

        # break loop if marginal gains
        push!(gain, (loglik - prev_loglik)/prev_loglik)
        if num_iter > 50
            if mean(gain[num_iter-20:num_iter]) < 1e-6
                break
            end
        end

        push!(LL, loglik)
        prev_loglik = loglik
        num_iter += 1
    end

    # Final run of the Kalman filter
    kf_output = kalman_filter(y_est, A, C, Q, R, Z0, V0)
    Zsmooth = transpose(kf_output[:Zsmooth])
    Xsmooth_std = Zsmooth[2:end,:] * transpose(C)
    Xsmooth = [std(skipmissing(row)) for row in eachrow(y_est)] |> x->
        reshape(x, 1, size(x)[1]) |> x->
        repeat(x, n_obs) |> x->
        x .* Xsmooth_std + ([mean(skipmissing(row)) for row in eachrow(y_est)] |> x-> reshape(x, 1, size(x)[1]) |> x-> repeat(x, n_obs))

    # final results
    nowcast = Dict(
        :Xsmooth_std => Xsmooth_std,
        :Xsmooth => Xsmooth,
        :Z => Zsmooth[2:end,:],
        :C => C,
        :R => R,
        :A => A,
        :Q => Q,
        :means => [mean(skipmissing(row)) for row in eachrow(y_est)] |> x-> reshape(x, 1, size(x)[1]),
        :stds => [std(skipmissing(row)) for row in eachrow(y_est)] |> x-> reshape(x, 1, size(x)[1]),
        :Z0 => Z0,
        :V0 => V0,
        :p => p,
        :num_iter => num_iter,
        :convergence => converged,
        :loglik => loglik,
        :LL => LL[2:end]
    )
    return nowcast
end

"outputs an output_dfm dict to disk (path = foldername)"
function export_dfm(;output_dfm::Dict, out_path::String)
    if !isdir(out_path)
        mkpath(out_path)
    end

    DataFrame(output_dfm[:Xsmooth]) |> x-> CSV.write("$out_path/X_smooth.csv", x) # :X_smooth
    DataFrame(output_dfm[:Xsmooth_std]) |> x-> CSV.write("$out_path/Xsmooth_std.csv", x) # :Xsmooth_std
    DataFrame(output_dfm[:Z]) |> x-> CSV.write("$out_path/Z.csv", x) # :Z
    DataFrame(output_dfm[:C]) |> x-> CSV.write("$out_path/C.csv", x) # :C
    DataFrame(output_dfm[:R]) |> x-> CSV.write("$out_path/R.csv", x) # :R
    DataFrame(output_dfm[:A]) |> x-> CSV.write("$out_path/A.csv", x) # :A
    DataFrame(output_dfm[:Q]) |> x-> CSV.write("$out_path/Q.csv", x) # :Q
    DataFrame(output_dfm[:means]) |> x-> CSV.write("$out_path/means.csv", x) # :means
    DataFrame(output_dfm[:stds]) |> x-> CSV.write("$out_path/stds.csv", x) # :stds
    DataFrame(Z0=output_dfm[:Z0]) |> x-> CSV.write("$out_path/Z0.csv", x) # :Z0
    DataFrame(output_dfm[:V0]) |> x-> CSV.write("$out_path/V0.csv", x) # :V0
    DataFrame(p=output_dfm[:p]) |> x-> CSV.write("$out_path/p.csv", x) # :p
    DataFrame(num_iter=output_dfm[:num_iter]) |> x-> CSV.write("$out_path/num_iter.csv", x) # :num_iter
    DataFrame(convergence=output_dfm[:convergence]) |> x-> CSV.write("$out_path/convergence.csv", x) # :convergence
    DataFrame(loglik=output_dfm[:loglik]) |> x-> CSV.write("$out_path/loglik.csv", x) # :loglik
    DataFrame(LL=output_dfm[:LL]) |> x-> CSV.write("$out_path/LL.csv", x) # :LL
end

"read an output_dfm directory to disk (folder name)"
function import_dfm(;path::String)::Dict
    tmp = Dict()
    tmp[:Xsmooth] = CSV.File("$path/X_smooth.csv") |> DataFrame!  |> Array # :X_smooth
    tmp[:Xsmooth_std] = CSV.File("$path/Xsmooth_std.csv") |> DataFrame!  |> Array # :Xsmooth_std
    tmp[:Z] = CSV.File("$path/Z.csv") |> DataFrame!  |> Array # :Z
    tmp[:C] = CSV.File("$path/C.csv") |> DataFrame!  |> Array # :C
    tmp[:R] = CSV.File("$path/R.csv") |> DataFrame!  |> Array # :R
    tmp[:A] = CSV.File("$path/A.csv") |> DataFrame!  |> Array # :A
    tmp[:Q] = CSV.File("$path/Q.csv") |> DataFrame!  |> Array # :Q
    tmp[:means] = CSV.File("$path/means.csv") |> DataFrame!  |> Array # :means
    tmp[:stds] = CSV.File("$path/stds.csv") |> DataFrame!  |> Array # :stds
    tmp[:Z0] = CSV.File("$path/Z0.csv") |> DataFrame! |> Array |> x-> reshape(x, size(x)[1]) # :Z0
    tmp[:V0] = CSV.File("$path/V0.csv") |> DataFrame! |> Array # :V0
    tmp[:p] = CSV.File("$path/p.csv") |> DataFrame!  |> x-> x[1,1] # :p
    tmp[:num_iter] = CSV.File("$path/num_iter.csv") |> DataFrame!  |> x-> x[1,1]  # :num_iter
    tmp[:convergence] = CSV.File("$path/convergence.csv") |> DataFrame!  |> x-> x[1,1] # :convergence
    tmp[:loglik] = CSV.File("$path/loglik.csv") |> DataFrame!  |> x-> x[1,1] # :loglik
    tmp[:LL] = CSV.File("$path/LL.csv") |> DataFrame! |> Array # :LL

    return tmp
end

"""
    get predictions from an already estimated DFM model for all series, all periods.
    parameters:
        data : DataFrame
            dataframe to generate predictions. Must include a date column, must contain the same columns the output_dfm was trained on.
        output_dfm : Dict
            output dictionary of the estimate_dfm function
        months_ahead : Int
            number of months ahead to forecast
        lag : Int
            number of lags for the kalman filter
    returns: DataFrame
      dataframe with all missing values filled + predictions
"""
function predict_dfm(data::DataFrame; output_dfm, months_ahead::Int, lag=0)
    orig_col_order = names(data)
    cp_data = copy(data)
    date_col = date_col_name(cp_data)
    dates = copy(cp_data[!, date_col])

    # put all quarterly variables to the end
    monthly_quarterly_array = gen_monthly_quarterly(dates, cp_data)
    cp_data = [cp_data[!, .!monthly_quarterly_array]  cp_data[!, monthly_quarterly_array]]

    Y = copy(cp_data[!, Not(date_col)])

    # add months
    last_date = cp_data[end, date_col]
    row = DataFrame(eachrow(Y)[1:months_ahead])
    row .= missing
    Y = vcat(Y, row)
    for i in 1:months_ahead
        push!(dates, last_date + Month(i))
    end

    predictions = kalman_filter_constparams(Y; output_dfm=output_dfm, lag=lag)[:X_smooth] |> DataFrame
    predictions[!, :date] = dates
    predictions =  predictions[!, [:date; Symbol.(names(predictions)[1:end-1])]]
    rename!(predictions, names(cp_data))
    predictions = predictions[!, orig_col_order]
    return predictions
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "SoapyRTLSDR"
version = v"0.3.3"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/pothosware/SoapyRTLSDR.git", "80c93fbe189def3ab68d47a6ad3f813f96d3cb99")
]

dependencies = [
    Dependency(PackageSpec(name="CompilerSupportLibraries_jll", uuid="e66e0078-7015-5450-92f7-15fbd957f2ae")),
    Dependency("librtlsdr_jll"; compat="0.6.0"),
    Dependency("soapysdr_jll"; compat="0.8.0")
]

# Bash recipe for building across all platforms
script = raw"""
cd SoapyRTLSDR
mkdir build
cd build
cmake -DCMAKE_INSTALL_PREFIX=${prefix} \
      -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
      -DCMAKE_BUILD_TYPE=Release \
      ..
make -j${nproc}
make install
if [[ "${target}" == *-apple-* ]]; then
    mv ${libdir}/SoapySDR/modules0.8/librtlsdrSupport.so  ${libdir}/SoapySDR/modules0.8/librtlsdrSupport.dylib
fi
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = filter!(p -> arch(p) != "armv6l", supported_platforms(;experimental=true))
platforms = expand_cxxstring_abis(platforms)

# The products that we will ensure are always built
products = Product[
    LibraryProduct("librtlsdrSupport", :librtlsdrSupport, ["lib/SoapySDR/modules0.8/"])
]

# Build the tarballs, and possibly a `build.jl` as well.
# gcc7 constraint from boost
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6")

max_stars_repo_name: staticfloat/BinaryBuilderRepos

max_issues_repo_name: staticfloat/BinaryBuilderRepos

max_forks_repo_name: staticfloat/BinaryBuilderRepos


text:
@everywhere function whoami()
    println(myid(), gethostname())
end

remotecall_fetch(whoami, 2)
remotecall_fetch(whoami, 4)


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "TopOpt_in_PETSc"
version = v"0.1"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/topopt/TopOpt_in_PETSc", "26eecbf3b1d0135956e0364d77c30e43e9bc3db2"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
# New makefiles added, the patches fix some weird include issues mostly.
# There is likely a better way to fix them, or upstream the fixes.
script = raw"""
cd TopOpt_in_PETSc
make -f ../Makefile libtopopt.${dlext} topopt${exeext}
install -Dvm 755 "topopt${exeext}" "${bindir}/topopt${exeext}"
install -Dvm 755 "libtopopt.${dlext}" "${libdir}/libtopopt.${dlext}"
install_license lesser.txt
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = expand_cxxstring_abis(supported_platforms(exclude=[Platform("i686", "windows")]))

# The products that we will ensure are always built
products = [
    LibraryProduct("libtopopt", :libtopopt),
    ExecutableProduct("topopt", :topopt)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("PETSc_jll"; compat="3.16.5"),
    Dependency(PackageSpec(name="CompilerSupportLibraries_jll", uuid="e66e0078-7015-5450-92f7-15fbd957f2ae"))
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6", preferred_gcc_version=v"9")

max_stars_repo_name: cmhamel/Yggdrasil

max_issues_repo_name: cmhamel/Yggdrasil

max_forks_repo_name: cmhamel/Yggdrasil


text:
# # Plots

## load packages
using Pkg.Artifacts
using PlotPlants
#------------------------------------------------------------------------------

# ## Plot shapes
# ### Ellipse
_fig,_axes = create_canvas("Ellipse");
_ax1 = _axes[1];
plot_ellipse!(_ax1, (0,0))
plot_ellipse!(_ax1, (12,0); width=20, height=8, angle=45, color="c",
              alpha=0.3);
plot_ellipse!(_ax1, (34,0); width=20, height=8, angle=15, edgecolor="r",
              facecolor="c", alpha=0.8);
_ax1.set_xlim(-11,45);
_ax1.set_ylim(-11,45);
_fig
#------------------------------------------------------------------------------

# ### Stoma
_fig,_axes = create_canvas("Stoma");
_ax1 = _axes[1];
plot_stoma!(_ax1, (0,0));
plot_stoma!(_ax1, (10,0); width=12, height=20, stoma=0.5, angle=15)
_ax1.set_xlim(-11,17);
_ax1.set_ylim(-11,17);
_fig
#------------------------------------------------------------------------------

# ## Preview global dataset
# ### lon×lat×3rd dimension
## load Artifacts.toml from PlotPlants
## leaf area index dataset, variable name: LAI
artifact_toml = joinpath(@__DIR__, "../../../Artifacts.toml");
data_hash     = artifact_hash("LAI_4X_1M_V1", artifact_toml);
data_path     = artifact_path(data_hash) * "/LAI_4X_1M_V1.nc";
ensure_artifact_installed("LAI_4X_1M_V1", artifact_toml);

_fig,_axes = create_canvas("Preview 3D dataset");
_ax1 = _axes[1];
preview_dataset!(_ax1, data_path, "LAI");
_fig
#------------------------------------------------------------------------------

# ### lon×lat dimension
## canopy height dataset, variable name: Band1
artifact_toml = joinpath(@__DIR__, "../../../Artifacts.toml");
data_hash     = artifact_hash("CH_20X_1Y_V1", artifact_toml);
data_path     = artifact_path(data_hash) * "/CH_20X_1Y_V1.nc";
ensure_artifact_installed("CH_20X_1Y_V1", artifact_toml);

_fig,_axes = create_canvas("Preview 2D dataset");
_ax1 = _axes[1];
preview_dataset!(_ax1, data_path, "Band1");
_fig
#------------------------------------------------------------------------------

# ### Non-supported case
artifact_toml = joinpath(@__DIR__, "../../../Artifacts.toml");
data_hash     = artifact_hash("CH_20X_1Y_V1", artifact_toml);
data_path     = artifact_path(data_hash) * "/CH_20X_1Y_V1.nc";
ensure_artifact_installed("CH_20X_1Y_V1", artifact_toml);

@info "Expecting a warning and thus a blank figure here!";
_fig,_axes = create_canvas("Expecting a warning");
_ax1 = _axes[1];
preview_dataset!(_ax1, data_path, 1);
_fig
#------------------------------------------------------------------------------

# ## Plot regressions
# ### Linear regressions
xs = collect(1:100);
y1 = xs .* 0.01 .+ rand(100);
y2 = xs .* 0.01 .+ rand(100) .+ 2;

_fig,_axes = create_canvas("Linear regression");
_ax1 = _axes[1];
_ax1.plot(xs, y1, "ro", alpha=0.2);
plot_line_regress!(_ax1, xs, y1);
_ax1.plot(xs, y2, "co", alpha=0.2);
plot_line_regress!(_ax1, xs, y2; color="c");
_fig
#------------------------------------------------------------------------------

# ### Set intercept = 0
xs = collect(1:100);
ys = xs .* 0.01 .+ rand(100) .+ 2;

_fig,_axes = create_canvas("Force intercept = 0");
_ax1 = _axes[1];
_ax1.plot(xs, ys, "ro", alpha=0.2);
plot_line_regress!(_ax1, xs, ys; intercept=false);
_fig
#------------------------------------------------------------------------------

# ### Plot confidence interval
xs = collect(1:100);
ys = xs .* 0.01 .+ rand(100) .+ 2;

_fig,_axes = create_canvas("Plot confidence interval");
_ax1 = _axes[1];
_ax1.plot(xs, ys, "ko", alpha=0.2);
plot_line_regress!(_ax1, xs, ys; interval=true);
_fig
#------------------------------------------------------------------------------

# ## Plot density
# ### Density scatter plot
xs = [rand(100); rand(100).* 0.5; rand(100).* 0.2; rand(100).* 0.1];
ys = [rand(100); rand(100).* 0.5; rand(100).* 0.2; rand(100).* 0.1];

_fig,_axes = create_canvas("Plot confidence interval");
_ax1 = _axes[1];
plot_density!(_ax1, xs, ys; cmap="viridis", markersize=10);
_fig
#------------------------------------------------------------------------------

# ### Set maximal density limit
xs = [rand(100); rand(100).* 0.5; rand(100).* 0.2; rand(100).* 0.1];
ys = [rand(100); rand(100).* 0.5; rand(100).* 0.2; rand(100).* 0.1];

_fig,_axes = create_canvas("Set maximal density limit");
_ax1 = _axes[1];
plot_density!(_ax1, xs, ys; cmap="viridis", markersize=10, dmax=1);
_fig
#------------------------------------------------------------------------------


text:
@testset "Square root Unscented Kalman filter" begin

    @testset "Covariance" begin
        weight_params = ScaledSetWeightingParameters(0.5, 2, 1)
        x = randn(5)
        PL_prior = randn(5, 5)
        P_prior = PL_prior * PL_prior'
        χ = KalmanFilters.calc_sigma_points(x, P_prior, weight_params)
        F = randn(5, 5)
        f(x) = F * x
        𝓨 = KalmanFilters.transform(f, χ)
        QL = randn(5, 5)
        Q = QL * QL'
        Q_chol = cholesky(Q)
        P = @inferred KalmanFilters.cov(𝓨, Q_chol)
        @test P.L * P.U ≈ KalmanFilters.cov(𝓨, Q)
    end

    @testset "Posterior covariance" begin
        weight_params = ScaledSetWeightingParameters(0.5, 2, 1)
        x = randn(5)
        PL = randn(5, 5)
        P = PL * PL'
        P_chol = cholesky(P)
        χ = KalmanFilters.calc_sigma_points(x, P, weight_params)
        F = randn(3, 5)
        h(x) = F * x
        𝓨 = KalmanFilters.transform(h, χ)
        y_est = KalmanFilters.mean(𝓨)
        unbiased_𝓨 = KalmanFilters.substract_mean(𝓨, y_est)
        RL = randn(3, 3)
        R = RL * RL'
        R_chol = cholesky(R)
        S = KalmanFilters.cov(unbiased_𝓨, R)
        S_chol = KalmanFilters.cov(unbiased_𝓨, R_chol)
        Pᵪᵧ = KalmanFilters.cov(χ, unbiased_𝓨)
        K = Pᵪᵧ / S_chol
        K_temp, P_post = @inferred KalmanFilters.calc_kalman_gain_and_posterior_covariance(P_chol, Pᵪᵧ, S_chol, [])
        @test K_temp ≈ K
        @test P_post.L * P_post.U ≈ KalmanFilters.calc_posterior_covariance(P, Pᵪᵧ, K, [])
    end

    @testset "Time update" begin

        x = randn(6)
        A = randn(6,6)
        P = A'A
        P_chol = cholesky(P)
        B = randn(6,6)
        Q = B'B
        Q_chol = cholesky(Q)
        F = randn(6,6)
        f(x) = F * x

        tu = @inferred time_update(x, P, F, Q)
        tu_chol = @inferred time_update(x, P_chol, f, Q_chol)
        @test @inferred(get_covariance(tu_chol)) ≈ get_covariance(tu)
        @test @inferred(get_state(tu_chol)) ≈ get_state(tu)

        f!(y, x) = mul!(y, F, x)
        tu_inter = @inferred SRUKFTUIntermediate(6)
        tu_chol_inplace = @inferred time_update!(tu_inter, x, P_chol, f!, Q_chol)
        @test @inferred(get_covariance(tu_chol_inplace)) ≈ get_covariance(tu)
        @test @inferred(get_state(tu_chol_inplace)) ≈ get_state(tu)
    end

    @testset "Measurement update" begin

        x = randn(6)
        A = randn(6,6)
        P = A'A
        P_chol = cholesky(P)
        B = randn(3,3)
        R = B'B
        y = randn(3)
        R_chol = cholesky(R)
        H = randn(3,6)
        h(x) = H * x

        mu = @inferred measurement_update(x, P, y, H, R)
        mu_chol = @inferred measurement_update(x, P_chol, y, h, R_chol)
        @test @inferred(get_covariance(mu_chol)) ≈ get_covariance(mu)
        @test @inferred(get_state(mu_chol)) ≈ get_state(mu)

        h!(y, x) = mul!(y, H, x)
        mu_inter = @inferred SRUKFMUIntermediate(6, 3)
        mu_chol_inplace = @inferred measurement_update!(mu_inter, x, P_chol, y, h!, R_chol)
        @test @inferred(get_covariance(mu_chol_inplace)) ≈ get_covariance(mu)
        @test @inferred(get_state(mu_chol_inplace)) ≈ get_state(mu)
    end
end


text:
struct GeoConst{T,N}
    adv::Advection{T,N}
    #    rho::T
    #    g::T
    #    buoyancy_freq_n::T
    odg_b::T
    #    hv_order::Int
    #    hv_val::T
    coefrsqk::NTuple{N,Array{Complex{T},N}}
    pfftbig::Union{Missing,PrepareFftBig{T}}
    function GeoConst(
        adv::Advection{T,N};
        #        rho::T = 1e3,
        #        g::T = 9.81,
        buoyancy_freq_n::T = 3*2*T(2pi)/(24*3600)*sin(T(pi)/4),
        odg_b::T = T(1e-3),
        #        hv_order::Int = 8,
        #        hv_val::T = T(0),
    ) where {T,N}
        @assert N == 2 "the number of dimension must be 2"
        sz = sizeall(adv)
        coefrsqk = ntuple(x -> zeros(Complex{T}, sz), N)
        pfftbig = if T != Float64
            PrepareFftBig(sizeall(adv), T; numdims = N, dims = ntuple(x -> x, N))
        else
            missing
        end
        veck = vec_k_fft.(adv.t_mesh)
        for ind in CartesianIndices(sz)
            tp = ntuple(i -> veck[i][ind.I[i]], N)
            k = sqrt(sum(tp .^ 2))
            v = k == 0 ? 0 : 1 / k
            for d = 1:N
                coefrsqk[d][ind] = v .* tp[d%2+1]
            end
        end
        coefrsqk[1] .*= im # /buoyancy_freq_n
        coefrsqk[2] .*= -im # /buoyancy_freq_n

        return new{T,N}(
            adv,
            #            rho,
            #            g,
            #            buoyancy_freq_n,
            odg_b,
            #            hv_order,
            #            hv_val,
            coefrsqk,
            pfftbig,
        )
    end
end

struct GeoVar{T,N} <: AbstractExtDataAdv
    gc::GeoConst{T,N}
end

getgeovar(adv::Advection{T,N}; kwargs...) where {T,N} =
    GeoVar{T,N}(GeoConst(adv; kwargs...))

function initdata!(geoc::GeoVar{T,N}, advd::AdvectionData{T,N}) where {T,N}

    mesh_x = advd.adv.t_mesh[1]
    mesh_y = advd.adv.t_mesh[2]
    lx = stop(mesh_x) - start(mesh_x)
    ly = stop(mesh_y) - start(mesh_y)
    x, y = points.(advd.adv.t_mesh)
    ee = 4
    σ = lx / 15 # Length scale close to the Rossby radius

    ## Spatial buoyancy field
    anticyclone(cx, cy) =
        (exp.(-ee .* (x .- cx) .^ 2 ./ 2σ^2) .* transpose(exp.(-(y .- cy) .^ 2 ./ 2σ^2)))

    # Warm anticyclones
    advd.data = anticyclone(lx / 4, ly / 4)
    advd.data .+= anticyclone(3lx / 4, ly / 4)

    # Cold cyclones
    advd.data .-= anticyclone(lx / 4, 3ly / 4)
    advd.data .-= anticyclone(3lx / 4, 3ly / 4)

    # Specify the amplitude of the buoyancy
    advd.data .*= geoc.gc.odg_b

end

function initcoef!(geoc::GeoVar{T,N}, advd::AdvectionData{T,N}) where {T,N}
    pfft = geoc.gc.pfftbig
    dt = getcur_t(advd)
    buf = fftgenall(pfft, advd.data)

    if ismissing(advd.bufcur)
        advd.bufcur = zeros(OpTuple{N,T}, sizeall(advd.adv))
    end

    advd.bufcur .=
        dt * OpTuple.(zip(ntuple(x -> real(ifftgenall(pfft, geoc.gc.coefrsqk[x] .* buf)), N)...))
end
@inline function getalpha(
    ::GeoVar{T,N},
    self::AdvectionData{T},
    indext::CartesianIndex, indbuf::CartesianIndex
) where {T,N,Nsp,Nv}
    st = getst(self)
    ind = CartesianIndex((indbuf.I...,indext.I...)[st.invp])
    return ntuple( x -> self.bufcur[ind][st.invp[x]], st.ndims)
end