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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
        el
    end
end

# ╔═╡ 0232943f-74f8-40ec-96a3-495d5942e897
if occursin("Sesam", Base.current_project())
	import Pkg
	# activate the shared project environment
	Pkg.activate(Base.current_project())
	# instantiate, i.e. make sure that all packages are downloaded
	Pkg.develop("MTKHelpers")
	Pkg.instantiate()
	using PlutoLinks: @revise
	using Sesam, MTKHelpers
	# currently using @revise enters an endless loop of reloading
	#@revise using Sesam
	#@revise using MTKHelpers
else
	using Sesam, MTKHelpers
end

# ╔═╡ d6f85fa6-97cc-40e9-841a-053563713993
begin
	using Plots, PlutoUI
	import PlutoUI: combine
	gr()
end;

# ╔═╡ 2c1e57e8-2069-4268-ab11-9b514054c9f6
using Suppressor

# ╔═╡ 01be1fe6-39f5-4ecb-9d9e-3c5bcf39c936
using ModelingToolkit, DifferentialEquations, DataFrames, Tables,Distributions,Chain

# ╔═╡ f1a94d8f-927e-40f3-9af6-667890a6e733
using DistributionFits

# ╔═╡ 032c42f2-5103-11eb-0dce-e7ec59924648
html"<button onclick='present()'>present</button>"

# ╔═╡ 827fbc3a-512d-11eb-209e-cd74ddc17bae
begin
	struct TwoColumn{A, B}
		left::A
		right::B
	end
	function Base.show(io, mime::MIME"text/html", tc::TwoColumn)
		write(io,
			"""

			""")
		show(io, mime, tc.left)
		write(io,
			"""

			""")
		show(io, mime, tc.right)
		write(io,
			"""

		""")
	end
end

# ╔═╡ 6e60edbd-82eb-4283-91ec-8fcb33a64ff7
md"## Testing local image"

# ╔═╡ a9ff79e4-f0dd-4c23-8755-379ab84475bf
pwd()

# ╔═╡ a392d095-19f6-495f-8c59-a93eace2b9a6
@suppress LocalResource("fig/logo.png", :width => 200)

# ╔═╡ 6382a73e-5102-11eb-1cfb-f192df63435a
md"""
# Sesam3 interactive exploration
**Thomas Wutzler 2022**
"""

# ╔═╡ 26a15193-a85f-443c-b007-2d14915e69f7
md"## Setup prior parameter distributions"

# ╔═╡ 078f11d4-9630-4e5c-8048-ef0c0232294f
md"Define the System components"

# ╔═╡ 36b322a5-ef4a-4697-90a8-720d8c224219
begin
	tspinup = 2000.0; tface=100.0
	#@named s = sesam3(use_seam_revenue=true)
	@named s = sesam3CN(use_seam_revenue=false)
	@named pl = plant_face(t1=0.0,t2=tface)
	@named sp = plant_sesam_system(s,pl)
end;

# ╔═╡ 351ebf17-102c-4497-aab0-1e0a3f55f799
md"Define parameters and their distributions"

# ╔═╡ 7c778bf7-de69-4c14-99b6-b76cdf06596c
begin
	p = pC = Dict(
	    s.ϵ_tvr => 0.45,   # carbon use efficiency of microbial tvr (part by predators
	        #which respire and corresponding amount of N must be mineralized)
	    s.κ_E => 0.8,     ##<< amount of recycling enzyme turnover by biomass (
	        # added to assimilable, i.e. uptake instead of R)
	    s.a_E => 0.001*365,   ##<< C biomass allocated to enzymes 1/day /microbial biomass
	    s.m => 0.005*365,    ##<< maintenance respiration rate   1/day /microbial biomass,
	    s.τ => 1/60*365,  ##<< biomass turnover rate (12 days)
	    s.k_L => 1.0,       ##<< 1/(x years)
	    #s.k_L => 5.0,       ##<< 1/(x years)   # formerly 1 year
	    s.k_R => 1/(40.0),        ##<< 1/(x years) # to demonstrate changes on short time scale
	    s.k_mN_L => 0.05 * 60, # enzyme half-saturation constant, in magnitude of enzymes *
	        # /yr enzyme turnover 60 times a year
		s.k_mN_R => 0.05 * 60, # enzyme half-saturation constant, in magnitude of enzymes *
		# /yr enzyme turnover 60 times a year
	    s.ϵ => 0.5,      ##<< carbon use efficiency for growth respiration
	    #i_L => t -> 1 - exp(-t),  # litter input
	    pl.i_L0 => 400.0,         # g/m2 input per year (half NPP)
	    #pl.β_Ni0 => 25.0,
	    pl.β_Ni0 => 30.0,
	    #pl.i_IN0 => 0.0,   ##<< input of mineral N,
	    pl.i_IN0 => 0.7,   ##<< 7kg/ha/yr
	    #pl.k_Lagr => 12/2, # above ground litter turnover of 2 month
		# P from plant model parameters not used in CN-Sesam soil model
		pl.β_Pi0 => Inf, #25*20, ## leaf litter N:P ~20(massratio Kang10)
		pl.i_IP0 => Inf, #0.65,   ##<< input of mineral P, weathering: Table3 mixed sedimentary rocks 0.65g/m2/yr Hartmann14 10.1016/j.chemgeo.2013.10.025
		pl.s_EP0 => Inf, # 0.5, # plant 1/20 of typical total microbial enzyme synthesis flux
		pl.u_PlantPmax0 => Inf,
		pl.k_PlantP0 => Inf,
	)
	pN = Dict(
	    s.i_BN => 0.4, ##<< potential immobilization flux rate
	    s.β_NEnz => 3.1,     # Sterner02: Protein (Fig. 2.2.), high N investment (low P) need
	    s.β_NB => 11.0,
	    s.l_N => 0.96,       #0.00262647*365     ##<< leaching rate of mineralN lN IN
	    #s.l_N => 0.0,
	    s.ν_N =>  0.9,     # microbial N use efficiency accounting for apparent
	    ## minceralization of N during uptake in heterogeneous soils
	    # for N uptake: take the defaults which take as much N as supplied by litter
	    # pl.u_PlantNmax0 => Inf32, # only constrained by k_PlantN in min function
	    # pl.k_PlantN0 => 10.57, #0.0289652*365     ##<< plant uptake rate first order of IN
	)
	p = merge(pC, pN)

	u0 = u0C = Dict(
	    s.B => 17,
	    #s.B => 34,
	    s.L => 100,
	    #s.L => 110,
	    s.R => 1100,
	    #s.R => 3250,
	    #s.cumresp => 0.0,
	    s.α_R => 0.5,
	)
	u0C[s.α_L] = 1.0 - u0C[s.α_R]
	u0N = Dict(
	    s.I_N => 0.04, ##<< inorganic pool gN/m2
	    s.L_N => u0[s.L]/p[pl.β_Ni0],
	    #s.R_N => u0[s.R]/calculate_β_NR_sesam3(p,s) #p[s.β_NB],
	    s.R_N => u0[s.R]/7.0
	    )
	u0 = merge(u0C, u0N)
end;

# ╔═╡ 8247cdc0-8b23-4b1f-a4ee-6e6ac09a2232
begin
	cols = (:par, :dType, :mode, :upper)
	parmsModeUpperRows = [
	     (s.β_NB, LogNormal, 8. , 16.),
	     (s.β_NEnz, LogNormal, 3.0 , 3.5),
	     (s.k_mN_L, LogNormal, 60*0.05 , 120*2.),
	     (s.k_mN_R, LogNormal, 60*0.05 , 120*2.),
	     (s.κ_E, LogNormal, 0.7 , 0.9),
	     (s.k_R, LogNormal, 1/100, 1/10),
	     (s.k_L, LogNormal, 12/18, 12/4), # tvr time between 18 months and 4 months
	     (s.a_E, LogNormal, 0.001*365 , 0.005*365),
	     (s.m, LogNormal, 0.005*365 , 0.02*365),
	     (s.τ, LogNormal, 1.0/60*365 , 1.0/5*365),
	     (s.ϵ, LogitNormal, 0.5 , 0.7),
	     (s.ϵ_tvr, LogitNormal, 0.3 , 0.9),
	     (s.i_BN, LogNormal, 0.4 , 4),
	     (s.l_N, LogNormal, 0.96 , 5*0.96),
	     (s.ν_N, LogNormal, 0.9 , 0.99),
	     #(:kIPlant, LogNormal, 10.57 , 20)
	]
	df_dist = rename!(DataFrame(columntable(parmsModeUpperRows)), collect(cols))
end;

# ╔═╡ ad1b457f-f98e-4556-83f2-477e2b857a07
#df_dist

# ╔═╡ 383cfe5b-922e-4151-9375-6093523d7f03
begin
	f1v = (dType, mode, upper) -> fit(dType, @qp_m(mode), @qp_uu(upper))
	transform!(df_dist, Cols(:dType,:mode,:upper) => ByRow(f1v) => :dist)
end;

# ╔═╡ 5c6d02b6-0aed-4ad3-843f-7e46367171cf
begin
	transform!(df_dist, :dist => ByRow(d -> quantile(d, 0.025)) => :lower)
	transform!(df_dist, :par => ByRow(pname -> p[pname]) => :default)
	select!(df_dist, Not(:dist), :dist)
end

# ╔═╡ d45a7ad5-4803-4846-ad2e-a9d4fd294413
md"""
## Modifying parameters

Select the parameters to explore with sliders below. (Below the plot)
"""

# ╔═╡ 1d6e347e-1a1b-4f2a-bc48-524c8d36f771
md"## Plotting the solution"

# ╔═╡ 7f1a89ce-1b84-45f0-af85-87b1050b18ba
begin
	# tstart_scrubb = @bind tstart NumberField(tstart_coarse.+(-20:20),default=tstart_coarse)
	# tend_scrubb = @bind tend NumberField(tend_coarse.+(-20:20),default=tend_coarse)
	tstart_scrubb = @bind tstart NumberField(-tspinup:tface+20,default=-2)
	tend_scrubb = @bind tend NumberField(-tspinup:tface+20,default=+10)
end;

# ╔═╡ 9f251cbd-777e-4904-814d-6cef789fc560
@bind b_settimes Button("Set times programmatically")

# ╔═╡ 4b08433e-b8d8-4342-9153-2e4dad20356e
md"""
tstart: $(tstart_scrubb), tend: $(tend_scrubb)
"""

# ╔═╡ fdcefaa1-7ebf-49a1-ade5-00302ddeb6c9
@bind vars_disp MultiCheckBox(vcat(states(sp),getproperty.(observed(sp),:lhs)))

# ╔═╡ ba4b67ab-b3fe-44f8-aeaf-925f76b005b5
md"""
Try modifying
- β_NB between 11.6 and 11.8
- s.ϵ_tvr from 0.5 to 0.509
"""

# ╔═╡ 10214a12-e0e9-40aa-8a65-e8a9661125f6
@bind pars_mod MultiCheckBox(df_dist.par)

# ╔═╡ 9fe4cc1e-59d7-43c1-b10c-764756f5b4e8
function parameter_input(pars_mod)
	return PlutoUI.combine() do Child
		inputs = [begin
			row = subset(df_dist, :par => ByRow(x -> symbol(x) == symbol(name)))[1,:]
			el = Scrubbable(range(row.lower, row.upper, length = 100);default=row.default, format = ".3")
			# el = NumberField(range(first(row.lower), first(row.upper), length = 30);default=first(row.mode))
			md""" $(name): $(
				Child(name, el)
			)"""
		end
			for name in pars_mod
		]
		md"""
		$(inputs)
		"""
	end
end

# ╔═╡ 854ab5c5-a572-41e8-90ee-bf8a63d7adf9
@bind popt parameter_input(string.(pars_mod))

# ╔═╡ 66dae3d2-9ba5-491b-a3cd-77fcfbb0ab78
ps = ProblemParSetter(sp, propertynames(popt))

# ╔═╡ a8b21bd2-3914-46b2-8b11-f9279f88e732
begin
	tspan_sim = (-tspinup,tface) # simulate 500 yrs spinup, increase at yr 20
	prob0 = prob = ODEProblem(sp, u0, tspan_sim, p) #2ms
	if length(popt) != 0
		prob = update_statepar(ps, values(popt), prob0)
		p_upd = label_par(ps,prob.p)
	end
	sol0t = sol = solve(prob)
end;

# ╔═╡ 2ddb2869-4917-4973-8b3f-51dce7eba302
length(vars_disp) != 0 && plot(sol0t, tspan = (tstart, tend), vars=vars_disp)

# ╔═╡ b565d593-1369-4efe-8704-90f1611526c4
propertynames(popt)

# ╔═╡ 61e14646-568d-49a9-8d5e-afe9555e29ff
popt

# ╔═╡ 0b578017-ed4f-4dbf-80c8-7c5492799d17
plot(sol0t, tspan = (tstart, tend), vars=[s.dec_R p_upd[symbol(s.ϵ_tvr)] * s.tvr_B + (1-p_upd[symbol(s.κ_E)])*s.tvr_Enz], legend=:bottomright); vline!([0], linestyle=:dash, color = :lightgray, label=nothing)

# ╔═╡ 2448b012-864d-4797-b768-33174f9fd08f
plot(sol0t, tspan = (tstart, tend), vars=[s.α_L,s.α_R], legend=:bottomright); vline!([0], linestyle=:dash, color = :lightgray, label=nothing)

# ╔═╡ 773273a8-70fc-4897-8174-8b06aff47013
plot(sol0t, tspan = (tstart, tend), vars=[s.lim_C, s.lim_N]); vline!([0], linestyle=:dash, color = :lightgray, label=nothing)

# ╔═╡ 0cbe8944-5ebd-4804-863d-e567c703b89f
plot(sol0t, tspan = (tstart, tend), vars=[s.I_N]); vline!([0], linestyle=:dash, color = :lightgray,label=nothing)

# ╔═╡ 1b5413de-3d32-4086-b80c-11c10ee87dac
plot(sol0t, tspan = (tstart, tend), vars=[s.Φ_NB,s.Φ_N-s.i_L/s.β_Ni]); hline!([0], linestyle=:dash, color=:lightgray,label=nothing); vline!([0], linestyle=:dash, color = :lightgray, label=nothing)

# ╔═╡ 255018b9-28c0-442c-8fa9-9698e7d7b513
md"## Thanks"

# ╔═╡ Cell order:
# ╠═0232943f-74f8-40ec-96a3-495d5942e897
# ╟─d6f85fa6-97cc-40e9-841a-053563713993
# ╟─032c42f2-5103-11eb-0dce-e7ec59924648
# ╟─827fbc3a-512d-11eb-209e-cd74ddc17bae
# ╠═6e60edbd-82eb-4283-91ec-8fcb33a64ff7
# ╠═a9ff79e4-f0dd-4c23-8755-379ab84475bf
# ╠═2c1e57e8-2069-4268-ab11-9b514054c9f6
# ╠═a392d095-19f6-495f-8c59-a93eace2b9a6
# ╟─6382a73e-5102-11eb-1cfb-f192df63435a
# ╟─26a15193-a85f-443c-b007-2d14915e69f7
# ╠═01be1fe6-39f5-4ecb-9d9e-3c5bcf39c936
# ╠═f1a94d8f-927e-40f3-9af6-667890a6e733
# ╟─078f11d4-9630-4e5c-8048-ef0c0232294f
# ╠═36b322a5-ef4a-4697-90a8-720d8c224219
# ╟─351ebf17-102c-4497-aab0-1e0a3f55f799
# ╟─7c778bf7-de69-4c14-99b6-b76cdf06596c
# ╠═8247cdc0-8b23-4b1f-a4ee-6e6ac09a2232
# ╠═ad1b457f-f98e-4556-83f2-477e2b857a07
# ╠═383cfe5b-922e-4151-9375-6093523d7f03
# ╟─5c6d02b6-0aed-4ad3-843f-7e46367171cf
# ╟─d45a7ad5-4803-4846-ad2e-a9d4fd294413
# ╟─a8b21bd2-3914-46b2-8b11-f9279f88e732
# ╟─66dae3d2-9ba5-491b-a3cd-77fcfbb0ab78
# ╠═b565d593-1369-4efe-8704-90f1611526c4
# ╟─1d6e347e-1a1b-4f2a-bc48-524c8d36f771
# ╟─7f1a89ce-1b84-45f0-af85-87b1050b18ba
# ╠═9f251cbd-777e-4904-814d-6cef789fc560
# ╠═4b08433e-b8d8-4342-9153-2e4dad20356e
# ╠═fdcefaa1-7ebf-49a1-ade5-00302ddeb6c9
# ╠═2ddb2869-4917-4973-8b3f-51dce7eba302
# ╟─ba4b67ab-b3fe-44f8-aeaf-925f76b005b5
# ╠═10214a12-e0e9-40aa-8a65-e8a9661125f6
# ╟─854ab5c5-a572-41e8-90ee-bf8a63d7adf9
# ╟─61e14646-568d-49a9-8d5e-afe9555e29ff
# ╟─9fe4cc1e-59d7-43c1-b10c-764756f5b4e8
# ╠═0b578017-ed4f-4dbf-80c8-7c5492799d17
# ╠═2448b012-864d-4797-b768-33174f9fd08f
# ╠═773273a8-70fc-4897-8174-8b06aff47013
# ╠═0cbe8944-5ebd-4804-863d-e567c703b89f
# ╠═1b5413de-3d32-4086-b80c-11c10ee87dac
# ╟─255018b9-28c0-442c-8fa9-9698e7d7b513


text:
export buildtrainingdata, loadtrainingdata, loaddemanddata, savetrainingdata

# Can call this with ssp5 or ssp32 (global constants)
# Since not all GADM level 0 regions are assigned to SSP regions, some regions are missed
# (mostly small island nations and the like). This captures 99.85% of global population.
function make_sspregionlookup(ssp)
    println("Lookup SSP R5.2 region names for each GADM level 0 region...")
    _, _, regionlist, _, _ = loadregions("Global_GADM0")
    sspregions = fill("", length(regionlist))
    for (sspregionname, gadmregionnames) in ssp
        for gname in gadmregionnames
            index = findfirst(isequal(Symbol(gname)), regionlist)
            if index == nothing
                error("Region $gname is not a GADM level 0 region.")
            end
            sspregions[index] = sspregionname
        end
    end
    return sspregions # Vector{String}(length numcountries, indexed by gadm country code)
end

function getnationalpopulation(scenarioyear)
    filename = in_datafolder("nationalpopulation_$scenarioyear.jld")
    natpop = isfile(filename) ? JLD.load(filename, "natpop") : savenationalpopulation(scenarioyear)
    return natpop
end

# load population dataset and sum globally by country
# Vector{Float64}(length numcountries, indexed by gadm country code)
function savenationalpopulation(scenarioyear)
    println("Calculating population in all GADM level 0 regions...")
    pop = JLD.load(in_datafolder("population_$scenarioyear.jld"), "population")   # persons per pixel
    regions, _, regionlist, lonrange, latrange = loadregions("Global_GADM0")
    natpop = zeros(length(regionlist))
    for j in latrange
        for i in lonrange
            reg = regions[i,j]
            if reg > 0
                natpop[reg] += pop[i,j]
            end
        end
    end
    JLD.save(in_datafolder("nationalpopulation_$scenarioyear.jld"), "natpop", natpop, compress=true)
    return natpop
end

# Get current national electricity demand from IEA online statistics (actually "domestic supply").
# https://www.iea.org/statistics/?country=MOROCCO&year=2016&category=Electricity&indicator=ElecGenByFuel&mode=table&dataTable=ELECTRICITYANDHEAT
function ieademand()
    println("Get current national electricity demand from IEA statistics...")
    iea = CSV.read(in_datafolder("ieademand_2016.csv"), DataFrame)      # GWh/year
    _, _, regionlist, _, _ = loadregions("Global_GADM0")
    nationaldemand = zeros(length(regionlist))
    for row in eachrow(iea)
        country, demand = row
        country == uppercase(country) && continue   # skip IEA aggregated regions in uppercase, e.g. "ASIA"
        f = findfirst(regionlist .== Symbol(country))
        f == nothing && error("IEA country $country missing in GADM regions.")
        nationaldemand[f] = demand/1000
    end
    return nationaldemand        # TWh/year
end

# Lookup electricity demand from a SSP database dataframe. Return a single dataframe row.
ssplookup(ssp, model, scen, region) =
    ssp[(ssp[:,:MODEL] .== model) .& (ssp[:,:SCENARIO] .== scen) .& (ssp[:,:REGION] .== region), :][1,:]

# Calculate regional multiplers for demand from 2016 to target year
function calcdemandmultipliers(sspscenario::String, year::Int)
    println("Calculate demand multipliers...")
    # First read electricity demand from the SSP database into a dataframe.
    ssp = CSV.read(in_datafolder("SSP v2 Final Energy - Electricity.csv"), DataFrame)

    # SSP scenarios also include radiative forcing targets, e.g. SSP2-34
    # (only for IAM energy system results, not underlying population & GDP scenario)
    # We'll hardcode the 3.4 W/m2 scenario variant for now and make it a configurable option later.
    scen = uppercase(sspscenario)   # e.g. "SSP2-34"

    # We'll take the average result from IMAGE and MESSAGE models.
    demandmult2050 = Dict{String, Float64}()
    for sspreg in keys(ssp5)
        reg = "R5.2$sspreg"
        image = ssplookup(ssp, "IMAGE", scen, reg)
        message = ssplookup(ssp, "MESSAGE-GLOBIOM", scen, reg)
        # Take the year 2020 to represent 2016, then adjust for this difference.
        image_mult = image[Symbol(year)] / image[Symbol(2020)]
        message_mult = message[Symbol(year)] / message[Symbol(2020)]
        demandmult2050[sspreg] = ((image_mult + message_mult)/2) ^ ((year-2016)/(year-2020))
    end
    return demandmult2050
end

function makeregionaldemanddata(gisregion, sspscenario::String, year::Int)
    res = 0.01          # resolution of auxiliary datasets [degrees per pixel]

    scenarioyear = "$(sspscenario[1:4])_$year"
    datasetinfo = Dict(:gisregion=>gisregion, :scenarioyear=>scenarioyear, :res=>res)
    regions, _, regionlist, _, pop, _, _, _, lonrange, latrange = read_datasets(datasetinfo)     # pop unit: people/grid cell

    gdp = JLD.load(in_datafolder("gdp_$(scenarioyear).jld"))["gdp"][lonrange,latrange]       # unit: USD(2010)/grid cell, PPP

    gadm0regions, _, gadm0regionlist, _, _ = loadregions("Global_GADM0")
    countrycodes = gadm0regions[lonrange,latrange]

    lats = (90-res/2:-res:-90+res/2)[latrange]
    cellarea = rastercellarea.(lats, res)
    popdens = pop ./ cellarea'

    natpop = getnationalpopulation(scenarioyear)            # people/GADMregion (vector, length numGADMcountries)
    demandpercapita = ieademand()./natpop * 1e6             # MWh/year/capita (vector, length numGADMcountries)
    ssp5region = make_sspregionlookup(ssp5)                 # SSP 5-region names (vector, length numGADMcountries)
    demandmult = calcdemandmultipliers(sspscenario, year)   # Dict: SSP region name => multiplier

    numreg = length(regionlist)
    regionaldemand = zeros(numreg)
    regionalpop = zeros(numreg)
    regionalgdp = zeros(numreg)

    println("Calculating annual electricity demand for model regions...")
    nlons, nlats = size(regions)
    updateprogress = Progress(nlats, 1)
    for j = 1:nlats
        for i = 1:nlons
            reg = regions[i,j]
            countrycode = countrycodes[i,j]
            if reg > 0 && reg != NOREGION && countrycode > 0
                sspreg = ssp5region[countrycode]
                if isempty(sspreg)
                    # error("Oops, no SSP region assigned to region $(gadm0regionlist[countrycode]).")
                    continue
                end
                regionaldemand[reg] += demandpercapita[countrycode] * pop[i,j]/1e6 * demandmult[sspreg]     # TWh/year
                regionalpop[reg] += pop[i,j]        # unit: people
                regionalgdp[reg] += gdp[i,j]        # unit: USD(2010)
            end
        end
        next!(updateprogress)
    end
    regionaldemandpercapita = regionaldemand ./ regionalpop * 1e6   # MWh/year/capita
    regionalgdppercapita = regionalgdp ./ regionalpop               # USD(2010)/capita

    return regionlist, regionaldemandpercapita, regionalgdppercapita
end

function buildtrainingdata(; gisregion="Europe8", sspscenario="ssp2-34", sspyear=2050, era_year=2018, numcenters=3, mindist=3.3)
    println("\nBuilding training data for $gisregion...")
    regionlist, demandpercapita, gdppercapita =
                    makeregionaldemanddata(gisregion, sspscenario, sspyear)
    scenarioyear = "$(sspscenario[1:4])_$sspyear"
    hours, temp_popcenters = GIStemp(gisregion, scenarioyear, era_year, numcenters, mindist)
    offsets, zone_maxpop, population = regional_timezone_offsets_Jan1(gisregion=gisregion, scenarioyear=scenarioyear, era_year=era_year)

    numreg, numhours = length(regionlist), length(hours)
    firsttime = ZonedDateTime.(hours[1], zone_maxpop)
    zonedtime = hcat(collect.([firsttime[i]:Hour(1):firsttime[i]+Hour(8759) for i = 1:numreg])...)[:]

    println("\nShifting hourly temperatures from UTC to local time...")
    temperature_top3_mean = dropdims(mean(temp_popcenters, dims=3), dims=3)
    quantiles = hcat([quantile(temp_popcenters[:,c,1], [0.05, 0.5, 0.95]) for c = 1:numreg]...)
    shifted_temperature_top3_mean = similar(temperature_top3_mean)
    shifted_temp1 = similar(temperature_top3_mean)
    for r = 1:numreg
        shifted_temperature_top3_mean[:,r] = circshift(temperature_top3_mean[:,r], round(Int, offsets[r]))
        shifted_temp1[:,r] = circshift(temp_popcenters[:,r,1], round(Int, offsets[r]))
    end

    # dataframe with hourly data
    df_time = DataFrame(
        localtime = DateTime.(zonedtime, Local),
        country = repeat(string.(regionlist), inner=numhours),
        temp_top3 = shifted_temperature_top3_mean[:],
        temp1 = shifted_temp1[:],
        localhour = hour.(zonedtime),
        month = month.(zonedtime),
        weekend01 = Int.(dayofweek.(zonedtime) .>= 6)
    )

    # sort by average monthly temperature (in popcenter 1), store rank in ranked_month
    df_monthlytemp = combine(groupby(df_time, [:country, :month]), :temp1 => mean => :temp_monthly) |>
            d -> sort!(d, [:country, :temp_monthly]) |>
            d -> insertcols!(d, 4, :ranked_month => repeat(1:12, outer=numreg))

    # dataframe with regional data
    df_reg = DataFrame(country=string.(regionlist), demandpercapita=demandpercapita, gdppercapita=gdppercapita,
                        temp1_qlow=quantiles[1,:], temp1_mean=quantiles[2,:], temp1_qhigh=quantiles[3,:])

    # join everything together
    df = innerjoin(df_time, df_monthlytemp, on=[:country, :month]) |>
                    d -> innerjoin(d, df_reg, on=:country)
    return df, offsets, population
end

function loadtrainingdata()
    filename = in_datafolder("syntheticdemand_trainingdata.csv")
    if !isfile(filename)
        savetrainingdata()
    end
    df_train = CSV.read(filename, DataFrame)
    offsets = CSV.read(in_datafolder("syntheticdemand_timezoneoffsets.csv"), DataFrame)[:,1]
    return df_train, offsets
end

function savetrainingdata(; numcenters=3, mindist=3.3)
    create_scenario_datasets("SSP2", 2020)
    println("\nCreating training dataset for synthetic demand...")
    println("(This requires ERA5 temperature data for the year 2015 and scenario datasets for SSP2 2020.)")
    df_train, offsets, _ = buildtrainingdata(gisregion="SyntheticDemandRegions",
                sspscenario="SSP2-34", sspyear=2020, era_year=2015,
                numcenters=numcenters, mindist=mindist)
    CSV.write(in_datafolder("syntheticdemand_timezoneoffsets.csv"), DataFrame(offsets=offsets))
    CSV.write(in_datafolder("syntheticdemand_trainingdata.csv"), df_train)
end

loaddemanddata() = CSV.read(in_datafolder("syntheticdemand_demanddata.csv"), DataFrame)

function regional_timezone_offsets_Jan1(; gisregion="Europe8", scenarioyear="ssp2_2050", era_year=2018)
    println("\nCalculating population-weighted regional time zone offsets...")
    regions, offshoreregions, regionlist, lonrange, latrange = loadregions(gisregion)
    tzindices, tznames = loadtimezones(lonrange, latrange)
    popscale = 1e5
    pop = JLD.load(in_datafolder("population_$scenarioyear.jld"), "population")[lonrange,latrange] ./ popscale   # scale down for better precision
    numreg = length(regionlist)
    numhours = 24*daysinyear(era_year)
    offsets = zeros(numreg)
    population = zeros(numreg)
    zone_maxpop = fill(TimeZone("Europe/London"), numreg)
    updateprogress = Progress(numreg, 1)
    for r = 1:numreg
        reg = (regions .== r)
        regindices = tzindices[reg]
        regpop = pop[reg]
        zoneindices = unique(regindices)
        weightedoffset = 0.0
        zones = TimeZone[]
        pops = Float64[]
        for idx in zoneindices
            tzname = tznames[idx]
            zone = tzname[1:3] == "Etc" ? TimeZone(tzname, TimeZones.Class(:LEGACY)) : TimeZone(tzname)
            push!(zones, zone)
            firsthour = ZonedDateTime(DateTime(era_year,1,1,0), zone)
            hours = firsthour : Hour(1) : firsthour + Hour(numhours-1)
            offset = tzoffset(hours[1])
            localpop = sum(regpop[regindices.==idx])
            push!(pops, localpop)
            weightedoffset += localpop * offset
        end
        _, i = findmax(pops)    # find time zone with the most population
        zone_maxpop[r] = zones[i]
        population[r] = sum(pops) * popscale    # scale up again
        offsets[r] = weightedoffset / sum(pops)
        next!(updateprogress)
    end

    return offsets, zone_maxpop, population
end

function saveVilledemand()
    vv = CSV.read(in_datafolder("syntheticdemand", "data", "df_model_features.csv"), DataFrame)
    hours = DateTime(2015, 1, 1, 0) : Hour(1) : DateTime(2015, 12, 31, 23)
    demand1D = hcat(DataFrame(hours=repeat(hours, outer=44)), select(vv, [:country, :demand_total_mwh]))
    demand2D = unstack(demand1D, :country, :demand_total_mwh)
    rename!(demand2D, Dict("Bosnia and Herz." => "Bosnia and Herzegovina",
                                "Czech Rep." => "Czech Republic",
                                "Korea" => "South Korea"))
    select!(demand2D, [:hours; sort(names(demand2D)[2:end])])

    dem2 = stack(demand2D, variable_name=:country, value_name=:demand_MW)
    select!(dem2, [3,1,2])
    meandemand = by(dem2, :country, meandemand = :demand_MW => mean)
    demand = join(dem2, meandemand, on=:country)

    insertcols!(demand, 5, normdemand=demand[:,:demand_MW]./demand[:,:meandemand])
    CSV.write(in_datafolder("syntheticdemand_demanddata.csv"), demand)
end

tzoffset(tz::FixedTimeZone) = tz.offset.std.value / 3600
tzoffset(tz::VariableTimeZone) = tzoffset(tz.transitions[end].zone)
tzoffset(dt::ZonedDateTime) = TimeZones.value(dt.zone.offset) / 3600

function timezonetest()
    # zones = TimeZone.(["Europe/London", "Europe/Stockholm", "Europe/Helsinki", "Europe/Minsk"])
    zones = TimeZone.(["Australia/Perth", "Australia/Adelaide", "Australia/Darwin", "Australia/Brisbane", "Australia/Sydney"])
    # all hours in 2015 for each zone in local time
    hours = hcat([ZonedDateTime(DateTime(2015,1,1,0), z) : Hour(1) : ZonedDateTime(DateTime(2015,12,31,23), z) for z in zones]...)
    # shift hours by the time zone offset to ensure simultaneous rows
    # take the offset of ZonedDateTime in hour 1 (in case hour 1 is DST), not the UTC time zone offset
    shifted = hcat([circshift(hours[:,i], round(Int, -tzoffset(hours[1,i]))) for i = 1:length(zones)]...)
    # convert to UTC and ensure all columns identical (except last 24 hours because of shifted hours from previous year)
    utctime = DateTime.(shifted, UTC)
    # @assert all([all(utctime[1:end-24,1] .== utctime[1:end-24,i]) for i=2:length(zones)])  # doesn't work for non-integer time offsets
    # show the local time in each zone around the shift to summer time (daylight savings) in the EU
    localtime = Time.(DateTime.(shifted, Local))
    # display(localtime[2085:2095,:])   # EU
    display(localtime[6615:6620,:])     # Australia
    return hours, shifted, utctime, localtime
end

ext: jl

max_stars_repo_name: henrikhodel/GlobalEnergyGIS

max_issues_repo_name: henrikhodel/GlobalEnergyGIS

max_forks_repo_name: henrikhodel/GlobalEnergyGIS


text:
function initialize!(p, ::Type{PyPlotType})
    PyPlot=p.context[:Plotter]
    PyPlot.PyObject(PyPlot.axes3D)# see https://github.com/JuliaPy/PyPlot.jl/issues/351
    if !haskey(p.context,:figure)
        res=p.context[:resolution]
        p.context[:figure]=PyPlot.figure(p.context[:fignumber],dpi=100)
        p.context[:figure].set_size_inches(res[1]/100,res[2]/100,forward=true)
        for ctx in p.subplots
            ctx[:figure]=p.context[:figure]
        end
    end
    if p.context[:clear]
        p.context[:figure].clf()
        p.context[:revealed]=false
    end
end

function save(fname,p,::Type{PyPlotType})
    p.context[:figure].savefig(fname)
end

function save(fname,scene,PyPlot,::Type{PyPlotType})
       isnothing(scene) ? nothing : scene.savefig(fname)
end

function reveal(p::GridVisualizer,::Type{PyPlotType})
    p.context[:revealed]=true
    p.Plotter.tight_layout()
    if !(isdefined(Main,:PlutoRunner))&& isinteractive()
        p.Plotter.pause(1.0e-10)
        p.Plotter.draw()
        sleep(1.0e-3)
    end
    p.context[:figure]
end

function reveal(ctx::SubVisualizer,TP::Type{PyPlotType})
    yield()
    if ctx[:show]||ctx[:reveal]
        reveal(ctx[:GridVisualizer],TP)
    end
end

#translate Julia attribute symbols to pyplot-speak
const mshapes=Dict(
    :dtriangle => "v",
    :utriangle => "^",
    :rtriangle => ">",
    :ltriangle => "^",
    :circle => "o",
    :square => "s",
    :cross => "+",
    :+ => "+",
    :xcross => "x",
    :x => "x",
    :diamond => "D",
    :star5 => "*",
    :pentagon => "p",
    :hexagon => "h",
)

const lstyles=Dict(
    :solid => "-",
    :dot => "dotted",
    :dash => "--",
    :dashdot => "-.",
    :dashdotdot => (0, (3, 1, 1, 1))
)

const leglocs=Dict(
    :none => "",
    :best => "best",
    :lt => "upper left",
    :ct => "upper center",
    :rt => "upper right",
    :lc => "center left",
    :cc => "center center",
    :rc => "center right",
    :lb => "lower left",
    :cb => "lower center",
    :rb => "lower right"
)

"""
$(SIGNATURES)
Return tridata to be splatted to PyPlot calls
"""
function tridata(grid)
    coord=grid[Coordinates]
    cellnodes=Matrix(grid[CellNodes])
    coord[1,:], coord[2,:],transpose(cellnodes.-1)
end

# Interfaces to Colors/Colorschemes
plaincolormap(ctx)=colorschemes[ctx[:colormap]].colors

### 1D grid
function gridplot!(ctx, TP::Type{PyPlotType}, ::Type{Val{1}}, grid)
    PyPlot=ctx[:Plotter]

    if !haskey(ctx,:ax)
        ctx[:ax]=ctx[:figure].add_subplot(ctx[:layout]...,ctx[:iplot])
    end
    if ctx[:clear]
        ctx[:ax].cla()
    end
    ax=ctx[:ax]
    fig=ctx[:figure]

    cellregions=grid[CellRegions]
    cellnodes=grid[CellNodes]
    coord=grid[Coordinates]
    ncellregions=grid[NumCellRegions]
    bfacenodes=grid[BFaceNodes]
    bfaceregions=grid[BFaceRegions]
    nbfaceregions=grid[NumBFaceRegions]
    ncellregions=grid[NumCellRegions]

    crflag=ones(Bool,ncellregions)
    brflag=ones(Bool,nbfaceregions)

    xmin=minimum(coord)
    xmax=maximum(coord)
    h=(xmax-xmin)/20.0
    #    ax.set_aspect(ctx[:aspect])
    ax.grid(true)
    ax.get_yaxis().set_ticks([])
    ax.set_ylim(-5*h,xmax-xmin)
    cmap=region_cmap(ncellregions)

    for icell=1:num_cells(grid)
        ireg=cellregions[icell]
        label = crflag[ireg] ? "c$(ireg)" : ""
        crflag[ireg]=false

        x1=coord[1,cellnodes[1,icell]]
        x2=coord[1,cellnodes[2,icell]]
        ax.plot([x1,x2],[0,0],linewidth=3.0,color=rgbtuple(cmap[cellregions[icell]]),label=label)
        ax.plot([x1,x1],[-h,h],linewidth=ctx[:linewidth],color="k",label="")
        ax.plot([x2,x2],[-h,h],linewidth=ctx[:linewidth],color="k",label="")
    end

    cmap=bregion_cmap(nbfaceregions)
    for ibface=1:num_bfaces(grid)
        ireg=bfaceregions[ibface]
        if ireg >0
            label = brflag[ireg] ? "b$(ireg)" : ""
            brflag[ireg]=false
            x1=coord[1,bfacenodes[1,ibface]]
            ax.plot([x1,x1],[-2*h,2*h],linewidth=3.0,color=rgbtuple(cmap[ireg]),label=label)
        end
    end
    if ctx[:legend]!=:none
        ax.legend(loc=leglocs[ctx[:legend]],ncol=5)
    end
    reveal(ctx,TP)
end

### 2D grid
function gridplot!(ctx, TP::Type{PyPlotType}, ::Type{Val{2}},grid)
    PyPlot=ctx[:Plotter]
    if !haskey(ctx,:ax)
        ctx[:ax]=ctx[:figure].add_subplot(ctx[:layout]...,ctx[:iplot])
    end
    if ctx[:clear]
        ctx[:ax].cla()
        ctx[:ax].set_title(ctx[:title])

        xlimits=ctx[:xlimits]
        ylimits=ctx[:ylimits]
        if xlimits[1]<xlimits[2]
            ctx[:ax].set_xlim(xlimits...)
        end
        if ylimits[1]<ylimits[2]
            ctx[:ax].set_ylim(ylimits...)
        end
    end
    ax=ctx[:ax]
    fig=ctx[:figure]
    cellregions=grid[CellRegions]
    cellnodes=grid[CellNodes]
    coord=grid[Coordinates]
    ncellregions=grid[NumCellRegions]
    nbfaceregions=grid[NumBFaceRegions]
    ncellregions=grid[NumCellRegions]
    if nbfaceregions>0
        bfacenodes=grid[BFaceNodes]
        bfaceregions=grid[BFaceRegions]
    end

    crflag=ones(Bool,ncellregions)
    brflag=ones(Bool,nbfaceregions)
    ax.set_aspect(ctx[:aspect])
    tridat=tridata(grid)
    cmap=region_cmap(ncellregions)
    cdata=ax.tripcolor(tridat...,
                       facecolors=grid[CellRegions],
                       cmap=PyPlot.ColorMap(cmap,length(cmap)),
                       vmin=1.0,
                       vmax=length(cmap)
                       )
    if ctx[:colorbar]==:horizontal
        cbar=fig.colorbar(cdata,ax=ax,ticks=collect(1:length(cmap)),orientation="horizontal")
    end

    if ctx[:colorbar]==:vertical
        cbar=fig.colorbar(cdata,ax=ax,ticks=collect(1:length(cmap)),orientation="vertical")
    end

    ax.triplot(tridat...,color="k",linewidth=ctx[:linewidth])

    if nbfaceregions>0
        cmap=bregion_cmap(nbfaceregions)
        # see https://gist.github.com/gizmaa/7214002
        xc=[coord[:,bfacenodes[1,i]] for i=1:num_sources(bfacenodes)]
        yc=[coord[:,bfacenodes[2,i]] for i=1:num_sources(bfacenodes)]
        rgb=[rgbtuple(cmap[bfaceregions[i]]) for i=1:length(bfaceregions)]
        ax.add_collection(PyPlot.matplotlib.collections.LineCollection(collect(zip(xc,yc)),colors=rgb,linewidth=3))
        for i=1:nbfaceregions
            ax.plot(coord[:,1], coord[:,1],label="$(i)", color=rgbtuple(cmap[i]))
        end
    end
    if ctx[:legend]!=:none
        ax.legend(loc=leglocs[ctx[:legend]])
    end
    reveal(ctx,TP)
end

### 3D Grid
function gridplot!(ctx, TP::Type{PyPlotType}, ::Type{Val{3}},grid)
    # See https://jakevdp.github.io/PythonDataScienceHandbook/04.12-three-dimensional-plotting.html

    PyPlot=ctx[:Plotter]
    if !haskey(ctx,:ax)
        ctx[:ax]=ctx[:figure].add_subplot(ctx[:layout]...,ctx[:iplot],projection="3d")
    end

    ax=ctx[:ax]

    if ctx[:clear]
        ctx[:ax].cla()
    end

    fig=ctx[:figure]

    nregions=num_cellregions(grid)
    nbregions=num_bfaceregions(grid)

    xyzmin=zeros(3)
    xyzmax=ones(3)
    coord=grid[Coordinates]
    @views for idim=1:3
        xyzmin[idim]=minimum(coord[idim,:])
        xyzmax[idim]=maximum(coord[idim,:])
    end

    ax.set_xlim3d(xyzmin[1],xyzmax[1])
    ax.set_ylim3d(xyzmin[2],xyzmax[2])
    ax.set_zlim3d(xyzmin[3],xyzmax[3])
    ax.view_init(ctx[:elev],ctx[:azim])

    cmap=region_cmap(nregions)
    bcmap=bregion_cmap(nbregions)

    xyzcut=[ctx[:xplanes][1],ctx[:yplanes][1],ctx[:zplanes][1]]

    if ctx[:interior]
        regpoints0,regfacets0=extract_visible_cells3D(grid,
                                                      xyzcut,
                                                      primepoints=hcat(xyzmin,xyzmax)
                                                      )
        regfacets=[reshape(reinterpret(Int32,regfacets0[i]),(3,length(regfacets0[i]))) for i=1:nregions]
        regpoints=[reshape(reinterpret(Float32,regpoints0[i]),(3,length(regpoints0[i]))) for i=1:nregions]

        for ireg=1:nregions
            if size(regfacets[ireg],2)>0
                ax.plot_trisurf(regpoints[ireg][1,:],regpoints[ireg][2,:],transpose(regfacets[ireg].-1),regpoints[ireg][3,:],
                                color=rgbtuple(cmap[ireg]),edgecolors=:black,linewidth=0.5)
            end
        end
    end

    bregpoints0,bregfacets0=extract_visible_bfaces3D(grid,
                                                     xyzcut,
                                                     primepoints=hcat(xyzmin,xyzmax)
                                                     )
    bregfacets=[reshape(reinterpret(Int32,bregfacets0[i]),(3,length(bregfacets0[i]))) for i=1:nbregions]
    bregpoints=[reshape(reinterpret(Float32,bregpoints0[i]),(3,length(bregpoints0[i]))) for i=1:nbregions]
    for ireg=1:nbregions
        if size(bregfacets[ireg],2)>0
            ax.plot_trisurf(bregpoints[ireg][1,:],bregpoints[ireg][2,:],transpose(bregfacets[ireg].-1),bregpoints[ireg][3,:],
                            color=rgbtuple(bcmap[ireg]),edgecolors=:black,linewidth=0.5)
        end
    end

    if ctx[:legend]!=:none
        ax.legend(loc=leglocs[ctx[:legend]])
    end
    reveal(ctx,TP)
end

### 1D Function
function scalarplot!(ctx, TP::Type{PyPlotType}, ::Type{Val{1}},grid, func::AbstractVector)
    PyPlot=ctx[:Plotter]
    if !haskey(ctx,:ax)
        ctx[:ax]=ctx[:figure].add_subplot(ctx[:layout]...,ctx[:iplot])
    end
    if ctx[:clear]
        ctx[:ax].cla()
        xlimits=ctx[:xlimits]
        ylimits=ctx[:limits]
        ax=ctx[:ax]

        if xlimits[1]<xlimits[2]
            ax.set_xlim(xlimits...)
        end
        if ylimits[1]<ylimits[2]
            ax.set_ylim(ylimits...)
        end
        ax.grid()
        ax.set_title(ctx[:title])
    end
    ax=ctx[:ax]
    fig=ctx[:figure]

    cellnodes=grid[CellNodes]
    coord=grid[Coordinates]

    pplot=ax.plot
    if ctx[:xscale]==:log
        if ctx[:yscale]==:log
            pplot=ax.loglog
        else
            pplot=ax.semilogx
        end
    end
    if ctx[:yscale]==:log
        if ctx[:xscale]==:log
            pplot=ax.loglog
        else
            pplot=ax.semilogy
        end
    end

    if ctx[:cellwise] # not checked
        for icell=1:num_cells(grid)
            i1=cellnodes[1,icell]
            i2=cellnodes[2,icell]
            x1=coord[1,i1]
            x2=coord[1,i2]
            if icell==1
                ax.plot([x1,x2],[func[i1],func[i2]],color=rgbtuple(ctx[:color]),label=ctx[:label])
            else
                ax.plot([x1,x2],[func[i1],func[i2]],color=rgbtuple(ctx[:color]))
            end
        end
    else
        if ctx[:markershape]==:none
            if ctx[:label]!==""
                pplot(coord[1,:],func,
                        linestyle=lstyles[ctx[:linestyle]],
                        color=rgbtuple(ctx[:color]),
                        linewidth=ctx[:linewidth],
                        label=ctx[:label])
            else
                pplot(coord[1,:],func,
                        linestyle=lstyles[ctx[:linestyle]],
                        linewidth=ctx[:linewidth],
                        color=rgbtuple(ctx[:color]))
            end
        else
            if ctx[:label]!==""
                pplot(coord[1,:],func,
                      linestyle=lstyles[ctx[:linestyle]],
                      color=rgbtuple(ctx[:color]),
                      label=ctx[:label],
                      marker=mshapes[ctx[:markershape]],
                      markevery=ctx[:markevery],
                      markersize=ctx[:markersize],
                      linewidth=ctx[:linewidth]
                      )
            else
                pplot(coord[1,:],func,
                      linestyle=lstyles[ctx[:linestyle]],
                      color=rgbtuple(ctx[:color]),
                      marker=mshapes[ctx[:markershape]],
                      markevery=ctx[:markevery],
                      markersize=ctx[:markersize],
                      linewidth=ctx[:linewidth]
                      )
            end
        end
        # points=[Point2f(coord[1,i],func[i]) for i=1:length(func)]
        # Hard to get this robust, as we need to get axislimits
        # mpoints=markerpoints(points,ctx[:markers],Diagonal([1,1]))
        # ampoints=reshape(reinterpret(Float32,mpoints),(2,length(mpoints)))
        # ax.scatter(ampoints[1,:], ampoints[2,:],color=ctx[:color],label="")

    end
    ax.grid(true)
    ax.set_xlabel(ctx[:xlabel])
    ax.set_ylabel(ctx[:ylabel])
    if ctx[:legend]!=:none
        ax.legend(loc=leglocs[ctx[:legend]])
    end

    reveal(ctx,TP)
end

### 2D Function
function scalarplot!(ctx, TP::Type{PyPlotType}, ::Type{Val{2}},grid, func)
    PyPlot=ctx[:Plotter]
    if !haskey(ctx,:ax)
        ctx[:ax]=ctx[:figure].add_subplot(ctx[:layout]...,ctx[:iplot])
    end
    if ctx[:clear]
        if haskey(ctx,:cbar)
            ctx[:cbar].remove()
        end
        ctx[:ax].remove()
        ctx[:ax]=ctx[:figure].add_subplot(ctx[:layout]...,ctx[:iplot])

        xlimits=ctx[:xlimits]
        ylimits=ctx[:ylimits]
        if xlimits[1]<xlimits[2]
            ctx[:ax].set_xlim(xlimits...)
        end
        if ylimits[1]<ylimits[2]
            ctx[:ax].set_ylim(ylimits...)
        end
    end

    ax=ctx[:ax]
    fig=ctx[:figure]
    ax.set_aspect(ctx[:aspect])
    ax.set_title(ctx[:title])

    isolines,crange=isolevels(ctx,func)
    colorlevels=collect(crange[1]:(crange[2]-crange[1])/(ctx[:colorlevels]-1):crange[2])

    if !haskey(ctx,:grid) || !seemingly_equal(ctx[:grid],grid)
        ctx[:grid]=grid
        ctx[:tridata]=tridata(grid)
    end
    cnt=ax.tricontourf(ctx[:tridata]...,func;levels=colorlevels,cmap=PyPlot.ColorMap(plaincolormap(ctx)))
    for c in cnt.collections
        c.set_edgecolor("face")
    end
    ax.tricontour(ctx[:tridata]...,func,colors="k",levels=isolines)

    if ctx[:colorbarticks] == :default
        ticks = isolines
    elseif isa(ctx[:colorbarticks],Number)
        ticks = collect(crange[1]:(crange[2]-crange[1])/(ctx[:colorbarticks]-1):crange[2])
    else
        ticks = ctx[:colorbarticks]
    end

    if ctx[:colorbar]==:horizontal
        ctx[:cbar]=fig.colorbar(cnt,ax=ax,ticks=ticks,boundaries=colorlevels, orientation="horizontal")
    end

    if ctx[:colorbar]==:vertical
        ctx[:cbar]=fig.colorbar(cnt,ax=ax,ticks=ticks,boundaries=colorlevels, orientation="vertical")
    end

    ax.set_xlabel(ctx[:xlabel])

    ax.set_ylabel(ctx[:ylabel])

    reveal(ctx,TP)
end

function scalarplot!(ctx, TP::Type{PyPlotType}, ::Type{Val{3}},grid,func)

    PyPlot=ctx[:Plotter]
    if !haskey(ctx,:ax)
        ctx[:ax]=ctx[:figure].add_subplot(ctx[:layout]...,ctx[:iplot],projection="3d")
    end
    ax=ctx[:ax]
    fig=ctx[:figure]

    nregions=num_cellregions(grid)
    nbregions=num_bfaceregions(grid)

    xyzmin=zeros(3)
    xyzmax=ones(3)
    coord=grid[Coordinates]
    @views for idim=1:3
        xyzmin[idim]=minimum(coord[idim,:])
        xyzmax[idim]=maximum(coord[idim,:])
    end
    xyzcut=[ctx[:xplanes],ctx[:yplanes],ctx[:zplanes]]

    levels,crange=isolevels(ctx,func)

    planes=makeplanes(xyzmin,xyzmax,ctx[:xplanes],ctx[:yplanes],ctx[:zplanes])

    ccoord0,faces0,values=marching_tetrahedra(grid,func,planes,levels)

    faces=reshape(reinterpret(Int32,faces0),(3,length(faces0)))
    ccoord=reshape(reinterpret(Float32,ccoord0),(3,length(ccoord0)))

    nfaces=size(faces,2)
    if nfaces>0
        colors=zeros(nfaces)
        for i=1:nfaces
            colors[i]=(values[faces[1,i]]+values[faces[2,i]]+values[faces[3,i]])/3
        end
        # thx, https://stackoverflow.com/a/24229480/8922290
        collec=ctx[:ax].plot_trisurf(ccoord[1,:],ccoord[2,:],transpose(faces.-1),ccoord[3,:],
                                     cmap=PyPlot.ColorMap(plaincolormap(ctx)),
                                     vmin=crange[1],
                                     vmax=crange[2])
        collec.set_array(colors)
        collec.autoscale()
    end

    ax.set_xlim3d(xyzmin[1],xyzmax[1])
    ax.set_ylim3d(xyzmin[2],xyzmax[2])
    ax.set_zlim3d(xyzmin[3],xyzmax[3])
    ax.view_init(ctx[:elev],ctx[:azim])

    ax.set_title(ctx[:title])

    if ctx[:legend]!=:none
        ax.legend(loc=leglocs[ctx[:legend]])
    end
    reveal(ctx,TP)
end

### 2D Vector
function vectorplot!(ctx, TP::Type{PyPlotType}, ::Type{Val{2}},grid, func)
    PyPlot=ctx[:Plotter]
    if !haskey(ctx,:ax)
        ctx[:ax]=ctx[:figure].add_subplot(ctx[:layout]...,ctx[:iplot])
    end
    if ctx[:clear]
        if haskey(ctx,:cbar)
            ctx[:cbar].remove()
        end
        ctx[:ax].remove()
        ctx[:ax]=ctx[:figure].add_subplot(ctx[:layout]...,ctx[:iplot])

        xlimits=ctx[:xlimits]
        ylimits=ctx[:ylimits]
        if xlimits[1]<xlimits[2]
            ctx[:ax].set_xlim(xlimits...)
        end
        if ylimits[1]<ylimits[2]
            ctx[:ax].set_ylim(ylimits...)
        end
    end

    ax=ctx[:ax]
    fig=ctx[:figure]
    ax.set_aspect(ctx[:aspect])
    ax.set_title(ctx[:title])

    rc,rv=vectorsample(grid,func,spacing=ctx[:spacing], offset=ctx[:offset])
    qc,qv=quiverdata(rc,rv;vscale=ctx[:vscale],vnormalize=ctx[:vnormalize])

    ax.quiver(qc[1,:], qc[2,:], qv[1,:],qv[2,:])
    ax.set_xlabel(ctx[:xlabel])
    ax.set_ylabel(ctx[:ylabel])

    reveal(ctx,TP)
end

### 2D stream
function streamplot!(ctx, TP::Type{PyPlotType}, ::Type{Val{2}},grid, func)
    PyPlot=ctx[:Plotter]
    if !haskey(ctx,:ax)
        ctx[:ax]=ctx[:figure].add_subplot(ctx[:layout]...,ctx[:iplot])
    end
    if ctx[:clear]
        if haskey(ctx,:cbar)
            ctx[:cbar].remove()
        end
        ctx[:ax].remove()
        ctx[:ax]=ctx[:figure].add_subplot(ctx[:layout]...,ctx[:iplot])
    end

    ax=ctx[:ax]
    fig=ctx[:figure]
    ax.set_aspect(ctx[:aspect])
    ax.set_title(ctx[:title])

    xlimits=ctx[:xlimits]
    ylimits=ctx[:ylimits]
    if xlimits[1]<xlimits[2]
        ax.set_xlim(xlimits...)
    end
    if ylimits[1]<ylimits[2]
        ax.set_ylim(ylimits...)
    end

    # thx https://discourse.julialang.org/t/meshgrid-function-in-julia/48679/4?u=j-fu
    function meshgrid(rc)
        nx=length(rc[1])
        ny=length(rc[2])
        xout=zeros(ny,nx)
        yout=zeros(ny,nx)
        for jx=1:nx
            for ix=1:ny
                xout[ix,jx]=rc[1][jx]
                yout[ix,jx]=rc[2][ix]
            end
        end
        xout,yout
    end

    rc,rv=vectorsample(grid,func,spacing=ctx[:spacing], offset=ctx[:offset])

    X,Y=meshgrid(rc)

    ax.streamplot(X,Y, rv[1,:,:,1],rv[2,:,:,1],color=ctx[:color])
    ax.set_xlabel(ctx[:xlabel])
    ax.set_ylabel(ctx[:ylabel])

    reveal(ctx,TP)
end


text:
import ValueHistories: push!

function checkstatus(::NoStatus, z)
    return false
end

function logextra(::NoStatus, z, override=false)
    return false
end

printstatusheader(::NoStatus) = nothing


text:
__precompile__()
module SpikeSorting

using Interpolations, DistributedArrays, Gtk.ShortNames, Cairo, MultivariateStats

using SharedArrays, LinearAlgebra

if VERSION > v"0.7-"
    const Void = Nothing
    const setproperty! = set_gtk_property!
    const getproperty = get_gtk_property
    const linspace(x,y,z) = range(x,stop=y,length=z)
    const indmin = argmin
    const method_exists = hasmethod
end

export Sorting,Spike,create_multi,output_buffer,benchmark, SortView

include("types.jl")
include("constants.jl")
include("cluster.jl")
include("gui/gui_types.jl")
include("gui/template_cluster.jl")
include("gui/threshold.jl")
include("gui/cluster_treeview.jl")
include("gui/multi_dim_view.jl")
include("gui/single_channel.jl")
include("gui/drawing_helpers.jl")
include("gui/isi_canvas.jl")
include("gui/offline_gui.jl")
include("detect.jl")
include("align.jl")
include("feature.jl")
include("reduction.jl")
include("run.jl")
include("threshold.jl")
include("benchmark.jl")

end


text:
import Base.+
using Test
import Test: AbstractTestSet, record, finish, get_testset_depth, get_testset, Broken, Pass, Fail, Error, TestSetException
include("types.jl")

@MaracasTestSet("DescribeTestSet")
@MaracasTestSet("SpecTestSet")
@MaracasTestSet("TestTestSet")
format_title(::Type{DescribeTestSet}, desc) = string(MARACAS_SETTING[:title], desc, MARACAS_SETTING[:default])
format_title(::Type{SpecTestSet}, desc) = string(MARACAS_SETTING[:spec], "[Spec] ", MARACAS_SETTING[:default], "it ", desc)
format_title(::Type{TestTestSet}, desc) = string(MARACAS_SETTING[:test], "[Test] ", MARACAS_SETTING[:default], desc)

ResultsCount(ts) = ResultsCount(0, 0, 0, 0)
total(count::ResultsCount)  = count.passes + count.fails + count.errors + count.broken
+(a::ResultsCount, b::ResultsCount) = ResultsCount(a.passes + b.passes, a.fails + b.fails, a.errors + b.errors, a.broken + b.broken)
+(a::ResultsCount, b::Pass) = ResultsCount(a.passes + 1, a.fails, a.errors, a.broken)
+(a::ResultsCount, b::Fail) = ResultsCount(a.passes, a.fails + 1, a.errors, a.broken)
+(a::ResultsCount, b::Error) = ResultsCount(a.passes, a.fails, a.errors + 1, a.broken)
+(a::ResultsCount, b::Broken) = ResultsCount(a.passes, a.fails, a.errors, a.broken + 1)
+(a::ResultsCount, b::AbstractTestSet) = (a + ResultsCount(b))

passes(count::Dict) = get(count, :passes, 0)
fails(count::Dict)  = get(count, :fails, 0)
errors(count::Dict) = get(count, :errors, 0)
broken(count::Dict) = get(count, :broken, 0)
total(count::Dict)  = passes(count) + fails(count) + errors(count) + broken(count)
# Backtrace utility functions copied from Test.jl
function ip_has_file_and_func(ip, file, funcs)
    return any(fr -> (occursin(file, string(fr.file)) && fr.func in funcs), StackTraces.lookup(ip))
end

function scrub_backtrace(bt)
    do_test_ind = findfirst(ip -> ip_has_file_and_func(ip, r"[Tt]est.jl$", (:do_test, :do_test_throws)), bt)
    if do_test_ind !== nothing && length(bt) > do_test_ind
        bt = bt[do_test_ind + 1:end]
    end
    name_ind = findfirst(ip -> ip_has_file_and_func(ip, r"Maracas.jl", (:describe, )), bt)
    if name_ind !== nothing && length(bt) != 0
        bt = bt[1:name_ind]
    end
    return bt
end

ResultsCount(ts::MaracasTestSet) = ts.count

passes(ts::MaracasTestSet) = ts.count.passes
fails(ts::MaracasTestSet)  = ts.count.fails
errors(ts::MaracasTestSet) = ts.count.errors
broken(ts::MaracasTestSet) = ts.count.broken
total(ts::MaracasTestSet)  = ts.count.passes + ts.count.fails + ts.count.errors + ts.count.broken

failed(ts::MaracasTestSet) = (errors(ts) > 0 || fails(ts) > 0 )

function hwidth(ts::MaracasTestSet)
    return Dict(
    :passes => passes(ts) > 0 ? 8 : 0,
    :fails => fails(ts) > 0 ? 8 : 0,
    :errors => errors(ts) > 0 ? 8 : 0,
    :broken => broken(ts) > 0 ? 8 : 0,
    :total => total(ts) > 0 ? 8 : length("No Tests")
    )
end

# For a broken result, simply store the result
function record(ts::MaracasTestSet, t::Broken)
    ts.count += t;
    push!(ts.results, t)
    return t
end
# For a passed result, do not store the result since it uses a lot of memory
function record(ts::MaracasTestSet, t::Pass)
    ts.count += t;
    return t
end
# For the other result types, immediately print the error message
# but do not terminate. Print a backtrace.
function record(ts::MaracasTestSet, t::Union{Fail, Error})
    ts.count += t;
    if myid() == 1
        printstyled(ts.description, ": ", bold=true)
        print(t)
        # don't print the backtrace for Errors because it gets printed in the show
        # method
        isa(t, Error) || Base.show_backtrace(stdout, scrub_backtrace(backtrace()))
        println()
    end
    push!(ts.results, t)
    t, isa(t, Error) || backtrace()
end

function record(ts::MaracasTestSet, t::MaracasTestSet)
    ts.count += t;
    ts.max_depth = t.max_depth + 1
    push!(ts.results, t)
end

print_test_errors(ts::MaracasTestSet) = map(print_test_errors, ts.results)
function print_test_errors(t::Union{Fail, Error})
    if myid() == 1
        println("Error in testset $(ts.description):")
        Base.show(stdout, t)
        println()
    end
end
print_test_errors(t) = nothing

function print_col_header(color::Symbol, title::AbstractString, result::Int)
    result == 0 && return
    printstyled(lpad(title, 8); bold=true, color=color)
end

function print_test_results(ts::MaracasTestSet)
    print_summary()
    print_col_header(MARACAS_SETTING[:pass], "Pass", passes(ts))
    print_col_header(MARACAS_SETTING[:error], "Fail", fails(ts))
    print_col_header(MARACAS_SETTING[:error], "Error", errors(ts))
    print_col_header(MARACAS_SETTING[:warn], "Broken", broken(ts))
    print_col_header(MARACAS_SETTING[:info], "Total", total(ts))
    println()
    # Recursively print a summary at every level
    print_counts(ts, 0, hwidth(ts))
end

const TESTSET_PRINT_ENABLE = Ref(true)

# Called at the end of a @testset, behaviour depends on whether
# this is a child of another testset, or the "root" testset
function finish(ts::MaracasTestSet)
    # If we are a nested test set, do not print a full summary
    # now - let the parent test set do the printing
    if get_testset_depth() != 0
        # Attach this test set to the parent test set
        parent_ts = get_testset()
        record(parent_ts, ts)
        return ts
    end

    if TESTSET_PRINT_ENABLE[]
        print_test_results(ts)
    end

    # Finally throw an error as we are the outermost test set
    if failed(ts)
        # Get all the error/failures and bring them along for the ride
        throw(TestSetException(passes(ts), fails(ts), errors(ts), broken(ts), filter_errors(ts)))
    end

    # return the testset so it is returned from the @testset macro
    ts
end

# Recursive function that fetches backtraces for any and all errors
# or failures the testset and its children encountered
function filter_errors(ts::MaracasTestSet)
    efs = []
    for t in ts.results
        if isa(t, MaracasTestSet)
            append!(efs, filter_errors(t))
        elseif isa(t, Union{Fail, Error})
            append!(efs, [t])
        end
    end
    efs
end

function print_result_column(color, result, width)
    if result > 0
        printstyled(lpad(string(result), width, " "), "  "; bold=true, color=color)
    elseif width > 0
        print(lpad(" ", width), "  ")
    end
end
# Recursive function that prints out the results at each level of
# the tree of test sets
function print_counts(ts::MaracasTestSet, depth, headers_width)
    print_title(ts.description, depth)
    print_result_column(MARACAS_SETTING[:pass], passes(ts), passes(headers_width))
    print_result_column(MARACAS_SETTING[:error], fails(ts), fails(headers_width))
    print_result_column(MARACAS_SETTING[:error], errors(ts), errors(headers_width))
    print_result_column(MARACAS_SETTING[:warn], broken(ts), broken(headers_width))

    subtotal = total(ts)
    result = subtotal == 0 ? "No tests" : string(subtotal)
    printstyled(lpad(result, total(headers_width), " "); bold=true, color=MARACAS_SETTING[:info])
    println()

    for t in ts.results
        print_counts(t, depth + 1, headers_width)
    end
end
print_counts(args...) = nothing

caesura(text, quantity) = string(text[1:(end + quantity - 4)], Base.text_colors[:red], "... ")

rm_spec_char(text) = replace(text, r"\e\[[0-9;]+m" => "")
function rpad_title(text)
    space_repeat = MARACAS_SETTING[:title_length] - length(rm_spec_char(text))
    return space_repeat > 0 ? string(text, " "^space_repeat) : caesura(text, space_repeat)
end

print_title(text, depth=0) = print(rpad_title(string("  "^depth, text)), "$(MARACAS_SETTING[:default])|")
print_summary() = print_title("$(MARACAS_SETTING[:bold])Test Summary:")


text:
# ---
# title: 1047. Remove All Adjacent Duplicates In String
# id: problem1047
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Stack
# link: <https://leetcode.com/problems/remove-all-adjacent-duplicates-in-string/description/>
# hidden: true
# ---
#
# Given a string `S` of lowercase letters, a _duplicate removal_ consists of
# choosing two adjacent and equal letters, and removing them.
#
# We repeatedly make duplicate removals on S until we no longer can.
#
# Return the final string after all such duplicate removals have been made.  It
# is guaranteed the answer is unique.
#
#
#
# **Example 1:**
#
#
#
#     Input: "abbaca"
#     Output: "ca"
#     Explanation:
#     For example, in "abbaca" we could remove "bb" since the letters are adjacent and equal, and this is the only possible move.  The result of this move is that the string is "aaca", of which only "aa" is possible, so the final string is "ca".
#
#
#
#
# **Note:**
#
#   1. `1 <= S.length <= 20000`
#   2. `S` consists only of English lowercase letters.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using NeutronTransport
import Gridap: DiscreteModelFromFile

jsonfile = joinpath(@__DIR__, "pincell.json")
geometry = DiscreteModelFromFile(jsonfile)

# number of azimuthal angles
nφ = 32

# azimuthal spacing
δ = 1e-2

# boundary conditions
bcs = BoundaryConditions(top=Reflective, bottom=Reflective, left=Reflective, right=Reflective)
# bcs = BoundaryConditions(top=Vaccum, bottom=Vaccum, left=Vaccum, right=Vaccum)
# bcs = BoundaryConditions(top=Periodic, bottom=Periodic, left=Periodic, right=Periodic)

# initialize track generator
tg = TrackGenerator(geometry, nφ, δ, bcs=bcs)

# perform ray tracing
trace!(tg)

# proceed to segmentation
segmentize!(tg)

# polar quadrature
pq = TabuchiYamamoto(6)

# materials
pin = CrossSections("pin", 2;
    νΣf = [1.86278e-2, 3.44137e-1],
    Σt  = [3.62022e-1, 5.72155e-1],
    Σs0 = [3.33748e-1  6.64881e-4; 0.0e-0 3.80898e-1]
)

cladding = CrossSections("cladding", 2;
    Σt  = [2.74144e-1, 2.80890e-1],
    Σs0 = [2.72377e-1 1.90838e-4; 0.0e-0 2.77230e-1]
)

water = CrossSections("water", 2;
    Σt  = [6.40711e-1, 1.69131e-0],
    Σs0 = [6.07382e-1 3.31316e-2; 0.0e-0 1.68428e-0]
)

xs = [pin, cladding, water]

# define the problem
prob = MoCProblem(tg, pq, xs)

# solve
sol = solve(prob, debug=true)

import Gridap: writevtk
import Gridap.Geometry: get_triangulation
trian = get_triangulation(tg.mesh.model)
writevtk(trian, "fluxes-pincell", cellfields=["g1" => sol(1), "g2" => sol(2)])


text:
struct InferenceTiming
    mi_info::Core.Compiler.Timings.InferenceFrameInfo
    inclusive_time::Float64
    exclusive_time::Float64
end
"""
    inclusive(frame)

Return the time spent inferring `frame` and its callees.
"""
inclusive(it::InferenceTiming) = it.inclusive_time
"""
    exclusive(frame)

Return the time spent inferring `frame`, not including the time needed for any of its callees.
"""
exclusive(it::InferenceTiming) = it.exclusive_time

struct InferenceTimingNode
    mi_timing::InferenceTiming
    start_time::Float64
    children::Vector{InferenceTimingNode}
    bt
    parent::InferenceTimingNode

    # Root constructor
    InferenceTimingNode(mi_timing::InferenceTiming, start_time, @nospecialize(bt)) =
        new(mi_timing, start_time, InferenceTimingNode[], bt)
    # Child constructor
    function InferenceTimingNode(mi_timing::InferenceTiming, start_time, @nospecialize(bt), parent::InferenceTimingNode)
        child = new(mi_timing, start_time, InferenceTimingNode[], bt, parent)
        push!(parent.children, child)
        return child
    end
end
inclusive(node::InferenceTimingNode) = inclusive(node.mi_timing)
exclusive(node::InferenceTimingNode) = exclusive(node.mi_timing)
InferenceTiming(node::InferenceTimingNode) = node.mi_timing

function InferenceTimingNode(t::Core.Compiler.Timings.Timing)
    ttree = timingtree(t)
    it, start_time, ttree_children = ttree::Tuple{InferenceTiming, Float64, Vector{Any}}
    root = InferenceTimingNode(it, start_time, t.bt)
    addchildren!(root, t, ttree_children)
    return root
end

# Compute inclusive times and store as a temporary tree.
# To allow InferenceTimingNode to be both bidirectional and immutable, we need to create parent node before the child nodes.
# However, each node stores its inclusive time, which can only be computed efficiently from the leaves up (children before parents).
# This performs the inclusive-time computation, storing the result as a "temporary tree" that can be used during
# InferenceTimingNode creation (see `addchildren!`).
function timingtree(t::Core.Compiler.Timings.Timing)
    time, start_time = t.time/10^9, t.start_time/10^9
    incl_time = time
    tchildren = []
    for child in t.children
        tchild = timingtree(child)
        push!(tchildren, tchild)
        incl_time += inclusive(tchild[1])
    end
    return (InferenceTiming(t.mi_info, incl_time, time), start_time, tchildren)
end

function addchildren!(parent::InferenceTimingNode, t::Core.Compiler.Timings.Timing, ttrees)
    for (child, ttree) in zip(t.children, ttrees)
        it, start_time, ttree_children = ttree::Tuple{InferenceTiming, Float64, Vector{Any}}
        node = InferenceTimingNode(it, start_time, child.bt, parent)
        addchildren!(node, child, ttree_children)
    end
end

function start_deep_timing()
    Core.Compiler.Timings.reset_timings()
    Core.Compiler.__set_measure_typeinf(true)
end
function stop_deep_timing()
    Core.Compiler.__set_measure_typeinf(false)
    Core.Compiler.Timings.close_current_timer()
end

function finish_snoopi_deep()
    return InferenceTimingNode(Core.Compiler.Timings._timings[1])
end

function _snoopi_deep(cmd::Expr)
    return quote
        start_deep_timing()
        try
            $(esc(cmd))
        finally
            stop_deep_timing()
        end
        finish_snoopi_deep()
    end
end

"""
    tinf = @snoopi_deep commands

Produce a profile of julia's type inference, recording the amount of time spent inferring
every `MethodInstance` processed while executing `commands`. Each fresh entrance to
type inference (whether executed directly in `commands` or because a call was made
by runtime-dispatch) also collects a backtrace so the caller can be identified.

`tinf` is a tree, each node containing data on a particular inference "frame" (the method,
argument-type specializations, parameters, and even any constant-propagated values).
Each reports the [`exclusive`](@ref)/[`inclusive`](@ref) times, where the exclusive
time corresponds to the time spent inferring this frame in and of itself, whereas
the inclusive time includes the time needed to infer all the callees of this frame.

The top-level node in this profile tree is `ROOT`. Uniquely, its exclusive time
corresponds to the time spent _not_ in julia's type inference (codegen, llvm_opt, runtime, etc).

There are many different ways of inspecting and using the data stored in `tinf`.
The simplest is to load the `AbstracTrees` package and display the tree with
`AbstractTrees.print_tree(tinf)`.
See also:  `flamegraph`, `flatten`, `inference_triggers`, `SnoopCompile.parcel`,
`runtime_inferencetime`.

# Example
```jldoctest; setup=:(using SnoopCompile), filter=r"([0-9]*\\.?[0-9]+([eE][-+]?[0-9]+)?/[0-9]*\\.?[0-9]+([eE][-+]?[0-9]+)?|\\d direct)"
julia> tinf = @snoopi_deep begin
           sort(rand(100))  # Evaluate some code and profile julia's type inference
       end
InferenceTimingNode: 0.110018224/0.131464476 on InferenceFrameInfo for Core.Compiler.Timings.ROOT() with 2 direct children
```

"""
macro snoopi_deep(cmd)
    return _snoopi_deep(cmd)
end

# These are okay to come at the top-level because we're only measuring inference, and
# inference results will be cached in a `.ji` file.
precompile(start_deep_timing, ())
precompile(stop_deep_timing, ())
precompile(finish_snoopi_deep, ())


text:
# function plotmesh(m :: TriMesh;
#                         linewidth :: Real = 1,
#                         marker :: String = "None",
#                         markersize :: Real = 10,
#                         color = :black)
#     p = plot()
#     for edge in eachcol(m.edge)
#         xy = m.point[:, edge]
#         plot!(p, xy[1, :], xy[2, :], linewidth=linewidth,
#             markersize=markersize, color=color, legend=false)
#     end
#     return p
# end
#

"""
    make_border_points(points, n, expansion, k = size(points, 2)) -> Matrix

Create a set of border points around the hull of a domain.

# Arguments
- `points`: points to form a hull around, as a matrix with a point in each column
- `n`: number of edge points
- `expansion`: factor to extend the hull beyond the provided points
- `k`: number of neighbors to use when determining hull. Default is all points, giving a convex hull; can be decreased to make the hull concave.

# Returns
- `Matrix`: coordinates of domain border
"""
function make_border_points(points, n, expansion, k=size(points, 2))
    hull = concave_hull(collect(eachcol(points)), k)
    cm = mean(points, dims=2)
    points2 = cm .+ 1.01 * (points .- cm)
    edge_points = reduce(hcat, [p for p in eachcol(points2) if ! in_hull(p, hull)])

    θ = [atan(y, x) for (x, y) in eachcol(edge_points .- cm)]
    x = edge_points[1, sortperm(θ)]
    y = edge_points[2, sortperm(θ)]
    θ = [θ; θ[1]]
    x = [x; x[1]]
    y = [y; y[1]]

    ds = vec(sqrt.(sum(diff([x y]', dims=2).^2, dims=1)))
    s = cumsum([0; ds])
    xinterp = LinearInterpolation(s, x)
    yinterp = LinearInterpolation(s, y)

    Δs = (s[end] - s[1]) / n
    ss = s[1]:Δs:s[end]
    border_points = [xinterp(ss) yinterp(ss)]'
    border_points = (border_points .- cm) * expansion .+ cm
end

"""
    inverse_distance_weights(dists::AbstractVector) -> Vector

Calculate the relative weights via inverse distance.

# Arguments
- `dists::AbstractVector`: vector of distances

# Returns
- `Vector`: relative weight of points at each distance
"""
function inverse_distance_weights(dists::AbstractVector)
    w = [d > 0 ? 1/d : 1/(d+eps()) for d in dists]
    return w / sum(w)
end

"""
    observation_matrix(mesh::TriMesh, points::AbstractMatrix) -> SparseMatricCSC

Create matrix of weights to linearly interpolate from nodes to observation locations.

# Arguments
- `mesh::TriMesh`: a mesh from TriangleMesh.jl
- `points::AbstractMatrix`: observation locations in columns

# Returns
- `SparseMatrixCSC`: sparse matrix with weights
"""
function observation_matrix(mesh::TriMesh, points::AbstractMatrix)
    npoint = size(points, 2)
    nmesh = size(mesh.point, 2)
    tree = KDTree(mesh.point)
    idx, dists = knn(tree, points, 3)

    ii = repeat(1:npoint, inner=3)
    jj = reduce(vcat, idx)
    ww = mapreduce(inverse_distance_weights, vcat, dists)
    return sparse(ii, jj, ww, npoint, nmesh)
end


text:
# The Computer Language Benchmarks Game
# https://salsa.debian.org/benchmarksgame-team/benchmarksgame/
#

const line_width = 60

const alu = string(
   "GGCCGGGCGCGGTGGCTCACGCCTGTAATCCCAGCACTTTGG",
   "GAGGCCGAGGCGGGCGGATCACCTGAGGTCAGGAGTTCGAGA",
   "CCAGCCTGGCCAACATGGTGAAACCCCGTCTCTACTAAAAAT",
   "ACAAAAATTAGCCGGGCGTGGTGGCGCGCGCCTGTAATCCCA",
   "GCTACTCGGGAGGCTGAGGCAGGAGAATCGCTTGAACCCGGG",
   "AGGCGGAGGTTGCAGTGAGCCGAGATCGCGCCACTGCACTCC",
   "AGCCTGGGCGACAGAGCGAGACTCCGTCTCAAAAA")

const iub1 = b"acgtBDHKMNRSVWY"
const iub2 = [0.27, 0.12, 0.12, 0.27, 0.02,0.02, 0.02, 0.02, 0.02, 0.02,0.02, 0.02, 0.02, 0.02, 0.02]

const homosapiens1 = b"acgt"
const homosapiens2 = [0.3029549426680, 0.1979883004921,0.1975473066391, 0.3015094502008]

const IM  = 139968.0
const IA  =   3877.0
const IC  =  29573.0

function gen_random()
    global rng_state = ((rng_state::Float64 * IA + IC) % IM)
end

function repeat_fasta(src, n)
    k = length(src)
    s = string(src, src, src[1:(n % k)])
    I = Iterators.cycle(src)
    col = 1
    count = 1
    c, state = iterate(I)
    print(c)
    while count < n
        col += 1
        c, state = iterate(I, state)
        print(c)
        if col == line_width
            println()
            col = 0
        end
        count += 1
    end
    println()
    return
end

function choose_char(cs)
    k = length(cs)
    r = gen_random() / IM
    r < cs[1] && return 1
    a = 1
    b = k
    while b > a + 1
        c = fld(a + b, 2)
        if r < cs[c]
            b = c
        else
            a = c
        end
    end
    return b
end

function random_fasta(symb, pr, n)
    cs = cumsum(pr)
    line = Vector{UInt8}(undef, line_width)
    k = n
    while k > 0
        m = min(k, line_width)
        resize!(line, m)
        for i = 1:m
            line[i] = symb[choose_char(cs)]
        end
        println(String(copy(line)))
        k -= line_width
    end
    return
end

rng_state = 42.0

function perf_fasta(n=25000000)
  println(">ONE Homo sapiens alu")
  repeat_fasta(alu, 2n)

  println(">TWO IUB ambiguity codes")
  random_fasta(iub1, iub2, 3n)
  println(">THREE Homo sapiens frequency")
  random_fasta(homosapiens1, homosapiens2, 5n)
end

n = parse(Int,ARGS[1])
perf_fasta(n)


text:
# ------- R 2 Julia Conversion ---------------
#function saveModels(d::Dict,fname::String)
#    j = JSON.json([ d])
#    open(fname,"w") do f
#        write(f,j)
#    end
#end

function expandM(di::Dict)
    d2=deepcopy(di)
    if d2[:modelName]==:occ
        d2[:dist] = Poisson()
        d2[:lnk] = LogLink()
        d2[:Buyer_Pos_P1_is1] = true
    end
    if d2[:modelName]==:dolocc
        d2[:dist] = Gamma()
        d2[:lnk] = LogLink()
        d2[:Buyer_Pos_P1_is1] = true
    end
    if d2[:modelName]==:pen
        d2[:dist] = Bernoulli()
        d2[:lnk] = LogitLink()
        d2[:Buyer_Pos_P1_is1] = false
    end
    return d2
end

#function Dict_Str2Sym(d::Dict,keys::Array{Any}=Any[])
#    keys = length(keys)==0 ? keys(d) : keys
#    dout=Dict()
#    for (k, v) in d #Dict(k=>d[k] for k in keys)
#        if k in keys
#            dout[Symbol(k)] = isString(v) ? Symbol(v) : isArray(v) ? [Symbol(vx) for vx in v] : isDict(v) ? Dict_Str2Sym(v) : v
#        else
#            dout[k] = v
#        end
#    end
#    dout
#end

function mergeDict(dnew::Dict, dold::OrderedDict)
    dout=deepcopy(dold)
    for (k,v) in dnew dout[k] = v end
    dout
end

function dict_Sym(d::Dict)
           kys = keys(d)
           dout=Dict()
           for (k, v) in d #Dict(k=>d[k] for k in kys)
              dout[Symbol(k)] = isString(v) ? Symbol(v) : isArray(v) ? [Symbol(vx) for vx in v] : isDict(v) ? Dict_Sym(v) : v
           end
           dout
 end

function readModels(fname::String)
    if isfile(fname)
        d=Dict()
        j = JSON.parsefile(fname)[1]
        #for (key, value) in dict
        for (k, v) in j
            d2=Dict()
            if Symbol(k)==:mocc
                d2[:dist] = Poisson()
                d2[:lnk] = LogLink()
                d2[:Buyer_Pos_P1_is1] = true
            end
            if Symbol(k)==:mdolocc
                d2[:dist] = Gamma()
                d2[:lnk] = LogLink()
                d2[:Buyer_Pos_P1_is1] = true
            end
            if Symbol(k)==:mpen
                d2[:dist] = Bernoulli()
                d2[:lnk] = LogitLink()
                d2[:Buyer_Pos_P1_is1] = false
            end

            if typeof(v) == Array{Any,1}
                d[Symbol(k)] = [Symbol(vx) for vx in v]
            else
                for (k1, v1) in v
                    if typeof(v1) in [ Array{Any,1} ]
                        v2=[Symbol(v1x) for v1x in v1]
                    else
                        v2=Symbol(v1)
                    end
                    d2[Symbol(k1)] = v2
                end
                d[Symbol(k)] = d2
            end
        end
        for (k,v) in d
            if k in [:dolocc, :pen, :occ]
                d[k] = expandM(v)
            end
        end
        return d
    else
        return "file not found!"
    end
end

function read_orig2(root::String)
    push!(Base.DL_LOAD_PATH, pwd())
    fread(fname::String)=unsafe_string(ccall((:xfread,"fread"), Cstring, (Cstring,), fname))
    #df_data = readtable(root*"/orig.csv",header=false);
    if isfile(root*"/origHead.csv")
        df_data = readtable(IOBuffer(fread(root*"/orig.csv")),header=false)
        df_h = readtable(root*"/origHead.csv",header=false);
        names!(df_data, convert(Array{Symbol}, df_h[:x1]) )
    else
        df_data = readtable(IOBuffer(fread(root*"/orig.csv")))
    end
    return df_data
end
""" C Code
#include <stdlib.h>
#include <string.h>
/*gcc fread.c -shared -fPIC -o fread.so
  push!(Base.DL_LOAD_PATH, pwd())
  x=ccall((:xmean,"/mapr/mapr04p/analytics0001/analytic_users/jmdl/Hormel_Rev_8_967/fread"),Float64,(Float64,Float64),2.0,5.0)

  x=ccall((:myFunction,"/mapr/mapr04p/analytics0001/analytic_users/jmdl/Hormel_Rev_8_967/fread"),String)

  push!(Base.DL_LOAD_PATH, pwd())
  x=ccall((:myFunction,"fread"), Cstring, (Cstring,), "tst.csv")
*/
/*const char* xfread()  / * http://www.fundza.com/c4serious/fileIO_reading_all/  */
/*   http://docs.julialang.org/en/stable/manual/calling-c-and-fortran-code/
 * Julia Code *
  x=ccall((:xfread,"fread"), Cstring, (Cstring,), "tst.csv")
  unsafe_string(x)
  fread(fname::String)=unsafe_string(ccall((:xfread,"fread"), Cstring, (Cstring,), fname)) */
const char* xfread(char* fname)
{   FILE    *infile; /* declare a file pointer */
    char    *buffer;
    long    numbytes;
    infile = fopen(fname, "r");   /* open an existing file for reading */
    if(infile == NULL)  /* quit if the file does not exist */
        return "NULL"; //return 1;
    fseek(infile, 0L, SEEK_END); /* Get the number of bytes */
    numbytes = ftell(infile);
    printf("The file called orig.csv s size %ld\n", numbytes);

    fseek(infile, 0L, SEEK_SET);        /* reset the file position indicator to * the beginning of the file */
    buffer = (char*)calloc(numbytes, sizeof(char));      /* grab sufficient memory for the  buffer to hold the text */
    if(buffer == NULL)  /* memory error */
        return "NULL"; //return 1;
    fread(buffer, sizeof(char), numbytes, infile);  /* copy all the text into the buffer */
    fclose(infile);
    //printf("The file called orig.csv contains this text\n\n%s", buffer); /* confirm we have read the file by outputing it to the console */
    //free(buffer);  /* free the memory we used for the buffer */ /*xxzreturn buffer;*/
    //return 0;
    return buffer;
}
double xmean(double a, double b) { return (a+b) / 2; }
const char* myFunction(char* fname) {  printf("working on %s\n\n",fname); return "done with My String"; }
//gcc fread.c -shared -fPIC -o fread.so; gcc fread.c; ./a.out // x=unsafe_string(ccall((:xfread,"fread"), Cstring, (Cstring,), "tst.csv"))  // dfd = readtable(IOBuffer( x ) )
void main()
{   double  numbytes;
    numbytes =  xmean(2.0,5.0);
    printf("xmean is :  %d\n", numbytes);

    printf("finished and returned %s\n\n", myFunction("tst.csv") )   ;
    //const char* szSomeString = xfread(); // fraught with problems
    //printf("%s", szSomeString);
    printf("file contents :  \n%s\n\n", xfread("tst.csv"));
}
/* using DataFrames
push!(Base.DL_LOAD_PATH, pwd())
root=pwd();
function loadDF()
    x=unsafe_string(ccall((:xfread,"fread"), Cstring, (Cstring,), "orig.csv"));
    dfd = readtable(IOBuffer(x))
    #df_data = readtable(root*"/orig.csv",header=false);
    df_h = readtable(root*"/origHead.csv",header=false);
    names!(dfd, convert(Array{Symbol}, df_h[:x1]) )
end
@time df_in=loadDF();
*/
"""

function read_orig(root::String=pwd(),useGC::Bool=true)
    if !useGC gc_enable(false); end
    df_data = readtable(root*"/orig.csv",header=false);
    df_h = readtable(root*"/origHead.csv",header=false);
    names!(df_data, convert(Array{Symbol}, df_h[:x1]) )
    if !useGC gc_enable(true); end
    return df_data
end

dict2json(d::Dict) = JSON.json(d)
savejson(fname::String, j::String) = open(fname,"w") do f write(f,j) end
readjson(fname::String) = JSON.parsefile(fname)

save_cfg(root::String, cfg::Dict) = savejson("cfg.json", dict2json(cfg))
read_cfg(root::String) = readjson(root*"/cfg.json")

read_dfx(root::String) = readtable(root*"/dfx.csv")
save_dfx(root::String,dfx::DataFrame) = writetable(root*"/dfx.csv", dfx)
read_modelsDict(root::String) = readModels(root*"/modelsDict.json")
function save_modelsDict(root::String, d::Dict)
    j = JSON.json([d])
    open(root*"/modelsDict.json","w") do f
        write(f,j)
    end
end
read_dfd(root::String) = readtable(root*"/dfd.csv",header=true);
save_dfd(root::String,dfd::DataFrame) = writetable(root*"/dfd.csv", dfd)

read_dfr(root::String) = readtable(root*"/scored.csv",header=true);
function save_dfr(root::String,dfr::DataFrame)
    writetable(root*"/scored.tmp", dfr)
    outfile = open(root*"/scored.csv", "w")
    open(root*"/scored.tmp") do filehandle
        for line in eachline(filehandle)
            write(outfile, replace(replace(line,"NA",""),"NaN","") )
        end
    end
    close(outfile)
    rm(root*"/scored.tmp")
end

function write2disk(idf::DataFrame,fname::AbstractString, repvalues::Bool=true)
    ofile=fname
    if repvalues
        writetable(ofile*"tmp", idf, separator = ',', header = true)
        outfile = open(ofile, "w")
        open(ofile*"tmp") do filehandle
            for line in eachline(filehandle)
                write(outfile, replace(line,"NA",""))
            end
        end
        close(outfile)
        rm(ofile*"tmp")
    else
        writetable(ofile, idf, separator = ',', header = true)
    end
    println("Complete !!!! - output written to ",ofile)
end

# ------- R 2 Julia Conversion - END ---------------

#include("trd.jl")
#function Base.lowercase(df::DataFrame)
function lowercase!(df::DataFrame)
    for (i,v) in enumerate(names(df))   if v != Symbol(lowercase(string(v)))  rename!(df,v,Symbol(lowercase(string(v)))) end end
end
function Base.lowercase(arr::Array{Symbol})    return  Symbol[Symbol(lowercase(string(s))) for s in arr] end
function Base.lowercase(arr::Array{String})    return  String[lowercase(string(s)) for s in arr] end
function Base.lowercase(arr::Array{String})    return  String[lowercase(string(s)) for s in arr] end
function Base.lowercase(s::Symbol) return Symbol(lowercase(string(s))) end
function Base.lowercase(d::OrderedDict)
    vout=OrderedDict()
    for (k, v) in d  vout[lowercase(k)] = typeof(v) in [Bool,Float64,Int64] ? v : lowercase(v) end
    return vout
end

function grep(v::String,vars::Array{Symbol})   #function grep(v::AbstractString) filter(x->ismatch(Regex("$(v)"), string(x)), vars) end
    filter(x->ismatch(Regex("$(v)"), string(x)), vars)
end
function grep1(v::String,vars::Array{Symbol})   #function grep(v::AbstractString) filter(x->ismatch(Regex("$(v)"), string(x)), vars) end
    lst = filter(x->ismatch(Regex("$(v)"), string(x)), vars)
    return length(lst)>0 ? lst[1] : nothing
end

function grep(varr::Array{String},vars::Array{Symbol})
    arr_out=Symbol[]
    for v in varr
        res = filter(x->ismatch(Regex("$(v)"), string(x)), vars)
        arr_out = vcat(arr_out, res)
    end
    return arr_out
end
#function grep(varr::Array{String},vars::Array{Symbol}) return grep(convert(Array{String},varr),vars) end  #grep(["M","pan"],names(df_in))

function df2dict(idf::DataFrame)
    dout=OrderedDict()
    for n in names(idf)
       dout[n] = idf[1,Symbol(n)]
    end
    return dout
end

function df2dict(row::DataFrameRow)
    dout=OrderedDict()
    for n in names(row)
       dout[n] = row[Symbol(n)]
    end
    return dout
end

"""
app.cfg
P2_Competitor = true
pvalue_lvl = 0.20
excludedKeys = key1,key2,key3
exposed_flag_var = new_exposed_flag
sigLevel = 0.2
random_demos =
random_campaigns =
dropvars = goop,snop
scoring_vars =
TotalModelsOnly=false
"""

function loadCFG(d::OrderedDict,fn::String)
    function cvt(v::String, t::DataType)
        if t in [Bool]
            return lowercase(v) =="true" ? true : lowercase(v) =="false" ? false : NA
        elseif t in [nothing]
            return NA
        elseif t in [Array{String,1}]
            #return AbstractString[string(x) split(v,",")]
            return String[string(x) for x in split(replace(v,"\"",""),",")]
        elseif t in [Array{Symbol,1}]
            return Symbol[Symbol(x) for x in split(replace(v,"\"",""),",")]
        elseif t in [Symbol]
            return Symbol(v)
        elseif t in [String]
            return string(v)
        elseif t in [Float64]
            return parse(Float64,v)
        else
            return "__ooops__ TYPE not found!"
        end
    end
    if isfile(fn)
        f = open(fn)
        for ln in eachline(f)
            ln=replace(ln,"\n","")
            sidx=search(ln,"=").start-1
            eidx=search(ln,"=").stop+1
            if sidx>0
                k=strip(ln[1:sidx])
                v=strip(ln[eidx:end])
                if (length(k)>0) & (length(v)>0)
                    println("LINE:: ~"*k*"~"*v*"~")
                    #println("DEF::       ______ ",d[Symbol(k)],"  ----> ",typeof(d[Symbol(k)]) )
                    newk=Symbol(k)
                    dv=get(d,newk,nothing)
                    if dv!=nothing
                        newv = cvt(v,typeof(dv))
                        #println("IN::        ______ ",Symbol(k),"  ----> ",ct, " ....... ", typeof(ct))
                        #println( d[Symbol(k)],"|",v," {",typeof(d[Symbol(k)]),"}  ----->  ",newk,"|",newv,"  {",typeof(newv),"}")
                        d[newk] = newv
                    else
                        println(newk," in app.cfg, but not defaults")
                        #d[newk] = v
                    end
                else
                    #println("CFG ERROR :: invalid length : "*"~"*k*"~"*v*"~",typeof(d[Symbol(k)]))
                end
            else
                if length(ln)>0   println("CFG ERROR :: invalid cfg parameter : ~",ln,"~") end
            end
        end
        close(f)
    else
        println("Using Default configuration")
    end
    return d
end

# Retype DataFrame Columns
# Ensure/Convert DataFrame Column Types
function retypeBools!(dfd::DataFrame,vars::Array{Symbol}=Symbol[])
    if length(vars)==0 vars=names(dfd) end
    for c in vars
        if typeof(dfd[c]) in [DataArray{Any,1}, DataArray{String,1}]  #DataArray{Bool,1}
            println(c," ~ ",typeof(dfd[c]))
            if  typeof(dfd[c]) == DataArray{Any,1}
                if unique(map(x-> typeof(x),dfd[c])) == [Bool]
                    println("Converting Col")
                    dfd[c] = convert(Array{Bool}, dfd[c])
                else
                    println("Type = Any; but not Bool")
                end
            else
                if true in unique(isna(dfd[c]))
                    println("String Col has NA's")
                else
                    ul=convert(Array{String},unique(map(x -> lowercase(x), dfd[c] )))
                    if (length(setdiff(ul,["true","false"]))==0) & (length(intersect(ul,["true","false"]))>0)
                        println("Mapping String Col")
                        dfd[c]=  map(x-> lowercase(x)=="true" ? true :  false , dfd[c])
                        dfd[c] = convert(Array{Bool}, dfd[c])
                    else
                        println("String col not Bool")
                    end
                end
            end
        end
    end
end

max_stars_repo_name: werbac/jlib

max_issues_repo_name: werbac/jlib

max_forks_repo_name: werbac/jlib


text:
using DataStructures
using Base.Test

l1 = nil()
@test length(l1) == 0

l2 = cons(1, l1)
@test length(l2) == 1
@test head(l2) == 1

l3 = list(2, 3)
@test is(typeof(l3), Cons{Int})
@test length(l3) == 2
@test head(l3) == 2
@test head(tail(l3)) == 3
@test collect(l3) == [2; 3]

l4 = cat(l1, l2, l3)
@test length(l4) == 3
@test collect(l4) == [1; 2; 3]

l5 = map((x) -> x*2, l4)
@test typeof(l5) == Cons{Int}
@test collect(l5) == [2; 4; 6]

l5b = map((x) -> "$x", l5)
@test typeof(l5b) == Cons{ASCIIString}
@test collect(l5b) == ["2"; "4"; "6"]

l6 = filter((x) -> x < 6, l5)
@test length(l6) == 2
@test l6.head == 2
@test l6.tail.head == 4

l7 = reverse(l6)
@test length(l7) == 2
@test l7.head == 4
@test l7.tail.head == 2

l8 = cat(l5b, l6)
@test collect(l8) == ["2"; "4"; "6"; 2; 4]

l9 = cat(list(1, 2), list(3.0, 4.0))
@test typeof(l9) == Cons{Real}
@test collect(l9) == [1; 2; 3.0; 4.0]


text:
@test_logs (:warn, r"only informs an external profiler to start") CUDA.Profile.start()
CUDA.Profile.stop()

@test_logs (:warn, r"only informs an external profiler to start") CUDA.@profile begin end


text:
using Parameters, JuMP, JuMP.Containers

# --- Utility ---

data(a::DenseAxisArray) = a.data
data(a::Any) = a

round_int(x::AbstractFloat) = Integer(round(x))
convert_int(x, ::Type{<:Integer}) = round_int(x)
convert_int(x, ::Type{Array{T, N}}) where T <: Integer where N = round_int.(x)
convert_int(x, ::Any) = x

"""The function queries values from the model to data type based on its field names. It extracts values from DenseAxisArray from its `data` field. Then, it converts the values to the corresponding field type. The function rounds integers before conversion because JuMP outputs integer variables as floats."""
function model_to_dtype(dtype::DataType, model::Model)
    fields = (
        value.(model[n]) |> data |> x -> convert_int(x, t)
        for (n, t) in zip(fieldnames(dtype), fieldtypes(dtype)))
    dtype(fields...)
end

# --- Model ---

"""ShelfSpaceAllocationModel type as JuMP.Model"""
const ShelfSpaceAllocationModel = Model

"""Specs"""
@with_kw struct Specs
    height_placement::Bool = true
    blocking::Bool = true
end

# TODO: I <: Integer, F <: AbstractFloat

"""Parameters"""
@with_kw struct Params
    # --- Sets and Subsets ---
    products::Array{Int, 1}
    shelves::Array{Int, 1}
    blocks::Array{Int, 1}
    modules::Array{Int, 1}
    P_b::Array{Array{Int, 1}, 1}
    S_m::Array{Array{Int, 1}, 1}
    # --- Parameters ---
    # Products
    N_p_min::Array{Float64, 1}
    N_p_max::Array{Float64, 1}
    G_p::Array{Float64, 1}
    R_p::Array{Float64, 1}
    D_p::Array{Float64, 1}
    L_p::Array{Float64, 1}
    W_p::Array{Float64, 1}
    H_p::Array{Float64, 1}
    M_p::Array{Float64, 1}
    SK_p::Array{Float64, 1}
    # Shelves
    M_s_min::Array{Float64, 1}
    M_s_max::Array{Float64, 1}
    W_s::Array{Float64, 1}
    H_s::Array{Float64, 1}
    L_s::Array{Int, 1}
    # Product-shelves
    P_ps::Array{Float64, 2}
    # Constants
    SL::Float64 = 0.0
    w1::Float64 = 0.5
    w2::Float64 = 10.0
    w3::Float64 = 0.1
end

"""Variables"""
@with_kw struct Variables
    # --- Basic Variables ---
    s_p::Array{Float64, 1}
    e_p::Array{Float64, 1}
    o_s::Array{Float64, 1}
    n_ps::Array{Int, 2}
    y_p::Array{Int, 1}
    # --- Blocking Variables ---
    b_bs::Array{Float64, 2}
    m_bm::Array{Float64, 2}
    z_bs::Array{Int, 2}
    z_bs_f::Array{Int, 2}
    z_bs_l::Array{Int, 2}
    x_bs::Array{Float64, 2}
    x_bm::Array{Float64, 2}
    w_bb::Array{Int, 2}
    v_bm::Array{Int, 2}
end

"""Objectives"""
@with_kw struct Objectives
    empty_shelf_space::Float64
    profit_loss::Float64
    height_placement_penalty::Float64
end

"""Variable values from model.

# Arguments
- `model::ShelfSpaceAllocationModel`
"""
Variables(model::ShelfSpaceAllocationModel) = model_to_dtype(Variables, model)

"""Objective values from model.

# Arguments
- `model::ShelfSpaceAllocationModel`
"""
Objectives(model::ShelfSpaceAllocationModel) = model_to_dtype(Objectives, model)

"""Mixed Integer Linear Program (MILP) formulation of the Shelf Space Allocation
Problem (SSAP).

# Arguments
- `parameters::Params`
- `specs::Specs`
"""
function ShelfSpaceAllocationModel(parameters::Params, specs::Specs)
    # Unpack parameters values
    @unpack products, shelves, blocks, modules, P_b, S_m, N_p_min, N_p_max,
            G_p, R_p, D_p, L_p, W_p, H_p, M_p, SK_p, M_s_min, M_s_max, W_s,
            H_s, L_s, P_ps, SL, w1, w2, w3 = parameters

    # Initialize the model
    model = ShelfSpaceAllocationModel()

    # --- Basic Variables ---
    @variable(model, s_p[products] ≥ 0)
    @variable(model, e_p[products] ≥ 0)
    @variable(model, o_s[shelves] ≥ 0)
    @variable(model, n_ps[products, shelves] ≥ 0, Int)
    @variable(model, y_p[products], Bin)

    # --- Block Variables ---
    @variable(model, b_bs[blocks, shelves] ≥ 0)
    @variable(model, m_bm[blocks, modules] ≥ 0)
    @variable(model, z_bs[blocks, shelves], Bin)
    @variable(model, z_bs_f[blocks, shelves], Bin)
    @variable(model, z_bs_l[blocks, shelves], Bin)
    @variable(model, x_bs[blocks, shelves] ≥ 0)
    @variable(model, x_bm[blocks, modules] ≥ 0)
    @variable(model, w_bb[blocks, blocks], Bin)
    @variable(model, v_bm[blocks, modules], Bin)

    # --- Height and weight constraints ---
    for p in products, s in shelves
        if (H_p[p] > H_s[s]) | (M_p[p] > M_s_max[s])
            fix(n_ps[p, s], 0, force=true)
        end
    end

    # --- Objective ---
    @expression(model, empty_shelf_space,
        sum(o_s[s] for s in shelves))
    @expression(model, profit_loss,
        sum(G_p[p] * e_p[p] for p in products))
    if specs.height_placement
        @expression(model, height_placement_penalty,
            sum(L_p[p] * L_s[s] * n_ps[p, s] for p in products for s in shelves))
    else
        @expression(model, height_placement_penalty, AffExpr(0.0))
    end
    @objective(model, Min,
        w1 * empty_shelf_space +
        w2 * profit_loss +
        w3 * height_placement_penalty
    )

    # --- Basic constraints ---
    @constraints(model, begin
        [p = products],
        s_p[p] ≤ sum(30 / R_p[p] * P_ps[p, s] * n_ps[p, s] for s in shelves)
        [p = products],
        s_p[p] ≤ D_p[p]
    end)
    @constraint(model, [p = products],
        s_p[p] + e_p[p] == D_p[p])
    @constraint(model, [p = products],
        sum(n_ps[p, s] for s in shelves) ≥ y_p[p])
    @constraints(model, begin
        [p = products],
        N_p_min[p] * y_p[p] ≤ sum(n_ps[p, s] for s in shelves)
        [p = products],
        sum(n_ps[p, s] for s in shelves) ≤ N_p_max[p] * y_p[p]
    end)
    @constraint(model, [s = shelves],
        sum(W_p[p] * n_ps[p, s] for p in products) + o_s[s] == W_s[s])

    # --- Block constraints ---
    if specs.blocking
        @constraint(model, [s = shelves, b = blocks],
            sum(W_p[p] * n_ps[p, s] for p in P_b[b]) ≤ b_bs[b, s])
        @constraint(model, [s = shelves],
            sum(b_bs[b, s] for b in blocks) ≤ W_s[s])
        @constraint(model, [b = blocks, m = modules, s = S_m[m]],
            b_bs[b, s] ≥ m_bm[b, m] - W_s[s] * (1 - z_bs[b, s]) - SL)
        @constraint(model, [b = blocks, m = modules, s = S_m[m]],
            b_bs[b, s] ≤ m_bm[b, m] + W_s[s] * (1 - z_bs[b, s]) + SL)
        @constraint(model, [b = blocks, s = shelves],
            b_bs[b, s] ≤ W_s[s] * z_bs[b, s])
        # ---
        @constraint(model, [b = blocks, s = 1:length(shelves)-1],
            z_bs_f[b, s+1] + z_bs[b, s] == z_bs[b, s+1] + z_bs_l[b, s])
        @constraint(model, [b = blocks],
            sum(z_bs_f[b, s] for s in shelves) ≤ 1)
        @constraint(model, [b = blocks],
            sum(z_bs_l[b, s] for s in shelves) ≤ 1)
        @constraint(model, [b = blocks],
            z_bs_f[b, 1] == z_bs[b, 1])
        @constraint(model, [b = blocks],
            z_bs_l[b, end] == z_bs[b, end])
        # ---
        @constraint(model, [b = blocks, s = shelves],
            sum(n_ps[p, s] for p in P_b[b]) ≥ z_bs[b, s])
        @constraint(model, [b = blocks, s = shelves, p = P_b[b]],
            n_ps[p, s] ≤ N_p_max[p] * z_bs[b, s])
        # ---
        @constraint(model, [b = blocks, b′ = filter(a->a≠b, blocks), s = shelves],
            x_bs[b, s] + W_s[s] * (1 - z_bs[b, s]) ≥ x_bs[b′, s] + b_bs[b, s] - W_s[s] * (1 - w_bb[b, b′]))
        @constraint(model, [b = blocks, b′ = filter(a->a≠b, blocks), s = shelves],
            x_bs[b′, s] + W_s[s] * (1 - z_bs[b′, s]) ≥ x_bs[b, s] + b_bs[b, s] - W_s[s] * w_bb[b, b′])
        @constraint(model, [b = blocks, m = modules, s = S_m[m]],
            x_bm[b, m] ≥ x_bs[b, s] - W_s[s] * (1 - z_bs[b, s]) - SL)
        @constraint(model, [b = blocks, m = modules, s = S_m[m]],
            x_bm[b, m] ≤ x_bs[b, s] + W_s[s] * (1 - z_bs[b, s]) + SL)
        @constraint(model, [b = blocks, s = shelves],
            x_bs[b, s] ≤ W_s[s] * z_bs[b, s])
        @constraint(model, [b = blocks, s = shelves],
            x_bs[b, s] + b_bs[b, s] ≤ W_s[s])
        # ---
        @constraint(model, [b = blocks, m = modules, s = S_m[m], p = P_b[b]],
            n_ps[p, s] ≤ N_p_max[p] * v_bm[b, m])
        @constraint(model, [b = blocks],
            sum(v_bm[b, m] for m in modules) ≤ 1)
    end

    return model
end


text:
using ApproxFun, SingularIntegralEquations, Test
import Base.MathConstants: γ, catalan

@testset "Convolution ProductFun" begin
    @testset "Chebyshev addition" begin
        d = Segment(-1.,1.)
        fK(x,y) = exp(abs2(y-x))
        FLR = LowRankFun(fK,Chebyshev(d),Chebyshev(d))
        FPF = ProductFun(fK,Chebyshev(d),Chebyshev(d))
        G = convolutionProductFun(fK,Chebyshev(d),Chebyshev(d))
        G1 = convolutionProductFun(fK,Ultraspherical(1,d),Chebyshev(d))
        G2 = convolutionProductFun(fK,Ultraspherical(1,d),ChebyshevDirichlet{1,1}(d))

        @test norm(pad(coefficients(FLR),size(coefficients(G),1),size(coefficients(G),2))-coefficients(G))<100000eps()
        @test norm(pad(coefficients(FPF),size(coefficients(G),1),size(coefficients(G),2))-coefficients(G))<100000eps()
        @test norm(fK(.123,.456)-G(.123,.456))≤10000eps()
        @test norm(G(.123,.456)-G1(.123,.456))≤10000eps()
        @test norm(G(.123,.456)-G2(.123,.456))≤10000eps()
    end

    @testset "Fourier on PeriodicSegment" begin
        f2 = Fun(θ->exp(sin(θ))+sin(cos(θ)),Fourier())
        FLR = LowRankFun((θ,ϕ)->f2(ϕ-θ),Fourier(),Fourier())
        FPF = ProductFun((θ,ϕ)->f2(ϕ-θ),Fourier(),Fourier())
        G = convolutionProductFun(f2,Fourier(),Fourier())

        @test norm(pad(coefficients(FLR),size(coefficients(G),1),size(coefficients(G),2))-coefficients(G))<400eps()
        @test norm(pad(coefficients(FPF),size(coefficients(G),1),size(coefficients(G),2))-coefficients(G))<400eps()
        @test norm(f2(.456-.123)-G(.123,.456))≤400eps()

        f2 = Fun(θ->1+sin(cos(θ)),CosSpace())
        G = convolutionProductFun(f2,Fourier(),Fourier())
        @test norm(f2(.456-.123)-G(.123,.456))≤400eps()

        f2 = Fun(θ->sin(sin(θ)),SinSpace())
        G = convolutionProductFun(f2,Fourier(),Fourier())
        @test norm(f2(.456-.123)-G(.123,.456))≤400eps()
    end

    @testset "Laurent on PeriodicSegment" begin
        f2 = Fun(θ->exp(sin(θ))+sin(cos(θ)),Laurent())
        G = convolutionProductFun(f2,Laurent(),Laurent())
        @test norm(f2(.456-.123)-G(.123,.456))≤400eps()

        f2 = Fun(θ->π+ℯ*exp(im*θ)+sqrt(2)*exp(im*2θ)+catalan*exp(im*3θ)+γ*exp(im*4θ),Taylor(PeriodicSegment()))
        G = convolutionProductFun(f2,Laurent(),Laurent())
        @test norm(f2(.456-.123)-G(.123,.456))≤100eps()

        f2 = Fun(θ->ℯ*exp(-im*θ)+sqrt(2)*exp(-im*2θ)+catalan*exp(-im*3θ)+γ*exp(-im*4θ),Hardy{false}(PeriodicSegment()))
        G = convolutionProductFun(f2,Laurent(),Laurent())
        @test norm(f2(.456-.123)-G(.123,.456))≤100eps()
    end
    @testset "Timing" begin
        d = Segment(-2.5,-.5)
        fK2(x,y) = besselj0(100(y-x))
        convolutionProductFun(fK2,Chebyshev(d),Chebyshev(d))
        @time G = convolutionProductFun(fK2,Chebyshev(d),Chebyshev(d))
        convolutionProductFun(fK2,Ultraspherical(1,d),ChebyshevDirichlet{1,1}(d))
        @time G = convolutionProductFun(fK2,Ultraspherical(1,d),ChebyshevDirichlet{1,1}(d))
        println("Chebyshev addition: Time should be ~0.01 seconds.")

        f2 = Fun(θ->besselj0(500*abs(2sin(θ/2))),CosSpace())
        convolutionProductFun(f2,Fourier(),Fourier())
        @time G = convolutionProductFun(f2,Fourier(),Fourier())
        println("CosSpace addition: Time should be ~0.05 seconds.")

        f2 = Fun(θ->besselj0(500*abs(2sin(θ/2))),Laurent())
        convolutionProductFun(f2,Laurent(),Laurent())
        @time G = convolutionProductFun(f2,Laurent(),Laurent())
        println("Laurent addition: Time should be ~0.14 seconds.")
    end
end


text:
"""
    function rr_blup(X, gt, y, par)
## Description
Estimate breeding values and locus effects with method RR-BLUP.
Where,
- `X`: is `n × p` incidence matrix for the non-genetic fixed effects.
- `gt`: are 0, 1, 2 genotypes.
- `y`: is a vector of phenotypes.
Argument `par` is a named tuple of parameters for the methods.
- `nit`: number of iterations, or chain length.
- `dflv`: hyperparameter (degree of freedom) for locus effect variance.
- `dfrv`: hyperparameter (degree of freedom) for residual variance.
- `vg`: used to derive hyper parameter (scale) for locus effect variance.
- `vr`: used to derive hyper parameter (scale) for residual variance.
- `window`: number of consecutive markers in a genomic window.

## Notes
This function is based on chapter 10, Bayesian methods applied to GWAS
by Fernando and Garrick of Gondro, van der Werf, and Heyes, Genome-wide
assocation studies and genomic prediction, 2013, Humana Press.
"""
function rr_blup(X, gt, y, par)
    num = (loci    = size(gt)[1],
           ID      = size(gt)[2],
           windows = Int(ceil(size(gt)[1] / par.window)),
           samples = Int(floor(par.n/par.ofrq))
           )
    p = mean(gt, dims = 2) ./ 2 # realized allele frequency
    Z = zeros(Float64, size(gt))
    copyto!(Z, gt)
    Z .-= 2p

    m2pq = mean(2 .* p .* (1 .- p))
    vq = par.var.genotype / num.loci / m2pq  # QTL variance
    # Scale parameters for the scaled inverse Χ² priors
    S² = (genotype = vq * (par.df.genotype - 2) / par.df.genotype,
          residue  = par.var.residue * (par.df.residue - 2) / par.df.residue)

    α, mα, μ, mμ = zeros(num.loci), zeros(num.loci), mean(y), 0
    mdfrq = zeros(num.loci)
    vg = zeros(num.samples)
    wvp = zeros(num.samples, num.windows) # window var prop.
    sc = 0                      # sample counts
    yc = y .- μ
    dχ² = Chisq(num.ID + par.df.residue)

    # the MCMC part
    for i in 1:par.n
        # sample residue variance
        ve = (yc'yc + par.df.residue * S².residue) / rand(dχ²)
        # sample intercept
        yc .+= μ
        rhs = sum(yc)

    end
end

function _test_ba(dat)
    T, y, b = begin
        if isfile(dat)
            @debug "Load formerly generated data"
            deserialize(dat)
        else
            @debug "Create data and share them with R"
            nlc, nid, μ, σ, ϵ = 5, 10, 0., .5, 1.0
            T, y, b = qsgp(nlc, nid, μ, σ, ϵ)
            dir = dirname(dat)
            open("$dir/T.csv", "w") do io
                for i in 1:nid
                    println(io, join(T[:, i], ", "))
                end
            end
            open("$dir/y.csv", "w") do io
                for x in y
                    println(io, x)
                end
            end
            open("$dir/b.csv", "w") do io
                for x in b
                    println(io, b)
                end
            end
            serialize(dat, (T = T, y = y, b = b))
            T, y, b
        end
    end

    par = (; Dict(:n => 10_000,
                  :df => (genotype = 4., residue = 4.),
                  :var => (genotype = 1., residue = 1.),
                  :window => 40, # windows size, n markers
                  :ofrq => 10     # output frequency
                  )...)
    nlc, nid = size(T)
    rr_blup(ones(nid), T, y, par)
end

function bayes_a()
end

function bayes_b(; nit = 10_000)
end

function bayes_c(; nit = 10_000)
end

function bayes_cπ(; nit = 10_000)
end