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

using Markdown
using InteractiveUtils

# ╔═╡ d948dc6e-2de1-11eb-19e7-cb3bb66353b6
md"# Immer"

# ╔═╡ 3a99e22d-42d6-4b2d-9381-022b41b0e852
md"### wrappath"

# ╔═╡ aad7ab32-eecf-4aad-883d-1c802cad6c0c
md"### =="

# ╔═╡ f658a72d-871d-49b3-9b73-7efedafbd7a6
md"### convert(::Type{Dict}, ::JSONPatch)"

# ╔═╡ daf9ec12-2de1-11eb-3a8d-59d9c2753134
md"## Diff"

# ╔═╡ db116c0a-2de1-11eb-2a56-872af797c547
abstract type JSONPatch end

# ╔═╡ 73631aea-5e93-4da2-a32d-649029660d4e
Patches = Array{JSONPatch, 1}

# ╔═╡ 0fd3e910-abcc-4421-9d0b-5cfb90034338
NoChanges = Patches()

# ╔═╡ bd0d46bb-3e58-4522-bae0-83eb799196c4
PatchPath = AbstractVector{Any}

# ╔═╡ 831d84a6-1c71-4e68-8c7c-27d9093a82c4
function wrappath(path::PatchPath, patches::Vector{JSONPatch})
	map(patches) do patch
		wrappath(path, patch)
	end
end

# ╔═╡ db2d8a3e-2de1-11eb-02b8-9ffbfaeff61c
struct AddPatch <: JSONPatch
	path::PatchPath
	value::Any
end

# ╔═╡ 2ad11c73-4691-4283-8f98-3d2a87926b99
function wrappath(path, patch::AddPatch)
	AddPatch([path..., patch.path...], patch.value)
end

# ╔═╡ f649f67c-aab0-4d35-a799-f398e5f3ecc4
function Base.:(==)(a::AddPatch, b::AddPatch)
	a.value == b.value && a.path == b.path
end

# ╔═╡ 230bafe2-aaa7-48f0-9fd1-b53956281684
function Base.convert(::Type{Dict}, patch::AddPatch)
	Dict(:op => :add, :path => patch.path, :value => patch.value)
end

# ╔═╡ 6d67f8a5-0e0c-4b6e-a267-96b34d580946
add_patch = AddPatch(["counter"], 10)

# ╔═╡ ffe9b3d9-8e35-4a31-bab2-8787a4140594
struct RemovePatch <: JSONPatch
	path::PatchPath
end

# ╔═╡ 5513ea3b-9498-426c-98cb-7dc23d32f72e
function wrappath(path, patch::RemovePatch)
	RemovePatch([path..., patch.path...])
end

# ╔═╡ 63087738-d70c-46f5-b072-21cd8953df35
function Base.:(==)(a::RemovePatch, b::RemovePatch)
	a.path == b.path
end

# ╔═╡ b48e2c08-a94a-4247-877d-949d92dde626
function Base.convert(::Type{Dict}, patch::RemovePatch)
	Dict(:op => :remove, :path => patch.path)
end

# ╔═╡ 56b28842-4a67-44d7-95e7-55d457a44fb1
remove_patch = RemovePatch(["counter"])

# ╔═╡ 894de8a7-2757-4d7a-a2be-1069fa872911
struct ReplacePatch <: JSONPatch
	path::PatchPath
	value::Any
end

# ╔═╡ 0c2d6da1-cad3-4c9f-93e9-922457083945
function wrappath(path, patch::ReplacePatch)
	ReplacePatch([path..., patch.path...], patch.value)
end

# ╔═╡ aa81974a-7254-45e0-9bfe-840c4793147f
function Base.:(==)(a::ReplacePatch, b::ReplacePatch)
	a.path == b.path && a.value == b.value
end

# ╔═╡ fafcb8b8-cde9-4f99-9bab-8128025953a4
function Base.convert(::Type{<:Dict}, patch::ReplacePatch)
	Dict(:op => :replace, :path => patch.path, :value => patch.value)
end

# ╔═╡ 921a130e-b028-4f91-b077-3bd79dcb6c6d
function Base.convert(::Type{JSONPatch}, patch_dict::Dict)
	if patch_dict[:op] == :add
		AddPatch(patch_dict[:path], patch_dict[:value])
	elseif patch_dict[:op] == :remove
		RemovePatch(patch_dict[:path])
	elseif patch_dict[:op] == :replace
		ReplacePatch(patch_dict[:path], patch_dict[:value])
	else
		throw("Unknown operation :$(patch_dict[:op]) in Dict to JSONPatch conversion")
	end
end

# ╔═╡ 07eeb122-6706-4544-a007-1c8d6581eec8
Base.convert(Dict, AddPatch([:x, :y], 10))

# ╔═╡ c59b30b9-f702-41f1-bb2e-1736c8cd5ede
Base.convert(Dict, RemovePatch([:x, :y]))

# ╔═╡ 7feeee3a-3aec-47ce-b8d7-74a0d9b0b381
Base.convert(Dict, ReplacePatch([:x, :y], 10))

# ╔═╡ f10e31c0-1d2c-4727-aba5-dd676a10041b
replace_patch = ReplacePatch(["counter"], 10)

# ╔═╡ 9a364714-edb1-4bca-9387-a8bbacccd10d
struct CopyPatch <: JSONPatch
	path::PatchPath
	from::PatchPath
end

# ╔═╡ 84c87031-7733-4d1f-aa90-f8ab71506251
function wrappath(path, patch::CopyPatch)
	CopyPatch([path..., patch.path...], patch.from)
end

# ╔═╡ 31188a03-76ba-40cf-a333-4d339ce37711
function Base.:(==)(a::CopyPatch, b::CopyPatch)
	a.path == b.path && a.from == b.from
end

# ╔═╡ 9321d3be-cb91-4406-9dc7-e5c38f7d377c
struct MovePatch <: JSONPatch
	path::PatchPath
	from::PatchPath
end

# ╔═╡ 8f265a33-3a2d-4508-9477-ca62e8ce3c12
function wrappath(path, patch::MovePatch)
	MovePatch([path..., patch.path...], patch.from)
end

# ╔═╡ 7524a9e8-1a6d-4851-b50e-19415f25a84b
function Base.:(==)(a::MovePatch, b::MovePatch)
	a.path == b.path && a.from == b.from
end

# ╔═╡ 5e360fcd-9943-4a17-9672-f1fded2f7e3a
function diff(old::T, new::T)::Patches where T
	if old == new
		NoChanges
	else
		[ReplacePatch([], new)]
	end
end

# ╔═╡ 9cbaaec2-709c-4769-886c-ec92b12c18bc
struct Deep{T} value::T end

# ╔═╡ db75df12-2de1-11eb-0726-d1995cebd382
function diff(old::Deep{T}, new::Deep{T})::Patches where T
	changes = JSONPatch[]
	for property in propertynames(old.value)
		for change in diff(getproperty(old.value, property), getproperty(new.value, property))
			push!(changes, wrappath([property], change))
		end
	end
	changes

	# changes = []
	# for property in fieldnames(T)
	# 	for change in diff(getfield(old.value, property), getfield(new.value, property))
	# 		push!(changes, wrappath([property], change))
	# 	end
	# end
	# changes
end

# ╔═╡ dbc7f97a-2de1-11eb-362f-055a734d1a9e
function diff(o1::AbstractDict, o2::AbstractDict)::Patches
	changes = JSONPatch[]
	for key in keys(o1) ∪ keys(o2)
		for change in diff(get(o1, key, nothing), get(o2, key, nothing))
			push!(changes, wrappath([key], change))
		end
	end
	changes
end

# ╔═╡ 67ade214-2de3-11eb-291d-135a397d629b
function diff(o1, o2)::Patches
	[ReplacePatch([], o2)]
end

# ╔═╡ b8c58aa4-c24d-48a3-b2a8-7c01d50a3349
function diff(o1::Nothing, o2)::Patches
	[AddPatch([], o2)]
end

# ╔═╡ 5ab390f9-3b0c-4978-9e21-2aaa61db2ce4
function diff(o1, o2::Nothing)::Patches
	[RemovePatch([])]
end

# ╔═╡ 09f53db0-21ae-490b-86b5-414eba403d57
function diff(o1::Nothing, o2::Nothing)::Patches
	NoChanges
end

# ╔═╡ c7de406d-ccfe-41cf-8388-6bd2d7c42d64
md"### Struct example"

# ╔═╡ b9cc11ae-394b-44b9-bfbe-541d7720ead0
struct Cell
	id
	code
	folded
end

# ╔═╡ c3c675be-9178-4176-afe0-30501786b72c
deep_diff(old::Cell, new::Cell) = diff(Deep(old), Deep(new))

# ╔═╡ 02585c72-1d92-4526-98c2-1ca07aad87a3
function direct_diff(old::Cell, new::Cell)
	changes = []
	if old.id ≠ new.id
		push!(changes, ReplacePatch([:id], new.id))
	end
	if old.code ≠ new.code
		push!(changes, ReplacePatch([:code], new.code))
	end
	if old.folded ≠ new.folded
		push!(changes, ReplacePatch([:folded], new.folded))
	end
end

# ╔═╡ 2d084dd1-240d-4443-a8a2-82ae6e0b8900
cell1 = Cell(1,2,3)

# ╔═╡ 3e05200f-071a-4ebe-b685-ff980f07cde7
cell2 = Cell(1,2,4)

# ╔═╡ 200516da-8cfb-42fe-a6b9-cb4730168923
celldict1 = Dict(:x => 1, :y => 2, :z => 3)

# ╔═╡ 76326e6c-b95a-4b2d-a78c-e283e5fadbe2
celldict2 = Dict(:x => 1, :y => 2, :z => 4)

# ╔═╡ dd312598-2de1-11eb-144c-f92ed6484f5d
md"## Update"

# ╔═╡ d2af2a4b-8982-4e43-9fd7-0ecfdfb70511
STRICT = false

# ╔═╡ 640663fc-06ba-491e-bd85-299514237651
begin
	function force_convert_key(::Dict{T,<:Any}, value::T) where T
		value
	end
	function force_convert_key(::Dict{T,<:Any}, value::Any) where T
		T(value)
	end
end

# ╔═╡ 48a45941-2489-4666-b4e5-88d3f82e5145
function getpath(value, path)
	if length(path) == 0
		return value
	end

	current = path[firstindex(path)]
	rest = path[firstindex(path)+1:end]
	if value isa AbstractDict
		key = force_convert_key(value, current)
		getpath(getindex(value, key), rest)
	else
		getpath(getproperty(value, Symbol(current)), rest)
	end
end

# ╔═╡ 752b2da3-ff24-4758-8843-186368069888
function update!(value, patches::Array{JSONPatch})
	for patch in patches
		update!(value, patch)
	end
	return value
end

# ╔═╡ 3e285076-1d97-4728-87cf-f71b22569e57
md"### update! AddPatch"

# ╔═╡ a11e4082-4ff4-4c1b-9c74-c8fa7dcceaa6
md"*Should throw in strict mode:*"

# ╔═╡ dd87ca7e-2de1-11eb-2ec3-d5721c32f192
function update!(value, patch::AddPatch)
	if length(patch.path) == 0
		throw("Impossible")
	else
		last = patch.path[end]
		rest = patch.path[firstindex(patch.path):end-1]
		subvalue = getpath(value, rest)
		if subvalue isa AbstractDict
			key = force_convert_key(subvalue, last)
			if STRICT
				@assert get(subvalue, key, nothing) === nothing
			end
			subvalue[key] = patch.value
		else
			key = Symbol(last)
			if STRICT
				@assert getproperty(subvalue, key) === nothing
			end
			setproperty!(subvalue, key, patch.value)
		end
	end
	return value
end

# ╔═╡ be6b6fc4-e12a-4cef-81d8-d5115fda50b7
md"### update! ReplacePatch"

# ╔═╡ f1dde1bd-3fa4-48b7-91ed-b2f98680fcc1
md"*Should throw in strict mode:*"

# ╔═╡ 6509d62e-77b6-499c-8dab-4a608e44720a
function update!(value, patch::ReplacePatch)
	if length(patch.path) == 0
		throw("Impossible")
	else
		last = patch.path[end]
		rest = patch.path[firstindex(patch.path):end-1]
		subvalue = getpath(value, rest)
		if subvalue isa AbstractDict
			key = force_convert_key(subvalue, last)
			if STRICT
				@assert get(subvalue, key, nothing) !== nothing
			end
			subvalue[key] = patch.value
		else
			key = Symbol(last)
			if STRICT
				@assert getproperty(subvalue, key) !== nothing
			end
			setproperty!(subvalue, key, patch.value)
		end
	end
	return value
end

# ╔═╡ f3ef354b-b480-4b48-8358-46dbf37e1d95
md"### update! RemovePatch"

# ╔═╡ df41caa7-f0fc-4b0d-ab3d-ebdab4804040
md"*Should throw in strict mode:*"

# ╔═╡ ddaf5b66-2de1-11eb-3348-b905b94a984b
function update!(value, patch::RemovePatch)
	if length(patch.path) == 0
		throw("Impossible")
	else
		last = patch.path[end]
		rest = patch.path[firstindex(patch.path):end-1]
		subvalue = getpath(value, rest)
		if subvalue isa AbstractDict
			key = force_convert_key(subvalue, last)
			if STRICT
				@assert get(subvalue, key, nothing) !== nothing
			end
			delete!(subvalue, key)
		else
			key = Symbol(last)
			if STRICT
				@assert getproperty(subvalue, key) !== nothing
			end
			setproperty!(subvalue, key, nothing)
		end
	end
	return value
end

# ╔═╡ e55d1cea-2de1-11eb-0d0e-c95009eedc34
md"## Testing"

# ╔═╡ e598832a-2de1-11eb-3831-371aa2e54828
abstract type TestResult end

# ╔═╡ e5b46afe-2de1-11eb-0de5-6d571c0fbbcf
const Code = Any

# ╔═╡ e5dbaf38-2de1-11eb-13a9-a994ac40bf9f
struct Pass <: TestResult
	expr::Code
end

# ╔═╡ e616c708-2de1-11eb-2e66-f972030a7ec5
abstract type Fail <: TestResult end

# ╔═╡ e6501fda-2de1-11eb-33ba-4bb34dc13d00
struct Wrong <: Fail
	expr::Code
	result
end

# ╔═╡ e66c8454-2de1-11eb-1d79-499e6873d0d2
struct Error <: Fail
	expr::Code
	error
end

# ╔═╡ e699ae9a-2de1-11eb-3ff0-c31222ac399e
function Base.show(io::IO, mime::MIME"text/html", value::Pass)
	show(io, mime, HTML("""

			$(value.expr)

	"""))
end

# ╔═╡ e6c17fae-2de1-11eb-1397-1b1cdfcc387c
function Base.show(io::IO, mime::MIME"text/html", value::Wrong)
	show(io, mime, HTML("""

			$(value.expr)

	"""))
end

# ╔═╡ e705bd90-2de1-11eb-3759-3d59a90e6e44
function Base.show(io::IO, mime::MIME"text/html", value::Error)
	show(io, mime, HTML("""

				$(value.expr)

Error: $(sprint(showerror, value.error))

	"""))
end

# ╔═╡ e748600a-2de1-11eb-24be-d5f0ecab8fa4
macro test(expr)
	quote
		expr_raw = $(expr |> QuoteNode)
		try
			result = $(esc(expr))
			if result == true
				Pass(expr_raw)
			else
				Wrong(expr_raw, result)
			end
		catch e
			Error(expr_raw, e)
		end

		# Base.@locals()
	end
end

# ╔═╡ 7b8ab89b-bf56-4ddf-b220-b4881f4a2050
@test Base.convert(JSONPatch, convert(Dict, add_patch)) == add_patch

# ╔═╡ 48ccd28a-060d-4214-9a39-f4c4e506d1aa
@test Base.convert(JSONPatch, convert(Dict, remove_patch)) == remove_patch

# ╔═╡ 34d86e02-dd34-4691-bb78-3023568a5d16
@test Base.convert(JSONPatch, convert(Dict, replace_patch)) == replace_patch

# ╔═╡ c3e4738f-4568-4910-a211-6a46a9d447ee
@test update!(Dict(:y => "x"), AddPatch([:x], "-")) == Dict(:y => "x", :x => "-")

# ╔═╡ 0f094932-10e5-40f9-a3fc-db27a85b4999
@test update!(Dict(:x => "x"), AddPatch([:x], "-")) == Dict(:x => "-")

# ╔═╡ a560fdca-ee12-469c-bda5-62d7203235b8
@test update!(Dict(:x => "x"), ReplacePatch([:x], "-")) == Dict(:x => "-")

# ╔═╡ 01e3417e-334e-4a8d-b086-4bddc42737b3
@test update!(Dict(:y => "x"), ReplacePatch([:x], "-")) == Dict(:x => "-", :y => "x")

# ╔═╡ 96a80a23-7c56-4c41-b489-15bc1c4e3700
@test update!(Dict(:x => "x"), RemovePatch([:x])) == Dict()

# ╔═╡ fac65755-2a2a-4a3c-b5a8-fc4f6d256754
@test update!(Dict(:y => "x"), RemovePatch([:x])) == Dict(:y => "x")

# ╔═╡ e8d0c98a-2de1-11eb-37b9-e1df3f5cfa25
md"## DisplayOnly"

# ╔═╡ e907d862-2de1-11eb-11a9-4b3ac37cb0f3
function displayonly(m::Module)
	if isdefined(m, :PlutoForceDisplay)
		return m.PlutoForceDisplay
	else
		isdefined(m, :PlutoRunner) && parentmodule(m) == Main
	end
end

# ╔═╡ e924a0be-2de1-11eb-2170-71d56e117af2
"""
	@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 displayonly(ex) displayonly(__module__) ? esc(ex) : nothing end

# ╔═╡ 664cd334-91c7-40dd-a2bf-0da720307cfc
@displayonly notebook1 = Dict(
	:x => 1,
	:y => 2,
)

# ╔═╡ e7e8d076-2de1-11eb-0214-8160bb81370a
@test notebook1 == deepcopy(notebook1)

# ╔═╡ b7fa5625-6178-4da8-a889-cd4f014f43ba
@displayonly notebook2 = Dict(
	:y => 4,
	:z => 5
)

# ╔═╡ dbdd1df0-2de1-11eb-152f-8d1af1ad02fe
@displayonly notebook1_to_notebook2 = diff(notebook1, notebook2)

# ╔═╡ e78b7408-2de1-11eb-2f1a-3f0783049c7d
@test to_jsonpatch(flatten(notebook1_to_notebook2)) == test_jsonpatch

# ╔═╡ e7c85c1c-2de1-11eb-1a2a-65f8f21e4a05
@test update(notebook1, from_jsonpatch(to_jsonpatch(notebook1_to_notebook2))) == notebook2

# ╔═╡ 595fdfd4-3960-4fbd-956c-509c4cf03473
@displayonly @test update!(deepcopy(notebook1), notebook1_to_notebook2) == notebook2

# ╔═╡ e9d2eba8-2de1-11eb-16bf-bd2a16537a97
@displayonly x = 2

# ╔═╡ ea45104e-2de1-11eb-3248-5dd833d350e4
@displayonly @test 1 + 1 == x

# ╔═╡ ea6650bc-2de1-11eb-3016-4542c5c333a5
@displayonly @test 1 + 1 + 1 == x

# ╔═╡ ea934d9c-2de1-11eb-3f1d-3b60465decde
@displayonly @test throw("Oh my god") == x

# ╔═╡ ee70e282-36d5-4772-8585-f50b9a67ca54
md"## Track"

# ╔═╡ 3ab3b8a0-19fe-424d-8857-604b1e805a26
xxxxxz = 4

# ╔═╡ a3e8fe70-cbf5-4758-a0f2-d329d138728c
function prettytime(time_ns::Number)
    suffices = ["ns", "μs", "ms", "s"]

	current_amount = time_ns
	suffix = ""
	for current_suffix in suffices
    	if current_amount >= 1000.0
        	current_amount = current_amount / 1000.0
		else
			suffix = current_suffix
			break
		end
	end

    # const roundedtime = time_ns.toFixed(time_ns >= 100.0 ? 0 : 1)
	roundedtime = if current_amount >= 100.0
		round(current_amount; digits=0)
	else
		round(current_amount; digits=1)
	end
    return "$(roundedtime) $(suffix)"
end

# ╔═╡ 0e1c6442-9040-49d9-b754-173583db7ba2
begin
Base.@kwdef struct Tracked
	expr
	value
	time
	bytes
	times_ran=1
	which=nothing
	code_info=nothing
end
function Base.show(io::IO, mime::MIME"text/html", value::Tracked)
	times_ran = if value.times_ran === 1
		""
	else
		""" ($(value.times_ran)×)"""
	end
	# method = sprint(show, MIME("text/plain"), value.which)
	code_info = if value.code_info ≠ nothing
		codelength = length(value.code_info.first.code)
		"$(codelength) frames in @code_typed"
	else
		""
	end
	color = if value.time > 1
		"red"
	elseif value.time > 0.001
		"orange"
	elseif value.time > 0.0001
		"blue"
	else
		"green"
	end

	show(io, mime, HTML("""

				$(value.expr)

					$(prettytime(value.time * 1e9))
					$(times_ran)

$(code_info)

	"""))
end
	Tracked
end

# ╔═╡ 7618aef7-1884-4e32-992d-0fd988e1ab20
macro track(expr)
	times_ran_expr = :(1)
	expr_to_show = expr
	if expr.head == :for
		@assert expr.args[1].head == :(=)
		times_ran_expr = expr.args[1].args[2]
		expr_to_show = expr.args[2].args[2]
	end

	Tracked # reference so that baby Pluto understands

	quote
		local times_ran = length($(esc(times_ran_expr)))
		local value, time, bytes = @timed $(esc(expr))

		local method = nothing
		local code_info = nothing
		try
			# Uhhh
			method = @which $(expr_to_show)
			code_info = @code_typed $(expr_to_show)
		catch nothing end
		Tracked(
			expr=$(QuoteNode(expr_to_show)),
			value=value,
			time=time,
			bytes=bytes,
			times_ran=times_ran,
			which=method,
			code_info=code_info
		)
	end
end

# ╔═╡ fa959806-3264-4dd5-9f94-ba369697689b
@displayonly @track for _ in 1:1000 direct_diff(cell2, cell1) end

# ╔═╡ a9088341-647c-4fe1-ab85-d7da049513ae
@displayonly @track for _ in 1:1000 diff(Deep(cell1), Deep(cell2)) end

# ╔═╡ 95ff676d-73c8-44cb-ac35-af94418737e9
@displayonly @track for _ in 1:100 diff(celldict1, celldict2) end

# ╔═╡ 1020b11e-1364-42bc-a91a-2e96f6228c42
@track for _ in 1:100 Base.isless(1, 2) end

# ╔═╡ f958a11c-27bd-47c5-b239-564fd082f014
@track for _ in 1:100 Base.isless(xxxxxz, 2) end

# ╔═╡ 71b95468-c934-43dd-833e-33683015597f
@track for _ in 1:100 Base.isless(1, 2) end

# ╔═╡ 1a26eed8-670c-43bf-9726-2db84b1afdab
@track sleep(0.1)

# ╔═╡ Cell order:
# ╟─d948dc6e-2de1-11eb-19e7-cb3bb66353b6
# ╟─3a99e22d-42d6-4b2d-9381-022b41b0e852
# ╟─831d84a6-1c71-4e68-8c7c-27d9093a82c4
# ╟─2ad11c73-4691-4283-8f98-3d2a87926b99
# ╟─5513ea3b-9498-426c-98cb-7dc23d32f72e
# ╟─0c2d6da1-cad3-4c9f-93e9-922457083945
# ╟─84c87031-7733-4d1f-aa90-f8ab71506251
# ╟─8f265a33-3a2d-4508-9477-ca62e8ce3c12
# ╟─aad7ab32-eecf-4aad-883d-1c802cad6c0c
# ╠═f649f67c-aab0-4d35-a799-f398e5f3ecc4
# ╠═63087738-d70c-46f5-b072-21cd8953df35
# ╠═aa81974a-7254-45e0-9bfe-840c4793147f
# ╠═31188a03-76ba-40cf-a333-4d339ce37711
# ╠═7524a9e8-1a6d-4851-b50e-19415f25a84b
# ╟─f658a72d-871d-49b3-9b73-7efedafbd7a6
# ╟─07eeb122-6706-4544-a007-1c8d6581eec8
# ╠═230bafe2-aaa7-48f0-9fd1-b53956281684
# ╟─c59b30b9-f702-41f1-bb2e-1736c8cd5ede
# ╠═b48e2c08-a94a-4247-877d-949d92dde626
# ╟─7feeee3a-3aec-47ce-b8d7-74a0d9b0b381
# ╠═fafcb8b8-cde9-4f99-9bab-8128025953a4
# ╠═921a130e-b028-4f91-b077-3bd79dcb6c6d
# ╟─6d67f8a5-0e0c-4b6e-a267-96b34d580946
# ╟─7b8ab89b-bf56-4ddf-b220-b4881f4a2050
# ╟─56b28842-4a67-44d7-95e7-55d457a44fb1
# ╟─48ccd28a-060d-4214-9a39-f4c4e506d1aa
# ╟─f10e31c0-1d2c-4727-aba5-dd676a10041b
# ╟─34d86e02-dd34-4691-bb78-3023568a5d16
# ╟─daf9ec12-2de1-11eb-3a8d-59d9c2753134
# ╟─0fd3e910-abcc-4421-9d0b-5cfb90034338
# ╠═db116c0a-2de1-11eb-2a56-872af797c547
# ╟─73631aea-5e93-4da2-a32d-649029660d4e
# ╠═bd0d46bb-3e58-4522-bae0-83eb799196c4
# ╠═db2d8a3e-2de1-11eb-02b8-9ffbfaeff61c
# ╠═ffe9b3d9-8e35-4a31-bab2-8787a4140594
# ╠═894de8a7-2757-4d7a-a2be-1069fa872911
# ╠═9a364714-edb1-4bca-9387-a8bbacccd10d
# ╠═9321d3be-cb91-4406-9dc7-e5c38f7d377c
# ╠═5e360fcd-9943-4a17-9672-f1fded2f7e3a
# ╠═9cbaaec2-709c-4769-886c-ec92b12c18bc
# ╠═db75df12-2de1-11eb-0726-d1995cebd382
# ╠═dbc7f97a-2de1-11eb-362f-055a734d1a9e
# ╠═67ade214-2de3-11eb-291d-135a397d629b
# ╠═b8c58aa4-c24d-48a3-b2a8-7c01d50a3349
# ╠═5ab390f9-3b0c-4978-9e21-2aaa61db2ce4
# ╠═09f53db0-21ae-490b-86b5-414eba403d57
# ╟─c7de406d-ccfe-41cf-8388-6bd2d7c42d64
# ╠═b9cc11ae-394b-44b9-bfbe-541d7720ead0
# ╠═c3c675be-9178-4176-afe0-30501786b72c
# ╠═02585c72-1d92-4526-98c2-1ca07aad87a3
# ╟─2d084dd1-240d-4443-a8a2-82ae6e0b8900
# ╟─3e05200f-071a-4ebe-b685-ff980f07cde7
# ╠═fa959806-3264-4dd5-9f94-ba369697689b
# ╟─a9088341-647c-4fe1-ab85-d7da049513ae
# ╟─200516da-8cfb-42fe-a6b9-cb4730168923
# ╟─76326e6c-b95a-4b2d-a78c-e283e5fadbe2
# ╟─95ff676d-73c8-44cb-ac35-af94418737e9
# ╠═664cd334-91c7-40dd-a2bf-0da720307cfc
# ╠═b7fa5625-6178-4da8-a889-cd4f014f43ba
# ╠═dbdd1df0-2de1-11eb-152f-8d1af1ad02fe
# ╠═595fdfd4-3960-4fbd-956c-509c4cf03473
# ╟─dd312598-2de1-11eb-144c-f92ed6484f5d
# ╠═d2af2a4b-8982-4e43-9fd7-0ecfdfb70511
# ╠═640663fc-06ba-491e-bd85-299514237651
# ╠═48a45941-2489-4666-b4e5-88d3f82e5145
# ╟─752b2da3-ff24-4758-8843-186368069888
# ╟─3e285076-1d97-4728-87cf-f71b22569e57
# ╟─c3e4738f-4568-4910-a211-6a46a9d447ee
# ╟─a11e4082-4ff4-4c1b-9c74-c8fa7dcceaa6
# ╟─0f094932-10e5-40f9-a3fc-db27a85b4999
# ╟─dd87ca7e-2de1-11eb-2ec3-d5721c32f192
# ╟─be6b6fc4-e12a-4cef-81d8-d5115fda50b7
# ╟─a560fdca-ee12-469c-bda5-62d7203235b8
# ╟─f1dde1bd-3fa4-48b7-91ed-b2f98680fcc1
# ╟─01e3417e-334e-4a8d-b086-4bddc42737b3
# ╟─6509d62e-77b6-499c-8dab-4a608e44720a
# ╟─f3ef354b-b480-4b48-8358-46dbf37e1d95
# ╟─96a80a23-7c56-4c41-b489-15bc1c4e3700
# ╟─df41caa7-f0fc-4b0d-ab3d-ebdab4804040
# ╟─fac65755-2a2a-4a3c-b5a8-fc4f6d256754
# ╟─ddaf5b66-2de1-11eb-3348-b905b94a984b
# ╟─e55d1cea-2de1-11eb-0d0e-c95009eedc34
# ╟─e598832a-2de1-11eb-3831-371aa2e54828
# ╟─e5b46afe-2de1-11eb-0de5-6d571c0fbbcf
# ╟─e5dbaf38-2de1-11eb-13a9-a994ac40bf9f
# ╟─e616c708-2de1-11eb-2e66-f972030a7ec5
# ╟─e6501fda-2de1-11eb-33ba-4bb34dc13d00
# ╟─e66c8454-2de1-11eb-1d79-499e6873d0d2
# ╟─e699ae9a-2de1-11eb-3ff0-c31222ac399e
# ╟─e6c17fae-2de1-11eb-1397-1b1cdfcc387c
# ╠═e705bd90-2de1-11eb-3759-3d59a90e6e44
# ╠═e748600a-2de1-11eb-24be-d5f0ecab8fa4
# ╟─e78b7408-2de1-11eb-2f1a-3f0783049c7d
# ╟─e7c85c1c-2de1-11eb-1a2a-65f8f21e4a05
# ╟─e7e8d076-2de1-11eb-0214-8160bb81370a
# ╟─e8d0c98a-2de1-11eb-37b9-e1df3f5cfa25
# ╟─e907d862-2de1-11eb-11a9-4b3ac37cb0f3
# ╟─e924a0be-2de1-11eb-2170-71d56e117af2
# ╟─e9d2eba8-2de1-11eb-16bf-bd2a16537a97
# ╟─ea45104e-2de1-11eb-3248-5dd833d350e4
# ╟─ea6650bc-2de1-11eb-3016-4542c5c333a5
# ╟─ea934d9c-2de1-11eb-3f1d-3b60465decde
# ╟─ee70e282-36d5-4772-8585-f50b9a67ca54
# ╟─0e1c6442-9040-49d9-b754-173583db7ba2
# ╟─7618aef7-1884-4e32-992d-0fd988e1ab20
# ╠═1020b11e-1364-42bc-a91a-2e96f6228c42
# ╠═3ab3b8a0-19fe-424d-8857-604b1e805a26
# ╠═f958a11c-27bd-47c5-b239-564fd082f014
# ╠═71b95468-c934-43dd-833e-33683015597f
# ╠═1a26eed8-670c-43bf-9726-2db84b1afdab
# ╟─a3e8fe70-cbf5-4758-a0f2-d329d138728c


text:
using MixedModelsSim
using DataFrames, Tables
using Test

@testset "simdat_crossed" begin
    subj_btwn = Dict(:age => ["O", "Y"],
                     :pet => ["cat", "dog"])
    item_btwn = Dict(:cond => ["A", "B"])
    both_win = Dict(:time => ["morning", "evening"])

    for n in [4, 8]
        subj_n = n
        item_n = n
        # since we don't inspect the dv here, we don't worry about specifying the RNG
        # this has the advantage of testing the default/implicit-RNG method
        dat = simdat_crossed(subj_n, item_n,
                             subj_btwn = subj_btwn,
                             item_btwn = item_btwn,
                             both_win = both_win,
                             subj_prefix="SS",
                             item_prefix="LL")

        @test Tables.isrowtable(dat)

        @test first(dat).subj[1:2] == "SS"
        @test first(dat).item[1:2] == "LL"

        @test length(dat) == 2 * subj_n * item_n
        @test length(unique(getproperty.(dat, :subj))) == subj_n
        @test length(unique(getproperty.(dat, :item))) == item_n
        @test Set(Tables.columnnames(first(dat))) == Set([:subj, :pet, :age, :item, :cond, :time, :dv])
    end

    @test_throws ArgumentError simdat_crossed(3, 1, subj_btwn = subj_btwn)
    @test_throws ArgumentError simdat_crossed(1, 3, item_btwn = item_btwn)

    # the nothingness
    dat = simdat_crossed(1, 1)
    @test Tables.isrowtable(dat)
    @test Set(Tables.columnnames(first(dat))) == Set([:subj, :item, :dv])
end

ext: jl


text:
include("DM1_1.jl")
include("DM2_1.jl")

function ReactiveGRASP(matrix,cost, n, m, ite, coupe,temps)

    #initialisation
    (p,nb_iteration,z_cumul,zBest,zWorst) = intialiser(matrix,cost, n, m, ite, coupe)

    evol_p=Float64[]
    t=time()

    while (time()-t <= temps)
        cpt = 1
        append!(evol_p,p)
        while (cpt <= ite)
            for i in 1:length(p)
                (SOL,z) = GRASP(cost, matrix, n, m, p[i])
                z_cumul[i] += z
            end

            cpt = cpt+1

            prob=rand(Float64)

            alpha_choisit = choix_alpha(p,prob)
            (SOL,z) = GRASP(cost, matrix, n, m, p[alpha_choisit])

            nb_iteration[alpha_choisit] += 1
            z_cumul[alpha_choisit] += z

            if z > zBest
                zBest = z
            end

            if z < zWorst
                zWorst = z
            end
        end

        recalcul_p!(p,z_cumul,zBest,zWorst,nb_iteration,evol_p)
        nb_iteration = ones(Int64, length(p))
        z_cumul = zeros(Int64, length(p))

    end
#    println("^^^^^^^^ Fin ReactiveGRASP ^^^^^^^^")
    println("zBest: ", zBest, "    zWorst: ", zWorst)

    return(evol_p)
end

#=======================================#

function intialiser(matrix,cost, n, m, ite, coupe)

    p = Float64[]
    for i in 1:coupe
        push!(p,i/coupe)
    end

    nb_iteration = ones(Int64, length(p))

    z_cumul = zeros(Int64, length(p))

    (SOL,z) =Glouton(cost, matrix, n, m)

    zBest = z

    zWorst = z

    return(p,nb_iteration,z_cumul,zBest,zWorst)
end

#=======================================#

function recalcul_p!(p,z_cumul,zBest,zWorst,nb_iteration,evol_p)
    q=zeros(Float64,length(p))

    somme_q=0

    for i in 1:length(p)

        moyenne=(z_cumul[i]/nb_iteration[i])

        q[i]=(moyenne-zWorst)/(zBest-zWorst)

        somme_q+=q[i]

    end

    for i in 1:length(q)
        p[i]=q[i]/somme_q
    end

    for i in 2:length(p)-1
        p[i]+=p[i-1]
    end
    p[length(p)]=1

end

#=======================================#

function choix_alpha(p,prob)
    for i in 1:length(p)
        if prob < p[i]
            return(i)
        elseif i==length(p)
            return(i)
        end
    end
end

max_stars_repo_path: Archives/before the merge/DM2_2.jl

max_issues_repo_path: Archives/before the merge/DM2_2.jl

max_forks_repo_path: Archives/before the merge/DM2_2.jl


text:
using Random, StatsBase
Random.seed!(1)

A = Set(['a','e','i','o','u'])
B = Set(['x','y','z'])
omega = 'a':'z'

N = 10^6

println("mcEst1 \t \tmcEst2")
for _ in 1:5
    mcEst1 = sum([in(sample(omega),A) || in(sample(omega),B) for _ in 1:N])/N
    mcEst2 = sum([in(sample(omega),union(A,B)) for _ in 1:N])/N
    println(mcEst1,"\t",mcEst2)
end

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
let
    @test AutomotiveDrivingModels.cyclic_shift_left!([1,2,3,4], 1, 4) == [2,3,4,1]
    @test AutomotiveDrivingModels.cyclic_shift_left!([1,2,3,4], 2, 4) == [3,4,1,2]
    @test AutomotiveDrivingModels.cyclic_shift_left!([1,2,3,4,5,6,7], 1, 4) == [2,3,4,1,5,6,7]
    @test AutomotiveDrivingModels.cyclic_shift_left!([1,2,3,4,5,6,7], 2, 4) == [3,4,1,2,5,6,7]

    roadway = get_test_roadway()
    trajdata = get_test_trajdata(roadway)
    scene = Scene()

    col = get_first_collision(get!(scene, trajdata, 1))
    @test col.is_colliding
    @test col.A == 1
    @test col.B == 2

    @test get_first_collision(get!(scene, trajdata, 2), CPAMemory()).is_colliding == true
end


text:
"""
    FluxVariabilitySummary

Stores summary information about the result of a flux variability analysis.
"""
struct FluxVariabilitySummary
    biomass_fluxes::Dict{String,Vector{Maybe{Float64}}}
    exchange_fluxes::Dict{String,Vector{Maybe{Float64}}}
end

"""
    FluxVariabilitySummary()

A default empty constructor for [`FluxVariabilitySummary`](@ref).
"""
function FluxVariabilitySummary()
    FluxVariabilitySummary(
        Dict{String,Vector{Maybe{Float64}}}(),
        Dict{String,Vector{Maybe{Float64}}}(),
    )
end

"""
    flux_variability_summary(flux_result::Tuple{Dict{String, Dict{String, Float64}}, Dict{String, Dict{String, Float64}}};
        exclude_exchanges = false,
        exchange_prefixes = _constants.exchange_prefixes,
        biomass_strings = _constants.biomass_strings,
        exclude_biomass = false,
        )::FluxVariabilitySummary

Summarize a dictionary of flux dictionaries obtained eg. from
[`flux_variability_analysis_dict`](@ref). The simplified summary representation
is useful for pretty-printing and easily showing the most important results.

Internally this function uses [`looks_like_biomass_reaction`](@ref) and
[`looks_like_exchange_reaction`](@ref). The corresponding keyword arguments are
passed to these functions. Use this if your model has an uncommon naming of
reactions.

# Example
```
julia> sol = flux_variability_analysis_dict(model, Gurobi.Optimizer; bounds = objective_bounds(0.99))
julia> flux_res = flux_variability_summary(sol)
Biomass                     Lower bound   Upper bound
  BIOMASS_Ecoli_core_w_GAM: 0.8652        0.8652
Exchange
  EX_h2o_e:                 28.34         28.34
  EX_co2_e:                 22.0377       22.0377
  EX_o2_e:                  -22.1815      -22.1815
  EX_h_e:                   17.3556       17.3556
  EX_glc__D_e:              -10.0         -10.0
  EX_nh4_e:                 -4.8448       -4.8448
  EX_pi_e:                  -3.2149       -3.2149
  EX_for_e:                 0.0           0.0
  ...                       ...           ...
```
"""
function flux_variability_summary(
    flux_result::Tuple{Dict{String,Dict{String,Float64}},Dict{String,Dict{String,Float64}}};
    exclude_exchanges = false,
    exchange_prefixes = _constants.exchange_prefixes,
    biomass_strings = _constants.biomass_strings,
    exclude_biomass = false,
)
    isnothing(flux_result) && return FluxVariabilitySummary()

    rxn_ids = keys(flux_result[1])
    ex_rxns = filter(
        x -> looks_like_exchange_reaction(
            x,
            exclude_biomass = exclude_biomass,
            biomass_strings = biomass_strings,
            exchange_prefixes = exchange_prefixes,
        ),
        rxn_ids,
    )
    bmasses = filter(
        x -> looks_like_biomass_reaction(
            x;
            exclude_exchanges = exclude_exchanges,
            exchange_prefixes = exchange_prefixes,
            biomass_strings = biomass_strings,
        ),
        rxn_ids,
    )

    biomass_fluxes = Dict{String,Vector{Maybe{Float64}}}()
    for rxn_id in bmasses
        lb = isnothing(flux_result[1][rxn_id]) ? nothing : flux_result[1][rxn_id][rxn_id]
        ub = isnothing(flux_result[2][rxn_id]) ? nothing : flux_result[2][rxn_id][rxn_id]
        biomass_fluxes[rxn_id] = [lb, ub]
    end

    ex_rxn_fluxes = Dict{String,Vector{Maybe{Float64}}}()
    for rxn_id in ex_rxns
        lb = isnothing(flux_result[1][rxn_id]) ? nothing : flux_result[1][rxn_id][rxn_id]
        ub = isnothing(flux_result[2][rxn_id]) ? nothing : flux_result[2][rxn_id][rxn_id]
        ex_rxn_fluxes[rxn_id] = [lb, ub]
    end

    return FluxVariabilitySummary(biomass_fluxes, ex_rxn_fluxes)
end


text:
using Oceananigans.Architectures: architecture, device_event
using Oceananigans.AbstractOperations: KernelFunctionOperation
using Oceananigans.BuoyancyModels: ∂z_b
using Oceananigans.Operators: ℑzᵃᵃᶜ

struct ConvectiveAdjustmentVerticalDiffusivity{TD, CK, CN, BK, BN} <: AbstractTurbulenceClosure{TD}
    convective_κz :: CK
    convective_νz :: CN
    background_κz :: BK
    background_νz :: BN

    function ConvectiveAdjustmentVerticalDiffusivity{TD}(convective_κz::CK,
                                                         convective_νz::CN,
                                                         background_κz::BK,
                                                         background_νz::BN) where {TD, CK, CN, BK, BN}

        return new{TD, CK, CN, BK, BN}(convective_κz, convective_νz, background_κz, background_νz)
    end
end

"""
    ConvectiveAdjustmentVerticalDiffusivity([FT=Float64;]
                                            convective_κz = 0,
                                            convective_νz = 0,
                                            background_κz = 0,
                                            background_νz = 0,
                                            time_discretization = VerticallyImplicitTimeDiscretization())

The one positional argument determines the floating point type of the free parameters
of `ConvectiveAdjustmentVerticalDiffusivity`. The default is `Float64`.

Keyword arguments
=================

* `convective_κz` : Vertical tracer diffusivity in regions with negative (unstable) buoyancy gradients. Either
                    a single number, function, array, field, or tuple of diffusivities for each tracer.

* `background_κz` : Vertical tracer diffusivity in regions with zero or positive (stable) buoyancy gradients.

* `convective_νz` : Vertical viscosity in regions with negative (unstable) buoyancy gradients. Either
                    a number, function, array, or field.

* `background_κz` : Vertical viscosity in regions with zero or positive (stable) buoyancy gradients.

* `time_discretization` : Either `ExplicitTimeDiscretization` or `VerticallyImplicitTimeDiscretization`.
"""
function ConvectiveAdjustmentVerticalDiffusivity(FT = Float64;
                                                 convective_κz = zero(FT),
                                                 convective_νz = zero(FT),
                                                 background_κz = zero(FT),
                                                 background_νz = zero(FT),
                                                 time_discretization::TD = VerticallyImplicitTimeDiscretization()) where TD

    return ConvectiveAdjustmentVerticalDiffusivity{TD}(convective_κz, convective_νz,
                                                       background_κz, background_νz)
end

const CAVD = ConvectiveAdjustmentVerticalDiffusivity

#####
##### Diffusivity field utilities
#####

with_tracers(tracers, closure::ConvectiveAdjustmentVerticalDiffusivity{TD}) where TD =
    ConvectiveAdjustmentVerticalDiffusivity{TD}(closure.convective_κz,
                                                closure.convective_νz,
                                                closure.background_κz,
                                                closure.background_νz)

# Note: computing diffusivities at cell centers for now.
function DiffusivityFields(arch, grid, tracer_names, bcs, closure::CAVD)
    #data = new_data(Bool, arch, grid, (Center, Center, Center))
    κ = Field(Center, Center, Center, arch, grid)
    ν = Field(Center, Center, Center, arch, grid)
    return (; κ, ν)
end

function calculate_diffusivities!(diffusivities, closure::CAVD, model)

    arch = model.architecture
    grid = model.grid
    tracers = model.tracers
    buoyancy = model.buoyancy

    event = launch!(arch, grid, :xyz,
                    #compute_stability!, diffusivities, grid, closure, tracers, buoyancy,
                    compute_convective_adjustment_diffusivities!, diffusivities, grid, closure, tracers, buoyancy,
                    dependencies = device_event(arch))

    wait(device(arch), event)

    return nothing
end

@inline is_unstableᶜᶜᶠ(i, j, k, grid, tracers, buoyancy) = ∂z_b(i, j, k, grid, buoyancy, tracers) < 0

@kernel function compute_convective_adjustment_diffusivities!(diffusivities, grid, closure, tracers, buoyancy)
    i, j, k, = @index(Global, NTuple)

    unstable_cell = is_unstableᶜᶜᶠ(i, j, k+1, grid, tracers, buoyancy) |
                    is_unstableᶜᶜᶠ(i, j, k,   grid, tracers, buoyancy)

    @inbounds diffusivities.κ[i, j, k] = ifelse(unstable_cell,
                                                closure.convective_κz,
                                                closure.background_κz)

    @inbounds diffusivities.ν[i, j, k] = ifelse(unstable_cell,
                                                closure.convective_νz,
                                                closure.background_νz)
end

#=
@kernel function compute_stability!(diffusivities, grid, closure, tracers, buoyancy)
    i, j, k, = @index(Global, NTuple)
    @inbounds diffusivities.unstable_buoyancy_gradient[i, j, k] = is_unstableᶜᶜᶠ(i, j, k, grid, tracers, buoyancy)
end
=#

#####
##### Fluxes
#####

const VITD = VerticallyImplicitTimeDiscretization
const ATD = AbstractTimeDiscretization

@inline viscous_flux_ux(i, j, k, grid, ::ATD, closure::CAVD, args...) = zero(eltype(grid))
@inline viscous_flux_uy(i, j, k, grid, ::ATD, closure::CAVD, args...) = zero(eltype(grid))
@inline viscous_flux_vx(i, j, k, grid, ::ATD, closure::CAVD, args...) = zero(eltype(grid))
@inline viscous_flux_vy(i, j, k, grid, ::ATD, closure::CAVD, args...) = zero(eltype(grid))
@inline viscous_flux_wx(i, j, k, grid, ::ATD, closure::CAVD, args...) = zero(eltype(grid))
@inline viscous_flux_wy(i, j, k, grid, ::ATD, closure::CAVD, args...) = zero(eltype(grid))

@inline diffusive_flux_x(i, j, k, grid, ::ATD, closure::CAVD, args...) = zero(eltype(grid))
@inline diffusive_flux_y(i, j, k, grid, ::ATD, closure::CAVD, args...) = zero(eltype(grid))

@inline viscous_flux_ux(i, j, k, grid, closure::CAVD, args...) = zero(eltype(grid))
@inline viscous_flux_uy(i, j, k, grid, closure::CAVD, args...) = zero(eltype(grid))
@inline viscous_flux_vx(i, j, k, grid, closure::CAVD, args...) = zero(eltype(grid))
@inline viscous_flux_vy(i, j, k, grid, closure::CAVD, args...) = zero(eltype(grid))
@inline viscous_flux_wx(i, j, k, grid, closure::CAVD, args...) = zero(eltype(grid))
@inline viscous_flux_wy(i, j, k, grid, closure::CAVD, args...) = zero(eltype(grid))

@inline diffusive_flux_x(i, j, k, grid, closure::CAVD, args...) = zero(eltype(grid))
@inline diffusive_flux_y(i, j, k, grid, closure::CAVD, args...) = zero(eltype(grid))

#####
##### Diffusivity
#####

const etd = ExplicitTimeDiscretization()

@inline z_boundary_adj(k, grid) = k == 1 | k == grid.Nz+1

@inline z_diffusivity(closure::CAVD, c_idx, diffusivities, args...) = diffusivities.κ

@inline function diffusive_flux_z(i, j, k, grid, closure::CAVD, c, tracer_index, clock, diffusivities, args...)
    κ = κᶜᶜᶠ(i, j, k, grid, clock, diffusivities.κ)
    return - κ * ∂zᵃᵃᶠ(i, j, k, grid, c)
end

@inline function diffusive_flux_z(i, j, k, grid::VerticallyBoundedGrid, ::VITD, closure::CAVD, args...)
    explicit_flux_z = diffusive_flux_z(i, j, k, grid, etd, closure, args...)
    return ifelse(z_boundary_adj(k, grid), explicit_flux_z, zero(eltype(grid)))
end

#####
##### Viscosity
#####

@inline z_viscosity(closure::CAVD, diffusivities, args...) = diffusivities.ν

@inline function viscous_flux_uz(i, j, k, grid::VerticallyBoundedGrid, ::VITD, closure::CAVD, args...)
    explicit_flux_z = viscous_flux_uz(i, j, k, grid, etd, closure, args...)
    return ifelse(z_boundary_adj(k, grid), explicit_flux_z, zero(eltype(grid)))
end

@inline function viscous_flux_vz(i, j, k, grid::VerticallyBoundedGrid, ::VITD, closure::CAVD, args...)
    explicit_flux_z = viscous_flux_vz(i, j, k, grid, etd, closure, args...)
    return ifelse(z_boundary_adj(k, grid), explicit_flux_z, zero(eltype(grid)))
end

@inline function viscous_flux_wz(i, j, k, grid::VerticallyBoundedGrid, ::VITD, closure::CAVD, args...)
    explicit_flux_z = viscous_flux_wz(i, j, k, grid, etd, closure, args...)
    return ifelse(z_boundary_adj(k, grid), explicit_flux_z, zero(eltype(grid)))
end

@inline function viscous_flux_uz(i, j, k, grid, closure::CAVD, clock, velocities, diffusivities, args...)
    ν = νᶠᶜᶠ(i, j, k, grid, clock, diffusivities.ν)
    return - ν * ∂zᵃᵃᶠ(i, j, k, grid, velocities.u)
end

@inline function viscous_flux_vz(i, j, k, grid, closure::CAVD, clock, velocities, diffusivities, args...)
    ν = νᶜᶠᶠ(i, j, k, grid, clock, diffusivities.ν)
    return - ν * ∂zᵃᵃᶠ(i, j, k, grid, velocities.v)
end

@inline function viscous_flux_wz(i, j, k, grid, closure::CAVD, clock, velocities, diffusivities, args...)
    ν = νᶜᶜᶜ(i, j, k, grid, clock, diffusivities.ν)
    return - ν * ∂zᵃᵃᶜ(i, j, k, grid, velocities.w)
end


text:
function arr_sumsqr{T <: Number}(x::Array{T})
r = zero(T)
for i = 1:length(x)
r = r + x[i] ^ 2
end
return r
end

max_stars_repo_path: Module 2/Chapter03/arr_sumsqr.jl

max_issues_repo_path: Module 2/Chapter03/arr_sumsqr.jl

max_forks_repo_path: Module 2/Chapter03/arr_sumsqr.jl


text:
using StringDistances, Unicode, Test, Random

@testset "Pairwise" begin

	TestStrings1 = ["", "abc", "bc", "kitten"]
	TestStrings2 = ["mew", "ab"]
	TestStrings1missing = ["", "abc", "bc", missing]
	TestStrings2missing = ["mew", missing]

	for d in [Jaro(), Levenshtein(), OptimalStringAlignement(), RatcliffObershelp(),
				QGram(2), Cosine(2), Jaccard(2), SorensenDice(2), Overlap(2)]

		R = pairwise(d, TestStrings1)

		@test size(R) == (4, 4)

		# No distance on the diagonal, since comparing strings to themselves
		@test R[1, 1] == 0.0
		@test R[2, 2] == 0.0
		@test R[3, 3] == 0.0
		@test R[4, 4] == 0.0

		# Since the distance might be NaN:
		equalorNaN(x, y) = (x == y) || (isnan(x) && isnan(y))

		# First row is comparing "" to the other strings, so:
		@test equalorNaN(R[1, 2], evaluate(d, "", "abc"))
		@test equalorNaN(R[1, 3], evaluate(d, "", "bc"))
		@test equalorNaN(R[1, 4], evaluate(d, "", "kitten"))

		# Second row is comparing "abc" to the other strings, so:
		@test equalorNaN(R[2, 3], evaluate(d, "abc", "bc"))
		@test equalorNaN(R[2, 4], evaluate(d, "abc", "kitten"))

		# Third row row is comparing "bc" to the other strings, so:
		@test equalorNaN(R[3, 4], evaluate(d, "bc", "kitten"))

		# Matrix is symmetric
		for i in 1:4
			for j in (i+1):4
				@test equalorNaN(R[i, j], R[j, i])
			end
		end

		# Test also the assymetric version
		R2 = pairwise(d, TestStrings1, TestStrings2)
		@test size(R2) == (4, 2)

		@test equalorNaN(R2[1, 1], evaluate(d, "", "mew"))
		@test equalorNaN(R2[1, 2], evaluate(d, "", "ab"))

		@test equalorNaN(R2[2, 1], evaluate(d, "abc", "mew"))
		@test equalorNaN(R2[2, 2], evaluate(d, "abc", "ab"))

		@test equalorNaN(R2[3, 1], evaluate(d, "bc", "mew"))
		@test equalorNaN(R2[3, 2], evaluate(d, "bc", "ab"))

		@test equalorNaN(R2[4, 1], evaluate(d, "kitten", "mew"))
		@test equalorNaN(R2[4, 2], evaluate(d, "kitten", "ab"))

		R3 = pairwise(d, TestStrings2, TestStrings1)
		@test size(R3) == (2, 4)

		for i in 1:length(TestStrings1)
			for j in 1:length(TestStrings2)
				@test equalorNaN(R2[i, j], R3[j, i])
			end
		end

		# Ensure same result if preprocessing for QGramDistances
		if d isa AbstractQGramDistance
			R4 = pairwise(d, TestStrings1; preprocess = true)
			@test typeof(R4) == typeof(R)
			@test size(R4) == size(R)
			for i in 1:size(R4, 1)
				for j in 1:size(R4, 2)
					@test equalorNaN(R4[i, j], R[i, j])
				end
			end
		end
		# ensures missing
		R5 = pairwise(d, TestStrings1missing; preprocess = true)
		@test eltype(R5) == Union{result_type(d, String, String), Missing}
	end
end


text:
"""
    inv(::AbstractSecondOrderTensor)
    inv(::AbstractSymmetricSecondOrderTensor)
    inv(::AbstractFourthOrderTensor)
    inv(::AbstractSymmetricFourthOrderTensor)

Compute the inverse of a tensor.

# Examples
```jldoctest
julia> x = rand(SecondOrderTensor{3})
3×3 Tensor{Tuple{3, 3}, Float64, 2, 9}:
 0.590845  0.460085  0.200586
 0.766797  0.794026  0.298614
 0.566237  0.854147  0.246837

julia> inv(x)
3×3 Tensor{Tuple{3, 3}, Float64, 2, 9}:
  19.7146   -19.2802    7.30384
   6.73809  -10.7687    7.55198
 -68.541     81.4917  -38.8361

julia> x ⋅ inv(x) ≈ one(I)
true
```
"""
function inv end

@inline function _inv(x::AbstractSquareTensor{1})
    typeof(x)(1/det(x))
end

@generated function _inv(x::AbstractSquareTensor{2})
    x_11 = getindex_expr(x, :x, 1, 1)
    x_21 = getindex_expr(x, :x, 2, 1)
    x_12 = getindex_expr(x, :x, 1, 2)
    x_22 = getindex_expr(x, :x, 2, 2)
    exps = [x_22, :(-$x_21), :(-$x_12), x_11]
    quote
        @_inline_meta
        @inbounds typeof(x)($(exps[indices(x)]...)) / det(x)
    end
end

@generated function _inv(x::AbstractSquareTensor{3})
    x_11 = getindex_expr(x, :x, 1, 1)
    x_21 = getindex_expr(x, :x, 2, 1)
    x_31 = getindex_expr(x, :x, 3, 1)
    x_12 = getindex_expr(x, :x, 1, 2)
    x_22 = getindex_expr(x, :x, 2, 2)
    x_32 = getindex_expr(x, :x, 3, 2)
    x_13 = getindex_expr(x, :x, 1, 3)
    x_23 = getindex_expr(x, :x, 2, 3)
    x_33 = getindex_expr(x, :x, 3, 3)
    exps = [:(($x_22*$x_33 - $x_23*$x_32) * dinv),
            :(($x_23*$x_31 - $x_21*$x_33) * dinv),
            :(($x_21*$x_32 - $x_22*$x_31) * dinv),
            :(($x_13*$x_32 - $x_12*$x_33) * dinv),
            :(($x_11*$x_33 - $x_13*$x_31) * dinv),
            :(($x_12*$x_31 - $x_11*$x_32) * dinv),
            :(($x_12*$x_23 - $x_13*$x_22) * dinv),
            :(($x_13*$x_21 - $x_11*$x_23) * dinv),
            :(($x_11*$x_22 - $x_12*$x_21) * dinv)]
    quote
        @_inline_meta
        dinv = 1 / det(x)
        @inbounds typeof(x)($(exps[indices(x)]...))
    end
end

@inline function _inv(x::AbstractSquareTensor{dim}) where {dim}
    typeof(x)(inv(SMatrix{dim, dim}(x)))
end

@generated function toblocks(x::Mat{dim, dim}) where {dim}
    m = dim ÷ 2
    n = dim - m
    inds = independent_indices(x)
    a = [:(Tuple(x)[$(inds[I])]) for I in CartesianIndices((1:m, 1:m))]
    b = [:(Tuple(x)[$(inds[I])]) for I in CartesianIndices((1:m, m+1:dim))]
    c = [:(Tuple(x)[$(inds[I])]) for I in CartesianIndices((m+1:dim, 1:m))]
    d = [:(Tuple(x)[$(inds[I])]) for I in CartesianIndices((m+1:dim, m+1:dim))]
    quote
        @_inline_meta
        A = Mat{$m,$m}($(a...))
        B = Mat{$m,$n}($(b...))
        C = Mat{$n,$m}($(c...))
        D = Mat{$n,$n}($(d...))
        A, B, C, D
    end
end

@generated function fromblocks(A::Mat{m,m}, B::Mat{m,n}, C::Mat{n,m}, D::Mat{n,n}) where {m, n}
    exps = Expr[]
    for j in 1:m
        for i in 1:m
            push!(exps, getindex_expr(A, :A, i, j))
        end
        for i in 1:n
            push!(exps, getindex_expr(C, :C, i, j))
        end
    end
    for j in 1:n
        for i in 1:m
            push!(exps, getindex_expr(B, :B, i, j))
        end
        for i in 1:n
            push!(exps, getindex_expr(D, :D, i, j))
        end
    end
    quote
        @_inline_meta
        @inbounds Mat{$(m+n), $(m+n)}($(exps...))
    end
end

# https://en.wikipedia.org/wiki/Block_matrix#Block_matrix_inversion
# https://core.ac.uk/download/pdf/193046446.pdf
# https://math.stackexchange.com/questions/411492/inverse-of-a-block-matrix-with-singular-diagonal-blocks
function _inv_with_blocks(x::Mat{dim, dim}) where {dim}
    xᵀ = x
    M = xᵀ ⋅ x

    A, B,
    C, D = toblocks(M)

    A⁻¹ = inv(A)
    A⁻¹B = A⁻¹ ⋅ B
    E = inv(D - C ⋅ A⁻¹B)
    A⁻¹BE = A⁻¹B ⋅ E
    CA⁻¹ = C ⋅ A⁻¹

    X = A⁻¹ + A⁻¹BE ⋅ CA⁻¹
    Y = -A⁻¹BE
    Z = -E ⋅ CA⁻¹
    W = E

    fromblocks(X, Y,
               Z, W) ⋅ xᵀ
end

@inline function _inv_with_blocks(x::Tensor{Tuple{@Symmetry{dim, dim}}}) where {dim}
    typeof(x)(_inv_with_blocks(convert(Mat{dim, dim}, x)))
end

@inline inv(x::AbstractSquareTensor) = _inv(x)

# fast inv for dim ≤ 10
@inline fastinv(x::AbstractSquareTensor{1, Float64}) = _inv(x)
@inline fastinv(x::AbstractSquareTensor{2, Float64}) = _inv(x)
@inline fastinv(x::AbstractSquareTensor{3, Float64}) = _inv(x)
@inline fastinv(x::AbstractSquareTensor{4, Float64}) = _inv(x)
@inline fastinv(x::AbstractSquareTensor{5, Float64}) = _inv_with_blocks(x)
@inline fastinv(x::AbstractSquareTensor{6, Float64}) = _inv_with_blocks(x)
@inline fastinv(x::AbstractSquareTensor{7, Float64}) = _inv_with_blocks(x)
@inline fastinv(x::AbstractSquareTensor{8, Float64}) = _inv_with_blocks(x)
@inline fastinv(x::AbstractSquareTensor{9, Float64}) = _inv_with_blocks(x)
@inline fastinv(x::AbstractSquareTensor{10, Float64}) = _inv_with_blocks(x)
@inline fastinv(x::AbstractSquareTensor{<: Any, Float64}) = _inv(x)

# don't use `voigt` or `mandel` for fast computations
@generated function inv(x::FourthOrderTensor{dim}) where {dim}
    L = dim * dim
    quote
        @_inline_meta
        M = Mat{$L, $L}(Tuple(x))
        FourthOrderTensor{dim}(Tuple(inv(M)))
    end
end

@generated function inv(x::SymmetricFourthOrderTensor{dim, T}) where {dim, T}
    S = Space(Symmetry(dim, dim))
    L = ncomponents(S)
    c = Vec{L, T}([i == j ? 1 : √2 for j in 1:dim for i in j:dim])
    coef = c ⊗ c
    quote
        @_inline_meta
        M = _map(*, Mat{$L, $L}(Tuple(x)), $coef)
        M⁻¹ = inv(M)
        SymmetricFourthOrderTensor{dim}(Tuple(_map(/, M⁻¹, $coef)))
    end
end

@inline function _solve(A::AbstractSquareTensor, b::AbstractVec)
    SA = SArray(A)
    Sb = SArray(b)
    Vec(Tuple(SA \ Sb))
end

@inline Base.:\(A::AbstractSquareTensor, b::AbstractVec) = _solve(A, b)

# fast solve for dim ≤ 10
@inline fastsolve(A::AbstractSquareTensor{1, Float64}, b::AbstractVec{1, Float64}) = inv(A) ⋅ b
@inline fastsolve(A::AbstractSquareTensor{2, Float64}, b::AbstractVec{2, Float64}) = inv(A) ⋅ b
@inline fastsolve(A::AbstractSquareTensor{3, Float64}, b::AbstractVec{3, Float64}) = inv(A) ⋅ b
@inline fastsolve(A::AbstractSquareTensor{4, Float64}, b::AbstractVec{4, Float64}) = inv(A) ⋅ b
@inline fastsolve(A::AbstractSquareTensor{5, Float64}, b::AbstractVec{5, Float64}) = inv(A) ⋅ b
@inline fastsolve(A::AbstractSquareTensor{6, Float64}, b::AbstractVec{6, Float64}) = inv(A) ⋅ b
@inline fastsolve(A::AbstractSquareTensor{7, Float64}, b::AbstractVec{7, Float64}) = inv(A) ⋅ b
@inline fastsolve(A::AbstractSquareTensor{8, Float64}, b::AbstractVec{8, Float64}) = inv(A) ⋅ b
@inline fastsolve(A::AbstractSquareTensor{9, Float64}, b::AbstractVec{9, Float64}) = inv(A) ⋅ b
@inline fastsolve(A::AbstractSquareTensor{10, Float64}, b::AbstractVec{10, Float64}) = inv(A) ⋅ b
@inline fastsolve(A::AbstractSquareTensor{dim, Float64}, b::AbstractVec{dim, Float64}) where {dim} = inv(A) ⋅ b


text:
using LinearAlgebra
using FastLapackInterface
using PolynomialMatrixEquations
#using SolveEyePlusMinusAkronB: EyePlusAtKronBWs, generalized_sylvester_solver!

using LinearAlgebra.BLAS

struct BackwardLinearRationalExpectationsWs
    algo::String
    endogenous_nbr::Int64
    exogenous_nbr::Int64
    exogenous_deterministic_nbr::Int64
    forward_indices::Vector{Int64}
    purely_forward_indices::Vector{Int64}
    current_indices::Vector{Int64}
    backward_indices::Vector{Int64}
    both_indices::Vector{Int64}
    static_indices::Vector{Int64}
    dynamic_indices::Vector{Int64}
    current_dynamic_indices::Vector{Int64}
    forward_indices_d::Vector{Int64}
    backward_indices_d::Vector{Int64}
    current_dynamic_indices_d::Vector{Int64}
    current_dynamic_indices_d_j::Vector{Int64}
    exogenous_indices::Vector{Int64}
    static_indices_j::Vector{Int64}
    static_nbr::Int64
    dynamic_nbr::Int64
    forward_nbr::Int64
    backward_nbr::Int64
    both_nbr::Int64
    current_nbr::Int64
    jacobian_static::Matrix{Float64}
    qr_ws::QrWs
    solver_ws::Union{GsSolverWs, CyclicReductionWs}
    a::Matrix{Float64}
    b::Matrix{Float64}
    c::Matrix{Float64}
    d::Matrix{Float64}
    e::Matrix{Float64}
    x::Matrix{Float64}
    ghx::Matrix{Float64}
    gx::Matrix{Float64}
    hx::Matrix{Float64}
    temp1::Matrix{Float64}
    temp2::Matrix{Float64}
    temp3::Matrix{Float64}
    temp4::Matrix{Float64}
    temp5::Matrix{Float64}
    temp6::Matrix{Float64}
    temp7::Matrix{Float64}
    temp8::Matrix{Float64}
    temp9::Matrix{Float64}
    b10::Matrix{Float64}
    b11::Matrix{Float64}
    icolsD::Vector{Int64}
    icolsE::Vector{Int64}
    jcolsD::Vector{Int64}
    jcolsE::Vector{Int64}
    colsUD::Vector{Int64}
    colsUE::Vector{Int64}
    AGplusB::Matrix{Float64}
    linsolve_static_ws::LinSolveWs
    AGplusB_linsolve_ws::LinSolveWs
    #    eye_plus_at_kron_b_ws::EyePlusAtKronBWs

    function BackwardLinearRationalExpectationsWs(algo::String,
                                          endogenous_nbr::Int64,
                                          exogenous_nbr::Int64,
                                          exogenous_deterministic_nbr::Int64,
                                          forward_indices::Vector{Int64},
                                          current_indices::Vector{Int64},
                                          backward_indices::Vector{Int64},
                                          both_indices::Vector{Int64},
                                          static_indices::Vector{Int64})
        static_nbr = length(static_indices)
        forward_nbr = length(forward_indices)
        backward_nbr = length(backward_indices)
        both_nbr = length(both_indices)
        current_nbr = length(current_indices)
        dynamic_nbr = endogenous_nbr - static_nbr
        dynamic_indices = setdiff(collect(1:endogenous_nbr), static_indices)
        current_dynamic_indices = setdiff(current_indices, static_indices)
        purely_forward_indices = setdiff(forward_indices, both_indices)
        forward_indices_d = findall(in(forward_indices), dynamic_indices)
        backward_indices_d = findall(in(backward_indices), dynamic_indices)
        current_dynamic_indices_d = findall(in(current_dynamic_indices), dynamic_indices)
        current_dynamic_indices_d_j = backward_nbr .+ findall(in(dynamic_indices), current_indices)
        exogenous_indices = backward_nbr + current_nbr + forward_nbr .+ (1:exogenous_nbr)
        if static_nbr > 0
            jacobian_static = Matrix{Float64}(undef, endogenous_nbr, static_nbr)
            qr_ws = QrWs(jacobian_static)
            static_indices_j = backward_nbr .+ [findfirst(isequal(x), current_indices) for x in static_indices]
#        else
#            jacobian_static = Matrix{Float64}(undef, 0,0)
#            qr_ws = QrWs(Matrix{Float64}(undef, 0,0))
        end
        if algo == "GS"
            de_order = forward_nbr + backward_nbr
            d = zeros(de_order, de_order)
            e = similar(d)
            solver_ws = GsSolverWs(d, e, backward_nbr)
            a = Matrix{Float64}(undef, 0, 0)
            b = similar(a)
            c = similar(a)
            x = similar(a)
        elseif algo == "CR"
            a = Matrix{Float64}(undef, dynamic_nbr, dynamic_nbr)
            b = similar(a)
            c = similar(a)
            x = similar(a)
            solver_ws = CyclicReductionWs(dynamic_nbr)
            d = Matrix{Float64}(undef, 0, 0)
            e = similar(d)
        end
        ghx = Matrix{Float64}(undef, endogenous_nbr, backward_nbr)
        gx = Matrix{Float64}(undef, forward_nbr, backward_nbr)
        hx = Matrix{Float64}(undef,  backward_nbr, backward_nbr)
        temp1 = Matrix{Float64}(undef, static_nbr, forward_nbr)
        temp2 = Matrix{Float64}(undef, static_nbr, backward_nbr)
        temp3 = Matrix{Float64}(undef, forward_nbr, backward_nbr)
        temp4 = Matrix{Float64}(undef, endogenous_nbr - static_nbr, backward_nbr)
        temp5 = Matrix{Float64}(undef, endogenous_nbr, exogenous_nbr)
        temp6 = Matrix{Float64}(undef, static_nbr, backward_nbr)
        temp7 = Matrix{Float64}(undef, endogenous_nbr, backward_nbr)
        temp8 = Matrix{Float64}(undef, endogenous_nbr, forward_nbr)
        temp9 = Matrix{Float64}(undef, endogenous_nbr, current_nbr)
        b10 = Matrix{Float64}(undef, static_nbr,static_nbr)
        b11 = Matrix{Float64}(undef, static_nbr, endogenous_nbr - static_nbr)
        current_backward_indices = findall(in(backward_indices), current_indices)
        current_forward_indices = findall(in(forward_indices), current_indices)
        # derivatives of current values of variables that are both
        # forward and backward are included in the D matrix
        k1 = findall(in(current_indices), backward_indices)
        icolsD = [k1;  backward_nbr .+ (1:forward_nbr)]
        jcolsD = [backward_nbr .+ current_backward_indices;
                  backward_nbr + current_nbr .+ (1:forward_nbr)]
        k2a = findall(in(purely_forward_indices), forward_indices)
        k2b = findall(in(purely_forward_indices), current_indices)
        icolsE = [1:backward_nbr; backward_nbr .+ k2a]
        jcolsE = [1:backward_nbr; backward_nbr .+ k2b]
        colsUD = findall(in(forward_indices), backward_indices)
        colsUE = backward_nbr .+ findall(in(backward_indices), forward_indices)
        linsolve_static_ws = LinSolveWs(static_nbr)
        AGplusB = Matrix{Float64}(undef, endogenous_nbr, endogenous_nbr)
        AGplusB_linsolve_ws = LinSolveWs(endogenous_nbr)
        #        if m.serially_correlated_exogenous
        #            eye_plus_at_kron_b_ws = EyePlusAtKronBWs(ma, mb, mc, 1)
        #        else
        #            eye_plus_at_kron_b_ws = EyePlusAtKronBWs(1, 1, 1, 1)
        # end
        new(algo, endogenous_nbr, exogenous_nbr, exogenous_deterministic_nbr,
            forward_indices, purely_forward_indices, current_indices,
            backward_indices, both_indices, static_indices, dynamic_indices,
            current_dynamic_indices, forward_indices_d, backward_indices_d,
            current_dynamic_indices_d, current_dynamic_indices_d_j,
            exogenous_indices, static_indices_j,
            static_nbr, dynamic_nbr, forward_nbr, backward_nbr,
            both_nbr, current_nbr, jacobian_static, qr_ws, solver_ws,
            a, b, c, d, e, x,
            ghx, gx, hx, temp1, temp2, temp3, temp4, temp5, temp6, temp7,
            temp8, temp9, b10, b11, icolsD, icolsE, jcolsD, jcolsE,
            colsUD, colsUE, AGplusB, linsolve_static_ws, AGplusB_linsolve_ws)
    end
end

struct BackwardLinearRationalExpectationsResults
    g1::Matrix{Float64}  # full approximation
    gs1::Matrix{Float64} # state transition matrices
    g1_1::SubArray # solution first order derivatives w.r. to state variables
     g1_2::SubArray # solution first order derivatives w.r. to current exogenous variables
#    g1_3::SubArray # solution first order derivatives w.r. to lagged exogenous variables

    function BackwardLinearRationalExpectationsResults(endogenous_nbr::Int64,
                                               exogenous_nbr::Int64,
                                               backward_nbr::Int64)
        nstate = backward_nbr + exogenous_nbr
        g1 =  zeros(endogenous_nbr,(nstate + 1))
        gs1 = zeros(backward_nbr,backward_nbr)
        g1_1 = view(g1, :, 1:backward_nbr)
        g1_2 = view(g1, :, backward_nbr .+ (1:exogenous_nbr))
#        g1_3 = view(g[1], :, backward_nbr + exogenous_nbr .+ lagged_exogenous_nbr)
#        new(g, gs, g1_1, g1_2, g1_3, AGplusB, AGplusB_linsolve_ws)
        new(g1, gs1, g1_1, g1_2)
    end
end

"""
Computes the solution for the static variables:
G_y,static = -B_s,s^{-1}(A_s*Gy,fwrd*Gs + B_s,d*Gy,dynamic + C_s)
"""
function add_static!(results::BackwardLinearRationalExpectationsResults,
                     jacobian::Matrix{Float64},
                     ws::BackwardLinearRationalExpectationsWs)
    # static rows are at the top of the QR transformed Jacobian matrix
    # B_s,s
    fill!(ws.b10, 0.0)
    @inbounds for i = 1:ws.static_nbr
        k = ws.static_indices_j[i]
        for j = 1:ws.static_nbr
            ws.b10[j, i] = jacobian[j, k]
        end
    end
    # B_s,d
    @inbounds for i = 1:length(ws.current_dynamic_indices)
        k1 = ws.current_dynamic_indices_d[i]
        k2 = ws.current_dynamic_indices_d_j[i]
        for j = 1:ws.static_nbr
            ws.b11[j, k1] = jacobian[j, k2]
        end
    end
    # C_s
    @inbounds for i = 1:ws.bacward_nbr
        k = ws.backward_nbr + ws.current_nbr + i
        for j = 1:ws.static_nbr
            ws.temp1[j, i] = jacobian[j, k]
        end
    end
    # A_s
    for i = 1:ws.backward_nbr
        for j = 1:ws.static_nbr
            ws.temp2[j, i] = jacobian[j, i]
        end
    end
    # Gy.fwrd
    @inbounds for i = 1:ws.backward_nbr
        for j = 1:ws.forward_nbr
            ws.temp3[j, i] = results.g1_1[ws.forward_indices[j], i]
        end
    end
    # Gy.dynamic
    @inbounds for i = 1:ws.backward_nbr
        for j = 1:ws.dynamic_nbr
            ws.temp4[j, i] = results.g1_1[ws.dynamic_indices[j], i]
        end
    end
    # ws.temp2 = B_s,d*Gy.dynamic + A_s
    mul!(ws.temp1, ws.b11, ws.temp4, 1.0, 1.0)
    # ws.temp6 = C_s*Gy.fwrd*gs1
    mul!(ws.temp6, ws.temp2, ws.temp3)
    mul!(ws.temp1, ws.temp6, results.gs1, -1.0, -1.0)
    # ws.temp3 = B_s,s\ws.temp1
    linsolve_core!(ws.b10, ws.temp1, ws.linsolve_static_ws)
    @inbounds for i = 1:ws.backward_nbr
        for j=1:ws.static_nbr
            results.g1[ws.static_indices[j],i] = ws.temp1[j,i]
        end
    end
end

function get_abc!(ws::LinearRationalExpectationsWs, jacobian::AbstractMatrix{Float64})
    fill!(ws.a, 0.0)
    fill!(ws.b, 0.0)
    fill!(ws.c, 0.0)

    col_s = ws.backward_nbr + ws.current_nbr + 1
    @inbounds for (i,j) in enumerate(ws.forward_indices_d)
        row_s = ws.static_nbr + 1
        for k = 1:ws.dynamic_nbr
            ws.a[k, j] = jacobian[row_s, col_s]
            row_s += 1
        end
        col_s += 1
    end

    @inbounds for (i,j) in enumerate(ws.current_dynamic_indices_d)
        col_s = ws.current_dynamic_indices_d_j[i]
        row_s = ws.static_nbr + 1
        for k = 1:ws.dynamic_nbr
            ws.b[k, j] = jacobian[row_s, col_s]
            row_s += 1
        end
    end

    @inbounds for (i,j) in enumerate(ws.backward_indices_d)
        row_s = ws.static_nbr + 1
        for k = 1:ws.dynamic_nbr
            ws.c[k, j] = jacobian[row_s, i]
            row_s += 1
        end
    end
end

function get_de!(ws::LinearRationalExpectationsWs, jacobian::AbstractMatrix{Float64})
    n1 = ws.backward_nbr + ws.forward_nbr - ws.both_nbr
    fill!(ws.d, 0.0)
    fill!(ws.e, 0.0)
    for i = 1:length(ws.jcolsD)
        k = ws.jcolsD[i]
        for j = 1:ws.dynamic_nbr
            m = ws.static_nbr + j
            ws.d[j, ws.icolsD[i]] = jacobian[m, k]
        end
    end
    for i = 1:length(ws.jcolsE)
        k = ws.jcolsE[i]
        for j = 1:ws.dynamic_nbr
            m = ws.static_nbr + j
            ws.e[j, ws.icolsE[i]] = -jacobian[m, k]
        end
    end
    for i = 1:ws.both_nbr
        k = ws.dynamic_nbr + i
        ws.d[k, ws.colsUD[i]] = 1.0
    end
    for i = 1:ws.both_nbr
        k = ws.dynamic_nbr + i
        ws.e[k, ws.colsUE[i]] = 1.0
    end
end

"""
Computes LU decomposition of A*G + B
"""
function make_lu_AGplusB!(AGplusB, A, G, B, ws)
    fill!(AGplusB, 0.0)
    for i = 1:ws.current_nbr
        k = ws.current_indices[i]
        for j = 1:ws.endogenous_nbr
            AGplusB[j, k] = B[j, i]
        end
    end
    for i = 1:ws.backward_nbr
        for j = 1:ws.forward_nbr
            ws.temp3[j, i] = G[ws.forward_indices[j], i]
        end
    end
    mul!(ws.temp7, A, ws.temp3)
    for i = 1:ws.backward_nbr
        k = ws.backward_indices[i]
        for j = 1:ws.endogenous_nbr
            AGplusB[j, k] += ws.temp7[j, i]
        end
    end
        #=
    vAGplusB = view(AGplusB, :, ws.current_indices)
    vAGplusB .= B
    ws.temp3 .= view(G, ws.forward_indices, :)
    mul!(ws.temp7, A, ws.temp3)
    vAGplusB = view(AGplusB, :, ws.backward_indices)
    vAGplusB .+= ws.temp7
        =#
    LinSolveAlgo.lu!(AGplusB, ws.AGplusB_linsolve_ws)
end

function solve_for_derivatives_with_respect_to_shocks!(results::LinearRationalExpectationsResults, jacobian::AbstractMatrix, ws::LinearRationalExpectationsWs)
    #=
    if model.lagged_exogenous_nbr > 0
        f6 = view(jacobian,:,model.i_lagged_exogenous)
        for i = 1:model.current_exogenous_nbr
            for j = 1:model.endo_nbr
                results.g1_3[i,j] = -f6[i,j]
            end
        end
        linsolve_core_no_lu!(results.f1g1plusf2, results.g1_3, ws)
    end
    =#
    if ws.exogenous_nbr > 0
        #   results.g1_2 .= .-view(jacobian, :, ws.exogenous_indices)
        for i = 1:ws.exogenous_nbr
            k = ws.exogenous_indices[i]
            for j = 1:ws.endogenous_nbr
                results.g1_2[j, i] = -jacobian[j, k]
            end
        end

#        if ws.serially_correlated_exogenous
            # TO BE DONE
        #        else
        linsolve_core_no_lu!(ws.AGplusB, results.g1_2, ws.AGplusB_linsolve_ws)
#        end
    end
end

function first_order_solver!(results::LinearRationalExpectationsResults,
                             algo::String,
                             jacobian::Matrix,
                             options,
                             ws::LinearRationalExpectationsWs)
    remove_static!(jacobian, ws)
    if algo == "CR"
        get_abc!(ws, jacobian)
        cyclic_reduction!(ws.x, ws.c, ws.b, ws.a, ws.solver_ws, options["cyclic_reduction"]["tol"], 100)
        for i = 1:ws.backward_nbr
            for j = 1:ws.backward_nbr
                results.gs1[j, i] = ws.x[ws.backward_indices_d[j], ws.backward_indices_d[i]]
            end
            for j = 1:(ws.endogenous_nbr - ws.static_nbr)
                results.g1[ws.dynamic_indices[j],i] = ws.x[j, ws.backward_indices_d[i]]
            end
        end
    elseif algo == "GS"
        get_de!(ws, jacobian)
        gs_solver!(ws.solver_ws, ws.d, ws.e, ws.backward_nbr, options["generalized_schur"]["criterium"])
        results.gs1 .= ws.solver_ws.g1
        for i = 1:ws.backward_nbr
            for j = 1:ws.backward_nbr
                x = ws.solver_ws.g1[j,i]
                results.g1[ws.backward_indices[j],i] = x
            end
            for j = 1:(ws.forward_nbr - ws.both_nbr)
                results.g1[ws.purely_forward_indices[j], i] =
                    ws.solver_ws.g2[ws.icolsE[ws.backward_nbr + j] - ws.backward_nbr, i]
            end
        end
    else
        error("Algorithm $algo not recognized")
    end
    if ws.static_nbr > 0
        add_static!(results, jacobian, ws)
    end
    #    A = view(jacobian, :, ws.backward_nbr + ws.current_nbr .+ (1:ws.forward_nbr))
    #    B = view(jacobian, :, ws.backward_nbr .+ ws.current_indices)
    for i = 1:ws.forward_nbr
        k = ws.backward_nbr + ws.current_nbr + i
        for j = 1:ws.endogenous_nbr
            ws.temp8[j, i] = jacobian[j, k]
        end
    end
    for i = 1:ws.current_nbr
        k = ws.backward_nbr + i
        for j = 1:ws.endogenous_nbr
            ws.temp9[j, i] = jacobian[j, k]
        end
    end
    make_lu_AGplusB!(ws.AGplusB, ws.temp8, results.g1_1, ws.temp9, ws)
    solve_for_derivatives_with_respect_to_shocks!(results, jacobian, ws)
end


text:
@testset "dirclustercoeffs" begin
    (A,xy,labels) = load_matrix_network_metadata("celegans")
    (cc, cccyc, ccmid, ccin, ccout, nf) = dirclustercoeffs(A, true, true)
    @test length(cc) == 202
    (maxval, maxind) = findmax(cc)
    @test maxind == 113
    (cc, cccyc, ccmid, ccin, ccout, nf) = dirclustercoeffs(A)
    @test length(cc) == 202
    (maxval, maxind) = findmax(cc)
    @test maxind == 113
    (cc, cccyc, ccmid, ccin, ccout, nf) = dirclustercoeffs(A, true)
    @test length(cc) == 202
    (maxval, maxind) = findmax(cc)
    @test maxind == 113
    # negative weights
    @test_throws Exception dirclustercoeffs(-sparse(1.0I,5,5))
end


text:
pillars(numpill, dist, width) =
    max(0, numpill-1)*100*dist +
    max(0, numpill-2)*width

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
# syntax: proto2
using ProtoBuf
import ProtoBuf.meta

mutable struct UninterpretedOption_NamePart <: ProtoType
    name_part::AbstractString
    is_extension::Bool
    UninterpretedOption_NamePart(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct UninterpretedOption_NamePart
const __req_UninterpretedOption_NamePart = Symbol[:name_part,:is_extension]
meta(t::Type{UninterpretedOption_NamePart}) = meta(t, __req_UninterpretedOption_NamePart, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

mutable struct UninterpretedOption <: ProtoType
    name::Base.Vector{UninterpretedOption_NamePart}
    identifier_value::AbstractString
    positive_int_value::UInt64
    negative_int_value::Int64
    double_value::Float64
    string_value::Array{UInt8,1}
    aggregate_value::AbstractString
    UninterpretedOption(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct UninterpretedOption
const __fnum_UninterpretedOption = Int[2,3,4,5,6,7,8]
meta(t::Type{UninterpretedOption}) = meta(t, ProtoBuf.DEF_REQ, __fnum_UninterpretedOption, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

struct __enum_FieldOptions_CType <: ProtoEnum
    STRING::Int32
    CORD::Int32
    STRING_PIECE::Int32
    __enum_FieldOptions_CType() = new(0,1,2)
end #struct __enum_FieldOptions_CType
const FieldOptions_CType = __enum_FieldOptions_CType()

struct __enum_FieldOptions_JSType <: ProtoEnum
    JS_NORMAL::Int32
    JS_STRING::Int32
    JS_NUMBER::Int32
    __enum_FieldOptions_JSType() = new(0,1,2)
end #struct __enum_FieldOptions_JSType
const FieldOptions_JSType = __enum_FieldOptions_JSType()

mutable struct FieldOptions <: ProtoType
    ctype::Int32
    packed::Bool
    jstype::Int32
    lazy::Bool
    deprecated::Bool
    weak::Bool
    uninterpreted_option::Base.Vector{UninterpretedOption}
    FieldOptions(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct FieldOptions
const __val_FieldOptions = Dict(:ctype => FieldOptions_CType.STRING, :jstype => FieldOptions_JSType.JS_NORMAL, :lazy => false, :deprecated => false, :weak => false)
const __fnum_FieldOptions = Int[1,2,6,5,3,10,999]
meta(t::Type{FieldOptions}) = meta(t, ProtoBuf.DEF_REQ, __fnum_FieldOptions, __val_FieldOptions, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

mutable struct MessageOptions <: ProtoType
    message_set_wire_format::Bool
    no_standard_descriptor_accessor::Bool
    deprecated::Bool
    map_entry::Bool
    uninterpreted_option::Base.Vector{UninterpretedOption}
    MessageOptions(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct MessageOptions
const __val_MessageOptions = Dict(:message_set_wire_format => false, :no_standard_descriptor_accessor => false, :deprecated => false)
const __fnum_MessageOptions = Int[1,2,3,7,999]
meta(t::Type{MessageOptions}) = meta(t, ProtoBuf.DEF_REQ, __fnum_MessageOptions, __val_MessageOptions, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

mutable struct EnumOptions <: ProtoType
    allow_alias::Bool
    deprecated::Bool
    uninterpreted_option::Base.Vector{UninterpretedOption}
    EnumOptions(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct EnumOptions
const __val_EnumOptions = Dict(:deprecated => false)
const __fnum_EnumOptions = Int[2,3,999]
meta(t::Type{EnumOptions}) = meta(t, ProtoBuf.DEF_REQ, __fnum_EnumOptions, __val_EnumOptions, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

mutable struct ExtensionRangeOptions <: ProtoType
    uninterpreted_option::Base.Vector{UninterpretedOption}
    ExtensionRangeOptions(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct ExtensionRangeOptions
const __fnum_ExtensionRangeOptions = Int[999]
meta(t::Type{ExtensionRangeOptions}) = meta(t, ProtoBuf.DEF_REQ, __fnum_ExtensionRangeOptions, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

struct __enum_MethodOptions_IdempotencyLevel <: ProtoEnum
    IDEMPOTENCY_UNKNOWN::Int32
    NO_SIDE_EFFECTS::Int32
    IDEMPOTENT::Int32
    __enum_MethodOptions_IdempotencyLevel() = new(0,1,2)
end #struct __enum_MethodOptions_IdempotencyLevel
const MethodOptions_IdempotencyLevel = __enum_MethodOptions_IdempotencyLevel()

mutable struct MethodOptions <: ProtoType
    deprecated::Bool
    idempotency_level::Int32
    uninterpreted_option::Base.Vector{UninterpretedOption}
    MethodOptions(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct MethodOptions
const __val_MethodOptions = Dict(:deprecated => false, :idempotency_level => MethodOptions_IdempotencyLevel.IDEMPOTENCY_UNKNOWN)
const __fnum_MethodOptions = Int[33,34,999]
meta(t::Type{MethodOptions}) = meta(t, ProtoBuf.DEF_REQ, __fnum_MethodOptions, __val_MethodOptions, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

struct __enum_FileOptions_OptimizeMode <: ProtoEnum
    SPEED::Int32
    CODE_SIZE::Int32
    LITE_RUNTIME::Int32
    __enum_FileOptions_OptimizeMode() = new(1,2,3)
end #struct __enum_FileOptions_OptimizeMode
const FileOptions_OptimizeMode = __enum_FileOptions_OptimizeMode()

mutable struct FileOptions <: ProtoType
    java_package::AbstractString
    java_outer_classname::AbstractString
    java_multiple_files::Bool
    java_generate_equals_and_hash::Bool
    java_string_check_utf8::Bool
    optimize_for::Int32
    go_package::AbstractString
    cc_generic_services::Bool
    java_generic_services::Bool
    py_generic_services::Bool
    php_generic_services::Bool
    deprecated::Bool
    cc_enable_arenas::Bool
    objc_class_prefix::AbstractString
    csharp_namespace::AbstractString
    swift_prefix::AbstractString
    php_class_prefix::AbstractString
    php_namespace::AbstractString
    php_metadata_namespace::AbstractString
    ruby_package::AbstractString
    uninterpreted_option::Base.Vector{UninterpretedOption}
    FileOptions(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct FileOptions
const __val_FileOptions = Dict(:java_multiple_files => false, :java_string_check_utf8 => false, :optimize_for => FileOptions_OptimizeMode.SPEED, :cc_generic_services => false, :java_generic_services => false, :py_generic_services => false, :php_generic_services => false, :deprecated => false, :cc_enable_arenas => false)
const __fnum_FileOptions = Int[1,8,10,20,27,9,11,16,17,18,42,23,31,36,37,39,40,41,44,45,999]
meta(t::Type{FileOptions}) = meta(t, ProtoBuf.DEF_REQ, __fnum_FileOptions, __val_FileOptions, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

mutable struct EnumValueOptions <: ProtoType
    deprecated::Bool
    uninterpreted_option::Base.Vector{UninterpretedOption}
    EnumValueOptions(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct EnumValueOptions
const __val_EnumValueOptions = Dict(:deprecated => false)
const __fnum_EnumValueOptions = Int[1,999]
meta(t::Type{EnumValueOptions}) = meta(t, ProtoBuf.DEF_REQ, __fnum_EnumValueOptions, __val_EnumValueOptions, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

mutable struct OneofOptions <: ProtoType
    uninterpreted_option::Base.Vector{UninterpretedOption}
    OneofOptions(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct OneofOptions
const __fnum_OneofOptions = Int[999]
meta(t::Type{OneofOptions}) = meta(t, ProtoBuf.DEF_REQ, __fnum_OneofOptions, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

mutable struct ServiceOptions <: ProtoType
    deprecated::Bool
    uninterpreted_option::Base.Vector{UninterpretedOption}
    ServiceOptions(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct ServiceOptions
const __val_ServiceOptions = Dict(:deprecated => false)
const __fnum_ServiceOptions = Int[33,999]
meta(t::Type{ServiceOptions}) = meta(t, ProtoBuf.DEF_REQ, __fnum_ServiceOptions, __val_ServiceOptions, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

struct __enum_FieldDescriptorProto_Type <: ProtoEnum
    TYPE_DOUBLE::Int32
    TYPE_FLOAT::Int32
    TYPE_INT64::Int32
    TYPE_UINT64::Int32
    TYPE_INT32::Int32
    TYPE_FIXED64::Int32
    TYPE_FIXED32::Int32
    TYPE_BOOL::Int32
    TYPE_STRING::Int32
    TYPE_GROUP::Int32
    TYPE_MESSAGE::Int32
    TYPE_BYTES::Int32
    TYPE_UINT32::Int32
    TYPE_ENUM::Int32
    TYPE_SFIXED32::Int32
    TYPE_SFIXED64::Int32
    TYPE_SINT32::Int32
    TYPE_SINT64::Int32
    __enum_FieldDescriptorProto_Type() = new(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18)
end #struct __enum_FieldDescriptorProto_Type
const FieldDescriptorProto_Type = __enum_FieldDescriptorProto_Type()

struct __enum_FieldDescriptorProto_Label <: ProtoEnum
    LABEL_OPTIONAL::Int32
    LABEL_REQUIRED::Int32
    LABEL_REPEATED::Int32
    __enum_FieldDescriptorProto_Label() = new(1,2,3)
end #struct __enum_FieldDescriptorProto_Label
const FieldDescriptorProto_Label = __enum_FieldDescriptorProto_Label()

mutable struct FieldDescriptorProto <: ProtoType
    name::AbstractString
    number::Int32
    label::Int32
    _type::Int32
    type_name::AbstractString
    extendee::AbstractString
    default_value::AbstractString
    oneof_index::Int32
    json_name::AbstractString
    options::FieldOptions
    FieldDescriptorProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct FieldDescriptorProto
const __fnum_FieldDescriptorProto = Int[1,3,4,5,6,2,7,9,10,8]
meta(t::Type{FieldDescriptorProto}) = meta(t, ProtoBuf.DEF_REQ, __fnum_FieldDescriptorProto, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

mutable struct DescriptorProto_ExtensionRange <: ProtoType
    start::Int32
    _end::Int32
    options::ExtensionRangeOptions
    DescriptorProto_ExtensionRange(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct DescriptorProto_ExtensionRange

mutable struct MethodDescriptorProto <: ProtoType
    name::AbstractString
    input_type::AbstractString
    output_type::AbstractString
    options::MethodOptions
    client_streaming::Bool
    server_streaming::Bool
    MethodDescriptorProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct MethodDescriptorProto
const __val_MethodDescriptorProto = Dict(:client_streaming => false, :server_streaming => false)
meta(t::Type{MethodDescriptorProto}) = meta(t, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, __val_MethodDescriptorProto, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

mutable struct EnumValueDescriptorProto <: ProtoType
    name::AbstractString
    number::Int32
    options::EnumValueOptions
    EnumValueDescriptorProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct EnumValueDescriptorProto

mutable struct EnumDescriptorProto_EnumReservedRange <: ProtoType
    start::Int32
    _end::Int32
    EnumDescriptorProto_EnumReservedRange(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct EnumDescriptorProto_EnumReservedRange

mutable struct EnumDescriptorProto <: ProtoType
    name::AbstractString
    value::Base.Vector{EnumValueDescriptorProto}
    options::EnumOptions
    reserved_range::Base.Vector{EnumDescriptorProto_EnumReservedRange}
    reserved_name::Base.Vector{AbstractString}
    EnumDescriptorProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct EnumDescriptorProto

mutable struct OneofDescriptorProto <: ProtoType
    name::AbstractString
    options::OneofOptions
    OneofDescriptorProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct OneofDescriptorProto

mutable struct DescriptorProto_ReservedRange <: ProtoType
    start::Int32
    _end::Int32
    DescriptorProto_ReservedRange(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct DescriptorProto_ReservedRange

mutable struct DescriptorProto <: ProtoType
    name::AbstractString
    field::Base.Vector{FieldDescriptorProto}
    extension::Base.Vector{FieldDescriptorProto}
    nested_type::Base.Vector{DescriptorProto}
    enum_type::Base.Vector{EnumDescriptorProto}
    extension_range::Base.Vector{DescriptorProto_ExtensionRange}
    oneof_decl::Base.Vector{OneofDescriptorProto}
    options::MessageOptions
    reserved_range::Base.Vector{DescriptorProto_ReservedRange}
    reserved_name::Base.Vector{AbstractString}
    DescriptorProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct DescriptorProto
const __fnum_DescriptorProto = Int[1,2,6,3,4,5,8,7,9,10]
meta(t::Type{DescriptorProto}) = meta(t, ProtoBuf.DEF_REQ, __fnum_DescriptorProto, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

mutable struct ServiceDescriptorProto <: ProtoType
    name::AbstractString
    method::Base.Vector{MethodDescriptorProto}
    options::ServiceOptions
    ServiceDescriptorProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct ServiceDescriptorProto

mutable struct SourceCodeInfo_Location <: ProtoType
    path::Base.Vector{Int32}
    span::Base.Vector{Int32}
    leading_comments::AbstractString
    trailing_comments::AbstractString
    leading_detached_comments::Base.Vector{AbstractString}
    SourceCodeInfo_Location(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct SourceCodeInfo_Location
const __fnum_SourceCodeInfo_Location = Int[1,2,3,4,6]
const __pack_SourceCodeInfo_Location = Symbol[:path,:span]
meta(t::Type{SourceCodeInfo_Location}) = meta(t, ProtoBuf.DEF_REQ, __fnum_SourceCodeInfo_Location, ProtoBuf.DEF_VAL, true, __pack_SourceCodeInfo_Location, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

mutable struct SourceCodeInfo <: ProtoType
    location::Base.Vector{SourceCodeInfo_Location}
    SourceCodeInfo(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct SourceCodeInfo

mutable struct FileDescriptorProto <: ProtoType
    name::AbstractString
    package::AbstractString
    dependency::Base.Vector{AbstractString}
    public_dependency::Base.Vector{Int32}
    weak_dependency::Base.Vector{Int32}
    message_type::Base.Vector{DescriptorProto}
    enum_type::Base.Vector{EnumDescriptorProto}
    service::Base.Vector{ServiceDescriptorProto}
    extension::Base.Vector{FieldDescriptorProto}
    options::FileOptions
    source_code_info::SourceCodeInfo
    syntax::AbstractString
    FileDescriptorProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct FileDescriptorProto
const __fnum_FileDescriptorProto = Int[1,2,3,10,11,4,5,6,7,8,9,12]
meta(t::Type{FileDescriptorProto}) = meta(t, ProtoBuf.DEF_REQ, __fnum_FileDescriptorProto, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

mutable struct FileDescriptorSet <: ProtoType
    file::Base.Vector{FileDescriptorProto}
    FileDescriptorSet(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct FileDescriptorSet

mutable struct GeneratedCodeInfo_Annotation <: ProtoType
    path::Base.Vector{Int32}
    source_file::AbstractString
    _begin::Int32
    _end::Int32
    GeneratedCodeInfo_Annotation(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct GeneratedCodeInfo_Annotation
const __pack_GeneratedCodeInfo_Annotation = Symbol[:path]
meta(t::Type{GeneratedCodeInfo_Annotation}) = meta(t, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, __pack_GeneratedCodeInfo_Annotation, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

mutable struct GeneratedCodeInfo <: ProtoType
    annotation::Base.Vector{GeneratedCodeInfo_Annotation}
    GeneratedCodeInfo(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct GeneratedCodeInfo

export FileDescriptorSet, FileDescriptorProto, DescriptorProto_ExtensionRange, DescriptorProto_ReservedRange, DescriptorProto, ExtensionRangeOptions, FieldDescriptorProto_Type, FieldDescriptorProto_Label, FieldDescriptorProto, OneofDescriptorProto, EnumDescriptorProto_EnumReservedRange, EnumDescriptorProto, EnumValueDescriptorProto, ServiceDescriptorProto, MethodDescriptorProto, FileOptions_OptimizeMode, FileOptions, MessageOptions, FieldOptions_CType, FieldOptions_JSType, FieldOptions, OneofOptions, EnumOptions, EnumValueOptions, ServiceOptions, MethodOptions_IdempotencyLevel, MethodOptions, UninterpretedOption_NamePart, UninterpretedOption, SourceCodeInfo_Location, SourceCodeInfo, GeneratedCodeInfo_Annotation, GeneratedCodeInfo


text:
module IdentityCellNumbers

using Gridap

export IdentityCellNumber
import Gridap: reindex
import Base: size
import Base: getindex

struct IdentityCellNumber{T} <: IndexCellValue{T,1}
  length::Int
end

function IdentityCellNumber(::Type{T},l::Integer) where T <: Integer
  IdentityCellNumber{T}(l)
end

function getindex(c::IdentityCellNumber{T},i::Integer) where T
  @assert i > 0
  @assert i <= c.length
  j::T = i
  j
end

size(c::IdentityCellNumber) = (c.length,)

function reindex(values::IndexCellValue, indices::IdentityCellNumber{T}) where T
  @assert length(values) == length(indices)
  values
end

function reindex(values::IndexCellNumber, indices::IdentityCellNumber{T}) where T
  @assert length(values) == length(indices)
  values
end

function reindex(values::IndexCellArray, indices::IdentityCellNumber{T}) where T
  @assert length(values) == length(indices)
  values
end

end # module


text:
# a cache for library handles

export HandleCache

struct HandleCache{K,V}
    active_handles::Set{Pair{K,V}}      # for debugging, and to prevent handle finalization
    idle_handles::Dict{K,Vector{V}}
    lock::ReentrantLock

    max_entries::Int

    function HandleCache{K,V}(max_entries::Int=32) where {K,V}
        return new{K,V}(Set{Pair{K,V}}(), Dict{K,Vector{V}}(), ReentrantLock(), max_entries)
    end
end

# remove a handle from the cache, or create a new one
function Base.pop!(f::Function, cache::HandleCache{K,V}, key) where {K,V}
    function check_cache(f::Function=()->nothing)
        try
            GC.enable_finalizers(false)
            lock(cache.lock) do
                handle = if !haskey(cache.idle_handles, key) || isempty(cache.idle_handles[key])
                    f()
                else
                    pop!(cache.idle_handles[key])
                end

                if handle !== nothing
                    push!(cache.active_handles, key=>handle)
                end

                return handle
            end
        finally
            GC.enable_finalizers(true)
        end
    end

    handle = check_cache()

    if handle === nothing
        # if we didn't find anything, perform a quick GC collection to free up old handles.
        GC.gc(false)

        handle = check_cache(f)
    end

    return handle::V
end

# put a handle in the cache, or destroy it if it doesn't fit
function Base.push!(f::Function, cache::HandleCache{K,V}, key::K, handle::V) where {K,V}
    # XXX: take this lock in a normal way once we have JuliaLang/julia#35689
    @spinlock cache.lock begin
        delete!(cache.active_handles, key=>handle)

        if haskey(cache.idle_handles, key)
            if length(cache.idle_handles[key]) > cache.max_entries
                f()
            else
                push!(cache.idle_handles[key], handle)
            end
        else
            cache.idle_handles[key] = [handle]
        end
    end
end


text:
using InteractionNetworkModels, Plots, Distributions
using StatsPlots, Plots.Measures, StatsBase

E = [[1,2,1,2],
    [1,2,1],
    [3,4,3],
    [3,4],
    [1,2]]
d = MatchingDist(FastLSP(100))
K_inner, K_outer = (DimensionRange(2,50), DimensionRange(1,50))
model = SIM(
    E, 4.4,
    d,
    1:50,
    K_inner, K_outer)

mcmc_sampler = SimMcmcInsertDelete(
    ν_ed=4, ν_td=1, β=0.7,
    len_dist=TrGeometric(0.9, 1, model.K_inner.u),
    burn_in=2000, lag=50
)

mcmc_sampler_sp = SimMcmcInsertDeleteSubpath(
    ν_ed=7, ν_td=1, β=0.7,
    len_dist=TrGeometric(0.9, 1, model.K_inner.u),
    burn_in=3000, lag=50
)

@time test = mcmc_sampler(
    model, desired_samples=500)

plot(test)
summaryplot(test)
test.sample

@time mcmc_out = mcmc_sampler(
    model,
    desired_samples=50,
    lag=500,
    burn_in=10000
    )

plot(mcmc_out)
summaryplot(mcmc_out)

data = mcmc_out.sample
E_prior = SIM(E, 0.1, model.dist, model.V, model.K_inner, model.K_outer)
γ_prior = Uniform(0.5,7.0)

target = SimPosterior(data, E_prior, γ_prior)

# Construct posterior sampler
posterior_sampler = SimIexInsertDelete(
    mcmc_sampler,
    len_dist=TrGeometric(0.9,K_inner.l,K_inner.u),
    ν_ed=1, ν_td=1,
    β=0.7, ε=0.3
)

# Mode Conditional
E_init, ind = sample_frechet_mean(target.data, target.dist)
d(E_init, E)

@time posterior_out = posterior_sampler(
    target,
    desired_samples=4000, lag=1, burn_in=0,
    S_init=E_init[3:3], γ_init=4.5
);

posterior_out.S_sample
plot(posterior_out, E)
plot(length.(posterior_out.S_sample))
findall(diff(length.(posterior_out.S_sample)).>0)

plot(posterior_out.suff_stat_trace)

S_est = deepcopy(posterior_out.S_sample[end])

mapreduce(x->d_lsp(x,S_est), +, data)
mapreduce(x->d_lsp(x,E), +, data)

plot(length.(data))

i=1
tmp1, C1 = d_lsp(E, data[i])
tmp2, C2 = d_lsp(data[i], E)
C1
C2

pairwise(LSP(), E, data[i])
pairwise(LSP(), data[i], E)
d_lsp(S_est, data[i])

using Hungarian
data[1]
LSP()(data[i][1], E[1])
LSP()(E[1], data[i][1])
FastLSP(100)(data[i][1], E[1])
FastLSP(100)(E[1], data[i][1])

E[1]
data[i][1]

[d_lsp(x,E) for x in data]

d_lsp(E, S_est)
d_lsp(S_est, E)
tmp, C = d_lsp(S_est, E)
C
S_est
E

d([[1]], [[2]])
d([[2]],[[1]])

d_lcs = MatchingDist(LCS())

d_lcs(E, S_est)
d_lcs(S_est, E)

map(x->d(x,E), data)
E
d

plot(map(x->d_lsp(E,x), data))
plot!(map(x->d_lsp(x,S_est), data))
S_est

print_matching(d, S, data[1])

print_map_est(posterior_out)

plot(posterior_out, E)

d([[1,1,1],[2,2]], [[2,1]])

# Dispersion conditional
E_fix = Eᵐ
posterior_out = posterior_sampler(
    target,
    E_fix,
    desired_samples=1000, lag=1, burn_in=0,
    init=4.6
)
plot(posterior_out)


text:
# Loads all files
module EffectiveWaves

include("export_import.jl")
# try using BlackBoxOptim end
# using IterTools
using RecipesBase
using Optim
using OffsetArrays
using ApproxFun
using LinearAlgebra
using WignerSymbols
using GSL

include("specialfunctions.jl")
include("particle.jl")
include("t-matrix.jl")

include("effective_waves/effective_waves_export.jl")
include("average_waves/average_waves_export.jl")
include("match_waves/match_waves.jl")
include("match_waves/match_arrays.jl")
include("match_waves/reflection.jl")

include("materials.jl")

include("plot/graphics.jl")
include("plot/plot.jl")

end # module