text:
module PatientAllocationFigures

using DataFrames
using CSV
using OrderedCollections: OrderedDict

using Distributions
using Dates

using Memoize
using ProgressMeter

using Gadfly
using Compose
import Cairo, Fontconfig

include("FiguresBase.jl")
using .FiguresBase

function compute_results(config, data, sent; use_rounding=false)
	N, T = data.N, data.T
	sent_null = zeros(Int, N, N, T)

	if use_rounding
		sent = round.(Int, sent)
	end

	model_title = titlecase(replace(config.experiment, "_" => " ")) * " Model"
	model_title = model_title == "Null Model" ? "No Transfers" : model_title

	active_patients = compute_active(sent, data.initial, data.discharged, data.admitted, data.los_dist, use_rounding=use_rounding)
	active_patients_nosent = compute_active(sent_null, data.initial, data.discharged, data.admitted, data.los_dist, use_rounding=use_rounding)

	overflow = compute_overflow(sent, data.initial, data.discharged, data.admitted, data.beds, data.los_dist, use_rounding=use_rounding)
	overflow_nosent = compute_overflow(sent_null, data.initial, data.discharged, data.admitted, data.beds, data.los_dist, use_rounding=use_rounding)

	load = compute_load(sent, data.initial, data.discharged, data.admitted, data.beds, data.los_dist, use_rounding=use_rounding)
	load_nosent = compute_load(sent_null, data.initial, data.discharged, data.admitted, data.beds, data.los_dist, use_rounding=use_rounding)

	capacity_byday = permutedims(repeat(data.capacity, 1, 1, data.T), (1,3,2))
	overflow_bycapacity = compute_overflow(sent, data.initial, data.discharged, data.admitted, capacity_byday, data.los_dist, use_rounding=use_rounding)
	overflow_bycapacity_nosent = compute_overflow(sent_null, data.initial, data.discharged, data.admitted, capacity_byday, data.los_dist, use_rounding=use_rounding)
	load_bycapacity = compute_load(sent, data.initial, data.discharged, data.admitted, capacity_byday, data.los_dist, use_rounding=use_rounding)
	load_bycapacity_nosent = compute_load(sent_null, data.initial, data.discharged, data.admitted, capacity_byday, data.los_dist, use_rounding=use_rounding)

	results_table = DataFrame(
		date = permutedims(repeat(data.date_range, 1, N), (2,1))[:],
		location = repeat(data.locations, 1, T)[:],
		active_sent = active_patients[:],
		active_nosent = active_patients_nosent[:],
		overflow_sent = overflow[:],
		overflow_nosent = overflow_nosent[:],
		load_sent = load[:],
		load_nosent = load_nosent[:],
	)

	results = (
		sent = sent,
		active = active_patients,
		active_nosent = active_patients_nosent,
		overflow = overflow,
		overflow_nosent = overflow_nosent,
		load = load,
		load_nosent = load_nosent,
		overflow_bycapacity = overflow_bycapacity,
		overflow_bycapacity_nosent = overflow_bycapacity_nosent,
		load_bycapacity = load_bycapacity,
		load_bycapacity_nosent = load_bycapacity_nosent,
		table = results_table,
		model_title = model_title,
	)
	return results
end;

function compute_active(_sent, _initial, _discharged, _admitted, los_dist; use_rounding=false)
	N, _, T = size(_sent)
	_L = 1.0 .- cdf.(los_dist, 0:T)
	_active = [(
		_initial[i]
		- sum(_discharged[i,1:t])
		+ sum(_L[t-t₁+1] * (
			_admitted[i,t₁]
			- sum(_sent[i,:,t₁])
			+ sum(_sent[:,i,t₁])
		) for t₁ in 1:t)
		+ sum(_sent[i,:,t])
	) for i in 1:N, t in 1:T]
	if use_rounding
		_active = round.(Int, _active)
	end
	return _active
end;

function compute_overflow(_sent, _initial, _discharged, _admitted, _cap, los_dist; use_rounding=false)
	return max.(0, compute_active(_sent, _initial, _discharged, _admitted, los_dist, use_rounding=use_rounding) .- _cap)
end;

function compute_load(_sent, _initial, _discharged, _admitted, _cap, los_dist; use_rounding=false)
	return compute_active(_sent, _initial, _discharged, _admitted, los_dist, use_rounding=use_rounding) ./ _cap
end;

function compute_results_block(config, data, sent; use_rounding=false)
	G, N, T = data.G, data.N, data.T
	sent_null = zeros(Int, G, N, N, T)

	if use_rounding
		sent = round.(Int, sent)
	end

	model_title = titlecase(replace(config.experiment, "_" => " ")) * " Model"
	model_title = model_title == "Null Model" ? "No Transfers" : model_title

	active_patients = compute_active_block(sent, data.initial, data.discharged, data.admitted, data.los_bygroup, data.bedtype_bygroup, data.group_transfer_graph, use_rounding=use_rounding)
	active_patients_nosent = compute_active_block(sent_null, data.initial, data.discharged, data.admitted, data.los_bygroup, data.bedtype_bygroup, data.group_transfer_graph, use_rounding=use_rounding)

	overflow = compute_overflow_block(sent, data.initial, data.discharged, data.admitted, data.beds, data.los_bygroup, data.bedtype_bygroup, data.group_transfer_graph, use_rounding=use_rounding)
	overflow_nosent = compute_overflow_block(sent_null, data.initial, data.discharged, data.admitted, data.beds, data.los_bygroup, data.bedtype_bygroup, data.group_transfer_graph, use_rounding=use_rounding)

	load = compute_load_block(sent, data.initial, data.discharged, data.admitted, data.beds, data.los_bygroup, data.bedtype_bygroup, data.group_transfer_graph, use_rounding=use_rounding)
	load_nosent = compute_load_block(sent_null, data.initial, data.discharged, data.admitted, data.beds, data.los_bygroup, data.bedtype_bygroup, data.group_transfer_graph, use_rounding=use_rounding)

	results_table_icu = DataFrame(
		date = permutedims(repeat(data.date_range, 1, N), (2,1))[:],
		location = repeat(data.locations, 1, T)[:],
		active_sent = active_patients.icu[:],
		active_nosent = active_patients_nosent.icu[:],
		overflow_sent = overflow.icu[:],
		overflow_nosent = overflow_nosent.icu[:],
		load_sent = load.icu[:],
		load_nosent = load_nosent.icu[:],
	)
	results_table_ward = DataFrame(
		date = permutedims(repeat(data.date_range, 1, N), (2,1))[:],
		location = repeat(data.locations, 1, T)[:],
		active_sent = active_patients.ward[:],
		active_nosent = active_patients_nosent.ward[:],
		overflow_sent = overflow.ward[:],
		overflow_nosent = overflow_nosent.ward[:],
		load_sent = load.ward[:],
		load_nosent = load_nosent.ward[:],
	)

	results_icu = (
		sent = sent[2,:,:,:],
		active = active_patients.icu,
		active_nosent = active_patients_nosent.icu,
		overflow = overflow.icu,
		overflow_nosent = overflow_nosent.icu,
		load = load.icu,
		load_nosent = load_nosent.icu,
		table = results_table_icu,
		model_title = model_title,
	)
	results_ward = (
		sent = sent[4,:,:,:],
		active = active_patients.ward,
		active_nosent = active_patients_nosent.ward,
		overflow = overflow.ward,
		overflow_nosent = overflow_nosent.ward,
		load = load.ward,
		load_nosent = load_nosent.ward,
		table = results_table_ward,
		model_title = model_title,
	)

	results = (icu = results_icu, ward = results_ward)
	return results
end;

function compute_active_block(_sent, _initial, _discharged, _admitted, los, bed_types, transfer_graph; use_rounding=false)
	G, N, T = size(_admitted)
	B = length(unique(bed_types))

	F = [1.0 .- cdf(los[g], t) for g in 1:G, t in 0:T-1]
	f = [pdf(los[g], t) for g in 1:G, t in 0:T-1]
	f = f ./ sum(f, dims=2)

	groups_bybedtype = [sort(findall(x -> x == b, bed_types)) for b in 1:B]

	@memoize q(g,i,t) = (
		_admitted[g,i,t] + tr(g,i,t) + sum(_sent[g,j,i,t] - _sent[g,i,j,t] for j in 1:N)
	)
	@memoize l(g,i,t) = (
		_discharged[g,i,t] + sum(f[g,t-t₁+1] * q(g,i,t₁) for t₁ in 1:t)
	)
	@memoize tr(g,i,t) = (
		any(transfer_graph[:,g]) ?
			sum(l(g₁,i,t) for g₁ in findall(transfer_graph[:,g]))
			: 0
	)
	active_bygroup = [(
		_initial[g,i] + sum(_sent[g,i,:,t]) + sum((q(g,i,t₁) - l(g,i,t₁)) for t₁ in 1:t)
	) for g in 1:G, i in 1:N, t in 1:T]

	active_bybedtype = [sum(active_bygroup[g,i,t] for g in groups_bybedtype[b]) for b in 1:B, i in 1:N, t in 1:T]

	if use_rounding
		active_bybedtype = round.(Int, active_bybedtype)
	end

	_active = (icu = active_bybedtype[1,:,:], ward = active_bybedtype[2,:,:])
	return _active
end;

function compute_overflow_block(_sent, _initial, _discharged, _admitted, _cap, los, bed_types, transfer_graph; use_rounding=false)
	_active = compute_active_block(_sent, _initial, _discharged, _admitted, los, bed_types, transfer_graph, use_rounding=use_rounding)
	return (
		icu = max.(0, _active.icu .- _cap[1,:]),
		ward = max.(0, _active.ward .- _cap[2,:]),
	)
end;

function compute_load_block(_sent, _initial, _discharged, _admitted, _cap, los, bed_types, transfer_graph; use_rounding=false)
	_active = compute_active_block(_sent, _initial, _discharged, _admitted, los, bed_types, transfer_graph, use_rounding=use_rounding)
	return (
		icu = _active.icu ./ _cap[1,:],
		ward = _active.ward ./ _cap[2,:],
	)
end;

function compute_metrics(config, data, results, bedtype=:all)
	N, T = data.N, data.T

	sent_out = sum(results.sent, dims=2)[:,1,:]
	sent_in = sum(results.sent, dims=1)[1,:,:]
	tfr = sent_out + sent_in

	o = sum(results.overflow[i,t] for i in 1:N, t in 1:T)
	o_null = sum(results.overflow_nosent[i,t] for i in 1:N, t in 1:T)
	o_red = (o_null - o) / o_null

	total_patients = sum(data.initial) + sum(data.admitted)

	nonzero(x) = filter(!=(0), x)

	metrics = [
		("experiment",                                     :config,    :str,  true,   config.experiment),
		("experiment date",                                :config,    :str,  false,  config.rundate),
		("region",                                         :config,    :str,  false,  config.region),
		("allocation level",                               :config,    :str,  false,  config.alloc_level),
		("bed type",                                       :config,    :str,  false,  bedtype),
		("start date",                                     :config,    :dt,   false,  data.start_date),
		("end date",                                       :config,    :dt,   false,  data.end_date),

		("number of days",                                 :data,      :int,  false,  T),
		("number of hospitals",                            :data,      :int,  false,  N),
		("number of hospital-days",                        :data,      :int,  false,  N*T),

		("total patients",                                 :data,      :int,  false,  sum(data.admitted) + sum(data.initial)),
		("total patient-days",                             :data,      :int,  false,  sum(results.active_nosent)),
		("peak active patients",                           :data,      :int,  false,  maximum(sum(results.active_nosent, dims=1))),

		("overflow",                                       :results,   :int,  true,   o),
		("overflow reduction",                             :results,   :pct,  true,   o_red),

		("median non-zero overflow",                       :results,   :int,  true,   iszero(results.overflow) ? 0 : median(nonzero(results.overflow[:]))),
		("mean non-zero overflow",                         :results,   :flt,  true,   iszero(results.overflow) ? 0 : mean(nonzero(results.overflow[:]))),
		("maximum overflow",                               :results,   :int,  true,   maximum(results.overflow)),

		("median load",                                    :results,   :flt,  true,   median(skipbad(results.load[:]))),
		("mean load",                                      :results,   :flt,  true,   mean(skipbad(results.load[:]))),
		("max load",                                       :results,   :flt,  true,   maximum(skipbad(results.load[:]))),

		("number of hospital-days with an overflow",       :results,   :int,  true,   sum(results.overflow .> 0)),
		("percent of hospital-days with an overflow",      :results,   :pct,  true,   mean(results.overflow .> 0)),
		("proportion of hospital-days with an overflow",   :results,   :prop, false,  mean(results.overflow .> 0)),

		("total patients transferred",                     :transfer,  :int,  true,   sum(results.sent)),
		("percent of patients transferred",                :transfer,  :pct,  true,   sum(results.sent) / total_patients),

		("percent of hospital-days with a transfer",       :transfer,  :pct,  true,   mean(tfr .> 0)),
		("percent of days with a transfer",                :transfer,  :pct,  true,   mean(sum(tfr, dims=1) .> 0)),
		("percent of hospitals with a transfer",           :transfer,  :pct,  true,   mean(sum(tfr, dims=2) .> 0)),

		("median non-zero transfer",                       :transfer,  :int,  true,   iszero(results.sent) ? 0 : median(nonzero(results.sent[:]))),
		("mean non-zero transfer",                         :transfer,  :flt,  true,   iszero(results.sent) ? 0 : mean(nonzero(results.sent[:]))),
		("max non-zero transfer",                          :transfer,  :int,  true,   iszero(results.sent) ? 0 : maximum(nonzero(results.sent[:]))),

		("number of transfers",                            :transfer,  :int,  false,  sum(results.sent .> 0)),
		("proportion of patients transferred",             :transfer,  :prop, false,  sum(results.sent) / total_patients),

		("max sent by hospital-day",                       :transfer,  :int,  false,  maximum(results.sent)),
		("max sent by day",                                :transfer,  :int,  false,  maximum(sum(results.sent, dims=[1,2]))),
		("max sent by hospital",                           :transfer,  :int,  false,  maximum(sum(results.sent, dims=[2,3]))),

		("proportion of hospital-days with a transfer",    :transfer,  :prop, false,  mean(tfr .> 0)),
		("proportion of days with a transfer",             :transfer,  :prop, false,  mean(sum(tfr, dims=1) .> 0)),
		("proportion of hospitals with a transfer",        :transfer,  :prop, false,  mean(sum(tfr, dims=2) .> 0)),
	]

	return metrics
end;

function plot_metrics(config, data, results; bedtype=:all, display=true, save=true)
	metrics = compute_metrics(config, data, results, bedtype)

	roundmetric(x) = round(x, digits=1)
	to_percent(x) = string(roundmetric(x*100))*"%"
	format_value(x, m_type) = begin
		if m_type == :pct
			return to_percent(x)
		elseif m_type == :flt
			return string(roundmetric(x))
		elseif m_type == :prop
			return string(round(x, digits=3))
		elseif m_type == :int && (x isa Float64 || x isa Float32)
			return round(Int, x)
		else
			return string(x)
		end
	end
	printmetric(m, x; m_type=:str, m_width=0) = begin
		x = format_value(x, m_type)
		m = m * ":"
		m = (length(m) < m_width+1) ? rpad(m, m_width+1) : m
		println("$(m) $(x)")
	end

	if display
		m_cat_prev = nothing
		m_width = maximum([length(m[1]) for m in metrics]) + 2
		for (m_name, m_cat, m_type, _, m_val) in metrics
			if (m_cat_prev != m_cat) println() end
			m_cat_prev = m_cat
			printmetric(m_name, m_val, m_type=m_type, m_width=m_width)
		end
	end

	all_metrics = DataFrame(
		metric=[m[1] for m in metrics],
		value=[format_value(m[5], m[3]) for m in metrics],
	)
	config_metrics = DataFrame(
		metric=[m[1] for m in metrics if m[2] == :config],
		value=[format_value(m[5], m[3]) for m in metrics if m[2] == :config],
	)
	compare_metrics = DataFrame(
		metric=[m[1] for m in metrics if m[4]],
		value=[format_value(m[5], m[3]) for m in metrics if m[4]],
	)

	if save
		bedtype_ext = (bedtype == :all) ? "" : (bedtype == :icu ? "_icu" : (bedtype == :ward ? "_ward" : ""))
		ext = "$(config.experiment)$(bedtype_ext)"

		metrics_dir = "$(config.results_basepath)/$(config.region_abbrev)/$(config.rundate)/$(config.experiment)/metrics/"
		if !isdir(metrics_dir) mkpath(metrics_dir) end

		figures_dir = "$(config.paperfigures_basepath)/$(config.region_abbrev)/$(config.experiment)/"
		if !isdir(figures_dir) mkpath(figures_dir) end

		all_metrics |> CSV.write(joinpath(metrics_dir, "metrics_all_$(ext).csv"))
		all_metrics_tex = metrics_totextable(all_metrics)
		write(joinpath(metrics_dir, "metrics_all_$(ext).tex"), all_metrics_tex)
		write(joinpath(figures_dir, "metrics_all$(bedtype_ext).tex"), all_metrics_tex)

		compare_metrics |> CSV.write(joinpath(metrics_dir, "metrics_$(ext).csv"))
		compare_metrics_tex = metrics_totextable(compare_metrics, header=true)
		write(joinpath(metrics_dir, "metrics_$(ext).tex"), compare_metrics_tex)
		write(joinpath(figures_dir, "metrics$(bedtype_ext).tex"), compare_metrics_tex)

		config_metrics |> CSV.write(joinpath(metrics_dir, "config_$(ext).csv"))
		config_metrics_tex = metrics_totextable(config_metrics)
		write(joinpath(metrics_dir, "config_$(ext).tex"), config_metrics_tex)
		write(joinpath(figures_dir, "config$(bedtype_ext).tex"), config_metrics_tex)
	end

	return
end;

function plot_estimates_total(config, data; bedtype=:all, display=true, save=true)
	active_computed = compute_active(
		zeros(Int, data.N, data.N, data.T),
		data.initial,
		data.discharged,
		data.admitted,
		data.los_dist,
	)
	colors = Scale.default_discrete_colors(3)[[1,3]]
	plt_active_computed = plot(
		layer(x=data.date_range, y=sum(data.active, dims=1)[:], Geom.point, Geom.line, style(default_color=colors[1])),
		layer(x=data.date_range, y=sum(active_computed, dims=1)[:], Geom.point, Geom.line, style(default_color=colors[2])),
		Guide.xlabel("Date"), Guide.ylabel("Active COVID Patients"), Guide.title("Total Active COVID Patients - Raw vs Computed"),
		Guide.manual_color_key("", ["Data", "Computed"], colors[[1,2]], shape=[Shape.square], size=[1.6mm]),
		Guide.xticks(ticks=DateTime(data.start_date):Day(14):DateTime(data.end_date)),
		Scale.x_continuous(labels=x -> Dates.format(x, dateformat"mm-dd")),
		Scale.y_continuous(format=:plain),
		Coord.cartesian(ymin=0),
		style(
			key_position=:none,
			minor_label_font_size=20px,
			major_label_font_size=25px,
			minor_label_font="CMU Serif",
			major_label_font="CMU Serif",
			key_label_font_size=16px,
			key_title_font_size=18px,
			key_label_font="CMU Serif",
			key_title_font="CMU Serif",
			background_color=colorant"white",
		),
	)
	if display
		plt_active_computed |> SVG(20cm, 10cm)
	end

	plt_admitted_computed = plot(
		layer(x=data.date_range, y=sum(data.admitted, dims=1)[:], Geom.point, Geom.line, style(default_color=colors[1])),
		Guide.xlabel("Date"), Guide.ylabel("Admitted Patients"),
		Guide.title("Total Admitted Patients - Estimated"),
		Guide.xticks(ticks=DateTime(data.start_date):Day(14):DateTime(data.end_date)),
		Scale.x_continuous(labels=x -> Dates.format(x, dateformat"mm-dd")),
		Scale.y_continuous(format=:plain),
		style(
			key_position=:none,
			minor_label_font_size=20px,
			major_label_font_size=25px,
			minor_label_font="CMU Serif",
			major_label_font="CMU Serif",
			background_color=colorant"white",
		),
	)
	plt_discharged_computed = plot(
		layer(x=data.date_range, y=sum(data.discharged, dims=1)[:], Geom.point, Geom.line, style(default_color=colors[2])),
		Guide.xlabel("Date"), Guide.ylabel("Discharged Patients"),
		Guide.title("Total Discharge of Initial Patients - Estimated"),
		Guide.xticks(ticks=DateTime(data.start_date):Day(14):DateTime(data.end_date)),
		Scale.x_continuous(labels=x -> Dates.format(x, dateformat"mm-dd")),
		Scale.y_continuous(format=:plain),
		style(
			key_position=:none,
			minor_label_font_size=20px,
			major_label_font_size=25px,
			minor_label_font="CMU Serif",
			major_label_font="CMU Serif",
			background_color=colorant"white",
		),
	)
	plt_data_computed = hstack(plt_admitted_computed, plt_discharged_computed)
	if display
		plt_data_computed |> SVG(32cm, 10cm)
	end

	bedtype_ext = (bedtype == :all) ? "" : (bedtype == :icu ? "_icu" : (bedtype == :ward ? "_ward" : ""))

	if save
		save_all(plt_active_computed, 20cm, 10cm, "active_computed"*bedtype_ext, config)
		save_all(plt_data_computed, 32cm, 10cm, "adm_dis_computed"*bedtype_ext, config)
	end
end;

function plot_overflow_distribution(config, data, results; bedtype=:all, display=true, save=true)
	colors = Scale.default_discrete_colors(4)[[1,3]]
	binwidth = 5
	overflow_nonzero = filter(x -> x > 0, results.overflow)
	overflow_nosent_nonzero = filter(x -> x > 0, results.overflow_nosent)
	nbins_sent = !isempty(overflow_nonzero) ? ceil(maximum(overflow_nonzero) / binwidth) : 1
	nbins_nosent = !isempty(overflow_nosent_nonzero) ? ceil(maximum(overflow_nosent_nonzero) / binwidth) : 1
	overflow_dist_plot = plot(
		layer(
			x=overflow_nonzero,
			Geom.histogram(
				bincount=nbins_sent,
				limits=(min=0, max=nbins_sent*binwidth)
			),
			style(default_color=colors[1], alphas=[0.7])
		),
		layer(
			x=overflow_nosent_nonzero,
			Geom.histogram(
				bincount=nbins_nosent,
				limits=(min=0, max=nbins_nosent*binwidth)
			),
			style(default_color=colors[2], alphas=[0.7])
		),
		Guide.title("Distribution of (Non-Zero) Overflow by Hospital-Day"),
		Guide.xlabel("Overflow"),
		Guide.ylabel("Count"),
		Guide.manual_color_key("", ["With Transfers", "Without Transfers"], colors, shape=[Shape.square], size=[1.6mm]),
		style(
			minor_label_font_size=18px,
			major_label_font_size=22px,
			minor_label_font="CMU Serif",
			major_label_font="CMU Serif",
			key_label_font_size=16px,
			key_title_font_size=18px,
			key_label_font="CMU Serif",
			key_title_font="CMU Serif",
			background_color=colorant"white",
		),
	)

	bedtype_ext = (bedtype == :all) ? "" : (bedtype == :icu ? "_icu" : (bedtype == :ward ? "_ward" : ""))

	if display
		overflow_dist_plot |> SVG(23cm, 10cm)
	end
	if save
		save_all(overflow_dist_plot, 23cm, 10cm, "overflow_dist"*bedtype_ext, config)
	end
end;

function plot_load_distribution(config, data, results; bedtype=:all, display=true, save=true)
	colors = Scale.default_discrete_colors(4)[[1,3]]
	binwidth = 0.2
	load_sent = skipbad(results.load)
	load_nosent = skipbad(results.load_nosent)
	nbins_sent = ceil(maximum(load_sent) / binwidth)
	nbins_nosent = ceil(maximum(load_nosent) / binwidth)
	load_dist_plot = plot(
		layer(
			x=load_sent,
			Geom.histogram(
				bincount=nbins_sent,
				limits=(min=0, max=nbins_sent*binwidth)
			),
			style(default_color=colors[1], alphas=[0.7])
		),
		layer(
			x=load_nosent,
			Geom.histogram(
				bincount=nbins_nosent,
				limits=(min=0, max=nbins_nosent*binwidth)
			),
			style(default_color=colors[2], alphas=[0.7])
		),
		layer(xintercept=[1.0], Geom.vline(color="red", size=1.0mm), order=10),
		Guide.title("Distribution of Load by Hospital-Day"),
		Guide.xlabel("Load"),
		Guide.ylabel("Count"),
		Guide.manual_color_key("", ["With Transfers", "Without Transfers"], colors, shape=[Shape.square], size=[1.6mm]),
		Coord.Cartesian(xmax=5.0),
		style(
			minor_label_font_size=18px,
			major_label_font_size=22px,
			minor_label_font="CMU Serif",
			major_label_font="CMU Serif",
			key_label_font_size=16px,
			key_title_font_size=18px,
			key_label_font="CMU Serif",
			key_title_font="CMU Serif",
			background_color=colorant"white",
		),
	)

	bedtype_ext = (bedtype == :all) ? "" : (bedtype == :icu ? "_icu" : (bedtype == :ward ? "_ward" : ""))

	if display
		load_dist_plot |> SVG(23cm, 10cm)
	end
	if save
		save_all(load_dist_plot, 23cm, 10cm, "load_dist"*bedtype_ext, config);
	end
end;

function plot_maxload_distribution(config, data, results; bedtype=:all, display=true, save=true)
	colors = Scale.default_discrete_colors(4)[[1,3]]
	binwidth = 0.2
	maxload_sent = skipbad(maximum(results.load, dims=2))
	maxload_nosent = skipbad(maximum(results.load_nosent, dims=2))
	nbins_sent = ceil(maximum(maxload_sent) / binwidth)
	nbins_nosent = ceil(maximum(maxload_nosent) / binwidth)
	maxload_dist_plot = plot(
		layer(
			x=maxload_sent,
			Geom.histogram(
				bincount=nbins_sent,
				limits=(min=0, max=nbins_sent*binwidth)
			),
			style(default_color=colors[1], alphas=[0.7])
		),
		layer(
			x=maxload_nosent,
			Geom.histogram(
				bincount=nbins_nosent,
				limits=(min=0, max=nbins_nosent*binwidth)
			),
			style(default_color=colors[2], alphas=[0.7])
		),
		layer(xintercept=[1.0], Geom.vline(color="red", size=1.0mm), order=10),
		Guide.title("Distribution of Maximum Load by Hospital"),
		Guide.xlabel("Maximum Load"),
		Guide.ylabel("Count"),
		Guide.manual_color_key("", ["With Transfers", "Without Transfers"], colors, shape=[Shape.square], size=[1.6mm]),
		style(
			minor_label_font_size=18px,
			major_label_font_size=22px,
			minor_label_font="CMU Serif",
			major_label_font="CMU Serif",
			key_label_font_size=16px,
			key_title_font_size=18px,
			key_label_font="CMU Serif",
			key_title_font="CMU Serif",
			background_color=colorant"white",
		),
	)

	bedtype_ext = (bedtype == :all) ? "" : (bedtype == :icu ? "_icu" : (bedtype == :ward ? "_ward" : ""))

	if display
		maxload_dist_plot |> SVG(23cm, 10cm)
	end
	if save
		save_all(maxload_dist_plot, 23cm, 10cm, "maxload_dist"*bedtype_ext, config)
	end
end;

function plot_active_total(config, data, results; bedtype=:all, display=true, save=true)
	if bedtype == :icu
		plot_title = "Total Active COVID Patients in ICU Beds in $(config.region)"
		fn = "active_total_icu"
	elseif bedtype == :ward
		plot_title = "Total Active COVID Patients in Ward Beds in $(config.region)"
		fn = "active_total_ward"
	else
		plot_title = "Total Active COVID Patients in $(config.region)"
		fn = "active_total"
	end

	overall_plt = Gadfly.plot(
		layer(
			x = data.date_range,
			y = sum(results.active_nosent, dims=1)[:],
			Geom.point, Geom.line,
		),
		layer(
			yintercept = [sum(data.beds)],
			Geom.hline(color="red", size=0.7mm),
		),
		Guide.xlabel("Date"), Guide.ylabel("Active Patients"),
		Coord.Cartesian(ymin=0),
		Guide.xticks(ticks=DateTime(data.start_date):Day(14):DateTime(data.end_date)),
		Scale.x_continuous(labels=x -> Dates.format(x, dateformat"mm-dd")),
		Scale.y_continuous(format=:plain),
		Guide.title(plot_title),
		style(
			key_position=:none,
			minor_label_font_size=20px,
			major_label_font_size=25px,
			minor_label_font="CMU Serif",
			major_label_font="CMU Serif",
			background_color=colorant"white",
		),
	)

	if display
		overall_plt |> SVG(20cm, 12cm)
	end
	if save
		save_all(overall_plt, 20cm, 12cm, fn, config)
	end
end;

function plot_active(config, data, results; add_title=true, columns=3, display=true, save=true, subset=nothing, bedtype=:all)
	if subset == nothing
		locs = enumerate(data.locations)
	else
		ind = [findfirst(==(h), data.locations) for h in subset]
		locs = zip(ind, subset)
	end

	active_plot_list = []
	colors = Scale.default_discrete_colors(5)[[1,4]]
	for (h_idx, h_name) in locs
		plot_title = length(h_name) > 43 ? strip(h_name[1:40])*"..." : h_name
		plt = plot(
			layer(x=data.date_range, y=results.active_nosent[h_idx,:], Geom.point, Geom.line, style(default_color=colors[1])),
			layer(x=data.date_range, y=results.active[h_idx,:], Geom.point, Geom.line, style(default_color=colors[2])),
			layer(yintercept=[data.beds[h_idx]], Geom.hline(color="red", size=0.7mm)),
			Guide.xticks(ticks=DateTime(data.start_date):Day(14):DateTime(data.end_date)),
			Scale.x_continuous(labels=x -> Dates.format(x, dateformat"mm-dd")),
			Coord.Cartesian(ymin=0),
			Guide.xlabel("Date"), Guide.ylabel("Active Patients"),
			Guide.title(plot_title),
			style(
				key_position=:none,
				minor_label_font_size=16px,
				major_label_font_size=18px,
				minor_label_font="CMU Serif",
				major_label_font="CMU Serif",
				background_color=colorant"white",
			),
		)
		push!(active_plot_list, plt)
	end

	@assert columns == 2 || columns == 3
	plot_width = Dict(2 => 16, 3 => 12)[columns]
	active_plots, active_plots_w, active_plots_h = plot_grid(active_plot_list, columns, plot_width=plot_width, return_plot=true, display=false)

	active_key_colors = vcat(colors, colorant"red")
	active_key = make_horizontal_key(["Without Transfers", "With Transfers", "Capacity"], active_key_colors, height=0.5cm)

	if bedtype == :all
		plot_title = "Active COVID Patients With and Without Transfers"
	elseif bedtype == :icu
		plot_title = "Active COVID Patients in ICU Beds With and Without Transfers"
	elseif bedtype == :ward
		plot_title = "Active COVID Patients in Ward Beds With and Without Transfers"
	end

	active_plots, dy1 = add_key_top(active_key, active_plots)
	active_plots, dy2 = add_title ? add_title_top(plot_title, active_plots) : (active_plots, 0cm)
	active_plots = fill_background(active_plots)

	active_plots_w = active_plots_w
	active_plots_h = active_plots_h + dy1 + dy2

	if display
		active_plots |> SVG(active_plots_w, active_plots_h)
	end

	if save
		out_name = isnothing(subset) ? "active" : "active_subset"
		out_name = (bedtype == :all) ? out_name : (bedtype == :icu ? out_name*"_icu" : (bedtype == :ward ? out_name*"_ward" : out_name))
		save_all(active_plots, active_plots_w, active_plots_h, out_name, config)
	end
end;

function plot_load(config, data, results; key_columns=3, display=true, save=true, subset=nothing, bedtype=:all, multiple_capacity=false)
	if !multiple_capacity
		return _plot_load(config, data, results; key_columns=key_columns, display=display, save=save, subset=subset, bedtype=bedtype, capacity_name=nothing)
	end

	C = size(data.capacity, 2)
	_results_table = deepcopy(results.table)
	for c in 1:C
		cap_abbrev = data.capacity_names_abbrev[c]
		cap_name = data.capacity_names_full[c]

		_results_table[!,:load_sent] = results.load_bycapacity[:,:,c][:]
		_results_table[!,:load_nosent] = results.load_bycapacity_nosent[:,:,c][:]
		_results = merge(results, (table = _results_table,))

		_plot_load(config, data, _results; key_columns=key_columns, display=display, save=save, subset=subset, bedtype=bedtype, capacity_name=cap_abbrev)
	end
end

function _plot_load(config, data, results; key_columns=3, display=true, save=true, subset=nothing, bedtype=:all, capacity_name=nothing)
	if subset == nothing
		results_table = results.table
		locations = data.locations
		nlocations = data.N
	else
		results_table = filter(row -> row.location in subset, results.table)
		locations = subset
		nlocations = length(subset)
	end

	## load plots ##
	load_plot_max = 3.0
	plt_nosent = Gadfly.plot(
		layer(
			results_table, x=:date, y=:load_nosent, color=:location,
			Geom.line,
		),
		Guide.xlabel("Date"), Guide.ylabel("Patient Load"),
		Guide.title("COVID Patient Load in $(config.region) - Without Transfers"),
		layer(ymin=[0.0], ymax=[1.0],           Geom.hband, alpha=[0.1], style(default_color=colorant"green"), order=-1),
		layer(ymin=[1.0], ymax=[load_plot_max], Geom.hband, alpha=[0.1], style(default_color=colorant"red"),   order=-1),
		Coord.cartesian(xmin=data.start_date+Day(1), xmax=data.end_date-Day(1), ymin=0, ymax=load_plot_max),
		Guide.xticks(ticks=DateTime(data.start_date):Day(14):DateTime(data.end_date)),
		Scale.x_continuous(labels=x -> Dates.format(x, dateformat"mm-dd")),
		style(
			key_position=:none,
			minor_label_font_size=20px,
			major_label_font_size=25px,
			minor_label_font="CMU Serif",
			major_label_font="CMU Serif",
		),
	)
	plt_sent = Gadfly.plot(
		layer(
			results_table, x=:date, y=:load_sent, color=:location,
			Geom.line,
		),
		Guide.xlabel("Date"), Guide.ylabel("Patient Load"),
		Guide.title("COVID Patient Load in $(config.region) - With Transfers"),
		layer(ymin=[0.0], ymax=[1.0],           Geom.hband, alpha=[0.1], style(default_color=colorant"green"), order=-1),
		layer(ymin=[1.0], ymax=[load_plot_max], Geom.hband, alpha=[0.1], style(default_color=colorant"red"),   order=-1),
		Coord.cartesian(xmin=data.start_date+Day(1), xmax=data.end_date-Day(1), ymin=0, ymax=load_plot_max),
		Guide.xticks(ticks=DateTime(data.start_date):Day(14):DateTime(data.end_date)),
		Scale.x_continuous(labels=x -> Dates.format(x, dateformat"mm-dd")),
		style(
			key_position=:none,
			minor_label_font_size=20px,
			major_label_font_size=25px,
			minor_label_font="CMU Serif",
			major_label_font="CMU Serif",
		),
	)

	## load key ##
	load_key_labels = map(h_name -> length(h_name) > 53 ? strip(h_name[1:50])*"..." : h_name, locations)
	load_key, load_key_h = make_key(load_key_labels, Scale.default_discrete_colors(nlocations), line_height=0.6cm, ncolumns=key_columns, shape=:circle)

	load_plts_withkey = compose(
		context(0,0,1,1),
		compose(
			context(0, 0, 1w, 1h - load_key_h),
			compose(context(0, 0, 0.5w, 1cy), render(plt_nosent)),
			compose(context(0.5w, 0, 0.5w, 1cy), render(plt_sent)),
		),
		compose(
			context(0, 1h - load_key_h, 1w, load_key_h),
			load_key,
		),
	)

	load_plts_withkey = fill_background(load_plts_withkey)

	if display
		load_plts_withkey |> SVG(40cm, 14cm + load_key_h)
	end

	if save
		out_name = isnothing(subset) ? "load_bylocation" : "load_bylocation_subset"
		out_name = (bedtype == :all) ? out_name : (bedtype == :icu ? out_name*"_icu" : (bedtype == :ward ? out_name*"_ward" : out_name))
		out_name = isnothing(capacity_name) ? out_name : (out_name * "_" * capacity_name)
		save_all(load_plts_withkey, 40cm, 14cm + load_key_h, out_name, config)
	end

	return
end;

function plot_sent_total(config, data, results; bedtype=:all, display=true, save=true)
	sent_plot = plot(
		x=data.date_range,
		y=sum(results.sent, dims=[1,2])[:],
		Geom.point, Geom.line,
		Guide.xticks(ticks=DateTime(data.start_date):Day(14):DateTime(data.end_date)),
		Scale.x_continuous(labels=x -> Dates.format(x, dateformat"mm-dd")),
		Coord.Cartesian(ymin=0),
		Guide.xlabel("Date"), Guide.ylabel("Patient Transfers"),
		Guide.title("Total Patient Transfers"),
		style(
			key_position=:none,
			minor_label_font_size=16px,
			major_label_font_size=18px,
			minor_label_font="CMU Serif",
			major_label_font="CMU Serif",
			background_color=colorant"white",
		),
	)

	bedtype_ext = (bedtype == :all) ? "" : (bedtype == :icu ? "_icu" : (bedtype == :ward ? "_ward" : ""))

	if display
		sent_plot |> SVG(20cm, 10cm)
	end
	if save
		save_all(sent_plot, 20cm, 10cm, "transfers_total"*bedtype_ext, config)
	end
end;

function plot_transfers(config, data, results; add_title=true, columns=3, display=true, save=true, subset=nothing, bedtype=:all)
	if subset == nothing
		locs = collect(enumerate(data.locations))
	else
		ind = [findfirst(==(h), data.locations) for h in subset]
		locs = zip(ind, subset)
	end

	transfer_plot_list = []
	colors = [colorant"seagreen4", colorant"firebrick1"]
	for (h_idx, h_name) in locs
		plot_title = length(h_name) > 43 ? strip(h_name[1:40])*"..." : h_name
		plt = plot(
			layer(x=data.date_range, y=sum(results.sent[:,h_idx,:], dims=1)[:], Geom.point, Geom.line, style(default_color=colors[1])),
			layer(x=data.date_range, y=sum(results.sent[h_idx,:,:], dims=1)[:], Geom.point, Geom.line, style(default_color=colors[2])),
			Guide.xticks(ticks=DateTime(data.start_date):Day(14):DateTime(data.end_date)),
			Scale.x_continuous(labels=x -> Dates.format(x, dateformat"mm-dd")),
			Coord.Cartesian(ymin=0),
			Guide.xlabel("Date"), Guide.ylabel("Patients"),
			Guide.title(plot_title),
			style(
				key_position=:none,
				minor_label_font_size=16px,
				major_label_font_size=18px,
				minor_label_font="CMU Serif",
				major_label_font="CMU Serif",
				background_color=colorant"white",
			),
		)
		push!(transfer_plot_list, plt)
	end

	@assert columns == 2 || columns == 3
	plot_width = Dict(2 => 16, 3 => 12)[columns]
	transfer_plot, plt_w, plt_h = plot_grid(transfer_plot_list, columns, plot_width=plot_width, return_plot=true, display=false)

	transfer_key = make_horizontal_key(["Received", "Sent"], colors, height=0.5cm)

	transfer_plot, dy1 = add_key_top(transfer_key, transfer_plot)
	transfer_plot, dy2 = add_title ? add_title_top("Patients Sent and Received by Hospital", transfer_plot) : (transfer_plot, 0cm)
	transfer_plot = fill_background(transfer_plot)

	transfer_plot_w = plt_w
	transfer_plot_h = plt_h + dy1 + dy2

	if display
		transfer_plot |> SVG(transfer_plot_w, transfer_plot_h)
	end

	if save
		out_name = isnothing(subset) ? "transfers" : "transfers_subset"
		out_name = (bedtype == :all) ? out_name : (bedtype == :icu ? out_name*"_icu" : (bedtype == :ward ? out_name*"_ward" : out_name))
		save_all(transfer_plot, transfer_plot_w, transfer_plot_h, out_name, config)
	end
end;

function sample_random_walk(n=30, maxchange=0.1)
	dist = Uniform(-maxchange,maxchange)
	walk = zeros(Float64, n)
	walk[1] = rand(Uniform(-1,1))
	for i in 2:n
		z = walk[i-1] + rand(dist)
		while z > 1 || z < -1
			z = walk[i-1] + rand(dist)
		end
		walk[i] = z
	end
	return walk
end;

function sample_forecast(data)
	forecast = zeros(Float64, data.N, data.T)
	for i in 1:data.N
		μ = data.admitted[i,:]
		d₁ = data.admitted_uncertainty[i,:]
		d₂ = data.admitted_uncertainty[i,:]

		walk = sample_random_walk(data.T, 0.5)
		patient_walk = μ .+ (max.(0, walk) .* d₁) .+ (min.(0, walk) .* d₂)

		forecast[i,:] = patient_walk
	end
	return forecast
end;

function plot_robust_overflow_distribution(config, data, results; display=true, save=true, bedtype=:all, debug=false)
	sent_null = zeros(Float64, data.N, data.N, data.T)
	_compute_active(_sent, _admitted) = compute_active(_sent, data.initial, data.discharged, _admitted, data.los_dist)
	_compute_total_overflow(_sent, _admitted) = sum(compute_overflow(_sent, data.initial, data.discharged, _admitted, data.beds, data.los_dist))

	n_samples = debug ? 20 : 300
	overflow_sampled_robust = Array{Float64, 1}(undef, n_samples)
	overflow_sampled_nonrobust = Array{Float64, 1}(undef, n_samples)
	overflow_sampled_null = Array{Float64, 1}(undef, n_samples)
	@showprogress for i in 1:n_samples
		_admitted = sample_forecast(data)
		overflow_sampled_robust[i] = _compute_total_overflow(results.sent_robust, _admitted)
		overflow_sampled_nonrobust[i] = _compute_total_overflow(results.sent_nonrobust, _admitted)
		overflow_sampled_null[i] = _compute_total_overflow(sent_null, _admitted)
	end

	binwidth = 100
	nbins_robust = ceil((maximum(overflow_sampled_robust) - minimum(overflow_sampled_robust)) / binwidth)
	nbins_nonrobust = ceil((maximum(overflow_sampled_robust) - minimum(overflow_sampled_robust)) / binwidth)
	nbins_null = ceil((maximum(overflow_sampled_robust) - minimum(overflow_sampled_robust)) / binwidth)

	colors = Scale.default_discrete_colors(4)[[4,1,3]]
	# colors = Scale.default_discrete_colors(4)[[4,2,3]]

	plt = plot(
		layer(
			x=overflow_sampled_robust,
			Geom.histogram(bincount=nbins_robust, density=true),
			style(
				default_color=colors[1],
				alphas=[0.6],
			),
		),
		layer(
			x=overflow_sampled_nonrobust,
			Geom.histogram(bincount=nbins_nonrobust, density=true),
			style(
				default_color=colors[2],
				alphas=[0.6],
			),
		),
		layer(
			x=overflow_sampled_null,
			Geom.histogram(bincount=nbins_null, density=true),
			style(
				default_color=colors[3],
				alphas=[0.6],
			),
		),
		Scale.x_continuous(format=:plain),
		Coord.Cartesian(xmin=0),
		Guide.manual_color_key(
			"Model", ["Robust", "Non-Robust", "No Transfers"], colors,
			shape=[Shape.square], size=[8px],
		),
		Guide.xlabel("Overflow (Patient-Days)"), Guide.ylabel("Frequency"),
		Guide.title("Overflow Distribution by Model"),
		style(
			key_position=:none,
			minor_label_font_size=16px,
			major_label_font_size=18px,
			key_title_font_size=16px,
			key_label_font_size=16px,
			minor_label_font="CMU Serif",
			major_label_font="CMU Serif",
			key_title_font="CMU Serif",
			key_label_font="CMU Serif",
			background_color=colorant"white",
		),
	)

	if display
		plt |> SVG(26cm, 10cm)
	end

	if save
		bedtype_ext = (bedtype == :all) ? "" : (bedtype == :icu ? "_icu" : (bedtype == :ward ? "_ward" : ""))
		save_all(plt, 26cm, 10cm, "overflow_distribution_robust"*bedtype_ext, config)
	end
end;

function plot_active_robust(config, data, results; display=true, save=true, columns=3, add_title=false, bedtype=:all, subset=nothing)
	if subset == nothing
		locs = collect(enumerate(data.locations))
	else
		ind = [findfirst(==(h), data.locations) for h in subset]
		locs = collect(zip(ind, subset))
	end

	sent_null = zeros(Float64, data.N, data.N, data.T)
	_compute_active(_sent, _admitted) = compute_active(_sent, data.initial, data.discharged, _admitted, data.los_dist)

	n_samples = 5
	sampled_admitted = [sample_forecast(data) for _ in 1:n_samples]
	active_sampled_robust = [_compute_active(results.sent_robust, sampled_admitted[i]) for i in 1:n_samples]
	active_sampled_nonrobust = [_compute_active(results.sent_nonrobust, sampled_admitted[i]) for i in 1:n_samples]
	active_sampled_null = [_compute_active(sent_null, sampled_admitted[i]) for i in 1:n_samples]

	colors = Scale.default_discrete_colors(4)[[4,1,3]]
	# colors = Scale.default_discrete_colors(4)[[4,2,3]]

	plots = []
	for (i, h_name) in locs
		sample_layers = []
		for j in 1:n_samples
			push!(sample_layers, layer(
				x = data.date_range,
				y = active_sampled_robust[j][i,:],
				Geom.line,
				style(default_color=colors[1]),
			))
			push!(sample_layers, layer(
				x = data.date_range,
				y = active_sampled_nonrobust[j][i,:],
				Geom.line,
				style(default_color=colors[2]),
			))
			push!(sample_layers, layer(
				x = data.date_range,
				y = active_sampled_null[j][i,:],
				Geom.line,
				style(default_color=colors[3]),
			))
		end

		plot_title = length(h_name) > 43 ? strip(h_name[1:40])*"..." : h_name
		plt = plot(
			sample_layers...,
			layer(yintercept=[data.beds[i]], Geom.hline(color="red", size=0.7mm)),
			Guide.xticks(ticks=DateTime(data.start_date):Day(14):DateTime(data.end_date)),
			Scale.x_continuous(labels=x -> Dates.format(x, dateformat"mm-dd")),
			Coord.Cartesian(ymin=0),
			Guide.xlabel("Date"), Guide.ylabel("Active Patients"),
			Guide.title(plot_title),
			style(
				key_position=:none,
				minor_label_font_size=16px,
				major_label_font_size=18px,
				minor_label_font="CMU Serif",
				major_label_font="CMU Serif",
				background_color=colorant"white",
			),
		)
		push!(plots, plt)
	end

	@assert columns == 2 || columns == 3
	plot_width = Dict(2 => 16, 3 => 12)[columns]
	active_plots, active_plots_w, active_plots_h = plot_grid(plots, columns, plot_width=plot_width, return_plot=true, display=false)

	active_key_colors = vcat(colors, colorant"red")
	active_key = make_horizontal_key(["Robust", "Non-Robust", "No Transfers", "Capacity"], active_key_colors, height=0.5cm)

	if bedtype == :all
		plot_title = "Active COVID Patients By Model"
	elseif bedtype == :icu
		plot_title = "Active COVID Patients in ICU Beds By Model"
	elseif bedtype == :ward
		plot_title = "Active COVID Patients in Ward Beds By Model"
	end

	active_plots, dy1 = add_key_top(active_key, active_plots)
	active_plots, dy2 = add_title ? add_title_top(plot_title, active_plots) : (active_plots, 0cm)
	active_plots = fill_background(active_plots)

	active_plots_w = active_plots_w
	active_plots_h = active_plots_h + dy1 + dy2

	if display
		active_plots |> SVG(active_plots_w, active_plots_h)
	end

	if save
		out_name = isnothing(subset) ? "active_robust" : "active_robust_subset"
		out_name = (bedtype == :all) ? out_name : (bedtype == :icu ? out_name*"_icu" : (bedtype == :ward ? out_name*"_ward" : out_name))
		save_all(active_plots, active_plots_w, active_plots_h, out_name, config)
	end
end;

function plot_active_samples_notransfers(config, data, results; display=true, save=true, columns=3, add_title=true, bedtype=:all, subset=nothing)
	if subset == nothing
		locs = collect(enumerate(data.locations))
	else
		ind = [findfirst(==(h), data.locations) for h in subset]
		locs = collect(zip(ind, subset))
	end

	sent_null = zeros(Float64, data.N, data.N, data.T)
	_compute_active(_admitted) = compute_active(sent_null, data.initial, data.discharged, _admitted, data.los_dist)

	active_null_lb = _compute_active(data.admitted - data.admitted_uncertainty)
	active_null_ub = _compute_active(data.admitted + data.admitted_uncertainty)

	n_samples = 5
	sampled_admitted = [sample_forecast(data) for _ in 1:n_samples-1]
	insert!(sampled_admitted, 1, data.admitted)
	sampled_active = [_compute_active(sampled_admitted[i]) for i in 1:n_samples]

	colors = Scale.default_discrete_colors(n_samples)

	plots = []
	for (i, h_name) in locs
		sample_layers = []
		for j in 1:n_samples
			push!(sample_layers, layer(
				x = data.date_range,
				y = sampled_active[j][i,:],
				Geom.line,
				style(default_color=colors[j]),
				order=10,
			))
		end

		plot_title = length(h_name) > 43 ? strip(h_name[1:40])*"..." : h_name
		plt = plot(
			layer(
				x = data.date_range,
				ymin = active_null_lb[i,:],
				ymax = active_null_ub[i,:],
				Geom.ribbon,
				style(alphas=[0.5]),
			),

			sample_layers...,
			layer(yintercept=[data.beds[i]], Geom.hline(color="red", size=0.7mm)),

			Guide.xticks(ticks=DateTime(data.start_date):Day(14):DateTime(data.end_date)),
			Scale.x_continuous(labels=x -> Dates.format(x, dateformat"mm-dd")),
			Coord.Cartesian(ymin=0),
			Guide.xlabel("Date"), Guide.ylabel("Active Patients"),
			Guide.title(plot_title),
			style(
				key_position=:none,
				minor_label_font_size=16px,
				major_label_font_size=18px,
				minor_label_font="CMU Serif",
				major_label_font="CMU Serif",
				background_color=colorant"white",
			),
		)
		push!(plots, plt)
	end

	@assert columns == 2 || columns == 3
	plot_width = Dict(2 => 16, 3 => 12)[columns]
	active_plots, active_plots_w, active_plots_h = plot_grid(plots, columns, plot_width=plot_width, return_plot=true, display=false)

	if bedtype == :all
		plot_title = "Active COVID Patients Sampled from the Uncertainty Set"
	elseif bedtype == :icu
		plot_title = "Active COVID Patients in ICU Beds Sampled from the Uncertainty Set"
	elseif bedtype == :ward
		plot_title = "Active COVID Patients in Ward Beds Sampled from the Uncertainty Set"
	end

	active_plots, dy1 = add_title ? add_title_top(plot_title, active_plots) : (active_plots, 0cm)
	active_plots = fill_background(active_plots)

	active_plots_w = active_plots_w
	active_plots_h = active_plots_h + dy1

	if display
		active_plots |> SVG(active_plots_w, active_plots_h)
	end

	if save
		out_name = isnothing(subset) ? "active_samples" : "active_samples_subset"
		out_name = (bedtype == :all) ? out_name : (bedtype == :icu ? out_name*"_icu" : (bedtype == :ward ? out_name*"_ward" : out_name))
		save_all(active_plots, active_plots_w, active_plots_h, out_name, config)
	end
end;

function plot_metrics_compare(shared_config, data, results_list; bedtype=:all, display_table=true, save=true)

	roundmetric(x) = round(x, digits=1)
	to_percent(x) = string(roundmetric(x*100))*"%"
	format_value(x, m_type) = begin
		if m_type == :pct
			return to_percent(x)
		elseif m_type == :flt
			return string(roundmetric(x))
		elseif m_type == :prop
			return string(round(x, digits=3))
		elseif m_type == :int && (x isa Float64 || x isa Float32)
			return round(Int, x)
		else
			return string(x)
		end
	end

	metrics_table = DataFrame(metric=[])
	for (exp_name, _results) in results_list
		_config = merge(shared_config, (experiment = exp_name,))
		_metrics = compute_metrics(_config, data, _results)
		_compare_metrics = DataFrame(
			metric=[m[1] for m in _metrics if m[4]],
			value=[format_value(m[5], m[3]) for m in _metrics if m[4]],
		)
		rename!(_compare_metrics, "value" => exp_name)
		metrics_table = outerjoin(metrics_table, _compare_metrics, on=:metric)
	end

	if display_table
		display(metrics_table)
	end

	if save
		bedtype_ext = (bedtype == :all) ? "" : (bedtype == :icu ? "_icu" : (bedtype == :ward ? "_ward" : ""))

		metrics_dir = "$(shared_config.results_basepath)/$(shared_config.region_abbrev)/$(shared_config.rundate)/compare/metrics/"
		if !isdir(metrics_dir) mkpath(metrics_dir) end

		out_path_base = joinpath(metrics_dir, "metrics_compare$(bedtype_ext)")
		metrics_table |> CSV.write(out_path_base * ".csv")
		write(out_path_base * ".tex", metrics_totextable(metrics_table, header=true))
	end

	return
end;

function plot_figures_list(config)
end;

end;

max_stars_repo_name: flixpar/covid-resource-allocation

max_issues_repo_name: flixpar/covid-resource-allocation

max_forks_repo_name: flixpar/covid-resource-allocation


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

mutable struct FeatureVectorizer_InputColumn <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function FeatureVectorizer_InputColumn(; kwargs...)
        obj = new(meta(FeatureVectorizer_InputColumn), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct FeatureVectorizer_InputColumn
const __meta_FeatureVectorizer_InputColumn = Ref{ProtoMeta}()
function meta(::Type{FeatureVectorizer_InputColumn})
    ProtoBuf.metalock() do
        if !isassigned(__meta_FeatureVectorizer_InputColumn)
            __meta_FeatureVectorizer_InputColumn[] = target = ProtoMeta(FeatureVectorizer_InputColumn)
            allflds = Pair{Symbol,Union{Type,String}}[:inputColumn => AbstractString, :inputDimensions => UInt64]
            meta(target, FeatureVectorizer_InputColumn, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_FeatureVectorizer_InputColumn[]
    end
end
function Base.getproperty(obj::FeatureVectorizer_InputColumn, name::Symbol)
    if name === :inputColumn
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :inputDimensions
        return (obj.__protobuf_jl_internal_values[name])::UInt64
    else
        getfield(obj, name)
    end
end

mutable struct FeatureVectorizer <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function FeatureVectorizer(; kwargs...)
        obj = new(meta(FeatureVectorizer), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct FeatureVectorizer
const __meta_FeatureVectorizer = Ref{ProtoMeta}()
function meta(::Type{FeatureVectorizer})
    ProtoBuf.metalock() do
        if !isassigned(__meta_FeatureVectorizer)
            __meta_FeatureVectorizer[] = target = ProtoMeta(FeatureVectorizer)
            allflds = Pair{Symbol,Union{Type,String}}[:inputList => Base.Vector{FeatureVectorizer_InputColumn}]
            meta(target, FeatureVectorizer, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_FeatureVectorizer[]
    end
end
function Base.getproperty(obj::FeatureVectorizer, name::Symbol)
    if name === :inputList
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{FeatureVectorizer_InputColumn}
    else
        getfield(obj, name)
    end
end

export FeatureVectorizer_InputColumn, FeatureVectorizer

ext: jl


text:
push!(LOAD_PATH, joinpath(homedir(),"GitHub/KenyaCoV/src"))
using Plots,Parameters,Distributions,DifferentialEquations,JLD2,DataFrames,StatsPlots,FileIO,MAT,RecursiveArrayTools
using Revise
import KenyaCoV
using LinearAlgebra:eigen
using Statistics: median, quantile

gr()#Plotting frontend

"""
Load age structured data, define initial state and declare the KenyaCoV problem
"""

u0,P,P_dest = KenyaCoV.model_ingredients_from_data("data/data_for_age_structuredmodel.jld2",
                                            "data/flight_numbers.csv",
                                            "data/projected_global_prevelance.csv")

@load "data/detection_rates_for_different_taus.jld2" d_0 d_01 d_025 d_05 d_1

#Initial infecteds

"""
Declare the treatment/isolation rates considered, and a callback that limits intervention
"""
treatment_rates = [(0.,1),(0.,0.5),(1/7.,1.),(1/7,0.5),(1/3.5,1.),(1/3.5,0.5)]

function isolation_limit(u,t,integrator) # Isolation can only continue until 1000 detected cases
  integrator.p.isolating_detecteds && sum(u[:,:,8]) > 1e3
end
function affect_isolation_limit!(integrator)
  integrator.p.τ = 0
  integrator.p.isolating_detecteds = false
end
cb_iso_limit = DiscreteCallback(isolation_limit,affect_isolation_limit!)

@load "data/agemixingmatrix_Kenya_all_types.jld2" M_Kenya M_Kenya_ho M_Kenya_other M_Kenya_school M_Kenya_work
reduced_treatment_rates = [(0.,1),(1/3.5,0.5)]

"""
SOCIAL DISTANCING SCENARIOS
"""

"""
SCENARIO 2 --- SOCIAL DISTANCING:
* No age-specific susceptibilties --- its disease difference
* Asymptomatics are 10% as infectious as symptomatics
"""
u0,P,P_dest = KenyaCoV.model_ingredients_from_data("data/data_for_age_structuredmodel.jld2",
                                            "data/flight_numbers.csv",
                                            "data/projected_global_prevelance.csv")
#Redistribute susceptibility JUST to rescale infectiousness so we get the correct R₀/r
P.χ = ones(KenyaCoV.n_a)
P.rel_detection_rate = d_01
P.dt = 0.25;
P.ext_inf_rate = 0.;
P.ϵ = 0.1
P.γ = 1/2.5
R₀_scale = KenyaCoV.calculate_R₀_scale(P)
P.χ = ones(KenyaCoV.n_a)/R₀_scale

P.σ = 1/rand(KenyaCoV.d_incubation)
P.β = 2.5*P.γ

#Choosing initial conditions - scenario 2
rand_chosen_R,eq_case_distribution,fifth_gen_case_distribution = KenyaCoV.calculate_R₀(P)
num_initial_detected = ceil(Int64,sum(fifth_gen_case_distribution.*d_01))
detected_age_group = argmax(fifth_gen_case_distribution.*d_01)
undetected_age_profile = round.(Int64,fifth_gen_case_distribution.*(1 .- d_01))

#Implement social distancing
KenyaCoV.calculate_R₀(P) #Before
P.M = M_Kenya_ho .+ M_Kenya_work
KenyaCoV.calculate_R₀(P)#after

u0[KenyaCoV.ind_nairobi_as,detected_age_group,4] = num_initial_detected
u0[KenyaCoV.ind_nairobi_as,:,3] = undetected_age_profile
prob = KenyaCoV.create_KenyaCoV_non_neg_prob(u0,(0.,365.),P)

results_2S = KenyaCoV.run_scenario(P,prob,1000,reduced_treatment_rates,cb_iso_limit)
@save joinpath(homedir(),"Github/KenyaCoVOutputs/results_2S.jld2") results_2S

println("Finished 2 social distancing")
"""
SCENARIO 3 --- Social distancing:
* No age-specific susceptibilties --- its disease difference
* Asymptomatics are 25% as infectious as symptomatics
"""
u0,P,P_dest = KenyaCoV.model_ingredients_from_data("data/data_for_age_structuredmodel.jld2",
                                            "data/flight_numbers.csv",
                                            "data/projected_global_prevelance.csv")
#Redistribute susceptibility JUST to rescale infectiousness so we get the correct R₀/r
P.χ = ones(KenyaCoV.n_a)
P.rel_detection_rate = d_025
P.dt = 0.25;
P.ext_inf_rate = 0.;
P.ϵ = 0.25
P.γ = 1/2.5
R₀_scale = KenyaCoV.calculate_R₀_scale(P)
P.χ = ones(KenyaCoV.n_a)/R₀_scale

P.σ = 1/rand(KenyaCoV.d_incubation)
P.β = 2.5*P.γ

#Choosing initial conditions - scenario 3
rand_chosen_R,eq_case_distribution,fifth_gen_case_distribution = KenyaCoV.calculate_R₀(P)
num_initial_detected = ceil(Int64,sum(fifth_gen_case_distribution.*d_025))
detected_age_group = argmax(fifth_gen_case_distribution.*d_025)
undetected_age_profile = round.(Int64,fifth_gen_case_distribution.*(1 .- d_025))

#Implement social distancing
KenyaCoV.calculate_R₀(P)
P.M = M_Kenya_ho .+ M_Kenya_work
KenyaCoV.calculate_R₀(P)

u0[KenyaCoV.ind_nairobi_as,detected_age_group,4] = num_initial_detected
u0[KenyaCoV.ind_nairobi_as,:,3] = undetected_age_profile
prob = KenyaCoV.create_KenyaCoV_non_neg_prob(u0,(0.,365.),P)

results_3S = KenyaCoV.run_scenario(P,prob,1000,reduced_treatment_rates,cb_iso_limit)
@save joinpath(homedir(),"Github/KenyaCoVOutputs/results_3S.jld2") results_3S

println("Finished 3 Social distancing")


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:3] <= 1)

                     @variable(m, -5.992552298143351 <= q <= 4.750344613921084)

                     add_NL_constraint(m, :(log(1 + exp(0.47930083070023954 + -0.2811056889538768*log(1 + exp(-0.30287094621097754 + 0.8508304098923269*$(x[1]) + 0.544182872577263*$(x[2]) + -0.2945971031198069*$(x[3]))) + -0.8918376604733185*log(1 + exp(0.6507439169566887 + 0.9833163152014732*$(x[1]) + -0.36891084266703134*$(x[2]) + -0.5752518327899794*$(x[3]))) + -0.8083246377220759*log(1 + exp(-0.1240101222843979 + -0.304115262711357*$(x[1]) + -0.19393989104697695*$(x[2]) + -0.50508829642924*$(x[3]))))) + log(1 + exp(-0.740550943643274 + -0.007926005977438866*log(1 + exp(-0.30287094621097754 + 0.8508304098923269*$(x[1]) + 0.544182872577263*$(x[2]) + -0.2945971031198069*$(x[3]))) + 0.328386397895855*log(1 + exp(0.6507439169566887 + 0.9833163152014732*$(x[1]) + -0.36891084266703134*$(x[2]) + -0.5752518327899794*$(x[3]))) + 0.6266373732109729*log(1 + exp(-0.1240101222843979 + -0.304115262711357*$(x[1]) + -0.19393989104697695*$(x[2]) + -0.50508829642924*$(x[3]))))) + log(1 + exp(-0.3322291763851517 + -0.3879865017411772*log(1 + exp(-0.30287094621097754 + 0.8508304098923269*$(x[1]) + 0.544182872577263*$(x[2]) + -0.2945971031198069*$(x[3]))) + 0.13951112377587105*log(1 + exp(0.6507439169566887 + 0.9833163152014732*$(x[1]) + -0.36891084266703134*$(x[2]) + -0.5752518327899794*$(x[3]))) + -0.16669374723112762*log(1 + exp(-0.1240101222843979 + -0.304115262711357*$(x[1]) + -0.19393989104697695*$(x[2]) + -0.50508829642924*$(x[3]))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
abstract type KNB <: DataOneV2 end

base_url(::KNB) = "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/"

export ArcticDataCenter, KnowledgeNetworkforBiocomplexity
include("ArcticDataCenter.jl")
include("KnowledgeNetworkforBiocomplexity.jl")

function get_urls(repo::KNB, page)
    urls = []
    links = eachmatch(sel"distribution online", page.root)
    url_links = [replace(text_only(i), "ecogrid://knb/"=>base_url(repo)) for i in links]
    push!(urls, url_links)
    urls
end

function website(repo::KNB, mainpage_url, mainpage)
    replace(mainpage_url, base_url(repo) => "http://dx.doi.org/")
end

function desc_(repo::KNB, mainpage)
    desc_ele = eachmatch(sel"description para", mainpage.root)
    text_only(first(desc_ele))
end

function authors_(repo::KNB, mainpage)
    author_ele = eachmatch(sel"individualName", mainpage.root)
    [text_only(i) for i in author_ele if occursin("creator", string(i.parent))]
end

function pub_date(repo::KNB, mainpage)
    date_ele = eachmatch(sel"pubDate", mainpage.root)
    Dates.Date(text_only(first(date_ele)), "yyyy-mm-dd")
end

function data_fullname(::KNB, mainpage)
    text_only(first(eachmatch(sel"description title", mainpage.root)))
end


text:
# maps:
map((x) -> x * 10, [1, 2, 3]) #=> [10, 20, 30]
map(*, [1, 2, 3], [4, 5, 6]) #=> [4, 10, 18]
cubes = map((x)-> x*x*x, [1:5]) #=> [1, 8, 27, 64, 125]

# do block:
map( x-> begin
           if x == 0 return 0
           elseif iseven(x) return 2
           elseif isodd(x) return 1
           end
         end,
    [-3:3]) #=> [1,2,1,0,1,2,1]

map([-3:3]) do x
    if x == 0 return 0
    elseif iseven(x) return 2
    elseif isodd(x) return 1
    end
end #=> [1,2,1,0,1,2,1]

# filter:
filter( n -> iseven(n), [1:10]) #=> [2, 4, 6, 8, 10]
filter( (n) -> n % 2 == 0, [1:10]) #=> [2, 4, 6, 8, 10]

# comprehensions:
arr = [n for n in 1:1000] #=> [1,2,...,1000]
arr2 = Float64[x^2 for x in 1:4] #=> [1.0,4.0,9.0,16.0]
cubes = [x*x*x for x in [1:5]] #=> [1, 8, 27, 64, 125]
table10 = [x*y for x=1:10, y=1:10] #=>
#=
#1   2   3   4   5   6   7   8   9   10
#2   4   6   8  10  12  14  16  18   20
#3   6   9  12  15  18  21  24  27   30
#4   8  12  16  20  24  28  32  36   40
#5  10  15  20  25  30  35  40  45   50
#6  12  18  24  30  36  42  48  54   60
#7  14  21  28  35  42  49  56  63   70
#8  16  24  32  40  48  56  64  72   80
#9  18  27  36  45  54  63  72  81   90
#10  20  30  40  50  60  70  80  90  100
=#
mat1 = [x + y for x in 1:2, y in 1:3]
#=> 2  3  4
#=> 3  4  5

arrany = { i * 2 for i = 1:5 } #=> [2,4,6,8,10]

max_stars_repo_path: Module 1/Chapter03/Old/map_filters.jl

max_issues_repo_path: Module 1/Chapter03/Old/map_filters.jl

max_forks_repo_path: Module 1/Chapter03/Old/map_filters.jl


text:
# -------  No shocks

# No H shocks
mutable struct HcShockSwitchesNone <: HcShockSwitches
    nColleges :: CollInt
end

mutable struct HcShockSetNone <: HcShockSet
    objId :: ObjectId
    switches :: HcShockSwitchesNone
end

struct HcShockNone <: HcShock end

ModelParams.has_pvector(o :: HcShockSetNone) = false;
shock_probs(o :: HcShockSetNone) = fill(0.0, n_colleges(o));
shock_probs(o :: HcShockSetNone, idx :: Integer) = 0.0;

has_h_shocks(dh :: HcShockSetNone) = false;
has_h_shocks(dh :: HcShockNone) = false;

function validate_hshocks(dh :: HcShockSetNone)
    return true
end

function validate_hshocks(dh :: HcShockNone)
    return true
end

function settings_table(dh :: HcShockSwitchesNone)
    return ["Human capital shocks"  "none"];
end

function make_h_shock_set(objId :: ObjectId, switches :: HcShockSwitchesNone)
    return HcShockSetNone(objId, switches)
end

function make_h_shock(dh :: HcShockSetNone, iCollege :: Integer)
    return HcShockNone()
end

function sim_h_shocks(dh :: HcShockNone, hStartM, hIdxEndM, hGridNextV, rng)
    return copy(hIdxEndM)
end

function pr_hprime(dh :: HcShockNone, hEndIdx :: Integer,
    hStart :: Real, hEnd :: Real, hGridNextV)
    idxV = [hEndIdx];
    prV = [1.0];
    return idxV, prV
end

# """
# 	$(SIGNATURES)

# Given value `vNextV` at start of `t+1` by `h[t+1]`, find value at end of `t` before shocks are realized.
# """
# value_before_h_shock(dh :: HcShockNone, hStart :: F1, hEnd :: F1,
#     hGridNextV :: AbstractVector{F1},  vNextV :: AbstractVector{F1}) where F1 =
#     copy(vNextV);

# ----------------


text:
global rng = srand()
global my_precision = 1e-17

function set_seed(seed)
    srand(rng, seed)
end

function set_precision(prec)
    my_precision = prec
end


text:
# Note: some of CodeTracking's functionality can only be tested by Revise

using CodeTracking
using Test, InteractiveUtils, REPL, LinearAlgebra, SparseArrays
# Note: ColorTypes needs to be installed, but note the intentional absence of `using ColorTypes`

using CodeTracking: line_is_decl

if !isempty(ARGS) && "revise" ∈ ARGS
    # For running tests with and without Revise
    using Revise
end

isdefined(Main, :Revise) ? Main.Revise.includet("script.jl") : include("script.jl")

@testset "CodeTracking.jl" begin
    m = first(methods(f1))
    file, line = whereis(m)
    scriptpath = normpath(joinpath(@__DIR__, "script.jl"))
    @test file == scriptpath
    @test line == (line_is_decl ? 2 : 4)
    trace = try
        call_throws()
    catch
        stacktrace(catch_backtrace())
    end
    @test whereis(trace[2]) == (scriptpath, 9)

    src, line = definition(String, m)
    @test src == chomp("""
    function f1(x, y)
        # A comment
        return x + y
    end
    """)
    @test line == 2
    @test code_string(f1, Tuple{Any,Any}) == src
    @test @code_string(f1(1, 2)) == src

    m = first(methods(f2))
    src, line = definition(String, m)
    @test src == "f2(x, y) = x + y"
    @test line == 14

    m = first(methods(throws))
    src, line = definition(String, m)
    @test startswith(src, "@noinline")
    @test line == 7

    m = first(methods(multilinesig))
    src, line = definition(String, m)
    @test startswith(src, "@inline")
    @test line == 16
    @test @code_string(multilinesig(1, "hi")) == src
    @test_throws ErrorException("no unique matching method found for the specified argument types") @code_string(multilinesig(1, 2))

    m = first(methods(f50))
    src, line = definition(String, m)
    @test occursin("100x", src)
    @test line == 22

    # Issue #81
    m = which(hasrettype, (Int,))
    src, line = definition(String, m)
    @test occursin("Float32", src)
    @test line == 43

    info = CodeTracking.PkgFiles(Base.PkgId(CodeTracking))
    @test Base.PkgId(info) === info.id
    @test CodeTracking.basedir(info) == dirname(@__DIR__)

    io = IOBuffer()
    show(io, info)
    str = String(take!(io))
    @test startswith(str, "PkgFiles(CodeTracking [da1fd8a2-8d9e-5ec2-8556-3022fb5608a2]):\n  basedir:")
    ioctx = IOContext(io, :compact=>true)
    show(ioctx, info)
    str = String(take!(io))
    @test match(r"PkgFiles\(CodeTracking, .*CodeTracking(\.jl)?, String\[\]\)", str) !== nothing

    @test pkgfiles("ColorTypes") === nothing
    @test_throws ErrorException pkgfiles("NotAPkg")

    # Test a method marked as missing
    m = @which sum(1:5)
    CodeTracking.method_info[m.sig] = missing
    @test whereis(m) == (CodeTracking.maybe_fix_path(String(m.file)), m.line)
    @test definition(m) === nothing

    # Test that definitions at the REPL work with `whereis`
    ex = Base.parse_input_line("replfunc(x) = 1"; filename="REPL[1]")
    eval(ex)
    m = first(methods(replfunc))
    @test whereis(m) == ("REPL[1]", 1)
    # Test with broken lookup
    oldlookup = CodeTracking.method_lookup_callback[]
    CodeTracking.method_lookup_callback[] = m -> error("oops")
    @test whereis(m) == ("REPL[1]", 1)
    CodeTracking.method_lookup_callback[] = oldlookup

    m = first(methods(Test.eval))
    @test occursin(Sys.STDLIB, whereis(m)[1])

    # https://github.com/JuliaDebug/JuliaInterpreter.jl/issues/150
    function f150()
        x = 1 + 1
        @info "hello"
    end
    m = first(methods(f150))
    src = Base.uncompressed_ast(m)
    idx = findfirst(lin -> String(lin.file) == @__FILE__, src.linetable)
    lin = src.linetable[idx]
    file, line = whereis(lin, m)
    @test endswith(file, String(lin.file))
    idx = findfirst(lin -> String(lin.file) != @__FILE__, src.linetable)
    lin = src.linetable[idx]
    file, line = whereis(lin, m)
    @test endswith(file, String(lin.file))

    # Issues raised in #48
    m = @which(sum([1]; dims=1))
    if !isdefined(Main, :Revise)
        def = definition(String, m)
        @test def === nothing || isa(def[1], AbstractString)
        def = definition(Expr, m)
        @test def === nothing || isa(def, Expr)
    else
        def = definition(String, m)
        @test isa(def[1], AbstractString)
        def = definition(Expr, m)
        @test isa(def, Expr)
    end

    # Issue #64
    B = Hermitian(hcat([one(BigFloat) + im]))
    m = @which cholesky(B)
    @test startswith(definition(String, m)[1], "cholesky")

    # Ensure that we don't error on difficult cases
    m = which(+, (AbstractSparseVector, AbstractSparseVector))  # defined inside an `@eval`
    d = definition(String, m)
    @test d === nothing || isa(d[1], AbstractString)

    # Check for existence of file
    id = Base.PkgId("__PackagePrecompilationStatementModule")   # not all Julia versions have this
    mod = try Base.root_module(id) catch nothing end
    if isa(mod, Module)
        m = first(methods(getfield(mod, :eval)))
        @test definition(String, m) === nothing
    end
end

@testset "With Revise" begin
    if isdefined(Main, :Revise)
        m = @which gcd(10, 20)
        sigs = signatures_at(Base.find_source_file(String(m.file)), m.line)
        @test !isempty(sigs)
        ex = @code_expr(gcd(10, 20))
        @test ex isa Expr
        body = ex.args[2]
        idx = findfirst(x -> isa(x, LineNumberNode), body.args)
        @test occursin(String(m.file), String(body.args[idx].file))
        @test ex == code_expr(gcd, Tuple{Int,Int})

        m = first(methods(edit))
        sigs = signatures_at(String(m.file), m.line)
        @test !isempty(sigs)
        sigs = signatures_at(Base.find_source_file(String(m.file)), m.line)
        @test !isempty(sigs)

        # issue #23
        @test !isempty(signatures_at("script.jl", 9))

        @test_throws ArgumentError signatures_at("nofile.jl", 1)

        if isdefined(Revise, :add_revise_deps)
            Revise.add_revise_deps()
            sigs = signatures_at(CodeTracking, "src/utils.jl", 5)
            @test length(sigs) == 1       # only isn't available on julia 1.0
            @test first(sigs) == Tuple{typeof(CodeTracking.checkname), Expr, Any}
            @test pkgfiles(CodeTracking).id == Base.PkgId(CodeTracking)
        end

        # REPL (test copied from Revise)
        if isdefined(Base, :active_repl)
            hp = Base.active_repl.interface.modes[1].hist
            fstr = "__fREPL__(x::Int16) = 0"
            histidx = length(hp.history) + 1 - hp.start_idx
            ex = Base.parse_input_line(fstr; filename="REPL[$histidx]")
            f = Core.eval(Main, ex)
            if ex.head === :toplevel
                ex = ex.args[end]
            end
            push!(hp.history, fstr)
            m = first(methods(f))
            @test definition(String, first(methods(f))) == (fstr, 1)
            @test !isempty(signatures_at(String(m.file), m.line))
            pop!(hp.history)
        elseif haskey(ENV, "CI")
            error("CI Revise tests must be run with -i")
        end
    end
end

(a_34)(x::T, y::T) where {T<:Integer} = no_op_err("&", T)
(b_34)(x::T, y::T) where {T<:Integer} = no_op_err("|", T)
c_34(x::T, y::T) where {T<:Integer} = no_op_err("xor", T)

(d_34)(x::T, y::T) where {T<:Number} = x === y
(e_34)(x::T, y::T) where {T<:Real} = no_op_err("<" , T)
(f_34)(x::T, y::T) where {T<:Real} = no_op_err("<=", T)
l = @__LINE__
@testset "#34 last character" begin
    def, line = definition(String, @which d_34(1, 2))
    @test line == l - 3
    @test def == "(d_34)(x::T, y::T) where {T<:Number} = x === y"
end

function g()
    Base.@_inline_meta
    print("hello")
end
@testset "inline macros" begin
    def, line = CodeTracking.definition(String, @which g())
    @test def == """
    function g()
        Base.@_inline_meta
        print("hello")
    end"""
end

@testset "kwargs methods" begin
    m = nothing
    for i in 1:30
        s = Symbol("#func_2nd_kwarg#$i")
        if isdefined(Main, s)
            m = @eval $s
        end
    end
    m === nothing && error("couldn't find keyword function")
    body, loc = CodeTracking.definition(String, first(methods(m)))
    @test loc == 28
    @test body == "func_2nd_kwarg(; kw=2) = true"
end

@testset "method extensions" begin
    body, _ = CodeTracking.definition(String, @which Foo.Bar.fit(1))
    @test body == """
    function Foo.Bar.fit(m)
        return m
    end"""
    body, _ = CodeTracking.definition(String, @which Foo.Bar.fit(1, 2))
    @test body == "Foo.Bar.fit(a, b) = a + b"
end


text:
using GaussianDistributions
using Test

include("gaussian.jl")
include("correct.jl")
include("bivariate.jl")


text:
"""
`Ensemble{TP, TF<:AbstractFloat, TI<:Integer, TD}`: A particle ensemble structure
designed for WE with bins.

### Fields

* `ξ̂` - particle positions after selection, before mutation
* `ξ` - particle positions after mutation
* `ω̂` - partice weights after selection, before mutation
* `ω` - partice weights after mutation
* `b̂` - particle bin after selection, before mutation
* `b` - particle bin after mutation
* `o` - number of offspring of the particle
* `d̂` - auxiliary data for each particle after selection, before mutation
* `d` - auxiliary data for each particle, after mutation
"""
struct Ensemble{TP, TF<:AbstractFloat, TI<:Integer, TD} <: AbstractEnsemble
   # positions of the walkers after selection, before mutation
   ξ̂::Vector{TP}
   # positions of the walkers after mutation
   ξ::Vector{TP}
   # weights of the walkers after selection, before mutation
   ω̂::Vector{TF}
   # weights of the walkers after mutation
   ω::Vector{TF}
   # category ("bin") type of each of the walkers after selection, before mutation
   b̂::Vector{TI}
   # category ("bin") type of each of the walkers after mutation
   b::Vector{TI}
   # number of offspring of each particle
   o::Vector{TI}
   # auxiliary data for each particle after selection, before mutation
   d̂::Vector{TD}
   # auxiliary data for each particle after mutation
   d::Vector{TD}
end

function Base.eltype(E::TE) where {TE<:Ensemble}
   return typeof((E.ξ̂[1], E.ξ[1], E.ω̂[1], E.ω[1],E.b̂[1], E.b[1],E.o[1], E.d̂[1], E.d[1]))
end

function Base.push!(E::TE, ξ̂, ξ, ω̂, ω, b̂, b, o, d̂, d) where {TE<:Ensemble}
   push!(E.ξ̂, ξ̂);
   push!(E.ξ, ξ);
   push!(E.ω̂, ω̂);
   push!(E.ω, ω);
   push!(E.b̂, b̂);
   push!(E.b, b);
   push!(E.o, o);
   push!(E.d̂, d̂);
   push!(E.d, d);
end

function Base.pop!(E::TE) where {TE<:Ensemble}
   ξ̂ = pop!(E.ξ̂);
   ξ = pop!(E.ξ);
   ω̂ = pop!(E.ω̂);
   ω = pop!(E.ω);
   b̂ = pop!(E.b̂);
   b = pop!(E.b);
   o = pop!(E.o);
   d̂ = pop!(E.d̂);
   d = pop!(E.d);
   return ξ̂, ξ, ω̂, ω, b̂, b, o, d̂, d
end

function Base.popfirst!(E::TE) where {TE<:Ensemble}
   ξ̂ = popfirst!(E.ξ̂);
   ξ = popfirst!(E.ξ);
   ω̂ = popfirst!(E.ω̂);
   ω = popfirst!(E.ω);
   b̂ = popfirst!(E.b̂);
   b = popfirst!(E.b);
   o = popfirst!(E.o);
   d̂ = popfirst!(E.d̂);
   d = popfirst!(E.d);
   return ξ̂, ξ, ω̂, ω, b̂, b, o, d̂, d
end

function Base.iterate(E::TE, state = 1) where {TE<:Ensemble}

   if state > length(E)
      return nothing
   end
   return (E.ξ̂[state],E.ξ[state], E.ω̂[state], E.ω[state], E.b̂[state], E.b[state], E.o[state], E.d̂[state], E.d[state]), state+1
end

# common ensmble functions

function Base.length(E::TE) where {TE<:AbstractEnsemble}
   return length(E.ξ)
end

function Base.isempty(E::TE) where {TE<:AbstractEnsemble}
   return isempty(E.ξ)
end


text:
module FenicsPy

using PyCall

export @pydef, @py_str

# must be explicitly imported to be extended
import Base: split, inv, transpose, div, diff, abs, sign, sqrt, exp,  cos, sin, tan, acos, asin, atan,  cosh, sinh, tanh, log, replace, adjoint, *, +, -, /, ^, ==, <<

import PyPlot: plot

##################################
#    Base : FeObject
##################################

abstract type
  FeObject
end

##################################
#    input
##################################

FeType = Union{Real, String, Array, Tuple, Dict, Function, FeObject}

to_pyarray(arr::Array) = [isa(a, FeObject) ? a.pyobject : a for a in arr]

to_pytuple(arr::Tuple) = Tuple(isa(a, FeObject) ? a.pyobject : a for a in arr)

to_pydict(arr::Dict) = Dict(i => (isa(arr[i], FeObject) ? arr[i].pyobject : arr[i]) for i in eachindex(arr))

function to_pytype(obj::FeType)
    if isa(obj, FeObject)
        obj.pyobject
    elseif isa(obj, Array)
        to_pyarray(obj)
    elseif isa(obj, Tuple)
        to_pytuple(obj)
    elseif isa(obj, Dict)
        to_pydict(obj)
    else
        obj
    end
end

function args_conv(args...; kwargs...)
    _args = [to_pytype(arg) for arg in args]
    _kwargs = Dict(idx =>  to_pytype(kwargs[idx]) for idx in eachindex(kwargs))
    _args, _kwargs
end

##################################
#    output
##################################

PyType = Union{Real, String, Array, Tuple, Dict, Function, PyObject, Nothing}

function to_feobject(obj::PyObject)
    clsname = match(r"(\w+)'>", string(obj.__class__))[1]
    clsname = get(class_dict, clsname, :FeAny)
    getproperty(FenicsPy, clsname)(obj)
end

to_fearray(arr::Array) = [isa(a, PyObject) ? to_feobject(a) : a for a in arr]

to_fetuple(arr::Tuple) = Tuple(isa(a, PyObject) ? to_feobject(a) : a for a in arr)

to_fedict(arr::Dict) = Dict(i => (isa(arr[i], PyObject) ? to_feobject(arr[i]) : arr[i]) for i in eachindex(arr))

function to_fetype(obj::PyType)
    if isa(obj, PyObject)
        to_feobject(obj)
    elseif isa(obj, Array)
        to_fearray(obj)
    elseif isa(obj, Tuple)
        to_fetuple(obj)
    elseif isa(obj, Dict)
        to_fedict(obj)
    else
        obj
    end
end

##################################
#    @pyfunc
#
# examples :
#   1  @pyfunc dolfin TestFunction
#   2  @pyfunc mshr generate_mesh  you_alias
##################################
macro pyfunc(_module::Symbol, name::Symbol, alias::Symbol=:nothing)
    alias = alias === :nothing ? name : alias
    name = string(name)
    esc(quote
        $(alias)(args::FeType...; kwargs...) = begin
            _args, _kwargs = args_conv(args...; kwargs...)
            obj = getproperty($_module, $name)(_args...; _kwargs...)
	    to_fetype(obj)
        end
    end)
end

export pyfunc

##################################
#    @pyclass
#
# examples :
#    1  @pyclass dolfin Point
#    2  @pyclass mshr Rectangle CSGGeometry
#    3  @pyclass dolfin Function FeObject FeFunction
##################################
const class_dict = Dict()
macro pyclass(_module::Symbol, name::Symbol, _base::Symbol=:FeObject, alias::Symbol=:Nothing)
    alias = alias === :Nothing ? name : alias
    impl = Symbol(alias, "Impl")
    name = string(name)
    class_dict[name] = alias
    esc(quote
        abstract type
            $alias <: $_base
        end

        struct $impl <: $alias
            pyobject::PyObject
        end

        ###############################
        #    Constructors ( copy )
        #
        # examples :
        #     dx = Measure(dolfin.dx)
        ###############################
        $(alias)(pyobject::PyObject) = $impl(pyobject)

        ###############################
        #    Constructors ( parameter )
        #
        # examples :
        #     p = Point(1.2, 2.3)
        ###############################
        $(alias)(args::FeType...; kwargs...) = begin
            _args, _kwargs = args_conv(args...; kwargs...)
            name = hasproperty($_module, $name) ? $name : $_base
            _obj = getproperty($_module, name)(_args...; _kwargs...)
            $impl(_obj)
        end

        ###############################
        # examples :
        #     p.x()    # = 1.2
        #     p.pyobject
        ###############################
        function Base.getproperty(_obj::$impl, _sym::Symbol)
            if _sym === :pyobject
	        return getfield(_obj, _sym)
            end
            o = getproperty(_obj.pyobject, _sym)
            if isa(o, PyObject) && match(r"(\w+)'>", string(o.__class__))[1] == "method"
                return (args::FeType...; kwargs...) ->  begin
                       _args, _kwargs = args_conv(args...; kwargs...)
		       to_fetype(o(_args...; _kwargs...))
                end
            end
            to_fetype(o)
        end

        function Base.setproperty!(_obj::$impl, _sym::Symbol, value)
            if _sym === :pyobject
                return
            end
            setproperty!(_obj.pyobject, _sym, to_pytype(value))
        end

        Base.getindex(o::$impl, key::Union{String, Symbol}) = to_fetype(get(o.pyobject, key))
        Base.getindex(o::$impl, idx::Int) = to_fetype(get(o.pyobject, idx-1))
        Base.setindex!(o::$impl, value, key::Union{String, Symbol}) = set!(o.pyobject, key, to_pytype(value))
        Base.setindex!(o::$impl, value, idx::Int) = set!(o.pyobject, idx-1, to_pytype(value))

        ###############################
        # https://docs.julialang.org/en/v1/manual/methods/#Function-like-objects
        #
        # examples :
        #    dx = Measure("dx", domain=mesh, subdomain_data=markers)
        #    dx(1)
        ###############################
        function (o::$impl)(args...; kwargs...)
            _args, _kwargs = args_conv(args...; kwargs...)
            to_fetype(o.pyobject(_args...; _kwargs...))
	end

        ###############################
        # examples :
        #     CellType.Type.quadrilateral
        #     UnitSquareMesh.create
        ###############################
        function Base.getproperty(_obj::Type{$alias}, _sym::Symbol)
            if hasproperty(_obj, _sym)
                return getfield(_obj, _sym)
            end
            name = hasproperty($_module, $name) ? $name : $_base
            pyclass = getproperty($_module, name)
            o = getproperty(pyclass, _sym)
            if isa(o, PyObject) && match(r"(\w+)'>", string(o.__class__))[1] == "method"
                return (args::FeType...; kwargs...) ->  begin
                       _args, _kwargs = args_conv(args...; kwargs...)
                       to_fetype(o(_args...; _kwargs...))
                end
            end
            to_fetype(o)
        end

    end)
end

export pyclass

##################################
#    python module : dolfin,  mshr, ufl
##################################

include("common.jl")
include("la.jl")
include("function.jl")
include("fem.jl")
include("mesh.jl")
include("geometry.jl")
include("operator.jl")
include("mshr.jl")
include("io.jl")

##################################
#    init
##################################

const dolfin = PyCall.PyNULL()
const ufl = PyCall.PyNULL()
const mshr = PyCall.PyNULL()

export dolfin, ufl, mshr

function __init__()

    copy!(dolfin, pyimport_conda("dolfin", "fenics-dolfin", "conda-forge"))
    copy!(mshr, pyimport_conda("mshr", "mshr", "conda-forge"))
    copy!(ufl, pyimport_conda("ufl", "fenics-ufl", "conda-forge"))

    global dx = Measure(dolfin.dx)
    global ds = Measure(dolfin.ds)
    global dS = Measure(dolfin.dS)
    global dP = Measure(dolfin.dP)

    global tetrahedron = Cell(dolfin.tetrahedron)
    global hexahedron = Cell(dolfin.hexahedron) #matplotlib cannot handle hexahedron elements
    global triangle = Cell(dolfin.triangle)
    global quadrilateral = Cell(dolfin.quadrilateral)

    # parameters["linear_algebra_backend"] = backendname  # list_linear_algebra_backends()
    global parameters = GlobalParameters(dolfin.parameters)

    global DOLFIN_EPS = dolfin.DOLFIN_EPS
end

export dx, ds, dS, dP,
       tetrahedron, hexahedron, triangle, quadrilateral,
       parameters,
       DOLFIN_EPS

end # FenicsPy


text:
module Booster

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Booster (Green)" => Ahorn.EntityPlacement(
        Maple.GreenBooster
    ),
    "Booster (Red)" => Ahorn.EntityPlacement(
        Maple.RedBooster
    )
)

function boosterSprite(entity::Maple.Booster)
    red = get(entity.data, "red", false)

    if red
        return "objects/booster/boosterRed00"

    else
        return "objects/booster/booster00"
    end
end

function Ahorn.selection(entity::Maple.Booster)
    x, y = Ahorn.position(entity)
    sprite = boosterSprite(entity)

    return Ahorn.getSpriteRectangle(sprite, x, y)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::Maple.Booster, room::Maple.Room)
    sprite = boosterSprite(entity)

    Ahorn.drawSprite(ctx, sprite, 0, 0)
end

end


text:
#=
    -------------------------------------------------------------------------
    Implements functionalities for setting up the Markov chain Monte Carlo
    algorithm. The main object is `MCMCSetup` and its members comprise of
    all objects that need to be passed to the `mcmc` function from `mcmc.jl`.
    The remaining routines are used to populate instances of `MCMCSetup` and
    verify its fields in a structured manner.
    --------------------------------------------------------------------------
=#

#===============================================================================
                                MCMC setup
===============================================================================#
"""
    mutable struct MCMCSetup

Implements functionalities for setting up the Markov chain Monte Carlo
algorithm. The main object is `MCMCSetup` and its members comprise of
all objects that need to be passed to the `mcmc` function from `mcmc.jl`.
The remaining routines are used to populate instances of `MCMCSetup` and
verify its fields in a structured manner.
"""
mutable struct MCMCSetup
    updates
    #readjust_params
    # add parameters for fusion, correlation matrix, historical acceptance rates
    function MCMCSetup(updates...)
        new([updates...])
    end
end

function set_readjust_params!(setup::MCMCSetup, readjust_params)
    setup.readjust_params = readjust_params
end

#===============================================================================
                                Diffusion setup
===============================================================================#

"""
    DiffusionSetup{ObsScheme} <: ModelSetup

Setup choices relevant to the path augmentation step to be passed to `mcmc` function from
`mcmc.jl`.

# Example:
```
         DiffusionSetup(P, P̃)
         DiffusionSetup(P, P̃, fptOrPartObs)
```
where P is the target process and P̃ is the linear approximation and possibly a vector of Booleans whether
the observations are first passage times.
"""
mutable struct DiffusionSetup{ObsScheme} <: ModelSetup
    setup_completion    # (internal) Indicates progress of setting-up DiffusionSetup
    P˟                  # Target diffusion law
    P̃                   # Vector of auxiliary diffusion laws
    adaptive_prop       # Object for adapting guided proposals
    skip_for_save       # When saving path, thin the grid by a factor of ...
    Ls                  # Vector of observation operators
    Σs                  # Vector of covariance matrices
    obs                 # Observations
    obs_times           # Recorded times of observations
    fpt                 # Vector with information about first-passage times
    dt                  # Granularity of the imputation grid
    τ                   # Time transformation for the imputation grid
    x0_prior            # Prior over the starting position
    x0_guess            # Guess for a position of a starting point
    Wnr                 # Definition of the driving Wiener law
    XX                  # Container for the path
    WW                  # Container for the driving noise
    P                   # Container for the guided proposals

    """
        DiffusionSetup(P˟, P̃, ::ObsScheme)

    Initialise `DiffusionSetup` with a given target law `P˟`, auxiliary laws `P̃` and
    observation scheme `ObsScheme`.
    """
    function DiffusionSetup(P˟, P̃, ::ObsScheme) where ObsScheme <: AbstractObsScheme
        new{ObsScheme}(Dict(:obs => false,
                            :imput => false,
                            :prior => false),
                       P˟, P̃, NoAdaptation(), 1)#NoBlocking(), ([], 0.1, NoChangePt()),
    end
end

adaptation_object(setup::DiffusionSetup) = deepcopy(setup.adaptive_prop)

"""
    set_observations!(setup::DiffusionSetup, Ls, Σs, obs, obs_times,
                      fpt=fill(nothing, length(obs)-1))

Store observations to `setup`. The observations follow the scheme `V=LX+η`,
where V are the observations in `obs`, observed at times given in `obs_times`,
`L` are the observation operators given in `Ls`, `X` is the unobserved,
underlying diffusion and `η` is a Gaussian noise with mean `0` and covariance
`Σ` with the last stored in `Σs`. `fpt` provides additional information in case
the nature of observations has to do with first-passage times
"""
function set_observations!(setup::DiffusionSetup, Ls, Σs, obs, obs_times,
                           fpt=fill(nothing, length(obs)-1))
    setup.Ls = Ls
    setup.Σs = Σs
    setup.obs = obs
    setup.obs_times = obs_times
    setup.fpt = fpt
    setup.setup_completion[:obs] = true
end

function incomplete_message(::Val{:obs})
    print("\nThe observations have not been set up. Please call ",
          "set_observations!() passing the data on observations and the type ",
          "of observation scheme.\n")
end

"""
    set_imputation_grid!(setup::DiffusionSetup, dt,
                         time_transf=(t₀,T) -> ((x) ->  t₀ + (x-t₀) * (2-(x-t₀)/(T-t₀))))

Define the imputation grid in `setup`. `dt` defines the granulatrity of the
imputation grid and `time_transf` defines a time transformation to use for
transforming equidistant grid.
"""
function set_imputation_grid!(setup::DiffusionSetup, dt,
                              time_transf=(t₀,T) -> ((x) ->  t₀ + (x-t₀) * (2-(x-t₀)/(T-t₀))))
    setup.dt = dt
    setup.τ = time_transf
    setup.setup_completion[:imput] = true
end

function incomplete_message(::Val{:imput})
    print("\nThe imputation grid has not been set up. Please call ",
          "set_imputation_grid!() passing the `delta-t` and the time ",
          "transformation.\n")
end

"""
    set_x0_prior!(setup::DiffusionSetup, x0_prior, x0_guess=nothing)

Store the priors over the starting point into the object `setup`.
"""
function set_x0_prior!(setup::DiffusionSetup, x0_prior, x0_guess=nothing)
    setup.x0_prior = x0_prior
    if x0_guess == nothing
        @assert typeof(x0_prior) <: KnownStartingPt
        x0_guess = start_pt(nothing, x0_prior)
    end
    setup.x0_guess = x0_guess
    setup.setup_completion[:prior] = true
end

function incomplete_message(::Val{:prior})
    print("\nThe priors have not been set up. Please call set_priors!() ",
          "passing the priors on the parameters and on the starting point.\n")
end

function set_auxiliary!(setup::DiffusionSetup; skip_for_save=nothing,
                        adaptive_prop=nothing)
    if skip_for_save !== nothing
        setup.skip_for_save = skip_for_save
    end
    if adaptive_prop !== nothing
        setup.adaptive_prop = adaptive_prop
    end
end

"""
    check_if_complete(setup::DiffusionSetup, labels=keys(setup.setup_completion))

Check if all the set-up steps listed in `labels` have been finalised in `setup`
"""
function check_if_complete(setup::DiffusionSetup,
                           labels=keys(setup.setup_completion))
    complete = true
    for label in labels
        if !setup.setup_completion[label]
            incomplete_message(Val(label))
            complete = false
        end
    end
    complete
end

"""
    determine_data_type(setup::DiffusionSetup)

Determine the data type of the containers with path and driving noise
"""
function determine_data_type(setup::DiffusionSetup)
    check_if_complete(setup, [:prior]) || throw(UndefRefError())
    x = setup.x0_guess
    drift = b(0.0, x, setup.P˟)
    vola = σ(0.0, x, setup.P˟)
    # @assert typeof(x) == typeof(drift) # maybe this assertion is too strong
    typeof(drift), typeof(vola' * drift)
end

"""
    prepare_containers!(setup::DiffusionSetup)

Set-up the containers for paths, driving noise, observations, observation noise
and observation operators, automatically choosing an appropriate data type based
on the return values of the drift and volatility coefficients
"""
function prepare_containers!(setup::DiffusionSetup)
    T, S = determine_data_type(setup)
    Wnr = Wiener{S}()
    TW = typeof(sample([0], Wnr))
    TX = typeof(SamplePath([], zeros(T, 0)))

    # TODO modify so that it will work with GPUArrays
    m = length(setup.obs)-1
    WW, XX = Vector{TW}(undef,m), Vector{TX}(undef,m)
    setup.Wnr, setup.WW, setup.XX = Wnr, WW, XX

    prepare_obs_containers!(T, setup)
end

"""
    prepare_obs_containers!(::Type{T}, setup::DiffusionSetup)

Define containers for the observations, observation noise and observation
operators based on the passed data-type
"""
function prepare_obs_containers!(::Type{T}, setup::DiffusionSetup) where T <: Number
    correct_data_type!(setup, T, T)
end

"""
    prepare_obs_containers!(::Type{T}, setup::DiffusionSetup)

Define containers for the observations, observation noise and observation
operators based on the passed data-type
"""
function prepare_obs_containers!(::Type{T}, setup::DiffusionSetup) where T <: Array
    correct_data_type!(setup, Array, Array)
end

"""
    prepare_obs_containers!(::Type{T}, setup::DiffusionSetup)

Define containers for the observations, observation noise and observation
operators based on the passed data-type
"""
function prepare_obs_containers!(::Type{T}, setup::DiffusionSetup) where T <:SArray
    f(x) = SMatrix{_dim(x)...}(x)

    g(x) = g(x, Val{ndims(x)}())
    g(x, ::Val{0}) = SVector{1}(x)
    g(x, ::Val{1}) = SVector{size(x)...}(x)
    g(x, ::Val{2}) = SMatrix{size(x)...}(x)

    correct_data_type!(setup, f, g)
end

"""
    correct_data_type!(setup::DiffusionSetup, f, g)

Transform the elements of observations, observation noise and observation
operators collections stored in `setup` according to functions `f` and `g`
"""
function correct_data_type!(setup::DiffusionSetup, f, g)
    setup.Σs = map(f, setup.Σs)
    setup.Ls = map(f, setup.Ls)
    setup.obs = map(g, setup.obs)
end

"""
    _dim(mat::T) where T <: Number

Size of a matrix corresponding to a scalar is (1,1)
"""
_dim(mat::T) where T <: Number = 1, 1

"""
    _dim(mat::T) where T <: Union{Array, SArray}

Size of a matrix corresponding to an element `mat`
"""
_dim(mat::T) where T <: Union{Array, SArray} = _dim(mat, Val{ndims(mat)}())

"""
    _dim(mat, ::Val{1})
Size of a matrix corresponding to a vector is (vector length, 1)
"""
_dim(mat, ::Val{1}) = size(mat)[1], 1

"""
    _dim(mat, ::Val{1})
Size of a matrix corresponding to a matrix `mat` is just size(mat)
"""
_dim(mat, ::Val{2}) = size(mat)

"""
    _dim(mat, ::T) where T
Size of a matrix corresponding to tensor with dimension larger than matrix is
undefined
"""
_dim(mat, ::T) where T = throw(ArgumentError())

function _build_time_grid(τ, dt, t0, T)
    num_pts = Int64(ceil((T-t0)/dt))+1
    tt = τ(t0, T).( range(t0, stop=T, length=num_pts) )
    tt
end

 #TODO remove this K, it's most likely not needed
"""
    find_proposal_law!(::Type{K}, setup::DiffusionSetup, solver, change_pt
                            ) where K

Initialise the object with proposal law and all the necessary containers needed
for the simulation of the guided proposals
"""
function find_proposal_law!(::Type{K}, setup::DiffusionSetup, solver, change_pt
                            ) where K
    xx, tt, P˟, P̃ = setup.obs, setup.obs_times, setup.P˟, setup.P̃
    Ls, Σs, dt, τ = setup.Ls, setup.Σs, setup.dt, setup.τ
    #NOTE remember to pre-allocate space for the M,L,mu containers for change_pt with blocking
    #change_pt = typeof(setup.change_pt)(get_change_pt(setup.blocking_params[3]))

    m = length(xx)-1
    t = _build_time_grid(τ, dt, tt[m], tt[m+1])
    P = GuidPropBridge(K, t, P˟, P̃[m], Ls[m], xx[m+1], Σs[m];
                       change_pt=change_pt, solver=solver)
    #P = Array{ContinuousTimeProcess,1}(undef,m)
    Ps = []
    push!(Ps, deepcopy(P))

    for i in m-1:-1:1
        t = _build_time_grid(τ, dt, tt[i], tt[i+1])
        P = GuidPropBridge(K, t, P˟, P̃[i], Ls[i], xx[i+1], Σs[i], P.H[1],
                           P.Hν[1], P.c[1]; change_pt=change_pt, solver=solver)
        push!(Ps, deepcopy(P))
    end
    Ps = [P for P in Ps]
    setup.P = reverse(Ps)
end

"""
    initialise(::Type{K}, setup::DiffusionSetup; verbose=false,
                        change_pt=NoChangePt())

Initialise the internal containers of `setup`. Check if all the necessary data
has been passed to `setup`
"""
function initialise!(::Type{K}, setup::DiffusionSetup, solver, verbose=false,
                     change_pt=NoChangePt()) where K
    verbose && print("Initialising MCMC setup...\nPreparing containers...\n")
    prepare_containers!(setup)
    verbose && print("Initialising proposal laws...\n")
    find_proposal_law!(K, setup, solver, change_pt) #TODO remove this K, it's most likely not needed
    #TODO initialise for computation of gradients
end

"""
    check_if_recompute_ODEs(Ps, updt_coord)

Utility function for checking if H,Hν,c need to be re-computed for a respective
parameter update
"""
function check_if_recompute_ODEs(Ps, updt_coord)
    any([any([uc in depends_on_params(P) for uc in updt_coord]) for P in Ps])
end

#NOTE code to adapt:
#θ = params(P[1].Target)
#ϑs = [[θ[j] for j in idx(uc)] for uc in updtCoord]
#result = [DiffResults.GradientResult(ϑ) for ϑ in ϑs]
#resultᵒ = [DiffResults.GradientResult(ϑ) for ϑ in ϑs]
#Q = eltype(result)


text:
using Documenter, BioAlignments

makedocs(
    format = :html,
    sitename = "BioAlignments.jl",
    pages = [
        "Home" => "index.md",
        "Alignment representation" => "alignments.md",
        "Pairwise alignment" => "pairalign.md",
        "IO" => [
            "SAM and BAM" => "hts-files.md"
        ],
        "API Reference" => "references.md"
    ],
    authors = "Kenta Sato, Ben J. Ward, The BioJulia Organisation and other contributors."
)
deploydocs(
    repo = "github.com/BioJulia/BioAlignments.jl.git",
    julia = "1.0",
    osname = "linux",
    target = "build",
    deps = nothing,
    make = nothing
)


text:
module GlorotUniform
    function GUMatrix(input_size::Int64, layer_size::Int64, shape::Any)
        return rand(-sqrt(6/(input_size+layer_size)):0.00001:sqrt(6/(input_size+layer_size)), shape)
    end
    export GUMatrix
end

# Source: https://keras.io/api/layers/initializers/#glorotuniform-class


text:
#=
Popis
  Funkcia vypocita hodnotu radiality centrality pre kazdy vrchol.

Parametre
  g : datova struktura grafu :: GenericGraph
  dist_key : nazov atributu, v ktorom je ulozene ohodnotenie hrany :: String
  normalize : parameter pre aplikovanie normalizacie :: Bool

Navratova hodnota
  radiality : slovnik s mapovanim (kluc=vrchol_grafu => hodnota=radiality) :: Dict{Int64, Float64}

Poznamky
  Implementacia tejto funkcie vyuziva schopnost poli, ktora sa nazyva list comprehension.
  Vdaka tejto vlastnosti je mozne vykonat zlozitejsie operacie, funkcie, podmienky priamo pri tvorbe pola.
=#

function centrality_radiality_list_comprehension(g::GenericGraph; dist_key::AbstractString="weight", normalize::Bool=true)

  dia = graph_diameter(g)

  if normalize
    return radiality = [v.index => ((sum([ x!= Inf ? dia+1-x : 0.0 for x=dijkstra_shortest_paths(g, AttributeEdgePropertyInspector{Float64}(dist_key), [v]).dists])-dia-1)/(num_vertices(g) - 1))/dia for v=vertices(g)]
  elseif !(normalize)
    return radiality = [v.index => ((sum([ x!= Inf ? dia+1-x : 0.0 for x=dijkstra_shortest_paths(g, AttributeEdgePropertyInspector{Float64}(dist_key), [v]).dists])-dia-1)/(num_vertices(g) - 1)) for v=vertices(g)]
  end

end

max_stars_repo_name: mpuk/GraphAnalysis

max_issues_repo_name: mpuk/GraphAnalysis

max_forks_repo_name: mpuk/GraphAnalysis


text:
export ProductDomain

immutable ProductDomain{D,T,dim} <: Domain{T,dim}
    domains::D
end

ProductDomain(d::Tuple) =
    ProductDomain{typeof(d),mapreduce(eltype,promote_type,d),mapreduce(dimension,+,d)}(d)

fromcanonical(d::BivariateDomain,x::Tuple)=fromcanonical(d,x...)
tocanonical(d::BivariateDomain,x::Tuple)=tocanonical(d,x...)

# product domains are their own canonical domain
for OP in (:fromcanonical,:tocanonical)
    @eval $OP(::ProductDomain,x,y)=(x,y)
end

ProductDomain(A,B)=ProductDomain((A,B))
*(A::ProductDomain,B::ProductDomain)=ProductDomain(tuple(A.domains...,B.domains...))
*(A::ProductDomain,B::Domain)=ProductDomain(tuple(A.domains...,B))
*(A::Domain,B::ProductDomain)=ProductDomain(tuple(A,B.domains...))
*(A::Domain,B::Domain)=ProductDomain(A,B)

Base.length(d::ProductDomain)=length(d.domains)
Base.transpose(d::ProductDomain)=ProductDomain(d[2],d[1])
Base.getindex(d::ProductDomain,k::Integer)=d.domains[k]
==(d1::ProductDomain,d2::ProductDomain)=d1.domains==d2.domains

Base.first(d::ProductDomain)=(first(d[1]),first(d[2]))

function pushappendpts!(ret,xx,pts)
    if isempty(pts)
        push!(ret,xx)
    else
        for x in pts[1]
            pushappendpts!(ret,(xx...,x),pts[2:end])
        end
    end
    ret
end

function checkpoints(d::ProductDomain)
    pts=map(checkpoints,d.domains)
    ret=Array(Tuple{map(eltype,d.domains)...},0)
    pushappendpts!(ret,(),pts)
    ret
end

function points(d::ProductDomain,n::Tuple)
    @assert length(d.domains) == length(n)
    pts=map(points,d.domains,n)
    ret=Array(Tuple{map(eltype,d.domains)...},0)
    pushappendpts!(ret,(),pts)
    ret
end

Base.reverse(d::ProductDomain)=ProductDomain(map(reverse,d.domains))

domainscompatible(a::ProductDomain,b::ProductDomain) =
                        length(a.domains)==length(b.domains) &&
                        all(map(domainscompatible,a.domains,b.domains))