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

using Markdown
using InteractiveUtils

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

# ╔═╡ 1a6e1853-6db1-4074-bce0-5f274351cece
md"""
We define a _diffing system_ for Julia `Dict`s, which is analogous to the diffing system of immer.js.

This notebook is part of Pluto's source code (included in `src/webserver/Dynamic.jl`).
"""

# ╔═╡ 49fc1f97-3b8f-4297-94e5-2e24c001d35c
md"""
## Example

Computing a diff:
"""

# ╔═╡ d8e73b90-24c5-4e50-830b-b1dbe6224c8e
dict_1 = Dict{String,Any}(
	"a" => 1,
	"b" => Dict(
		"c" => [3,4],
		"d" => 99,
	),
	"e" => "hello!"
);

# ╔═╡ 19646596-b35b-44fa-bfcf-891f9ffb748c
dict_2 = Dict{String,Any}(
	"a" => 1,
	"b" => Dict(
		"c" => [3,4,5],
		"d" => 99,
		"🏝" => "👍",
	),
);

# ╔═╡ 9d2c07d9-16a9-4b9f-a375-2adb6e5b907a
md"""
Applying a set of patches:
"""

# ╔═╡ 336bfd4f-8a8e-4a2d-be08-ee48d6a9f747
md"""
## JSONPatch objects
"""

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

# ╔═╡ bd0d46bb-3e58-4522-bae0-83eb799196c4
PatchPath = Vector

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

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

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

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

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

# ╔═╡ 73631aea-5e93-4da2-a32d-649029660d4e
const Patches = Vector{JSONPatch}

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

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

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

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

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

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

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

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

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

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

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

# ╔═╡ 921a130e-b028-4f91-b077-3bd79dcb6c6d
function Base.convert(::Type{JSONPatch}, patch_dict::Dict)
	op = patch_dict["op"]
	if op == "add"
		AddPatch(patch_dict["path"], patch_dict["value"])
	elseif op == "remove"
		RemovePatch(patch_dict["path"])
	elseif op == "replace"
		ReplacePatch(patch_dict["path"], patch_dict["value"])
	else
		throw(ArgumentError("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))

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

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

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

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

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

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

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

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

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

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

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

# ╔═╡ 0b50f6b2-8e85-4565-9f04-f99c913b4592
const use_triple_equals_for_arrays = Ref(false)

# ╔═╡ 59e94cb2-c2f9-4f6c-9562-45e8c15931af
function diff(old::T, new::T) where T <: AbstractArray
	if use_triple_equals_for_arrays[] ?
		((old === new) || isequal(old, new)) :
		isequal(old, new)
		NoChanges
	else
		JSONPatch[ReplacePatch([], new)]
	end
end

# ╔═╡ 5e360fcd-9943-4a17-9672-f1fded2f7e3a
function diff(old::T, new::T) where T
	if old === missing && new !== missing
		JSONPatch[ReplacePatch([], new)]
	elseif old !== missing && new === missing
		JSONPatch[RemovePatch([])]
	elseif old === new || isequal(old, new)
		NoChanges
	else
		JSONPatch[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}) 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)
	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

	# same as above but faster:

	for (key1, val1) in o1
		for change in diff(val1, get(o2, key1, nothing))
			push!(changes, wrappath([key1], change))
		end
	end
	for (key2, val2) in o2
		if !haskey(o1, key2)
			for change in diff(nothing, val2)
				push!(changes, wrappath([key2], change))
			end
		end
	end
	changes
end

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

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

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

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

# ╔═╡ 7ca087b8-73ac-49ea-9c5a-2971f0da491f
example_patches = diff(dict_1, dict_2)

# ╔═╡ 59b46bfe-da74-43af-9c11-cb0bdb2c13a2
md"""
### Dict example
"""

# ╔═╡ 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)

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

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

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

# ╔═╡ 3924953f-787a-4912-b6ee-9c9d3030f0f0
md"""
### Large Dict example 1
"""

# ╔═╡ 80689881-1b7e-49b2-af97-9e3ab639d006
big_array = rand(UInt8, 1_000_000)

# ╔═╡ fd22b6af-5fd2-428a-8291-53e223ea692c
big_string = repeat('a', 1_000_000);

# ╔═╡ bcd5059b-b0d2-49d8-a756-92349aa56aca
large_dict_1 = Dict{String,Any}(
	"cell_$(i)" => Dict{String,Any}(
		"x" => 1,
		"y" => big_array,
		"z" => big_string,
	)
	for i in 1:10
);

# ╔═╡ e7fd6bab-c114-4f3e-b9ad-1af2d1147770
begin
	large_dict_2 = Dict{String,Any}(
		"cell_$(i)" => Dict{String,Any}(
			"x" => 1,
			"y" => big_array,
			"z" => big_string,
		)
		for i in 1:10
	)
	large_dict_2["cell_5"]["y"] = [2,20]
	delete!(large_dict_2, "cell_2")
	large_dict_2["hello"] = Dict("a" => 1, "b" => 2)
	large_dict_2
end;

# ╔═╡ 43c36ab7-e9ac-450a-8abe-435412f2be1d
diff(large_dict_1, large_dict_2)

# ╔═╡ 1cf22fe6-4b58-4220-87a1-d7a18410b4e8
md"""
With `===` comparison for arrays:
"""

# ╔═╡ ffb01ab4-e2e3-4fa4-8c0b-093d2899a536
md"""
### Large Dict example 2
"""

# ╔═╡ 8188de75-ae6e-48aa-9495-111fd27ffd26
many_items_1 = Dict{String,Any}(
	"cell_$(i)" => Dict{String,Any}(
		"x" => 1,
		"y" => [2,3],
		"z" => "four",
	)
	for i in 1:100
)

# ╔═╡ d807195e-ba27-4015-92a7-c9294d458d47
begin
	many_items_2 = deepcopy(many_items_1)
	many_items_2["cell_5"]["y"][2] = 20
	delete!(many_items_2, "cell_2")
	many_items_2["hello"] = Dict("a" => 1, "b" => 2)
	many_items_2
end

# ╔═╡ 2e91a1a2-469c-4123-a0d7-3dcc49715738
diff(many_items_1, many_items_2)

# ╔═╡ b8061c1b-dd03-4cd1-b275-90359ae2bb39
fairly_equal(a,b) = all(x -> x in b, a) && all(y -> y in a, b)

# ╔═╡ aeab3363-08ba-47c2-bd33-04a004ed72c4
diff(many_items_1, many_items_1)

# ╔═╡ 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
	changes
end

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

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

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

# ╔═╡ d2af2a4b-8982-4e43-9fd7-0ecfdfb70511
const strict_applypatch = Ref(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 applypatch!(value, patches::Array{JSONPatch})
	for patch in patches
		applypatch!(value, patch)
	end
	return value
end

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

# ╔═╡ dd87ca7e-2de1-11eb-2ec3-d5721c32f192
function applypatch!(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_applypatch[]
				@assert get(subvalue, key, nothing) === nothing
			end
			subvalue[key] = patch.value
		else
			key = Symbol(last)
			if strict_applypatch[]
				@assert getproperty(subvalue, key) === nothing
			end
			setproperty!(subvalue, key, patch.value)
		end
	end
	return value
end

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

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

# ╔═╡ 6509d62e-77b6-499c-8dab-4a608e44720a
function applypatch!(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_applypatch[]
				@assert get(subvalue, key, nothing) !== nothing
			end
			subvalue[key] = patch.value
		else
			key = Symbol(last)
			if strict_applypatch[]
				@assert getproperty(subvalue, key) !== nothing
			end
			setproperty!(subvalue, key, patch.value)
		end
	end
	return value
end

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

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

# ╔═╡ ddaf5b66-2de1-11eb-3348-b905b94a984b
function applypatch!(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_applypatch[]
				@assert get(subvalue, key, nothing) !== nothing
			end
			delete!(subvalue, key)
		else
			key = Symbol(last)
			if strict_applypatch[]
				@assert getproperty(subvalue, key) !== nothing
			end
			setproperty!(subvalue, key, nothing)
		end
	end
	return value
end

# ╔═╡ e65d483a-4c13-49ba-bff1-1d54de78f534
let
	dict_1_copy = deepcopy(dict_1)
	applypatch!(dict_1_copy, example_patches)
end

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

# ╔═╡ 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

# ╔═╡ b05fcb88-3781-45d0-9f24-e88c339a72e5
macro test2(expr)
	quote nothing end
end

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

# ╔═╡ e907d862-2de1-11eb-11a9-4b3ac37cb0f3
function skip_as_script(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 skip_as_script(ex) skip_as_script(__module__) ? esc(ex) : nothing end

# ╔═╡ c2c2b057-a88f-4cc6-ada4-fc55ac29931e
"The opposite of `@skip_as_script`"
macro only_as_script(ex) skip_as_script(__module__) ? nothing : esc(ex) end

# ╔═╡ e748600a-2de1-11eb-24be-d5f0ecab8fa4
# Only define this in Pluto - assume we are `using Test` otherwise
begin
	@skip_as_script 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
	# Do nothing inside pluto (so we don't need to have Test as dependency)
	# test/Firebasey is `using Test` before including this file
	@only_as_script ((@isdefined Test) ? nothing : macro test(expr) quote nothing end end)
end

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

# ╔═╡ 2983f6d4-c1ca-4b66-a2d3-f858b0df2b4c
@test fairly_equal(diff(large_dict_1, large_dict_2), [
	ReplacePatch(["cell_5","y"], [2,20]),
	RemovePatch(["cell_2"]),
	AddPatch(["hello"], Dict("b" => 2, "a" => 1)),
])

# ╔═╡ fdc427f0-dfe8-4114-beca-48fc15434534
@test isempty(diff(many_items_1, many_items_1))

# ╔═╡ 61b81430-d26e-493c-96da-b6818e58c882
@test fairly_equal(diff(many_items_1, many_items_2), [
	ReplacePatch(["cell_5","y"], [2,20]),
	RemovePatch(["cell_2"]),
	AddPatch(["hello"], Dict("b" => 2, "a" => 1)),
])

# ╔═╡ 62de3e79-4b4e-41df-8020-769c3c255c3e
@test isempty(diff(many_items_1, many_items_1))

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

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

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

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

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

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

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

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

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

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

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

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

# ╔═╡ 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 / value.times_ran))
					$(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

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

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

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

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

# ╔═╡ 8c069015-d922-4c60-9340-8d65c80b1a06
@skip_as_script @track for _ in 1:1000 diff(large_dict_1, large_dict_1) end

# ╔═╡ bc9a0822-1088-4ee7-8c79-98e06fd50f11
@skip_as_script @track for _ in 1:1000 diff(large_dict_1, large_dict_2) end

# ╔═╡ ddf1090c-5239-41df-ae4d-70aeb3a75f2b
@skip_as_script let
	old = use_triple_equals_for_arrays[]
	use_triple_equals_for_arrays[] = true

	result = @track for _ in 1:1000 diff(large_dict_1, large_dict_1) end

	use_triple_equals_for_arrays[] = old
	result
end

# ╔═╡ 88009db3-f40e-4fd0-942a-c7f4a7eecb5a
@skip_as_script let
	old = use_triple_equals_for_arrays[]
	use_triple_equals_for_arrays[] = true

	result = @track for _ in 1:1000 diff(large_dict_1, large_dict_2) end

	use_triple_equals_for_arrays[] = old
	result
end

# ╔═╡ c287009f-e864-45d2-a4d0-a525c988a6e0
@skip_as_script @track for _ in 1:1000 diff(many_items_1, many_items_1) end

# ╔═╡ 67a1ae27-f7df-4f84-8809-1cc6a9bcd1ce
@skip_as_script @track for _ in 1:1000 diff(many_items_1, many_items_2) end

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

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

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

# ╔═╡ Cell order:
# ╟─d948dc6e-2de1-11eb-19e7-cb3bb66353b6
# ╟─1a6e1853-6db1-4074-bce0-5f274351cece
# ╟─49fc1f97-3b8f-4297-94e5-2e24c001d35c
# ╠═d8e73b90-24c5-4e50-830b-b1dbe6224c8e
# ╠═19646596-b35b-44fa-bfcf-891f9ffb748c
# ╠═7ca087b8-73ac-49ea-9c5a-2971f0da491f
# ╟─9d2c07d9-16a9-4b9f-a375-2adb6e5b907a
# ╠═e65d483a-4c13-49ba-bff1-1d54de78f534
# ╟─336bfd4f-8a8e-4a2d-be08-ee48d6a9f747
# ╠═db116c0a-2de1-11eb-2a56-872af797c547
# ╠═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
# ╠═73631aea-5e93-4da2-a32d-649029660d4e
# ╠═0fd3e910-abcc-4421-9d0b-5cfb90034338
# ╟─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
# ╠═230bafe2-aaa7-48f0-9fd1-b53956281684
# ╟─07eeb122-6706-4544-a007-1c8d6581eec8
# ╠═b48e2c08-a94a-4247-877d-949d92dde626
# ╟─c59b30b9-f702-41f1-bb2e-1736c8cd5ede
# ╠═fafcb8b8-cde9-4f99-9bab-8128025953a4
# ╟─7feeee3a-3aec-47ce-b8d7-74a0d9b0b381
# ╠═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
# ╟─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
# ╟─daf9ec12-2de1-11eb-3a8d-59d9c2753134
# ╠═0b50f6b2-8e85-4565-9f04-f99c913b4592
# ╠═59e94cb2-c2f9-4f6c-9562-45e8c15931af
# ╠═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
# ╟─59b46bfe-da74-43af-9c11-cb0bdb2c13a2
# ╟─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
# ╟─3924953f-787a-4912-b6ee-9c9d3030f0f0
# ╠═80689881-1b7e-49b2-af97-9e3ab639d006
# ╠═fd22b6af-5fd2-428a-8291-53e223ea692c
# ╠═bcd5059b-b0d2-49d8-a756-92349aa56aca
# ╠═e7fd6bab-c114-4f3e-b9ad-1af2d1147770
# ╠═43c36ab7-e9ac-450a-8abe-435412f2be1d
# ╟─2983f6d4-c1ca-4b66-a2d3-f858b0df2b4c
# ╟─fdc427f0-dfe8-4114-beca-48fc15434534
# ╟─8c069015-d922-4c60-9340-8d65c80b1a06
# ╟─bc9a0822-1088-4ee7-8c79-98e06fd50f11
# ╟─1cf22fe6-4b58-4220-87a1-d7a18410b4e8
# ╟─ddf1090c-5239-41df-ae4d-70aeb3a75f2b
# ╟─88009db3-f40e-4fd0-942a-c7f4a7eecb5a
# ╟─ffb01ab4-e2e3-4fa4-8c0b-093d2899a536
# ╠═8188de75-ae6e-48aa-9495-111fd27ffd26
# ╠═d807195e-ba27-4015-92a7-c9294d458d47
# ╠═2e91a1a2-469c-4123-a0d7-3dcc49715738
# ╟─61b81430-d26e-493c-96da-b6818e58c882
# ╠═b8061c1b-dd03-4cd1-b275-90359ae2bb39
# ╠═aeab3363-08ba-47c2-bd33-04a004ed72c4
# ╟─62de3e79-4b4e-41df-8020-769c3c255c3e
# ╟─c287009f-e864-45d2-a4d0-a525c988a6e0
# ╟─67a1ae27-f7df-4f84-8809-1cc6a9bcd1ce
# ╟─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
# ╟─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
# ╠═dd87ca7e-2de1-11eb-2ec3-d5721c32f192
# ╟─c3e4738f-4568-4910-a211-6a46a9d447ee
# ╟─a11e4082-4ff4-4c1b-9c74-c8fa7dcceaa6
# ╟─0f094932-10e5-40f9-a3fc-db27a85b4999
# ╟─be6b6fc4-e12a-4cef-81d8-d5115fda50b7
# ╠═6509d62e-77b6-499c-8dab-4a608e44720a
# ╟─a560fdca-ee12-469c-bda5-62d7203235b8
# ╟─f1dde1bd-3fa4-48b7-91ed-b2f98680fcc1
# ╟─01e3417e-334e-4a8d-b086-4bddc42737b3
# ╟─f3ef354b-b480-4b48-8358-46dbf37e1d95
# ╠═ddaf5b66-2de1-11eb-3348-b905b94a984b
# ╟─96a80a23-7c56-4c41-b489-15bc1c4e3700
# ╟─df41caa7-f0fc-4b0d-ab3d-ebdab4804040
# ╟─fac65755-2a2a-4a3c-b5a8-fc4f6d256754
# ╟─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
# ╠═b05fcb88-3781-45d0-9f24-e88c339a72e5
# ╠═e7e8d076-2de1-11eb-0214-8160bb81370a
# ╟─e8d0c98a-2de1-11eb-37b9-e1df3f5cfa25
# ╠═e907d862-2de1-11eb-11a9-4b3ac37cb0f3
# ╠═e924a0be-2de1-11eb-2170-71d56e117af2
# ╠═c2c2b057-a88f-4cc6-ada4-fc55ac29931e
# ╠═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
# ╠═1a26eed8-670c-43bf-9726-2db84b1afdab
# ╠═0e1c6442-9040-49d9-b754-173583db7ba2
# ╠═7618aef7-1884-4e32-992d-0fd988e1ab20
# ╠═a3e8fe70-cbf5-4758-a0f2-d329d138728c


text:
batchsize(batch::Union{Tuple,NamedTuple}) = batchsize(batch[1])
batchsize(batch::Dict) = batchsize(batch[first(keys(batch))])
batchsize(batch::AbstractArray{T, N}) where {T, N} = size(batch, N)

unbatch(batch) = collect(obsslices(batch))

obsslices(batch) =
    (obsslice(batch, i) for i in 1:batchsize(batch))

function obsslice(batch::AbstractArray{T, N}, i) where {T, N}
    return view(batch, [(:) for _ in 1:N-1]..., i)
end

obsslice(batch::AbstractVector, i) = batch[i]

function obsslice(batch::Tuple, i)
    return Tuple(obsslice(batch[j], i) for j in 1:length(batch))
end

function obsslice(batch::NamedTuple, i)
    return (; zip(keys(batch), obsslice(values(batch), i))...)
end

function obsslice(batch::Dict, i)
    return Dict(k => obsslice(v, i) for (k, v) in batch)
end

ext: jl


text:
module LadeGCM

using DifferentialEquations
using Sundials
using CSV
using DataFrames
using Interpolations

"""Houses all constants for the model. See: [`constants`](@ref)."""
struct Constants
    cₐ₀::Float64 #PgC, Pre-industrial atmospheric carbon
    cₜ₀::Float64 #PgC, Pre-industrial soil and vegetation carbon
    cₘ₀::Float64 #PgC, Pre-industrial ocean mixed layer carbon
    λ::Float64 #K, Climate sensitivity (TCR)
    τ::Float64 #yr, Climate lag
    D::Float64 #yr⁻¹, Atmosphere–ocean mixed layer CO₂ equilibration rate
    r::Float64 #Revelle (buffer) factor
    DT::Float64 #% K⁻¹, Solubility temperature effect
    B₀::Float64 #PgC yr⁻¹, Pre-industrial biological pump
    BT::Float64 #% K⁻¹, Temperature dependence of biological pump
    w₀::Float64 #yr⁻¹, Solubility pump rate
    wT::Float64 #% K⁻¹, Weakening of overturning circulation with climate change
    QR::Float64 #Terrestrial respiration temperature dependence
    NPP₀::Float64 #PgC yr⁻¹, Pre-industrial NPP
    KC::Float64 #Fertilisation effect
end

"""
    constants(<keyword arguments>)

A constructor to generate all required constants for the model,
any of which can be overridden via a keyword argument.

# Arguments
- `cₐ₀::Float64=589.0`: PgC, Pre-industrial atmospheric carbon
- `cₜ₀::Float64=1875.0`: PgC, Pre-industrial soil and vegetation carbon
- `cₘ₀::Float64=900.0`: PgC, Pre-industrial ocean mixed layer carbon
- `λ::Float64=1.8`: K, Climate sensitivity (TCR)
- `τ::Float64=4.0`: yr, Climate lag
- `D::Float64=1.0`: yr⁻¹, Atmosphere–ocean mixed layer CO₂ equilibration rate
- `r::Float64=12.5`: Revelle (buffer) factor
- `DT::Float64=0.0423`: 4.23 % K⁻¹, Solubility temperature effect
- `B₀::Float64=13.0`: PgC yr⁻¹, Pre-industrial biological pump
- `BT::Float64=0.032`: 3.2 % K⁻¹, Temperature dependence of biological pump
- `w₀::Float64=0.1`: yr⁻¹, Solubility pump rate
- `wT::Float64=0.1`: 10 % K⁻¹, Weakening of overturning circulation with climate change
- `QR::Float64=1.72`: Terrestrial respiration temperature dependence
- `NPP₀::Float64=55.0`: PgC yr⁻¹, Pre-industrial NPP
- `KC::Float64=0.3`: Fertilisation effect

# Examples

Obtaining the default set of constants

```jldoctest
julia> constants()
LadeGCM.Constants(589.0, 1875.0, 900.0, 1.8, 4.0, 1.0, 12.5, 0.0423, 13.0, 0.032, 0.1, 0.1, 1.72, 55.0, 0.3)
```

Override ``c_{a0}`` and ``KC`` with custom values.

```jldoctest
julia> constants(cₐ₀=600.0, KC=0.5)
LadeGCM.Constants(600.0, 1875.0, 900.0, 1.8, 4.0, 1.0, 12.5, 0.0423, 13.0, 0.032, 0.1, 0.1, 1.72, 55.0, 0.5)
```

"""
function constants(;
    cₐ₀::Float64 = 589., #PgC, Ciais et al. (2013)
    cₜ₀::Float64 = 1875., #PgC, 1325 Pg C of soil organic carbon in top metre of soil Köchy et al. (2015) plus mid range of vegetation carbon estimate by the Ciais et al. (2013)
    cₘ₀::Float64 = 900., #PgC, Ciais et al. (2013)
    λ::Float64 = 1.8, #K, Multi-model mean transient climate response Flato et al. (2013)
    τ::Float64 = 4., #yr, Calculations on ocean heat uptake, the primary cause of climate lag, indicate a response time (e-folding time) of 4 years for timescales up to centuries, before deep ocean heat uptake dominates at millennial timescales Gregory et al. (2015). This result is consistent with simulations that indicate that maximum warming after a CO₂ pulse is reached after only a decade Ricke and Caldeira (2014) and with results fromimpulse response model experiments Joos et al. (2013)
    D::Float64 = 1., #yr⁻¹, Timescale of approximately 1 year, although highly spatially dependent Jones et al. (2014)
    r::Float64 = 12.5, #Williams et al. (2016)
    DT::Float64 = 0.0423, #4.23 % K⁻¹, Takahashi et al. (1993); Ciais et al. (2013, p. 498)
    B₀::Float64 = 13., #PgC yr⁻¹, Ciais et al. (2013)
    BT::Float64 = 0.032, #3.2 % K⁻¹, 12 % decrease Bopp et al., (2013, Fig. 9b) after approximately 3.7 K climate change Collins et al. (2013)
    w₀::Float64 = 0.1, #yr⁻¹, Disolved inorganic carbon flux rate from ocean mixed layer divided by DIC stock in mixed layer Ciais et al. (2013)
    wT::Float64 = 0.1, #10 % K⁻¹, Approximate fit to values reported by Collins et al. (2013, p. 1095)
    QR::Float64 = 1.72, #Raich et al. (2002); Xu and Shang (2016). Based on soil respiration, which contributes the majority of terrestrial ecosystem respiration.
    NPP₀::Float64 = 55., #PgC yr⁻¹, Wieder et al. (2015); Sitch et al. (2015)
    KC::Float64 = 0.3) #Estimated by substituting recent NPP ≈ 60 PgC yr⁻¹ Wieder et al. (2015); Sitch et al. (2015) and recent terrestrial carbon stocks, cₜ ≈ cₜ₀ + 240 Ciais et al. (2013), into Eq. (1). Alexandrov et al. (2003) found that values between 0.3 and 0.4 are compatible with results from a process-based global NPP model.
    Constants(cₐ₀, cₜ₀, cₘ₀, λ, τ, D, r, DT, B₀, BT, w₀, wT, QR, NPP₀, KC)
end

export constants

# Pathways
"""Representation of a *Representative Concentration Pathway*."""
abstract type Pathway end
"""RCP3PD Representation."""
struct RCP3PDPathway <: Pathway end
"""RCP45 Representation."""
struct RCP45Pathway <: Pathway end
"""RCP6 Representation."""
struct RCP6Pathway <: Pathway end
"""RCP85 Representation."""
struct RCP85Pathway <: Pathway end

"""Alias for the [`RCP3PDPathway`](@ref) constructor."""
RCP3PD = RCP3PDPathway()
"""Alias for the [`RCP45Pathway`](@ref) constructor."""
RCP45 = RCP45Pathway()
"""Alias for the [`RCP6Pathway`](@ref) constructor."""
RCP6 = RCP6Pathway()
"""Alias for the [`RCP85Pathway`](@ref) constructor."""
RCP85 = RCP85Pathway()

Base.show(io::IO, s::RCP3PDPathway) = print(io, "+2.6 W/m² (Peak & Decline) Representative Concentration Pathway")
Base.show(io::IO, s::RCP45Pathway) = print(io, "+4.5 W/m² Representative Concentration Pathway")
Base.show(io::IO, s::RCP6Pathway) = print(io, "+6.0 W/m² Representative Concentration Pathway")
Base.show(io::IO, s::RCP85Pathway) = print(io, "+8.5 W/m² Representative Concentration Pathway")
Base.summary(s::RCP3PDPathway) = "RCP3PD"
Base.summary(s::RCP45Pathway) = "RCP45"
Base.summary(s::RCP6Pathway) = "RCP6"
Base.summary(s::RCP85Pathway) = "RCP85"

export RCP3PD, RCP45, RCP6, RCP85

"""
    load_pathway_data(rcp)

Loads csv file for a given RCP scenario into a dataframe for processing.
The information in these files comes from [Meinshausen et al. (2011)](https://doi.org/10.1007/s10584-011-0156-z), which were generated using MAGICC6.
"""
function load_pathway_data(rcp::P) where P <: Pathway
    data = CSV.read(joinpath(@__DIR__, "..", "input", string(summary(rcp), "_EMISSIONS.csv")), header=37, datarow=38, types=[Int64; repeat([Float64]; outer=[39])]);
    tmpnames = names(data);
    tmpnames[1] = :Year;
    names!(data, tmpnames);
end

"""
    (E, LUC) = generate_emission_parameters(rcp)

For a given concentration pathway, generate continuous functions
for fossil fuel emissions `E(t)` and land use emissions `LUC(t)`.

Data comes from files on disk and is linearly interpolated to
provide the continuous output.
"""
function generate_emission_parameters(rcp::P) where P <: Pathway
    data = load_pathway_data(rcp);
    Years = (data[:Year],);
    E = interpolate(Years, data[:FossilCO2], Gridded(Linear()));
    LUC = interpolate(Years, data[:OtherCO2], Gridded(Linear()));
    E, LUC
end

"""
    calculate(rcp; [c=constants(), tspan=(1765., 2100.), solve_args...])

Entry point to the module. Will construct all required parameters and solve the DAE.

The `tspan` here is set to mimic the Lade *et al*. paper, but the emission parameters
in each of the RCP files project out to 2500; meaning the 1765–2500 range is valid here.

Any flags you wish to pass to the solver can be added to the end of the call.

# Examples
```julia
results_6 = calculate(RCP6);
```

```julia
results_45_grid = calculate(RCP45, reltol=1e-10, abstol=1e-10, saveat=1765:0.1:2100);
```

"""
function calculate(rcp::P; #Which scenario are we solving for?
    c::Constants = constants(), #Input conditions.
    tspan = (1765., 2100.), #Time span of calculation.
    solve_args...) where P <: Pathway #Any additional arguments to pass to the solver, like tolerances or timestep alterations.

    #Get emission data
    (E, LUC) = generate_emission_parameters(rcp);

    #Construct initial conditions and parameters for DAE solving
    #             cₜ     cₘ           cₛ          ΔT    cₐ
    u₀ = [c.cₜ₀, c.cₘ₀, c.cₐ₀+c.cₜ₀+c.cₘ₀, 0.0, c.cₐ₀];
    du₀ = similar(u₀);
    diff_vars = [true,true,true,true,false];
    params = [c.NPP₀, c.KC, c.cₐ₀, c.QR, c.cₜ₀, c.D, c.cₘ₀, c.r, c.DT, c.w₀, c.wT, c.B₀, c.BT, c.τ, c.λ];

    """A helper used for solving the DAE."""
    function system(out, du, u, p, t)
        (NPP₀, KC, cₐ₀, QR, cₜ₀, D, cₘ₀, r, DT, w₀, wT, B₀, BT, τ, λ) = p
        cₜ = u[1]
        cₘ = u[2]
        cₛ = u[3]
        ΔT = u[4]
        cₐ = u[5]
        #Helpers
        K(ca,DeltaT) = ((1+KC*log(ca/cₐ₀))/(QR^(DeltaT/10.)))*cₜ₀
        ph(cm,DeltaT) = cₐ₀*(cm/cₘ₀)^r*(1/(1-DT*DeltaT))
        S(cm,DeltaT) = w₀*(1-wT*DeltaT)*(cm-cₘ₀)
        B(DeltaT) = B₀*(1-BT*DeltaT)
        out[1] = (NPP₀/cₜ₀)*QR^(ΔT/10.)*(K(cₐ,ΔT)-cₜ)-LUC(t) - du[1] #cₜ
        out[2] = ((D*cₘ₀)/(r*cₐ₀))*(cₐ-ph(cₘ,ΔT))-S(cₘ,ΔT)-B(ΔT)+B₀ - du[2] #cₘ
        out[3] = E(t) - S(cₘ,ΔT) - (B(ΔT)-B₀) - du[3] #cₛ
        out[4] = (1/τ)*((λ/log(2))*log(cₐ/cₐ₀)-ΔT) - du[4] #ΔT
        out[5] = cₐ + cₜ + cₘ - cₛ #equivalence (cₐ)
    end

    prob = DAEProblem(system,du₀,u₀,tspan,params,differential_vars=diff_vars);
    #Solve the DAE.
    sol = solve(prob,IDA();solve_args...);

    results(sol, c)

end

export calculate

"""Storage for all model output. See [`results`](@ref)."""
struct Results
    cₜ::Array{Float64,1} #PgC, Terrestrial (soil and vegetation) carbon concentration
    cₘ::Array{Float64,1} #PgC, Ocean mixed later carbon concentration
    cₛ::Array{Float64,1} #PgC, System (total) carbon concentration (Atmospheric + Terrestrial + Ocean mixed layer)
    cₐ::Array{Float64,1} #PgC, Atmospheric carbon concentration
    ΔT::Array{Float64,1} #K, Surface temperature relative to per-industrial ΔT ≡ T - T₀
    Δcₜ::Array{Float64,1} #PgC yr⁻¹, Terrestrial carbon stock changes
    Δcₘ::Array{Float64,1} #PgC yr⁻¹, Ocean mixed layer carbon stock changes
    ΔcM::Array{Float64,1} #PgC yr⁻¹, Total ocean carbon stock changes (Mixed layer + deep ocean sink)
    Δcₛ::Array{Float64,1} #PgC yr⁻¹, System carbon stock changes
    Δcₐ::Array{Float64,1} #PgC yr⁻¹, Atmospheric carbon stock changes
    year::Array{Float64,1} #yr, Corresponding time to results
end

"""
    results(sol, c)

Collects all required outputs from the DAE solution.
"""
function results(sol::S, c::Constants) where S <: DiffEqBase.DAESolution
    u = hcat(sol(sol.t, Val{0})...)';
    du = hcat(sol(sol.t, Val{1})...)';

    cₜ = u[:,1];
    cₘ = u[:,2];
    cₛ = u[:,3];
    ΔT = u[:,4];
    cₐ = u[:,5];

    Δcₜ = du[:,1];
    Δcₘ = du[:,2];
    Δcₛ = du[:,3];
    Δcₐ = Δcₛ - Δcₜ - Δcₘ;
    ΔcM = [Δcₘ[i] + sum.(c.w₀*(1-c.wT*ΔT[i])*(cₘ[i]-c.cₘ₀)+c.B₀*(1-c.BT*ΔT[i])-c.B₀) for i in 1:length(sol.t)];

    Results(cₜ, cₘ, cₛ, cₐ, ΔT, Δcₜ, Δcₘ, ΔcM, Δcₛ, Δcₐ, sol.t)
end

end

max_stars_repo_name: Libbum/LadeGCM

max_issues_repo_name: Libbum/LadeGCM

max_forks_repo_name: Libbum/LadeGCM


text:
# Script to run the four formulations -- Bertsimas and Mišić, Belloni et al., McBride and Zufryden, and the Utility formulation --
# in their integer and relaxed forms, for synthetic data, with cardinality constraints (sum x_i <= b).
# This script produces part of the results needed for Tables 5 and 6 of the paper.
#
# The script takes two command line arguments that specify the "rep range". The MIOexpdata1 synthetic instances
# were generated in such a way that for each n and numPermutations combination, there are 100 instances; only the
# first 20 are used in Bertsimas and Mišić. To run the first twenty, one would type:
#
# > julia MIOexpdata1_Comparison_Constraints.jl 1 20
#
# in the terminal.

using JuMP, Gurobi, Stats, Distributions

include("../optimalPLD_code/optimalPLD_createPath.jl")
include("../optimalPLD_code/optimalPLD_solveSAO_constraints.jl")
include("../optimalPLD_code/optimalPLD_solveSAO_BFSS_constraints.jl")
include("../optimalPLD_code/optimalPLD_solveSAO_McBZ_constraints.jl")
include("../optimalPLD_code/optimalPLD_solveSAO_Utility1_constraints.jl")

srand(1204)
machine = 2;

n_vec = [20, 50, 100]
numPerm_vec = [100,200,500,1000]

numReps = 100;

counter = 0;

LR = parse(ARGS[1]);
UR = parse(ARGS[2]);

expdirpath = optimalPLD_createPath(machine, "exp_dir", "MIOexpdata1_constrained_comparison_reps$(LR)to$(UR)");
mkdir(expdirpath)

b_vec = [5; 10]; # test sum x_i <= b for b = 5, 10.

MIP_timeLimit = 1*3600; # Impose a time limit of 1 hour on every MIO formulation
isRelaxation = false;

data_prefix = "MIOexpdata1"

outcsvfilepath = string(expdirpath, "out.csv");
outcsvhandle = open(outcsvfilepath, "w")
row_string = string("n", ",",
					"numPermutations", ",",
					"rep", ",",
					"b_RHS", ",",
					"sumxi", ",",
					"BM_expRevenue", ",",
					"BM_upperBound", ",",
					"BM_solutiontime", ",",
					"BM_MIPgap", ",",
					"BM_LP_expRevenue", ",",
					"BM_LP_upperBound", ",",
					"BM_LP_solutiontime", ",",
					"BM_LP_MIPgap", ",",
					"BM_LP_isIntegral", ",",
					"BFSS_expRevenue", ",",
					"BFSS_upperBound", ",",
					"BFSS_solutiontime", ",",
					"BFSS_MIPgap", ",",
					"BFSS_LP_expRevenue", ",",
					"BFSS_LP_upperBound", ",",
					"BFSS_LP_solutiontime", ",",
					"BFSS_LP_MIPgap", ",",
					"BFSS_LP_isIntegral", ",",
					"Utility1_expRevenue", ",",
					"Utility1_upperBound", ",",
					"Utility1_solutiontime", ",",
					"Utility1_MIPgap", ",",
					"Utility1_LP_expRevenue", ",",
					"Utility1_LP_upperBound", ",",
					"Utility1_LP_solutiontime", ",",
					"Utility1_LP_MIPgap", ",",
					"Utility1_LP_isIntegral", ",",
					"McBZ_expRevenue", ",",
					"McBZ_upperBound", ",",
					"McBZ_solutiontime", ",",
					"McBZ_MIPgap", ",",
					"McBZ_LP_expRevenue", ",",
					"McBZ_LP_upperBound", ",",
					"McBZ_LP_solutiontime", ",",
					"McBZ_LP_MIPgap", ",",
					"McBZ_LP_isIntegral" )
row_string = string(row_string, "\n");

print(outcsvhandle, row_string);
flush(outcsvhandle);

rep_range = LR:UR

counter = 0;
for nind = 1:length(n_vec)
	n = n_vec[nind]
	for nPind = 1:length(numPerm_vec)
		numPermutations = numPerm_vec[nPind];

		println("**** Reading from data directory...")
		data_name = string(data_prefix,"_neq",n,"_Keq",numPermutations)
		data_dir = optimalPLD_createPath(machine, "data_dir", data_name)
		revenues_mat = readcsv(string(data_dir, "revenues_mat.csv"))
		lambda_mat = readcsv(string(data_dir, "lambda_mat.csv"))
		orderings_mat = readcsv(string(data_dir, "orderings_mat.csv"))

		orderings_mat = convert(Array{Int64}, orderings_mat);

		for rep in rep_range
			revenues = revenues_mat[rep,:][:]
			lambda = lambda_mat[rep,:] #[1:numPermutations]
			orderings = orderings_mat[numPermutations*(rep-1)+1:rep*numPermutations,:]

			for b_feas in b_vec

				if (b_feas >= n)
					continue;
				end

				A_feas = ones(1,n);

				x_initial = [];

				println("**** Solving n = $n, nP = $numPermutations, rep = $rep ****")
				println("**** Solving Bertsimas and Misic MIO ****")
				counter = counter + 1;
				isRelaxation = false;
				x_val, y_val, BM_expRevenue, BM_upperBound, BM_solutiontime, BM_MIPgap = optimalPLD_solveSAO_constraints( n, revenues, numPermutations, lambda, orderings, MIP_timeLimit, isRelaxation, x_initial, A_feas, b_feas)

				sumxi = sum(x_val);

				println("**** Solving Bertsimas and Misic LO ****")
				isRelaxation = true;
				BM_LP_x_val, BM_LP_y_val, BM_LP_expRevenue, BM_LP_upperBound, BM_LP_solutiontime, BM_LP_MIPgap = optimalPLD_solveSAO_constraints( n, revenues, numPermutations, lambda, orderings, MIP_timeLimit, isRelaxation, x_initial,  A_feas, b_feas )
				BM_LP_isIntegral = isinteger(BM_LP_x_val) ? 1 : 0;

				println("**** Solving Belloni et al. MIO ****")
				isRelaxation = false;
				x_val, y_val, BFSS_expRevenue, BFSS_upperBound, BFSS_solutiontime, BFSS_MIPgap = optimalPLD_solveSAO_BFSS_constraints( n, revenues, numPermutations, lambda, orderings, MIP_timeLimit, isRelaxation, x_initial, A_feas, b_feas )

				println("**** Solving Belloni et al. LO ****")
				isRelaxation = true;
				BFSS_LP_x_val, BFSS_y_val, BFSS_LP_expRevenue, BFSS_LP_upperBound, BFSS_LP_solutiontime, BFSS_LP_MIPgap = optimalPLD_solveSAO_BFSS_constraints( n, revenues, numPermutations, lambda, orderings, MIP_timeLimit, isRelaxation, x_initial, A_feas, b_feas )
				BFSS_LP_isIntegral = isinteger(BFSS_LP_x_val) ? 1 : 0;

				utilities = zeros(numPermutations, n+1)
				for k=1:numPermutations
					utilities[k,:] = n - (sortperm(orderings[k,:][:]) - 1)
				end

				println("**** Solving McBride & Zufryden MIO ****")
				isRelaxation = false;
				x_val, y_val, McBZ_expRevenue, McBZ_upperBound, McBZ_solutiontime, McBZ_MIPgap = optimalPLD_solveSAO_McBZ_constraints( n, revenues, numPermutations, lambda, utilities, MIP_timeLimit, isRelaxation, x_initial, A_feas, b_feas )

				println("**** Solving McBride Zufryden LO ****")
				isRelaxation = true;
				McBZ_LP_x_val, McBZ_LP_y_val, McBZ_LP_expRevenue, McBZ_LP_upperBound, McBZ_LP_solutiontime, McBZ_LP_MIPgap = optimalPLD_solveSAO_McBZ_constraints( n, revenues, numPermutations, lambda, utilities, MIP_timeLimit, isRelaxation, x_initial, A_feas, b_feas )
				McBZ_LP_isIntegral = isinteger(McBZ_LP_x_val) ? 1 : 0;

				println("**** Solving Utility1 MIO ****")
				isRelaxation = false;
				x_val, y_val, Utility1_expRevenue, Utility1_upperBound, Utility1_solutiontime, Utility1_MIPgap = optimalPLD_solveSAO_Utility1_constraints( n, revenues, numPermutations, lambda, utilities, MIP_timeLimit, isRelaxation, x_initial, A_feas, b_feas )

				println("**** Solving Utility1 LO ****")
				isRelaxation = true;
				Utility1_LP_x_val, Utility1_LP_y_val, Utility1_LP_expRevenue, Utility1_LP_upperBound, Utility1_LP_solutiontime, Utility1_LP_MIPgap = optimalPLD_solveSAO_Utility1_constraints( n, revenues, numPermutations, lambda, utilities, MIP_timeLimit, isRelaxation, x_initial, A_feas, b_feas )
				Utility1_LP_isIntegral = isinteger(Utility1_LP_x_val) ? 1 : 0;

				row_string = string(n, ",",
						numPermutations, ",",
						rep, ",",
						b_feas, ",",
						sumxi, ",",
						BM_expRevenue, ",",
						BM_upperBound, ",",
						BM_solutiontime, ",",
						BM_MIPgap, ",",
						BM_LP_expRevenue, ",",
						BM_LP_upperBound, ",",
						BM_LP_solutiontime, ",",
						BM_LP_MIPgap, ",",
						BM_LP_isIntegral, ",",
						BFSS_expRevenue, ",",
						BFSS_upperBound, ",",
						BFSS_solutiontime, ",",
						BFSS_MIPgap, ",",
						BFSS_LP_expRevenue, ",",
						BFSS_LP_upperBound, ",",
						BFSS_LP_solutiontime, ",",
						BFSS_LP_MIPgap, ",",
						BFSS_LP_isIntegral, ",",
						Utility1_expRevenue, ",",
						Utility1_upperBound, ",",
						Utility1_solutiontime, ",",
						Utility1_MIPgap, ",",
						Utility1_LP_expRevenue, ",",
						Utility1_LP_upperBound, ",",
						Utility1_LP_solutiontime, ",",
						Utility1_LP_MIPgap, ",",
						Utility1_LP_isIntegral, ",",
						McBZ_expRevenue, ",",
						McBZ_upperBound, ",",
						McBZ_solutiontime, ",",
						McBZ_MIPgap, ",",
						McBZ_LP_expRevenue, ",",
						McBZ_LP_upperBound, ",",
						McBZ_LP_solutiontime, ",",
						McBZ_LP_MIPgap, ",",
						McBZ_LP_isIntegral )
				row_string = string(row_string, "\n");

				print(outcsvhandle, row_string);
				flush(outcsvhandle);
			end
		end
	end
end

close(outcsvhandle)

max_stars_repo_name: vvmisic/EFCPLO

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_name: vvmisic/EFCPLO

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_name: vvmisic/EFCPLO

max_forks_repo_licenses: MIT Unlicense


text:
@testset "Generic" begin

    let P = 4, Q = 3, rng = MersenneTwister(123456), N = 100

        # Utility for generating square matrices, vectors, non-square matrices and scalars.
        mPP, mQQ = ()->randn(rng, P, P), ()->randn(rng, Q, Q)
        mPQ, mQP = ()->randn(rng, P, Q), ()->randn(rng, Q, P)
        mPQQP = ()->randn(rng, P * Q, P * Q)
        v, sc = ()->randn(rng, P), ()->randn(rng)
        psd = ()->(A = randn(rng, P, P); transpose(A) * A + 1e-3I)

        # Test all of the unary sensitivities.
        @test check_errs(N, unary_ȲD(-)..., mPQ, mPQ, mPQ)
        @test check_errs(N, unary_ȲD(trace)..., sc, mPP, mPP)
        @test check_errs(N, unary_ȲD(inv)..., mPP, mPP, mPP)
        @test check_errs(N, unary_ȲD(det)..., sc, mPP, mPP)
        @test check_errs(N, unary_ȲD(logdet)..., sc, psd, psd)
        @test check_errs(N, unary_ȲD(transpose)..., mQP, mPQ, mPQ)
        @test check_errs(N, unary_ȲD(ctranspose)..., mQP, mPQ, mPQ)
        @test check_errs(N, unary_ȲD(vecnorm)..., sc, mPQ, mPQ)
        @test check_errs(N, unary_ȲD(vecnorm)..., sc, sc, sc)

        # Test all of the binary sensitivities.
        @test check_errs(N, binary_ȲD(*, 1, mQP)..., mPP, mPQ, mPQ)
        @test check_errs(N, binary_ȲD(*, 2, mPQ)..., mPP, mQP, mQP)
        @test check_errs(N, binary_ȲD(At_mul_B, 1, mQP)..., mPP, mQP, mQP)
        @test check_errs(N, binary_ȲD(At_mul_B, 2, mQP)..., mPP, mQP, mQP)
        @test check_errs(N, binary_ȲD(A_mul_Bt, 1, mPQ)..., mPP, mPQ, mPQ)
        @test check_errs(N, binary_ȲD(A_mul_Bt, 2, mPQ)..., mPP, mPQ, mPQ)
        @test check_errs(N, binary_ȲD(At_mul_Bt, 1, mPQ)..., mPP, mQP, mQP)
        @test check_errs(N, binary_ȲD(At_mul_Bt, 2, mQP)..., mPP, mPQ, mPQ)
        @test check_errs(N, binary_ȲD(Ac_mul_B, 1, mQP)..., mPP, mQP, mQP)
        @test check_errs(N, binary_ȲD(Ac_mul_B, 2, mQP)..., mPP, mQP, mQP)
        @test check_errs(N, binary_ȲD(A_mul_Bc, 1, mPQ)..., mPP, mPQ, mPQ)
        @test check_errs(N, binary_ȲD(A_mul_Bc, 2, mPQ)..., mPP, mPQ, mPQ)
        @test check_errs(N, binary_ȲD(Ac_mul_Bc, 1, mPQ)..., mPP, mQP, mQP)
        @test check_errs(N, binary_ȲD(Ac_mul_Bc, 2, mQP)..., mPP, mPQ, mPQ)
        @test check_errs(N, binary_ȲD(/, 1, mQQ)..., mPQ, mPQ, mPQ)
        @test check_errs(N, binary_ȲD(/, 2, mPQ)..., mPQ, mQQ, mQQ)
        @test check_errs(N, binary_ȲD(At_rdiv_B, 1, mQQ)..., mPQ, mQP, mQP)
        @test check_errs(N, binary_ȲD(At_rdiv_B, 2, mQP)..., mPQ, mQQ, mQQ)
        @test check_errs(N, binary_ȲD(A_rdiv_Bt, 1, mQQ)..., mPQ, mPQ, mPQ)
        @test check_errs(N, binary_ȲD(A_rdiv_Bt, 2, mPQ)..., mPQ, mQQ, mQQ)
        @test check_errs(N, binary_ȲD(At_rdiv_Bt, 1, mQQ)..., mPQ, mQP, mQP)
        @test check_errs(N, binary_ȲD(At_rdiv_Bt, 2, mQP)..., mPQ, mQQ, mQQ)
        @test check_errs(N, binary_ȲD(Ac_rdiv_B, 1, mQQ)..., mPQ, mQP, mQP)
        @test check_errs(N, binary_ȲD(Ac_rdiv_B, 2, mQP)..., mPQ, mQQ, mQQ)
        @test check_errs(N, binary_ȲD(A_rdiv_Bc, 1, mQQ)..., mPQ, mPQ, mPQ)
        @test check_errs(N, binary_ȲD(A_rdiv_Bc, 2, mPQ)..., mPQ, mQQ, mQQ)
        @test check_errs(N, binary_ȲD(Ac_rdiv_Bc, 1, mQQ)..., mPQ, mQP, mQP)
        @test check_errs(N, binary_ȲD(Ac_rdiv_Bc, 2, mQP)..., mPQ, mQQ, mQQ)
        @test check_errs(N, binary_ȲD(\, 1, mQP)..., mQP, mQQ, mQQ)
        @test check_errs(N, binary_ȲD(\, 2, mQQ)..., mQP, mQP, mQP)
        @test check_errs(N, binary_ȲD(At_ldiv_B, 1, mQP)..., mQP, mQQ, mQQ)
        @test check_errs(N, binary_ȲD(At_ldiv_B, 2, mQQ)..., mQP, mQP, mQP)
        @test check_errs(N, binary_ȲD(A_ldiv_Bt, 1, mPQ)..., mQP, mQQ, mQQ)
        @test check_errs(N, binary_ȲD(A_ldiv_Bt, 2, mQQ)..., mQP, mPQ, mPQ)
        @test check_errs(N, binary_ȲD(At_ldiv_Bt, 1, mPQ)..., mQP, mQQ, mQQ)
        @test check_errs(N, binary_ȲD(At_ldiv_Bt, 2, mQQ)..., mQP, mPQ, mPQ)
        @test check_errs(N, binary_ȲD(Ac_ldiv_B, 1, mQP)..., mQP, mQQ, mQQ)
        @test check_errs(N, binary_ȲD(Ac_ldiv_B, 2, mQQ)..., mQP, mQP, mQP)
        @test check_errs(N, binary_ȲD(A_ldiv_Bc, 1, mPQ)..., mQP, mQQ, mQQ)
        @test check_errs(N, binary_ȲD(A_ldiv_Bc, 2, mQQ)..., mQP, mPQ, mPQ)
        @test check_errs(N, binary_ȲD(Ac_ldiv_Bc, 1, mPQ)..., mQP, mQQ, mQQ)
        @test check_errs(N, binary_ȲD(Ac_ldiv_Bc, 2, mQQ)..., mQP, mPQ, mPQ)
        @test check_errs(N, binary_ȲD(vecnorm, 1, sc)..., sc, mPQ, mPQ)
        @test check_errs(N, binary_ȲD(vecnorm, 2, mPQ)..., sc, sc, sc)
        @test check_errs(N, binary_ȲD(vecnorm, 1, sc)..., sc, sc, sc)
        @test check_errs(N, binary_ȲD(vecnorm, 2, sc)..., sc, sc, sc)
        @test check_errs(N, binary_ȲD(kron, 1, mQP)..., mPQQP, mPQ, mPQ)
        @test check_errs(N, binary_ȲD(kron, 2, mPQ)..., mPQQP, mQP, mQP)
        @test check_errs(N, binary_ȲD_inplace(kron, 1, mQP, mPQ())..., mPQQP, mPQ, mPQ)
        @test check_errs(N, binary_ȲD_inplace(kron, 2, mPQ, mQP())..., mPQQP, mQP, mQP)
    end
end


text:
import Dates
println(Dates.now())

include(joinpath(@__DIR__, "run_defs.jl"))

using Distributed
using Random
using Distributions
using BSON
using DelimitedFiles

# Simulation name.
simname = basename("$(@__DIR__)")
println(simname); flush(stdout)

# NOTE: write to scratchdir
function outdir_suffix(pmiss, phi, zind)
  return "pmiss$(pmiss)-phi$(phi)-zind$(zind)"
end

# read Commandline arguments.
RESULTS_DIR_PREFIX = ARGS[1]
AWS_BUCKET_PREFIX = ARGS[2]
PROPMISSINGSCALE = parse(Float64, ARGS[3])
PHI = parse(Float64, ARGS[4])
ZIND = parse(Int, ARGS[5])
ISTEST = parse(Int, ARGS[6]) > 0  # this argument should be 0 or 1 (for test)

# Create settings dictionary.
settings = let
  Nfac = 2000
  temperatures = [1.0, 1.003, 1.006, 1.01]
  Z = let
    path_to_Z = joinpath(@__DIR__, "Z$(ZIND).txt")
    readdlm(path_to_Z, Int, comments=true)
  end
  W = readdlm(joinpath(@__DIR__, "W.txt"), comments=true)
  _outdir_suffix = outdir_suffix(PROPMISSINGSCALE, PHI, ZIND)
  results_dir = "$(RESULTS_DIR_PREFIX)/$(_outdir_suffix)"
  aws_bucket = "$(AWS_BUCKET_PREFIX)/$(_outdir_suffix)"
  if ISTEST
    nsamps = 10
    nburn = 10
  else
    nsamps = 3000
    nburn = 6000
  end

  Dict(:Nfac => Nfac,
       :N => [1, 1] * Nfac,
       :thin_samps => 1,
       :Lmcmc => Dict(0 => 3, 1 => 3),
       :Kmcmc => 10,
       :pthin => 5,
       :dataseed => 1,
       :mcmcseed => 1,
       :temperatures => temperatures,
       :ntemps => length(temperatures),
       :Z => Z,
       :W => W,
       :results_dir => results_dir,
       :aws_bucket => aws_bucket,
       :nsamps => nsamps,
       :nburn => nburn,
       :propmissingscale => PROPMISSINGSCALE,
       :phi => PHI,
       :batchprop => 0.05,
       :dataseed => 1,
       :mcmcseed => 1)
end

# Set number of cores
ncores = settings[:ntemps]
addprocs(ncores)
println("nprocs: $(nprocs())")
println("nworkers: $(nworkers())")
println("system nproc: $(parse(Int, read(`nproc`, String)))")

# Import packages on all workers
@everywhere include("imports.jl")  # Load on all nodes.

println("Starting job ..."); flush(stdout)
println(Dates.now())

sim_and_run(settings)

println("DONE with job!"); flush(stdout)
println(Dates.now())


text:
using {{{PKG}}}
using BenchmarkTools

SUITE = BenchmarkGroup()

# Write your benchmarks here.


text:
using BDF, Compat.Test, HDF5, JLD

#test BDF read and write
origFilePath = "Newtest17-256.bdf"
bdfHeader = readBDFHeader(origFilePath)
dats, evtTab, trigs, statusChan = readBDF(origFilePath)

save("test_data.jld", "bdfHeader", bdfHeader, "EEG", dats,
     "evtTab", evtTab, "trigs", trigs)

headerKeys = collect(keys(bdfHeader))
evtTabKeys = ["idx", "dur", "code"]

h5Out = h5open("Newtest17-256_data.h5", "w")
g = g_create(h5Out, "data")
g["EEG"] = dats
g["trigs"] = trigs
for k in evtTabKeys
    g[k] = evtTab[k]
end
## for k in headerKeys
##     g[k] = bdfHeader[k]
## end
close(h5Out)

#fooData = h5read("Newtest17-256_data.h5", "data")


text:
using SequentialMonteCarlo
using RNGPool
using SMCExamples.MVLinearGaussian
using StaticArrays
import SMCExamples.Particles.MVFloat64Particle
using Test

setRNGs(0)

function mvlgtest(d::Int64, nthreads::Int64)
  model, theta, ys, ko = MVLinearGaussian.defaultMVLGModel(5, d)

  numParticles = 2^16

  smcio = SMCIO{model.particle, model.pScratch}(numParticles, model.maxn,
    nthreads, false)
  smc!(model, smcio)

  @test smcio.logZhats ≈ ko.logZhats atol=0.1*d
end

function testmvlgcsmc(nthreads::Int64, essThreshold::Float64)
  model, theta, ys, ko = MVLinearGaussian.defaultMVLGModel(2, 10)
  nsamples = 2^12

  smcio = SMCIO{model.particle, model.pScratch}(256, model.maxn, nthreads, true,
    essThreshold)

  v = Vector{MVFloat64Particle{2}}(undef, 10)
  for p = 1:10
    v[p] = MVFloat64Particle{2}()
    v[p].x .= zeros(MVector{2, Float64})
  end

  meanEstimates = Vector{MVector{2,Float64}}(undef, 10)
  for p = 1:10
    meanEstimates[p] = zeros(MVector{2, Float64})
  end
  for i = 1:nsamples
    csmc!(model, smcio, v, v)
    for p = 1:10
      meanEstimates[p] .+= v[p].x
    end
  end
  meanEstimates ./= nsamples

  testapproxequal((x -> x[1]).(meanEstimates), (x -> x[1]).(ko.smoothingMeans),
    0.1, false)
  testapproxequal((x -> x[2]).(meanEstimates), (x -> x[2]).(ko.smoothingMeans),
    0.1, false)
end

@time @testset "MV Linear Gaussian: SMC" begin
  mvlgtest(1, 1)
  mvlgtest(1, Threads.nthreads())
  mvlgtest(2, 1)
  mvlgtest(2, Threads.nthreads())
end

@time @testset "MV Linear Gaussian: Serial cSMC" begin
  testmvlgcsmc(1, 2.0)
  testmvlgcsmc(1, 0.5)
end

@time @testset "MV Linear Gaussian: Parallel cSMC" begin
  nt = min(Threads.nthreads(), 4)
  testmvlgcsmc(nt, 2.0)
  testmvlgcsmc(nt, 0.5)
end


text:
@testset "Colormaps" begin
    @test length(colormap("RdBu", 100)) == 100

    # test ability to add to previously shown array of colors - issue #328
    a = [colorant"white", colorant"red"]
    show(stdout, MIME"image/svg+xml"(), a)
    push!(a, colorant"blue")
    @test length(a) == 3
end


text:
using Distributions
using Random
using StaticArrays
using Test
using DocStringExtensions
using LinearAlgebra

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

import CLIMA.DGmethods: boundary_state!
import CLIMA.Atmos: atmos_boundary_state!, atmos_boundary_flux_diffusive!, flux_diffusive!
import CLIMA.DGmethods.NumericalFluxes: boundary_flux_diffusive!

# -------------------- Surface Driven Bubble ----------------- #
# Rising thermals driven by a prescribed surface heat flux.
# 1) Boundary Conditions:
#       Laterally periodic with no flow penetration through top
#       and bottom wall boundaries.
#       Momentum: No flow penetration [NoFluxBC()]
#       Energy:   Spatially varying non-zero heat flux up to time t₁
# 2) Domain: 1250m × 1250m × 1000m
# Configuration defaults are in `src/Driver/Configurations.jl`

"""
  SurfaceDrivenBubbleBC <: BoundaryCondition
Y ≡ state vars
Σ ≡ diffusive vars
A ≡ auxiliary vars
X⁺ and X⁻ refer to exterior, interior faces
X₁ refers to the first interior node

# Fields
$(DocStringExtensions.FIELDS)
"""
struct SurfaceDrivenBubbleBC{FT} <: BoundaryCondition
  "Prescribed MSEF Magnitude `[W/m^2]`"
  F₀::FT
  "Time Cutoff `[s]`"
  t₁::FT
  "Plume wavelength scaling"
  x₀::FT
end
function atmos_boundary_state!(nf::Union{NumericalFluxNonDiffusive,NumericalFluxGradient},
                               bc::SurfaceDrivenBubbleBC,
                               m::AtmosModel,
                               Y⁺::Vars, A⁺::Vars,
                               n⁻,
                               Y⁻::Vars, A⁻::Vars,
                               bctype, t,_...)
  # Use default NoFluxBC()
  atmos_boundary_state!(nf, NoFluxBC(), m,
                        Y⁺, A⁺, n⁻,
                        Y⁻, A⁻,
                        bctype, t)
end
function atmos_boundary_flux_diffusive!(nf::CentralNumericalFluxDiffusive,
                                        bc::SurfaceDrivenBubbleBC,
                                        m::AtmosModel,
                                        F,
                                        Y⁺::Vars, Σ⁺::Vars, A⁺::Vars,
                                        n⁻,
                                        Y⁻::Vars, Σ⁻::Vars, A⁻::Vars,
                                        bctype, t, Y₁⁻, Σ₁⁻, A₁⁻)
  # Working precision
  FT = eltype(Y⁻)
  # Assign vertical unit vector and coordinates
  k̂  = vertical_unit_vector(m.orientation, A⁻)
  x = A⁻.coord[1]
  y = A⁻.coord[2]
  # Unpack fields
  t₁    = bc.t₁
  x₀    = bc.x₀
  F₀ =  t < t₁ ? bc.F₀ : -zero(FT)
  # Apply boundary condition per face (1 == bottom wall)
  if bctype != 1
    atmos_boundary_flux_diffusive!(nf, NoFluxBC(), m, F,
                                   Y⁺, Σ⁺, A⁺,
                                   n⁻,
                                   Y⁻, Σ⁻, A⁻,
                                   bctype, t,
                                   Y₁⁻, Σ₁⁻, A₁⁻)
  else
    atmos_boundary_state!(nf, NoFluxBC(), m,
                          Y⁺, Σ⁺, A⁺,
                          n⁻,
                          Y⁻, Σ⁻, A⁻,
                          bctype, t)
    MSEF    = F₀ * (cospi(2*x/x₀))^2 * (cospi(2*y/x₀))^2
    ∇h_tot⁺ = MSEF * k̂
    _,τ⁺ = turbulence_tensors(m.turbulence, Y⁺, Σ⁺, A⁺, t)
    flux_diffusive!(m, F, Y⁺, τ⁺, ∇h_tot⁺)
  end
end

"""
  Surface Driven Thermal Bubble
"""
function init_surfacebubble!(bl, state, aux, (x,y,z), t)
  FT            = eltype(state)
  R_gas::FT     = R_d
  c_p::FT       = cp_d
  c_v::FT       = cv_d
  γ::FT         = c_p / c_v
  p0::FT        = MSLP

  xc::FT        = 1250
  yc::FT        = 1250
  zc::FT        = 1250
  θ_ref::FT     = 300
  Δθ::FT        = 0

  #Perturbed state:
  θ            = θ_ref + Δθ # potential temperature
  π_exner      = FT(1) - grav / (c_p * θ) * z # exner pressure
  ρ            = p0 / (R_gas * θ) * (π_exner)^ (c_v / R_gas) # density
  P            = p0 * (R_gas * (ρ * θ) / p0) ^(c_p/c_v) # pressure (absolute)
  T            = P / (ρ * R_gas) # temperature
  ρu           = SVector(FT(0),FT(0),FT(0))
  # energy definitions
  e_kin        = FT(0)
  e_pot        = grav * z
  ρe_tot       = ρ * total_energy(e_kin, e_pot, T)
  state.ρ      = ρ
  state.ρu     = ρu
  state.ρe     = ρe_tot
  state.moisture.ρq_tot = FT(0)
end

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

  # Boundary conditions
  # Heat Flux Peak Magnitude
  F₀ = FT(100)
  # Time [s] at which `heater` turns off
  t₁ = FT(500)
  # Plume wavelength scaling
  x₀ = xmax
  bc = SurfaceDrivenBubbleBC{FT}(F₀, t₁, x₀)

  C_smag = FT(0.23)

  imex_solver = CLIMA.DefaultSolverType()
  explicit_solver = CLIMA.ExplicitSolverType(solver_method=LSRK144NiegemannDiehlBusch)

  model = AtmosModel{FT}(AtmosLESConfiguration;
                         turbulence=SmagorinskyLilly{FT}(C_smag),
                         source=(Gravity(),),
                         boundarycondition=bc,
                         moisture=EquilMoist(),
                         init_state=init_surfacebubble!)
  config = CLIMA.Atmos_LES_Configuration("SurfaceDrivenBubble",
                                   N, resolution, xmax, ymax, zmax,
                                   init_surfacebubble!,
                                   solver_type=explicit_solver,
                                   model=model)
  return config
end

function main()
  CLIMA.init()
  FT = Float64
  # DG polynomial order
  N = 4
  # Domain resolution and size
  Δh = FT(50)
  Δv = FT(50)
  resolution = (Δh, Δh, Δv)
  xmax = 2000
  ymax = 2000
  zmax = 2000
  t0 = FT(0)
  timeend = FT(2000)

  CFL_max = FT(0.4)

  driver_config = config_surfacebubble(FT, N, resolution, xmax, ymax, zmax)
  solver_config = CLIMA.setup_solver(t0, timeend, Courant_number=CFL_max, driver_config, forcecpu=true)

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

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

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

main()

max_stars_repo_name: favba/CLIMA

max_issues_repo_name: favba/CLIMA

max_forks_repo_name: favba/CLIMA


text:
import JSON

export
    toquery,
    @slackmethod,
    getresponsetype,
    method_name

isset(t::Any) = true
isset{T}(t::Nullable{T}) = !isnull(t)

getfieldvalue(t::Any) = t
getfieldvalue{T}(t::Nullable{T}) = get(t)

function toquery{T}(t::T)
    query_vars = Dict()
    for name in fieldnames(T)
        field = getfield(t, name)
        if isset(field)
            v = getfieldvalue(field)
            query_vars[string(name)] = v
        end
    end
    query_vars
end

macro slackmethod(req_type::Symbol, method_name::AbstractString, req_block::Expr, resp_block::Expr)
    local resp_type = Symbol(string(req_type) * "Response")

    quote
        immutable $req_type
            token::Token
            $(req_block.args...)
        end

        immutable $resp_type
            $(resp_block.args...)
        end

        $(esc(:getresponsetype))(::Type{$(esc(req_type))}) = $(esc(resp_type))
        $(esc(:method_name))(::Type{$(esc(req_type))}) = $method_name
    end
end


text:
function train_surrogate_model(sModelsProblem::SModelsProblem; verbose::Bool=false, saveToDisk::Bool=false, robust::Bool=true)

  starting_date = date_now()

  #initialization:
  #--------------
  # regression model to predict the continuous output:
  clfr = MLPRegressor(solver="lbfgs", alpha=1e-5, hidden_layer_sizes = (2), random_state=1)
  # classifier to predict convergence of the model:
  clf = MLPClassifier(solver="lbfgs", alpha=1e-5, hidden_layer_sizes = (2), random_state=1)
  # maximum absolute value percentage error regressor:
  max_abs_per_regressor = Inf
  # accuracy score for the classifier:
  accuracy_classifier = Inf
  # size of the train sample:
  size_trainSample = 0

  if verbose == true
    info("Starting date = $(starting_date)")
    info("nworkers = $(nworkers())")
  end

  # Initialization objects on process 1 only:
  #-----------------------------------------
  # To store parameter values
  XDense = Array{Float64}(0,0)
  # To store output values
  YDense = Array{Float64}(0,0)
  # To store convergence status:
  YConvergenceDense = Array{Int64}(0)

  # Generate random draws from the parameter space
  listGrids = []
  for i=1:sModelsProblem.options.nbBatches
    push!(listGrids, create_grid_stochastic(sModelsProblem.lowerBound, sModelsProblem.upperBound, sModelsProblem.options.batchSize, gridType = sModelsProblem.options.gridType))
  end

  counterBatches = 0

  for i=1:sModelsProblem.options.nbBatches

    counterBatches+=1

    if verbose == true
      info("nbBatches = $(i)")
      info("nb function evaluations = $(i*sModelsProblem.options.batchSize)")
    end

    # Objects distributed on the cluster
    #-----------------------------------
    # row dimension = moving along observations
    # column dimension = moving along one observation
    # Model's input:
    #---------------
    # not going to be modified
    #-------------------------
    XDistributed = distribute(listGrids[i], procs = workers(), dist = [nworkers(), 1])

    # Model's output
    #---------------
    # are going to be modified
    #-------------------------
    YDistributed = distribute(zeros(sModelsProblem.options.batchSize, sModelsProblem.dimY),  procs = workers(), dist = [nworkers(), 1])
    YConvergenceDistributed = distribute(zeros(Int64, sModelsProblem.options.batchSize),  procs = workers(), dist = [nworkers()])

    # Split the work on the cluster
    #------------------------------
    refs = Array{Any,1}(nworkers())
    @sync for (wIndex, wID) in enumerate(workers())
      # Each worker modifies its part of XDistributed and YDistributed
      @async refs[wIndex] = @spawnat wID evaluateModel!(sModelsProblem, localpart(XDistributed), localpart(YDistributed), localpart(YConvergenceDistributed))
    end

    # synchronization barrier
    #------------------------
    for (wIndex, wID) in enumerate(workers())
      while isready(refs[wIndex]) == false
        if verbose == true
          info("Waiting for worker $(wID).")
        end
      end
    end

    # Concatenate results
    #---------------------
    # First iteration of the loop is a special case
    if XDense == Array{Float64}(0,0)
      XDense = listGrids[i]
    else
      XDense = vcat(XDense, listGrids[i])
    end

    if YDense == Array{Float64}(0,0)
      YDense = convert(Array, YDistributed)
    else
      YDense = vcat(YDense, convert(Array, YDistributed))
    end

    if YConvergenceDense == Array{Int64}(0)
      YConvergenceDense = convert(Array, YConvergenceDistributed)
    else
      YConvergenceDense = vcat(YConvergenceDense, convert(Array, YConvergenceDistributed))
    end

    # Split between train and test samples
    #--------------------------------------
    XTrain, XTest, yTrain, yTest, yConvergenceTrain, yConvergenceTest = split_train_test(XDense, YDense, YConvergenceDense, testRatio = false)

    # If robust is true, train the regression on parameter values for which
    # convergence was reached
    #---------------------------------------------------------------------------
    # Check that contains at least 2 observations
    if sum(yConvergenceTrain .== sModelsProblem.options.convergenceFlag) > 2
      XTrainRobust = XTrain[yConvergenceTrain .== sModelsProblem.options.convergenceFlag]
      XTestRobust = XTest[yConvergenceTest .== sModelsProblem.options.convergenceFlag]
      yTrainRobust = yTrain[yConvergenceTrain .== sModelsProblem.options.convergenceFlag]
      yTestRobust = yTest[yConvergenceTest .== sModelsProblem.options.convergenceFlag]
    else
      info("No observations for which the model has converged.")
      info("Drawing more observations before training")
      continue
    end

    if robust == false
      size_trainSample = size(XTrain, 1)
    else
      size_trainSample = size(XTrainRobust, 1)
    end

    # Scale the data
    #---------------
    sModelsProblem.scaler = StandardScaler()
    sModelsProblem.scalerRobust = StandardScaler()

    fit!(sModelsProblem.scaler, XTrain)
    fit!(sModelsProblem.scalerRobust, XTrainRobust)

    XTrainScaled = transform(sModelsProblem.scaler, XTrain)
    XTestScaled = transform(sModelsProblem.scaler, XTest)

    XTrainScaledRobust =  transform(sModelsProblem.scalerRobust, XTrainRobust)
    XTestScaledRobust = transform(sModelsProblem.scalerRobust, XTestRobust)

    #---------------------------------------------------------------------------
    #---------------------------------------------------------------------------
    # I.Train another model to predict convergence versus not convergence of the
    # model.
    #---------------------------------------------------------------------------
    #---------------------------------------------------------------------------
    if verbose == true
      info("Training the classifier")
    end

    d = Dict(:hidden_layer_sizes => ((10), (20), (30), (40),
                                    (10, 10), (10, 20), (10, 30), (10, 40),
                                    (20, 10), (20, 20), (20, 30), (20, 40),
                                    (30, 10), (30, 20), (30, 30), (30, 40),
                                    (40, 10), (40, 20), (40, 30), (40, 40),
                                    (10,10,10), (20,20,20), (30,30,30), (40,40,40)))
    bestAccuracyScore = 0
    bestParam = d[:hidden_layer_sizes][1]

    # Loop over hidden layers configuration
    # Looking fot the hidden layer configuration that minimizes the max error in the test set
    for (kIndex, kValue) in enumerate(d[:hidden_layer_sizes])

        println(kValue)
        clf = MLPClassifier(solver="lbfgs", alpha=1e-5, hidden_layer_sizes = kValue, random_state=1)
        fit!(clf, XTrainScaled, yConvergenceTrain)

        yPredicted = round.(Int64, predict(clf, XTestScaled))
        accuracyScore = accuracy_score(yConvergenceTest, yPredicted)

        # Maximize the accuracy:
        if accuracyScore > bestAccuracyScore
            bestAccuracyScore = accuracyScore
            bestParam = kValue
        end

    end

    if verbose == true
      info("Best hyper-parameters for classifier: $(bestParam)")
    end

    # training classifier using the "best" hidden layer
    #-----------------------------------------------
    clf = MLPRegressor(solver="lbfgs", alpha=1e-5, hidden_layer_sizes=bestParam, random_state=1)
    fit!(clf, XTrainScaled, yConvergenceTrain)

    # Train set
    #----------
    yPredicted = round(Int64, predict(clf, XTrainScaled))

    if verbose == true
      info("Classifier")
      info("Accuracy Score Train Set = $(accuracy_score(yConvergenceTrain, yPredicted))")
    end

    # Test set
    #---------
    yPredicted = round(Int64, predict(clf, XTestScaled))
    accuracy_classifier = accuracy_score(yConvergenceTest, yPredicted)

    if verbose == true
      info("Classifier")
      info("Accuracy Score Test Set = $(accuracy_classifier)")
    end

    #---------------------------------------------------------------------------
    #---------------------------------------------------------------------------
    # II. Train the model to predict the continuous output
    # Use "manual" cross-validation on train versus test samples
    #---------------------------------------------------------------------------
    #---------------------------------------------------------------------------
    if verbose == true
      info("Training the regression model")
      info("Robust = $(robust)")
    end

    # If robust, train the model only on

    # TODO: automatic way of doing that + adapt to the number of inputs and output
    # Low tech: There is certainly a better way of doing that
    d = Dict(:hidden_layer_sizes => ((10), (20), (30), (40),
                                    (10, 10), (10, 20), (10, 30), (10, 40),
                                    (20, 10), (20, 20), (20, 30), (20, 40),
                                    (30, 10), (30, 20), (30, 30), (30, 40),
                                    (40, 10), (40, 20), (40, 30), (40, 40),
                                    (10,10,10), (20,20,20), (30,30,30), (40,40,40)))
    bestMaxPerError = Inf
    bestParam = d[:hidden_layer_sizes][1]

    # Loop over hidden layers configuration
    # Looking fot the hidden layer configuration that minimizes the max error in the test set
    for (kIndex, kValue) in enumerate(d[:hidden_layer_sizes])

        println(kValue)
        clfr = MLPRegressor(solver="lbfgs", alpha=1e-5, hidden_layer_sizes = kValue, random_state=1)

        if robust == false

          fit!(clfr, XTrainScaled, yTrain)
          yPredicted = predict(clfr, XTestScaled)
          maxPerErr = calculate_maximum_abs_per_error(yTest, yPredicted)

        # If robust = true, use only the data for which convergence was reached
        #-----------------------------------------------------------------------
        else

          fit!(clfr, XTrainScaledRobust, yTrainRobust)
          yPredicted = predict(clfr, XTestScaledRobust)
          maxPerErr = calculate_maximum_abs_per_error(yTestRobust, yPredicted)

        end

        if maxPerErr < bestMaxPerError
            bestMaxPerError = maxPerErr
            bestParam = kValue
        end

    end

    if verbose == true
      info("Best hyper-parameters for regression model: $(bestParam)")
    end

    clfr = MLPRegressor(solver="lbfgs", alpha=1e-5, hidden_layer_sizes=bestParam, random_state=1)
    # training network using the "best" hidden layer
    #-----------------------------------------------
    if robust == false
      # Train Set:
      #-----------
      fit!(clfr, XTrainScaled, yTrain)
      yPredicted = predict(clfr, XTrainScaled) # Train set
      mean_per_regressorTrain = calculate_mean_per_error(yTest, yPredicted)
      median_per_regressorTrain = calculate_median_abs_per_error(yTest, yPredicted)
      max_abs_per_regressorTrain = calculate_maximum_abs_per_error(yTest, yPredicted)
      # Test Set:
      #-----------
      yPredicted = predict(clfr, XTestScaled)
      mean_per_regressor = calculate_mean_per_error(yTest, yPredicted)
      median_per_regressor = calculate_median_abs_per_error(yTest, yPredicted)
      max_abs_per_regressor = calculate_maximum_abs_per_error(yTest, yPredicted)

    # If robust = true, use only the data for which convergence was reached
    #-----------------------------------------------------------------------
    else
      # Train Set:
      #-----------
      fit!(clfr, XTrainScaledRobust, yTrainRobust)
      yPredicted = predict(clfr, XTrainScaledRobust) # Train set
      mean_per_regressorTrain = calculate_mean_per_error(yTestRobust, yPredicted)
      median_per_regressorTrain = calculate_median_abs_per_error(yTestRobust, yPredicted)
      max_abs_per_regressorTrain = calculate_maximum_abs_per_error(yTestRobust, yPredicted)
      # Test Set:
      #-----------
      yPredicted = predict(clfr, XTestScaledRobust)
      mean_per_regressor = calculate_mean_per_error(yTestRobust, yPredicted)
      median_per_regressor = calculate_median_abs_per_error(yTestRobust, yPredicted)
      max_abs_per_regressor = calculate_maximum_abs_per_error(yTestRobust, yPredicted)
    end

    if verbose == true
      info("Regressor:")
      info("Train set:")
      info("Mean Percentage Error Train Set = $(mean_per_regressorTrain)")
      info("Median Percentage Error Train Set = $(median_per_regressorTrain)")
      info("Maximum Abs Percentage Error Train Set = $(max_abs_per_regressorTrain)")
      info("Test set:")
      info("Mean Percentage Error Test Set = $(mean_per_regressor)")
      info("Median Percentage Error Test Set = $(median_per_regressor)")
      info("Maximum Abs Percentage Error Test Set = $(max_abs_per_regressor)")
      info("Train sample size = $(size_trainSample)")
    end

    # If requested, save the model
    #-----------------------------
    if saveToDisk == true

      if verbose == true
        info("Saving files to disk")
      end

      # regressor:
      JLD.save("clfr_$(starting_date).jld", "clfr", clfr)
      # classifier:
      JLD.save("clf_$(starting_date).jld", "clf", clf)
      # scaler:
      JLD.save("scaler_$(starting_date).jld", "scaler", sModelsProblem.scaler)
      # scaler robust:
      JLD.save("scalerRobust_$(starting_date).jld", "scalerRobust", sModelsProblem.scalerRobust)
      # data:
      JLD.save("XDense_$(starting_date).jld", "XDense", XDense)
      JLD.save("YDense_$(starting_date).jld", "YDense", YDense)
      JLD.save("YConvergenceDense_$(starting_date).jld", "YConvergenceDense", YConvergenceDense)

    end

    #-----------------------------------------
    # Test model accuracy
    # If good enough, stop and save the model
    # Else, move the next batch
    #-----------------------------------------
    # Test on the regressor:
    #-----------------------
    if max_abs_per_regressor < sModelsProblem.options.desiredMaxPerErrorRegressor

      if verbose == true
        info("Current max percentage error regressor = $(max_abs_per_regressor)")
        info("Desired max percentage error regressor = $(sModelsProblem.options.desiredMaxPerErrorRegressor)")
        info("Size of the train sample = $(size_trainSample)")
      end

      # Test on the classifier
      #-----------------------
      if accuracy_classifier < sModelsProblem.options.desiredAccuracyClassifier

        if verbose == true
          info("Current accuracy score classifier = $(accuracy_classifier)")
          info("Desired accuracy score classifier = $(sModelsProblem.options.desiredAccuracyClassifier)")
        end

        # Test on sample size:
        #---------------------
        if size_trainSample > sModelsProblem.options.desiredMinObs

          if verbose == true
            info("Minimum number of points in the train sample reached. Training succesfull.")
          end

          sModelsProblem.trainingSuccessful = true
          break

        end

      end

    end

  end # for i=1:nbBatches

  # If training was not succesfull
  #-------------------------------
  if counterBatches == sModelsProblem.options.nbBatches

    if verbose == true
      info("----------------------------------------------------")
      info("Max number of iterations reached. Training failure.")
      info("Size of the train sample = $(size_trainSample)")
      info("Max abs percentage error regressor = $(max_abs_per_regressor)")
      info("Accuracy score classifier = $(accuracy_classifier)")
      info("----------------------------------------------------")
    end

    sModelsProblem.trainingSuccessful = false
  end

  return clfr, clf

end

function evaluateModel!(sModelsProblem::SModelsProblem, XDistributed::Array, YDistributed::Array, YConvergenceDistributed::Array)

  # Loop over local part of XDistributed
  #-------------------------------------
  for i=1:size(XDistributed, 1)

    # set penalty value by default:
    YDistributed[i,:] = ones(sModelsProblem.dimY)*sModelsProblem.options.penaltyValue
    YConvergenceDistributed[i] = sModelsProblem.options.nonConvergenceFlag

    try
      # if model is evaluated successfully, penalty values are overwritten
      #-------------------------------------------------------------------
      YDistributed[i,:] = sModelsProblem.modelFunction(XDistributed[i,:])
      YConvergenceDistributed[i] = sModelsProblem.options.convergenceFlag
    catch errF
      info("$(errF)")
    end

  end

end


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

module PairwiseListMatrices

  using IndexedArrays
  using Plots
  using Requires

  import Base: size, getindex, setindex!, length, eltype,
  similar, copy, abs, full, zeros, ones, convert,
  transpose, transpose!, ctranspose, ctranspose!, diag, svd,
  triu!, triu, writedlm, writecsv, join,
  -, +, .*, .-, .+, .-, ./, *, /,
  # faster
  mean, sum, varm, var, std, sqrt

  export @iterateupper,
  @iteratelist,
  @iteratediag,

  PairwiseListMatrix, hasdiagonal, lengthlist, ij2k,
  getlist, getdiag, labels, labels!, getlabel, setlabel!,
  sum_nodiag, mean_nodiag, zscore!, zscore,
  from_table, to_table,

  protovis

  include("macros.jl")
  include("pairwiselistmatrix.jl")
  include("protovis.jl")

  @require DataFrames include(joinpath("DataFrames","dataframe.jl"))

end