text:
include( "uses.jl")

const FORM_EXTRA =Dict(
	"border-bottom"=>"1px dashed #aaaaaa",
	"margin-bottom"=>"5px",
	"margin-top"=>"5px")

const PREAMBLE = """
Regardless of politics, I imagine most people an agree on the outlines of what
our taxation and benefit system should do. Give some support to people who haven't got much,
(depending who you ask, maybe take that support away from better off people); as you
get better off, pay a bit (or a lot) more in tax. But that's not always what the systems
we have do. Here's a simulation showing how the rules of the tax and benefit system
interact to make it almost impossible for households to make sensible choices about
whether it's worth taking a job, working a bit more, or going for a promotion.  Using the
controls to the left, you can see from the charts how the choices between what you earn and what
you end up with vary for different family times (renting or owning, with and without children, high or low earners).
"""

const INFO = """
The graph shows how net the net income of a household (i.e. after deducting taxes and optionally housing costs, then adding benefits)
varies with gross earnings, under the new [Universal Credit](https://www.gov.uk/universal-credit)
benefit system and the
['legacy' Tax Credit system](https://www.gov.uk/topic/benefits-credits/tax-credits) which UC is gradually replacing. You can also switch to an 'economist's
view' of the labour/leisure choice. Move your mouse over the chart to see a detailed
breakdown of what's happening at each point.

For simplicity, we assume:

* only one person in the household works;
* the adults are of working age;
* there are no other sources of income, for example savings income or disability benefits; and
* the working person faces a fixed hourly wage and can work any hours at that wage.

You'll see that the graphs are complicated enough even with those simplifications.
This is especially so for families with large housing costs
or several children, for whom the [Benefit Cap](https://www.gov.uk/benefit-cap) bites.

See [here](https://stb.virtual-worlds.scot/bc-intro.html) for more on the ideas behind budget constraints.

* Created with [Julia](https://julialang.org/) | [Dash](https://dash-julia.plotly.com/) | [Plotly](https://plotly.com/julia/) | [Budget Constraint Generator](https://github.com/grahamstark/BudgetConstraints.jl);
* Part of the [Scottish Tax Benefit Model](https://github.com/grahamstark/ScottishTaxBenefitModel.jl);
* Open Source software released under the [MIT Licence](https://github.com/grahamstark/Visualisations.jl/blob/main/LICENSE). [Source Code](https://github.com/grahamstark/Visualisations.jl).
"""

"""
Generate a pair of budget constraints (as Dataframes) for the given household.
"""
function getbc(
	hh  :: Household,
	sys :: TaxBenefitSystem,
	wage :: Real,
	settings :: Settings )::Tuple
	defroute = settings.means_tested_routing

	settings.means_tested_routing = lmt_full
	lbc = BCCalcs.makebc( hh, sys, settings, wage )

	settings.means_tested_routing = uc_full
	ubc = BCCalcs.makebc( hh, sys, settings, wage )

	settings.means_tested_routing = defroute
    (lbc,ubc)
end

"""
Retrieve one of the model's example households & overwrite a few fields
to make things simpler.
"""
function get_hh(
	tenure    :: AbstractString,
	bedrooms  :: Integer,
	hcost     :: Real,
	marrstat    :: AbstractString,
	chu5      :: Integer,
	ch5p      :: Integer ) :: Household
	hh = get_example( single_hh )
	head = get_head(hh)
	head.age = 30
	sp = get_spouse(hh)
	enable!(head) # clear dla stuff from example
	hh.tenure = if tenure == "private"
		Private_Rented_Unfurnished
	elseif tenure == "council"
		Council_Rented
	elseif tenure == "owner"
		Mortgaged_Or_Shared
	else
		@assert false "$tenure not recognised"
	end
	hh.bedrooms = bedrooms
	hh.other_housing_charges = hh.water_and_sewerage = 0
	if hh.tenure == Mortgaged_Or_Shared
		hh.mortgage_payment = hcost
		hh.mortgage_interest = hcost
		hh.gross_rent = 0
	else
		hh.mortgage_payment = 0
		hh.mortgage_interest = 0
		hh.gross_rent = hcost
	end
	if marrstat == "couple"
		sex = head.sex == Male ? Female : Male # hetero ..
		add_spouse!( hh, 30, sex )
		sp = get_spouse(hh)
		enable!(sp)
		set_wage!( sp, 0, 10 )
	end
	age = 0
	for ch in 1:chu5
		sex = ch % 1 == 0 ? Male : Female
		age += 1
		add_child!( hh, age, sex )
	end
	age = 7
	for ch in 1:ch5p
		sex = ch % 1 == 0 ? Male : Female
		age += 1
		add_child!( hh, age, sex )
	end
	set_wage!( head, 0, 10 )
	for (pid,pers) in hh.people
		# println( "age=$(pers.age) empstat=$(pers.employment_status) " )
		empty!( pers.income )
		empty!( pers.assets )
	end
	return hh
end

const MAX_HRS = 80

function gross_to_leisure!( bc :: DataFrame, wage::Real)
	gross = bc[:,:gross]
	l = MAX_HRS .- (gross ./ wage)
	w = MAX_HRS .- l
	# println(l)
	bc.leisure = l
	bc.ls = w
	# some nice way
	p = 0
	for r in eachrow(bc)
		p += 1
		if r.leisure < 0
			break
		end
	end
	bc = bc[1:p,:]
end

"""
Plot two budget constraints (contained in dataframes) - legacy & universal credit.
"""
function econ_bcplot( lbc:: DataFrame, ubc :: DataFrame, wage :: Real, ytitle :: String )
	# legacy
	gross_to_leisure!(lbc, wage )
	bl = scatter(
           lbc,
		   x=:leisure,
		   y=:net,
           mode="line",
		   name="Legacy Benefits",
		   text=:simplelabel,
		   hoverinfo="text"
       )
	# uc
	gross_to_leisure!(ubc, wage )
	bu = scatter(
		ubc,
		x=:leisure,
		y=:net,
		mode="line",
		name="Universal Credit",
		text=:simplelabel,
		hoverinfo="text"
		# ,xaxis="x2"
	)
	layout = Layout(
		title="Budget Constraint: Legacy Benefits vs Universal Credit",
        xaxis_title="Leisure (hours p.w.)",
        yaxis_title=ytitle,
		xaxis_range=[0, MAX_HRS ],
		yaxis_range=[0, 1_500],
		#=
		xaxis2=attr(

            title="Hours of work",
			# titlefont_color="blue",
            overlaying="x",
			side="bottom",
			position=0.15,
			anchor="free"

        ),
		=#

		legend=attr(x=0.01, y=0.95),
		width=700,
		height=700)
	p = PlotlyJS.Plot( [bl, bu], layout)
	# (typeof(p))
	return p
end

"""
Plot two budget constraints (contained in dataframes) - legacy & universal credit.
"""
function bcplot( lbc:: DataFrame, ubc :: DataFrame, ytitle :: String )
	# legacy
	bl = scatter(
           lbc,
		   x=:gross,
		   y=:net,
           mode="line",
		   name="Legacy Benefits",
		   text=:simplelabel,
		   hoverinfo="text"
       )
	# uc
	bu = scatter(
		ubc,
		x=:gross,
		y=:net,
		mode="line",
		name="Universal Credit",
		text=:simplelabel,
		hoverinfo="text"
	)
	# 45% line
	gn = scatter(y=[0,1500], x=[0,1500], showlegend=false, name="")

	gn["marker"] = Dict(:color => "#ccc",
                        :line => Dict(:color=> "#ccc",
                        :width=> 0.5))
	layout = Layout(
		title="Budget Constraint: Legacy Benefits vs Universal Credit",
        xaxis_title="Household Gross Earnings £p.w.",
        yaxis_title=ytitle,
		xaxis_range=[0, 1_200],
		yaxis_range=[0, 1_200],
		legend=attr(x=0.01, y=0.95),
		width=700,
		height=700)
	p = PlotlyJS.Plot( [gn, bl, bu], layout)
	# (typeof(p))
	return p
end

"""
Create the whole plot for the named household.
"""
function doplot(
	wage :: Real,
	tenure :: AbstractString,
	bedrooms::Integer,
	hcost::Real,
	marrstat::AbstractString,
	chu5::Integer,
	ch5p::Integer,
	view :: AbstractString,
	target_str :: AbstractString )
	hh = get_hh( tenure, bedrooms, hcost, marrstat, chu5, ch5p )
	# println(to_md_table(hh))
	settings = Settings()
	ytitle = ""
	target = nothing
	if target_str == "ahc_hh"
		target = ahc_hh
		ytitle = "Household Net Income After Housing Costs £p.w."
	elseif target_str == "bhc_hh"
		target = bhc_hh
		ytitle = "Household Net Income Before Housing Costs £p.w."
	elseif target_str == "total_bens"
		ytitle = "Total State Benefits received £pw"
		target = total_bens
	elseif target_str == "total_taxes"
		ytitle = "Total Income Tax and NI Payments £pw"
		target = total_taxes
	end # forgot how to cover
	settings.target_bc_income = target
	lbc, ubc = getbc( hh, sys, wage, settings )
	if view == "l_vs_l"
		figure=econ_bcplot( lbc, ubc, wage, ytitle )
	else
		figure = bcplot( lbc, ubc, ytitle )
	end
	return figure
end

"""
Create the block of sliders and radios on the LHS
"""
function get_input_block()
	return dbc_form(

		[
		dbc_row([
			dbc_col(
				dbc_label("Wage (£ per hour)"; html_for="wage"), width=3
			),
			dbc_col(
				dcc_slider(
					id = "wage",
					min = 1,
					max = 40,
					marks = Dict([Symbol("$v") => Symbol("$v") for v in 0:10:100]),
					value = 10.0,
					step = 1
				)) #
		], style=FORM_EXTRA),
		dbc_row([
			dbc_col(
				dbc_label("Tenure:"; html_for="tenure"), width=3
			),
			dbc_col(
				dbc_radioitems(
					id = "tenure",
					options = [(value = "owner", label= "Owner Occupier"),
					           (value = "private", label="Private Rented"),
							   (value = "council", label="Council/Housing Association")],
					value = "private"
				)
			)
		], style=FORM_EXTRA),
		dbc_row([
			dbc_col(
				dbc_label("Bedrooms:"; html_for="bedrooms"), width=3
			),
			dbc_col(
				dcc_slider(
					id = "bedrooms",
					min = 1,
					max = 5,
					marks = Dict([Symbol("$v") => Symbol("$v") for v in 1:5]),
					value = 2,
					step = 1
    			)
			)
		], style=FORM_EXTRA), # row
		dbc_row([
			dbc_col(
				dbc_label("Housing Costs £pw:"; html_for="hcost"), width=3
			),
			dbc_col(
				dcc_slider(
					id = "hcost",
					min = 0,
					max = 300,
					marks = Dict([Symbol("$v") => Symbol("$v") for v in 0:100:300]),
					value = 100,
					step = 1
    			)
			)
		], style=FORM_EXTRA), # row
		dbc_row([
			dbc_col(
				dbc_label("Adults:"; html_for="marrstat"), width=3
			),
			dbc_col(
				dbc_radioitems(
					id = "marrstat",
					options = [(value = "single", label= "Single"),
							   (value="couple", label="Couple")],
					value = "single"
				)
			)
		], style=FORM_EXTRA), # row
		dbc_row([
			dbc_col(
				dbc_label("Children aged under 5:"; html_for="chu5"), width=3
			),
			dbc_col(
				dcc_slider(
					id = "chu5",
					min = 0,
					max = 5,
					marks = Dict([Symbol("$v") => Symbol("$v") for v in 0:5]),
					value = 0,
					step = 1
    			)
			)
		], style=FORM_EXTRA),
		dbc_row([
			dbc_col(
				dbc_label("Children aged 5+:"; html_for="ch5p"), width=3
			),
			dbc_col(
				dcc_slider(
					id = "ch5p",
					min = 0,
					max = 8,
					marks = Dict([Symbol("$v") => Symbol("$v") for v in 0:8]),
					value = 0,
					step = 1
    			)
			)
		], style=FORM_EXTRA), # row
		dbc_row([
			dbc_col(
				dbc_label("Graph type:"; html_for="view"), width=3
			),
			dbc_col(
				dbc_radioitems(
					id = "view",
					options = [(value = "g_vs_n", label= "Gross vs Net Income"),
					           (value = "l_vs_l", label="Labour/Leisure")],
					value = "g_vs_n")
			)
		], style=FORM_EXTRA), # row
		dbc_row([
			dbc_col(
				dbc_label("Income Measure:"; html_for="target"), width=3
			),
			dbc_col(
				dbc_radioitems(
					id = "target",
					options = [(value = "ahc_hh", label= "After Housing Costs"),
					           (value = "bhc_hh", label= "Before Housing Costs"),
							   (value = "total_bens", label="Total Benefits Received"),
							   (value = "total_taxes", label="Total Taxes Paid") ],
					value = "ahc_hh")
			)
		], style=FORM_EXTRA)
	])
end

"""
an experimental html table. Not actually used.
"""
function generate_table(df :: DataFrame)
	rows = []
	for r in eachrow(df)
		push!(rows,
			html_tr([
				html_td(r.gross),
				html_td(r.net),
				html_td(dcc_markdown(r.label)),
				html_td(dcc_markdown(r.label_p1))]
			)
		)
	end
    t = html_table([
        html_thead(html_tr([html_th("Gross"), html_th("Net"), html_th("Breakdown"), html_th("Breakdown+1p")])),
        html_tbody(rows)])
end

sys = load_file( joinpath( Definitions.MODEL_PARAMS_DIR, "sys_2021_22.jl" ))
load_file!( sys, joinpath( Definitions.MODEL_PARAMS_DIR, "sys_2021-uplift-removed.jl"))
weeklyise!( sys )

app = dash(external_stylesheets=[dbc_themes.UNITED],
	url_base_pathname="/bcd/")
# BOOTSTRAP|SIMPLEX|MINTY|COSMO|SANDSTONE|UNITED|SLATE|SOLAR|UNITED|
app.layout = dbc_container(fluid=true, className="p-5") do
	html_title( "Scotland's Kinkiest Families: Household Budget Constraints")
	html_h1("Scotland's Kinkiest Families: Household Budget Constraints"),
	dbc_row([
		dbc_col( dcc_markdown( PREAMBLE ), width=10)
	]),
	dbc_row([
    	dbc_col(get_input_block(), width=4),
	    dbc_col( dcc_graph( id = "bc-1" ))
		]
	),
	dbc_row([
		dbc_col( dcc_markdown( INFO ), width=10)
	])
end

callback!(
    app,
    Output("bc-1", "figure"),
	# Input( "famchoice", "value"),
	Input( "wage", "value"),
	Input( "tenure", "value"),
	Input( "bedrooms", "value"),
	Input( "hcost", "value" ),
	Input( "marrstat", "value"),
	Input( "chu5", "value"),
	Input( "ch5p", "value"),
	Input( "view", "value"),
	Input( "target", "value")) do wage, tenure, bedrooms, hcost, marrstat, chu5, ch5p, view, target
		return doplot( wage, tenure, bedrooms, hcost, marrstat, chu5, ch5p, view, target )
	end

run_server(app, "0.0.0.0", debug=true )


text:
module Vofinit

using libvofi_jll

export getcelltype, getcc

const cintegrand = Ref{Ptr{Cvoid}}(0)

function integrand(xyz, thunk)
    func = unsafe_pointer_to_objref(thunk)
    x = unsafe_load(xyz, 1)
    y = unsafe_load(xyz, 2)
    z = unsafe_load(xyz, 3)
    func(x, y, z)
end

function getcelltype(func, x0, h0, ndim0 = Cint(3))
    @ccall libvofi.vofi_get_cell_type(cintegrand[]::Ptr{Cvoid},
                                      func::Any,
                                      x0::Ref{Cdouble},
                                      h0::Ref{Cdouble},
                                      ndim0::Cint)::Cint
end

function getcc(func, x0, h0, xex, ndim0 = Cint(3);
               nex = Cint.((0, 0)),
               npt = Cint.((4, 4, 4, 4)),
               nvis = Cint.((0, 0)))
    @ccall libvofi.vofi_get_cc(cintegrand[]::Ptr{Cvoid},
                               func::Any,
                               x0::Ref{Cdouble},
                               h0::Ref{Cdouble},
                               xex::Ref{Cdouble},
                               nex::Ref{Cint},
                               npt::Ref{Cint},
                               nvis::Ref{Cint},
                               ndim0::Cint)::Cdouble
end

function __init__()
    cintegrand[] = @cfunction(integrand, Cdouble, (Ptr{Cdouble}, Ptr{Cvoid}))
end

end

ext: jl


text:
sparse_options = (; conic_solver = sparse_hypatia)
noext_options = (sparse_options..., use_extended_form = false)
nosubp_options = (sparse_options..., solve_subproblems = false) # for SOS2
noconic_options = (nosubp_options..., solve_relaxation = false)

insts = OrderedDict()

insts["test"] = [
    ((1, 1, 3, true),),
    ((1, 1, 3, false),),
    ((2, 3, 4, true),),
    ((2, 3, 4, false),),
    ((3, 2, 5, true),),
    ((3, 2, 5, false), noext_options),
    ((3, 2, 5, true), nosubp_options),
    ((3, 2, 5, false), noconic_options),
    ((1, 2, 3, true, true, SOS2(), 16), nosubp_options),
    ((2, 2, 4, true, true, SOS2(), 16), noconic_options),
    ((2, 2, 4, true, true, CCBounded(), 16), sparse_options),
    ((2, 2, 4, true, true, CCBounded(), 16), noconic_options),
    ((3, 3, 6, true, true, LogIBBounded(), 32), nosubp_options),
    ((3, 3, 6, true, true, LogIBBounded(), 32), noconic_options),
]

# convex instances
function modulardesign_insts(max_n::Int, options::NamedTuple)
    return [((6, 5, n, true), options) for n in vcat(20, 20:20:max_n)]
end

insts["nat"] = modulardesign_insts(220, sparse_options)
insts["noext"] = modulardesign_insts(220, noext_options)

# nonconvex instances
function modulardesign_insts(pwl::PWLSOS2, max_n::Int)
    return [((3, 3, n, true, true, pwl, 512), noconic_options) for n in vcat(4, 4:2:max_n)]
end

insts["sos2"] = modulardesign_insts(SOS2(), 30)
insts["logib"] = modulardesign_insts(LogIBBounded(), 30)
insts["cc"] = modulardesign_insts(CCBounded(), 18)

return (ModularDesign, insts)


text:
using AutoGrad
using Distributions
include("../utils/nn_utils.jl")
include("mixModel_grad_helper.jl")

function init_params(in_size, hidden_size, latent_size, n_comps, std=0.01)
  # encoder params
  w_encoder = append!([std*randn(in_size, hidden_size), std*randn(hidden_size, n_comps)], [std*randn(hidden_size, latent_size) for i=1:2*n_comps])
  b_encoder = append!([zeros(1, hidden_size), zeros(1, n_comps)], [zeros(1, latent_size) for i=1:2*n_comps])

  # decoder params
  w_decoder = Any[std*randn(latent_size, hidden_size), std*randn(hidden_size, in_size)]
  b_decoder = Any[zeros(1, hidden_size), zeros(1, in_size)]

  return Dict("w_encoder"=>w_encoder, "b_encoder"=>b_encoder, "w_decoder"=>w_decoder, "b_decoder"=>b_decoder)
end

# Splitting encoder and decoder fprops--will make grad easier
function fprop_encoder(params, x)
  n_components = size(params["w_encoder"][2], 2)
  latent_size = size(params["w_encoder"][3], 2)
  batchSize = size(x,1)

  h1 = relu(x*params["w_encoder"][1] .+ params["b_encoder"][1])
  return h1
end

function fprop_z(params, h1)
  # compute mixture parameters (all components -- needed for KL term)
  log_mix_weights = h1*params["w_encoder"][2] .+ params["b_encoder"][2]
  posterior = Dict("weights"=>softmax(log_mix_weights), "mu"=>[zeros(batchSize, latent_size) for i=1:n_components], "sigma"=>[zeros(batchSize, latent_size) for i=1:n_components])
  for k in 1:n_components
    posterior["mu"][k] += h1*params["w_encoder"][2*k+1] .+ params["b_encoder"][2*k+1]
    posterior["sigma"][k] += exp( h1*params["w_encoder"][2*k+2] .+ params["b_encoder"][2*k+2] )
  end

  # sample component index
  vals, idxs = findmax(log_mix_weights + rand(Gumbel(0,1), size(log_mix_weights)),2)
  comp_idxs = ind2sub(size(log_mix_weights), vec(idxs))[2]

  # sample latent variables from selected component distributions
  z = randn(batchSize, latent_size)
  for idx in 1:batchSize
    z[idx,:] = posterior["sigma"][comp_idxs[idx]][idx,:] .* z[idx,:] +  posterior["mu"][comp_idxs[idx]][idx,:]
  end

  return z, posterior
end

function fprop_decoder(params)

  h2 = relu(params["z"]*params["w_decoder"][1] .+ params["b_decoder"][1])
  x_recon = sigmoid(h2*params["w_decoder"][2] .+ params["b_decoder"][2])

  # clip the reconstruction to prevent NaN
  x_recon = min(max(x_recon, .000001), .999999)

  return x_recon
end

function fprop(params, x)
  z, posterior = fprop_encoder(params, x)
  x_recon = fprop_decoder(Dict("w_decoder"=>params["w_decoder"], "b_decoder"=>params["b_decoder"], "z"=>z))
  return x_recon, posterior, z
end

### KLD Terms ###
function gaussPDF(x, mu, sigma)
  return prod(1./sqrt(2*3.1416*sigma.^2) .* exp(-.5*(x - mu).^2 ./ sigma.^2), 2)
end

function gaussMixPDF(x, params)
  s = 0.
  for idx in 1:size(params["weights"],2)
    s += params["weights"][:,idx] .* gaussPDF(x, params["mu"][idx], params["sigma"][idx])
  end
  return s
end

function computeKLD(prior, posterior, x=None)
  # need MC approximation
  return log(gaussMixPDF(x, posterior))-log(gaussMixPDF(x, prior))
end

### Gradient Functions ###

function decoder_nll(params, x)
  # prop through decoder
  x_recon = fprop_decoder(params)

  # likelihood term
  expected_nll = sum(-x .* log(x_recon) - (1-x) .* log(1-x_recon), 2)

  return sum(expected_nll)
end

# get decoder grad via AutoGrad
nll_grad = grad(decoder_nll)
# ? does this work
kl_grad = grad(computeKLD)
# get grads for first weights
encoder_grad = grad(fprop_encoder)

function elbo_grad(params, x, priorParams)
  # get dimensions
  in_size, hidden_size = size(params["w_encoder"][1])
  n_components = size(params["w_encoder"][2],2)
  latent_size = size(params["w_encoder"][3],2)

  # init gradient dictionary
  grads = init_params(in_size, hidden_size, latent_size, n_components, 0.)

  # get derivatives of ELBO wrt decoder params and z
  h1 = fprop_encoder(params, x)
  z, posterior = fprop_z(params, h1)
  decoder_grads = nll_grad(Dict("w_decoder"=>params["w_decoder"], "b_decoder"=>params["b_decoder"], "z"=>z), x)
  kl_grads = kl_grad(Dict("w_encoder"=>params["w_encoder"], "b_encoder"=>params["b_encoder"], "z"=>z), x)

  # get derivatives of z wrt pi
  weight_grads = []
  for j in 1:n_components
    append!(weight_grads, compute_dh_dpi(z, params["weights"], posterior["mu"], posterior["cov"], j, 100))
  end

  # we have grads for dz/dpi now need grads of z wrt mu, sigma

  return grads
end


text:
module ExtendedFESpacesTests

using Test
using Gridap.Arrays
using Gridap.TensorValues
using Gridap.Algebra
using Gridap.ReferenceFEs
using Gridap.Geometry
using Gridap.FESpaces
using Gridap.FESpaces: ExtendedVector
using Gridap.Integration
using Gridap.Fields
using Gridap.CellData

oldcell_to_cell = [1,2,-1,3,-2,-3,-4]
cell_to_oldcell = [1,2,4]
void_to_oldcell = [3,5,6,7]
void_to_val = Float64[-10,-20,-30,-40]
cell_to_val = Float64[10,20,30]
a = ExtendedVector(
  void_to_val,
  cell_to_val,
  oldcell_to_cell,
  void_to_oldcell,
  cell_to_oldcell)

r = [10.0, 20.0, -10.0, 30.0, -20.0, -30.0, -40.0]
test_array(a,r)

n = 10
mesh = (n,n)
domain = (0,1,0,1) .- 1
order = 1
model = CartesianDiscreteModel(domain, mesh)

trian = Triangulation(model)

const R = 0.7

function is_in(coords)
  n = length(coords)
  x = (1/n)*sum(coords)
  d = x[1]^2 + x[2]^2 - R^2
  d < 0
end

oldcell_to_coods = get_cell_coordinates(trian)
oldcell_to_is_in = collect1d(apply(is_in,oldcell_to_coods))

incell_to_cell = findall(oldcell_to_is_in)
outcell_to_cell = findall(collect(Bool, .! oldcell_to_is_in))

model_in = DiscreteModel(model,incell_to_cell)

V = TestFESpace(
  model=model_in,valuetype=Float64,reffe=:Lagrangian,
  order=1,conformity=:H1,dof_space=:physical)
@test isa(V,ExtendedFESpace)

model_out = DiscreteModel(model,outcell_to_cell)

trian_in = RestrictedTriangulation(trian, incell_to_cell)

trian_Γ = InterfaceTriangulation(model_in,model_out)

quad_in = CellQuadrature(trian_in,2*order)

quad_Γ = CellQuadrature(trian_Γ,2*order)

reffes = [LagrangianRefFE(VectorValue{2,Float64},get_polytope(p),order) for p in get_reffes(trian_in)]

V_in = DiscontinuousFESpace(reffes,trian_in)

V = ExtendedFESpace(V_in, trian_in)
test_single_field_fe_space(V)

U = TrialFESpace(V)
test_single_field_fe_space(U)

u(x) = VectorValue(x[1]+x[2], x[1])

uh = interpolate(u,U)

uh_in = restrict(uh,trian_in)

uh_Γ = restrict(uh,trian_Γ)

t_in = AffineFETerm( (u,v) -> v⊙u, (v) -> v⊙u, trian_in, quad_in)
op_in = AffineFEOperator(U,V,t_in)

quad = CellQuadrature(trian,2*order)

t_Ω = AffineFETerm( (u,v) -> v⊙u, (v) -> v⊙u, trian, quad)
op_Ω = AffineFEOperator(U,V,t_Ω)

@test get_vector(op_in) ≈ get_vector(op_Ω)

t_Γ = AffineFETerm( (u,v) -> jump(v)⊙jump(u) + jump(ε(v))⊙jump(ε(u)), (v) -> jump(v)⊙u, trian_Γ, quad_Γ)
op_Γ = AffineFEOperator(U,V,t_Γ)

q_in = get_coordinates(quad_in)
collect(evaluate(uh_in,q_in))

q = get_coordinates(quad)
collect(evaluate(uh,q))

q_Γ = get_coordinates(quad_Γ)
collect(evaluate(jump(uh_Γ),q_Γ))
collect(evaluate(jump(∇(uh_Γ)),q_Γ))
collect(evaluate(jump(ε(uh_Γ)),q_Γ))

V = TestFESpace(model=model_in,valuetype=Float64,reffe=:Lagrangian,order=2,conformity=:H1)
@test isa(V,ExtendedFESpace)

V = TestFESpace(model=model_in,valuetype=Float64,reffe=:Lagrangian,order=2,conformity=:H1,constraint=:zeromean)
uh = FEFunction(V,rand(num_free_dofs(V)))
uh_in = restrict(uh,trian_in)
@test sum(integrate(uh_in,trian_in,quad_in)) + 1 ≈ 1

V = TestFESpace(model=model,valuetype=Float64,reffe=:Lagrangian,order=2,conformity=:H1)
@test !isa(V,ExtendedFESpace)

V = TestFESpace(triangulation=trian_in,valuetype=Float64,reffe=:Lagrangian,order=2,conformity=:L2)
@test isa(V,ExtendedFESpace)

V = TestFESpace(triangulation=trian_in,valuetype=Float64,reffe=:Lagrangian,order=2,conformity=:L2,constraint=:zeromean)
uh = FEFunction(V,rand(num_free_dofs(V)))
uh_in = restrict(uh,trian_in)
@test sum(integrate(uh_in,trian_in,quad_in)) + 1 ≈ 1

V = TestFESpace(triangulation=trian,valuetype=Float64,reffe=:Lagrangian,order=2,conformity=:L2)
@test !isa(V,ExtendedFESpace)

#using Gridap.Visualization
#writevtk(trian,"trian",cellfields=["uh"=>uh])
#writevtk(trian_in,"trian_in",cellfields=["uh"=>uh_in])

V_in = TestFESpace(model=model_in,valuetype=Float64,reffe=:Lagrangian,order=2,conformity=:H1)
V = TestFESpace(model=model,valuetype=Float64,reffe=:Lagrangian,order=2,conformity=:H1)

vh_in = interpolate(V_in) do x
    x[1]
end
vh_in = interpolate(vh_in, V_in)
vh = interpolate(vh_in, V)

#using Gridap.Visualization
#writevtk(trian,"trian",cellfields=["vh"=>vh,"vh_in"=>vh_in])

V_in = TestFESpace(
  model=model_in,valuetype=Float64,reffe=:Lagrangian,
  order=1,conformity=:H1,dof_space=:physical)

V_out = TestFESpace(
  model=model_out,valuetype=Float64,reffe=:Lagrangian,
  order=2,conformity=:H1,dof_space=:physical)

cell_axes = Fields.create_array_of_blocked_axes(get_cell_axes(V_in),get_cell_axes(V_out))
#using BenchmarkTools
#cache = array_cache(cell_axes)
#@btime getindex!($cache,$cell_axes,1)

end # module


text:
function grow_minus_clusters(xnew,ynew,kappanew,seeds_index,seednum,Lx,Ly)

xseeds = xnew[seeds_index]
yseeds = ynew[seeds_index]
size_check = size(xseeds)
if  size_check[2] == 1
    xseeds = xseeds'
    yseeds = yseeds'
end

xc = xseeds[seednum,:]
yc = yseeds[seednum,:]
nc = 2; # initial cluster charge is 2

# This defines the current vortex indices associated with the cluster
cluster_index = seeds_index[seednum,:]
num_added_to_cluster = nc # seeds are added to the cluster

xp = @. xnew[kappanew == 1.0]
yp = @. ynew[kappanew == 1.0]
xn = @. xnew[kappanew == -1.0]
yn = @. ynew[kappanew == -1.0]
failed_attempts = 0

while (num_added_to_cluster > 0 && failed_attempts < 10)
    # Remove the vortices that are already in the cluster from consideration
        xtemp = xn
        ytemp = yn
        for j = 1:nc
            Ij = ( (xc[j] == xn) & (yc[j] == yn) )
            xtemp[Ij] = NaN
            ytemp[Ij] = NaN

            # Work out the distances to all other vortices
            xdiff = @. xtemp - xc'
            ydiff = @. ytemp - yc'
        end

    # incorporate periodicity
    xdiff,ydiff = periodic_distances(xdiff,ydiff,Lx,Ly)
    rij = hypot.(xdiff,ydiff)

    # Nearest same-sign vortex to each vortex, numbered by x
    rmin,candidate_index = minimum(rij)
    # candidate_index = unique(candidate_index)

   ## Find the nearest opposites to the cluster
    xdiff3 = @. xp - xc'
    ydiff3 = @. yp - yc'
    rij3 = hypot.(xdiff3,ydiff3)

    r_min3,_ = minimum(rij3)

    # candidate_index = candidate_index(rmin < r_min3);

    ## Find the nearest opposite signed vortex of each candidate
    number_of_candidates = length(candidate_index)

    x_candidates = xn[candidate_index]
    y_candidates = yn[candidate_index]
    xdiff2 = @. xp - x_candidates'
    ydiff2 = @. yp - y_candidates'
    xdiff2,ydiff2 = periodic_distances(xdiff2,ydiff2,Lx,Ly)

    rij2 = hypot.(xdiff2,ydiff2)
    rmin_opposite,_ = minimum(rij2)

    add_to_cluster = ( (rmin < rmin_opposite) & (rmin < r_min3) )

    add_to_cluster = ( (rmin < rmin_opposite) && (rmin < r_min3) )
    num_added_to_cluster = length(unique(candidate_index[add_to_cluster]))
    cluster_index = [cluster_index unique(candidate_index[add_to_cluster])] ##ok

    xc = xnew[cluster_index]'
    yc = ynew[cluster_index]'
    nc = length(xc)
#     disp(['found ' num2str(number_of_candidates) ' candidates'])
#     disp(['added ' num2str(num_added_to_cluster) ' candidates to cluster'])
end
xc = xc';
yc = yc';
cluster_index = cluster_index';
return xc,yc,cluster_index
end


text:
### utils.jl
#
# Copyright (C) 2016 Mosè Giordano.
#
# Maintainer: Mosè Giordano <mose AT gnu DOT org>
# Keywords: uncertainty, error propagation, physics
#
# This file is a part of Measurements.jl.
#
# License is MIT "Expat".
#
### Commentary:
#
# This file defines some utility functions.
#
### Code:

export stdscore, weightedmean, value, uncertainty

# Standard Score
"""
    stdscore(measure::Measurement, expected_value::Real) -> standard_score

Gives the standard score between a measure, with uncertainty, and its expected
value (that may or may not have the uncertainty):

    (measure.val - expected_value)/measure.err
"""
stdscore(a::Measurement, b::Real) = (a.val - b)/(a.err)

# Weighted Average with Inverse-Variance Weighting
"""
    weightedmean(iterable) -> measurement(weighted_mean, standard_deviation)

Return the weighted mean of measurements listed in `iterable`, using
inverse-variance weighting.
"""
function weightedmean(iterable)
    v = [el.val for el in iterable]
    w = [inv(el.err)^2 for el in iterable]
    invsumw = inv(sum(w))
    return measurement(dot(v, w)*invsumw, sqrt(invsumw))
end

# Derivative and Gradient
derivative{F<:AbstractFloat, T<:AbstractFloat}(a::Measurement{F},
                                               tag::Tuple{T, T, Float64}) =
                                                   get(a.der, tag, zero(F))

"""
    derivative(x::Measurement, y::Measurement)

Return the value of the partial derivative of `x` with respect to the
independent measurement `y`, calculated on the nominal value of `y`.  Return
`0.0` if `x` does not depend on `y`.

Use `Measurements.gradient` to calculate the gradient of `x` with respect to an
arrays of independent measurements.
"""
derivative(a::Measurement, b::Measurement) =
    derivative(a, (b.val, b.err, b.tag))

"""
    gradient(x::Measurement, [y::AbstractArray{Measurement}])

Return the gradient of `x` with respect to the independent measurements
contained in `y` array, each calculated on the nominal value of the
corresponding measurement.  Return `0.0` for the variables on which `x` does not
depend.

Use `Measurements.derivative` to calculate the partial derivative of `x` with
respect to a single independent measurement.
"""
function gradient{F<:AbstractFloat, T<:Measurement}(a::Measurement{F},
                                                    b::AbstractArray{T})
    out = similar(b, F)
    for i in eachindex(b)
        out[i] = derivative(a, b[i])
    end
    return out
end

# value and uncertainty
for (f, field) in ((:value, :val), (:uncertainty, :err))
    @eval begin
        ($f)(a::Measurement) = a.$field
        ($f){T<:AbstractFloat}(a::Complex{Measurement{T}}) =
            complex(($f)(a.re), ($f)(a.im))

        function ($f){T<:AbstractFloat}(A::AbstractArray{Measurement{T}})
            out = similar(A, T)
            for i in eachindex(A)
                out[i] = ($f)(A[i])
            end
            return out
        end

        function ($f){T<:AbstractFloat}(A::AbstractArray{Complex{Measurement{T}}})
            out = similar(A, Complex{T})
            for i in eachindex(A)
                out[i] = ($f)(A[i])
            end
            return out
        end
    end
end

"""
    value(x::Measurement)
    value(x::Complex{Measurement})
    value(x::AbstractArray{Measurement})
    value(x::AbstractArray{Complex{Measurement}})

Return the nominal value of measurement `x`.  If `x` is an array, return an
array of the same length as `x` with the nominal values of each elements of `x`.
"""
value

"""
    uncertainty(x::Measurement)
    uncertainty(x::Complex{Measurement})
    uncertainty(x::AbstractArray{Measurement})
    uncertainty(x::AbstractArray{Complex{Measurement}})

Return the uncertainty of measurement `x`.  If `x` is an array, return an array
of the same length as `x` with the uncertainty of each elements of `x`.
"""
uncertainty


text:
using Metatheory

failme = @theory p begin
    p ≠ ¬p
    :foo == ¬:foo
    :foo --> :bazoo
    :bazoo --> :wazoo
end

g = EGraph(:foo)
report = saturate!(g, failme)
println(report)
@test report.reason === :contradiction
# @test !(@areequal failme foo wazoo)


text:
# This is for the PriorityQueue
using DataStructures

struct Leaf
    weight::Int64
    key::Char
end

struct Branch
    right::Union{Leaf, Branch}
    left::Union{Leaf, Branch}
    weight::Int64
end

const Node = Union{Leaf, Branch}
isbranch(branch::Branch) = true
isbranch(other::T) where {T} = false

function codebook_recurse!(leaf::Leaf, code::String,
                          dict::Dict{Char,String})
    dict[leaf.key] = code
end

function codebook_recurse!(branch::Branch, code::String,
                          dict::Dict{Char,String})
    codebook_recurse!(branch.left, string(code, "1"), dict)
    codebook_recurse!(branch.right, string(code, "0"), dict)
end

# This will depth-first search through the tree
# to create bitstrings for each character.
# Note: Any depth-first search method will work
# This outputs encoding Dict to be used for encoding
function create_codebook(n::Node)
    codebook = Dict{Char,String}()
    codebook_recurse!(n, "", codebook)
    return codebook
end

# This outputs huffman tree to generate dictionary for encoding
function create_tree(phrase::String)

    # creating weights
    weights = PriorityQueue()
    for i in phrase
        temp_string = string(i)
        if (haskey(weights, temp_string))
            weights[temp_string] += 1
        else
            weights[temp_string] = 1
        end
    end

    # Creating all nodes to iterate through
    nodes = PriorityQueue{Node, Int64}()
    while(length(weights) > 0)
        weight = peek(weights)[2]
        key = dequeue!(weights)[1]
        temp_node = Leaf(weight, key)
        enqueue!(nodes, temp_node, weight)
    end

    while(length(nodes) > 1)
        node1 = dequeue!(nodes)
        node2 = dequeue!(nodes)
        temp_node = Branch(node1, node2, node1.weight + node2.weight)
        enqueue!(nodes, temp_node, temp_node.weight)
    end

    huffman_tree = dequeue!(nodes)
    return huffman_tree

end

function encode(codebook::Dict{Char, String}, phrase::String)
    final_bitstring = ""
    for i in phrase
        final_bitstring = final_bitstring * codebook[i]
    end

    return final_bitstring
end

function decode(huffman_tree::Node, bitstring::String)
    current = huffman_tree
    final_string = ""
    for i in bitstring
        if (i == '1')
            current = current.left
        else
            current = current.right
        end
        if (!isbranch(current))
            final_string = final_string * string(current.key)
            current = huffman_tree
        end
    end

    return final_string
end

function two_pass_huffman(phrase::String)
    huffman_tree = create_tree(phrase)
    codebook = create_codebook(huffman_tree)
    println(codebook)
    bitstring = encode(codebook, phrase)
    final_string = decode(huffman_tree, bitstring)
    println(bitstring)
    println(final_string)
end

two_pass_huffman("bibbity bobbity")


text:
# ---
# title: 1200. Minimum Absolute Difference
# id: problem1200
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Array
# link: <https://leetcode.com/problems/minimum-absolute-difference/description/>
# hidden: true
# ---
#
# Given an array of **distinct**  integers `arr`, find all pairs of elements
# with the minimum absolute difference of any two elements.
#
# Return a list of pairs in ascending order(with respect to pairs), each pair
# `[a, b]` follows
#
#   * `a, b` are from `arr`
#   * `a < b`
#   * `b - a` equals to the minimum absolute difference of any two elements in `arr`
#
#
#
# **Example 1:**
#
#
#
#     Input: arr = [4,2,1,3]
#     Output: [[1,2],[2,3],[3,4]]
#     Explanation: The minimum absolute difference is 1. List all pairs with difference equal to 1 in ascending order.
#
# **Example 2:**
#
#
#
#     Input: arr = [1,3,6,10,15]
#     Output: [[1,3]]
#
#
# **Example 3:**
#
#
#
#     Input: arr = [3,8,-10,23,19,-4,-14,27]
#     Output: [[-14,-10],[19,23],[23,27]]
#
#
#
#
# **Constraints:**
#
#   * `2 <= arr.length <= 10^5`
#   * `-10^6 <= arr[i] <= 10^6`
#
#
## @lc code=start
using LeetCode

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


text:
#==============================================================================#
# AWSSES.jl
#
# SES API. See http://aws.amazon.com/documentation/ses/
#
# Copyright Sam O'Connor 2016 - All rights reserved
#==============================================================================#

__precompile__()

module AWSSES

export ses_send, ses_send_raw, ses_send_attachments

using AWSCore
using SymDict
using Retry

function ses(aws::AWSConfig, query)

    query["AWSAccessKeyId"] = aws[:creds].secret_key

    if aws[:creds].token != ""
        query["SecurityToken"] = aws[:creds].token
    end

    do_request(post_request(aws, "email", "", query))
end

function ses_send(aws::AWSConfig; to="", from="", subject="", body="")

    ses(aws, Dict(
        "Action" => "SendEmail",
        "Source" => from,
        "Destination.ToAddresses.member.1" => to,
        "Message.Subject.Data" => subject,
        "Message.Body.Text.Data" => body))
end

function ses_send_raw(aws::AWSConfig; to="", from="", raw="")

    ses(aws, Dict(
        "Action" => "SendRawEmail",
        "Source" => from,
        "Destination.ToAddresses.member.1" => to,
        "RawMessage.Data" => base64encode(raw)))
end

#    e.g. attachments=[("bar.txt", "text/plain", "bar\n"),
#                      ("foo.txt", "text/plain", "foo\n")])

function  ses_send_attachments(aws::AWSConfig;
                               to="", from="", subject="", body="",
                               attachments=[])

    ses_send_raw(
        aws,
        to=to,
        from=from,
        raw = AWSCore.mime_multipart(
            """
            To: $to
            From: $from
            Subject: $subject
            """,
            [("", "text/plain", body), attachments...]))
end

end # module AWSSES

#==============================================================================#
# End of file.
#==============================================================================#


text:
export IndPolyhedral

abstract type IndPolyhedral <: ProximableFunction end

is_convex(::IndPolyhedral) = true
is_set(::IndPolyhedral) = true

"""
**Indicator of a polyhedral set**

    IndPolyhedral([l,] A, [u, xmin, xmax])

Returns the indicator function of the polyhedral set:
```math
S = \\{ x : x_\\min \\leq x \\leq x_\\max, l \\leq Ax \\leq u \\}.
```
Matrix `A` is a mandatory argument; when any of the bounds is not provided,
it is assumed to be (plus or minus) infinity.
"""
function IndPolyhedral(args...; solver=:osqp)
    if solver == :osqp
        IndPolyhedralOSQP(args...)
    else
        error("unknown solver")
    end
end

# including concrete types

include("indPolyhedralOSQP.jl")


text:
using RefractiveIndex
using Documenter

DocMeta.setdocmeta!(RefractiveIndex, :DocTestSetup, :(using RefractiveIndex); recursive=true)

makedocs(;
    modules=[RefractiveIndex],
    authors="Alex Ames <alexander.m.ames@gmail.com> and contributors",
    repo="https://github.com/stillyslalom/RefractiveIndex.jl/blob/{commit}{path}#{line}",
    sitename="RefractiveIndex.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://stillyslalom.github.io/RefractiveIndex.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/stillyslalom/RefractiveIndex.jl",
)


text:
import Statistics: mean

import Base: show, IO
function show_events_array(io::IO, v::Vector{Event}, i1=1, l=length(v))
    i = i1
    if l >= i
        while true
            show(io, v[i])
            i += 1
            i > l && break

            print(io, ", ")
        end
    end
end

function show(io::IO, evtset::EventSet)
    limited = get(io, :limit, false)::Bool

    print(io, "EventSet: [")
    if limited && length(evtset) > 20
        f, l = 1, length(evtset)
        show_events_array(io, evtset.events, f, f+9)
        print(io, "  …  ")
        show_events_array(io, evtset.events, l-9, l)
    else
        show_events_array(io, evtset.events)
    end
    print(io, "]")
end

function show(io::IO, ::MIME"text/plain", stats::EventValues)
    print(io, "EventValues:")
    for (e,v) in zip(stats.events, stats.vals)
        print(io, "\n  ", e, " = ", v)
        print_shadow(io, e, v, stats)
    end
    print(io, "\n  runtime = ", stats.time, " nsecs")
end

function show(io::IO, stats::EventValues)
    print(io, "EventValues")
    if !get(io, :compact, false)
        for (e,v) in zip(stats.events, stats.vals)
            print(io, " ", e, "=", v)
        end
    end
    print(io, "\n  runtime = ", stats.time, " nsecs")
end

function show(io::IO, ::MIME"text/plain", stats::EventStats)
    print(io, "EventStats:")
    for (e,v) in zip(stats.events, eachcol(stats.samples))
        print(io, "\n  ", e, " = ", v)
        print_shadow(io, e, v, stats)
    end
    print(io, "\n  runtime = ", stats.time, " nsecs")
end

function show(io::IO, stats::EventStats)
    print(io, "EventStats")
    if !get(io, :compact, false)
        for (e,v) in zip(stats.events, eachcol(stats.samples))
            print(io, " ", e, "=", v)
        end
    end
    print(io, "\n  runtime = ", stats.time, " nsecs")
end

function has_event(f::Function, evt::Event, stats::EventValues)
    idx = findfirst(isequal(evt), stats.events)
    if idx !== nothing
        f(stats.vals[idx])
    end
end

function has_event(f::Function, evt::Event, stats::EventStats)
    idx = findfirst(isequal(evt), stats.events)
    if idx !== nothing
        f(stats.samples[:, idx])
    end
end

const COLOR_NONE = ""
const COLOR_RESET = "\033[m"
const COLOR_BOLD = "\033[1m"
const COLOR_RED	= "\033[31m"
const COLOR_GREEN = "\033[32m"
const COLOR_YELLOW = "\033[33m"
const COLOR_BLUE = "\033[34m"
const COLOR_MAGENTA ="\033[35m"
const COLOR_CYAN = "\033[36m"
const COLOR_BG_RED = "\033[41m"

pct_color(ratio) = if ratio > 20
    COLOR_RED
elseif ratio > 10
    COLOR_MAGENTA
elseif ratio > 5
    COLOR_YELLOW
else
    COLOR_NONE
end

function ratio(num::Union{Counts, AbstractVector{Counts}}, den::Union{Counts, AbstractVector{Counts}})
    round(mean(num./ den), digits=3)
end

function percentage(num::Union{Counts, AbstractVector{Counts}}, den::Union{Counts, AbstractVector{Counts}})
    round(mean(num./ den) * 100, digits=0)
end

function print_shadow(io::IO, evt::Event, value::Union{Counts, AbstractVector{Counts}}, stats::Union{EventStats, EventValues})
    if evt == TOT_INS
        has_event(TOT_CYC, stats) do cycles
            print(io, " # $(ratio(value, cycles)) insn per cycle")
        end

        has_event(RES_STL, stats) do stalled_cycles
            print(io, " # $(ratio(stalled_cycles, value)) stalled cycles per insn")
        end
    elseif evt == TOT_CYC
        print(io, " # $(ratio(value, stats.time)) Ghz")

        has_event(TOT_INS, stats) do cycles
            print(io, " # $(ratio(value, cycles)) cycles per insn")
        end
    elseif evt == BR_MSP
        has_event(BR_INS, stats) do branches
            pct = percentage(value, branches)
            print(io, " # $(pct_color(pct))$pct%$COLOR_RESET of all branches")
        end
    elseif evt == L1_LDM || evt == L2_LDM || evt == L3_LDM
        has_event(LD_INS, stats) do loads
            pct = percentage(value, loads)
            print(io, " # $(pct_color(pct))$pct%$COLOR_RESET of all loads")
        end
    elseif evt == L1_STM || evt == L2_STM || evt == L3_STM
        has_event(SR_INS, stats) do stores
            pct = percentage(value, stores)
            print(io, " # $(pct_color(pct))$pct%$COLOR_RESET of all stores")
        end
    elseif evt == L1_TCM || evt == L1_TCH
        has_event(L1_TCA, stats) do accesses
            pct = percentage(value, accesses)
            color = pct_color(evt == L1_TCM ? pct : 1.0 - pct)
            print(io, " # $color$pct%$COLOR_RESET of all L1 cache refs")
        end
    elseif evt == L1_DCM || evt == L1_DCA || evt == L1_DCW || evt == L1_DCR
        has_event(L1_TCA, stats) do accesses
            pct = percentage(value, accesses)
            print(io, " # $pct% of all L1 cache refs")
        end

        if evt == L1_DCM
            has_event(L1_DCA, stats) do accesses
                pct = percentage(value, accesses)
                print(io, " # $(pct_color(pct))$pct%$COLOR_RESET of data L1 cache refs")
            end
        end
    elseif evt == L2_TCM || evt == L2_TCH
        has_event(L2_TCA, stats) do accesses
            pct = percentage(value, accesses)
            color = pct_color(evt == L2_TCM ? pct : 1.0 - pct)
            print(io, " # $color$pct%$COLOR_RESET of all L2 cache refs")
        end
    elseif evt == L2_DCM || evt == L2_DCA || evt == L2_DCW || evt == L2_DCR
        has_event(L2_TCA, stats) do accesses
            pct = percentage(value, accesses)
            print(io, " # $pct% of all L2 cache refs")
        end

        if evt == L2_DCM
            has_event(L2_DCA, stats) do accesses
                pct = percentage(value, accesses)
                print(io, " # $(pct_color(pct))$pct%$COLOR_RESET of data L2 cache refs")
            end
        end
    elseif evt == L3_TCM || evt == L3_TCH
        has_event(L3_TCA, stats) do accesses
            pct = percentage(value, accesses)
            color = pct_color(evt == L3_TCM ? pct : 1.0 - pct)
            print(io, " # $color$pct%$COLOR_RESET of all L3 cache refs")
        end
    elseif evt == L3_DCM || evt == L3_DCA || evt == L3_DCW || evt == L3_DCR
        has_event(L3_TCA, stats) do accesses
            pct = percentage(value, accesses)
            print(io, " # $pct% of all L3 cache refs")
        end

        if evt == L3_DCM
            has_event(L3_DCA, stats) do accesses
                pct = percentage(value, accesses)
                print(io, " # $(pct_color(pct))$pct%$COLOR_RESET of data L3 cache refs")
            end
        end
    elseif evt == TLB_DM
        has_event(LST_INS, stats) do lst
            pct = percentage(value, lst)
            print(io, " # $(pct_color(pct))$pct%$COLOR_RESET of all load/stores")
        end
    elseif evt in (FMA_INS, FP_INS, BR_INS, VEC_INS, FAD_INS, SR_INS, LD_INS, INT_INS, LST_INS, SYC_INS, FML_INS, FDV_INS, FSQ_INS, FNV_INS)
        has_event(TOT_INS, stats) do insn
            pct = percentage(value, insn)
            print(io, " # $pct% of all instructions")
        end

        unit = 'G'
        r = ratio(value, stats.time)
        if r < 1.
            unit = 'M'
            r *= 1e3
        elseif r < 1e-6
            unit = 'K'
            r *= 1e6
        end
        print(io, " # $r $unit/sec ")
    elseif evt in (STL_ICY, FUL_ICY, STL_CCY, FUL_CCY, MEM_SCY, MEM_RCY, MEM_WCY)
        has_event(TOT_CYC, stats) do cycles
            pct = percentage(value, cycles)
            print(io, " # $pct% of all cycles")
        end
    elseif evt == SP_OPS || evt == DP_OPS
        has_event(FP_INS, stats) do insn
            pct = percentage(value, insn)
            print(io, " # $pct% of all fp instructions")
        end

        has_event(TOT_INS, stats) do insn
            pct = percentage(value, insn)
            print(io, " # $pct% of all instructions")
        end

        has_event(evt == SP_OPS ? VEC_SP : VEC_DP, stats) do vec
            x = ratio(value, vec)
            print(io, " # $(x)x vectorized")
        end
    elseif evt == VEC_SP || evt == VEC_DP
        has_event(VEC_INS, stats) do insn
            pct = percentage(value, insn)
            print(io, " # $pct% of all vec instructions")
        end

        has_event(TOT_INS, stats) do insn
            pct = percentage(value, insn)
            print(io, " # $pct% of all instructions")
        end
    elseif evt in (MEM_RCY, MEM_SCY, MEM_WCY, RES_STL)
        has_event(TOT_CYC, stats) do cycles
            pct = percentage(value, cycles)
            color = if pct > 50
                COLOR_RED
            elseif pct > 30
                COLOR_MAGENTA
            elseif pct > 10
                COLOR_YELLOW
            else
                COLOR_NONE
            end
            print(io, " # $color$pct%$COLOR_RESET of all cycles")
        end
    else
        unit = 'G'
        r = ratio(value, stats.time)
        if r < 1.
            unit = 'M'
            r *= 1e3
        elseif r < 1e-6
            unit = 'K'
            r *= 1e6
        end

        print(io, " # $r $unit/sec ")
    end

    if evt in (L1_DCM, L1_DCA, L1_DCW, L1_DCR, L2_DCM, L2_DCA, L2_DCW, L2_DCR, L3_DCM, L3_DCA, L3_DCW, L3_DCR)
        has_event(LST_INS, stats) do lst
            pct = percentage(value, lst)
            print(io, " # $pct of all load/stores")
        end
    end
end


text:
struct OutputFlags
    write_volt_maps::Bool
    write_cur_maps::Bool
    write_cum_cur_map_only::Bool
    write_max_cur_maps::Bool
    set_null_currents_to_nodata::Bool
    set_null_voltages_to_nodata::Bool
    compress_grids::Bool
    log_transform_maps::Bool
end

function compute_3col(resistances::Matrix{T}) where {T}
    fp = deleteat!(resistances[:,1], 1)
    l = length(fp)
    r3col = zeros(T, div(l * (l-1), 2), 3)
    k = 1
    for i = 1:l
        for j = i+1:l
            r3col[k,1] = fp[i]
            r3col[k,2] = fp[j]
            r3col[k,3] = resistances[j+1,i+1]
            k += 1
        end
    end
    r3col
end

function write_cur_maps(name, output, component_data, finitegrounds, flags, cfg)

    # Get desired data
    G = component_data.matrix
    voltages = flags.is_advanced ? output : output.voltages
    cc = component_data.cc
    nodemap = component_data.local_nodemap
    hbmeta = component_data.hbmeta
    cellmap = component_data.cellmap

    # Flags
    log_transform = flags.outputflags.log_transform_maps
    set_null_currents_to_nodata =
        flags.outputflags.set_null_currents_to_nodata
    write_max_cur_maps = flags.outputflags.write_max_cur_maps
    write_cum_cur_map_only = flags.outputflags.write_cum_cur_map_only

    node_currents, branch_currents = _create_current_maps(G, voltages, finitegrounds, cfg, nodemap = nodemap, hbmeta = hbmeta)

    if !flags.is_raster

        # Branch currents
        branch_currents_array = _convert_to_3col(branch_currents, cc)

        # TODO: implement cumulative maps for netowrk mode

        # Accumulate branch currents
        # cum_branch_curr[mycsid()] .+= branch_currents_array[:,3]

        # Node currents
        node_currents_array = _append_name_to_node_currents(node_currents, cc)

        # Accumulate node currents
        # cum_node_curr[mycsid()] .+= node_currents_array[:,2]

        # !write_cum_cur_map_only &&
        write_currents(node_currents_array, branch_currents_array, name, cfg)
    else

        cmap = node_currents
        cum_curr = output.cum.cum_curr
        max_curr = output.cum.max_curr

        # Process the current map
        process_grid!(cmap, cellmap, hbmeta, log_transform = log_transform,
                            set_null_to_nodata = set_null_currents_to_nodata)

        # Accumulate by default
        cum_curr[mycsid()] .+= cmap

        # Max current if user asks for it
        if write_max_cur_maps
            max_curr[mycsid()] .= max.(max_curr[mycsid()], cmap)
        end

        # Write current maps
        !write_cum_cur_map_only &&
                        write_grid(cmap, name, cfg, hbmeta)

   end
end

function write_currents(node_curr_arr, branch_curr_arr, name, cfg)
   pref = split(cfg["output_file"], ".out")[1]
   writedlm("$(pref)_node_currents$(name).txt", node_curr_arr, '\t')
   writedlm("$(pref)_branch_currents$(name).txt", branch_curr_arr, '\t')
end

_append_name_to_node_currents(node_currents, cc) = [cc node_currents]

function _convert_to_3col(branch_currents, cc)

    l = length(branch_currents.nzval)
    graph = zeros(l, 3)
    graph[:,3] = branch_currents.nzval

    # Inspired by show method for sparse matrices in Julia Base
    b = branch_currents

    k = 1
    for i = 1:size(b, 1)
        for j in nzrange(b, i)
            row = b.rowval[j]
            graph[k,1] = cc[row]
            graph[k,2] = cc[i]
            k += 1
        end
    end

    graph
end

function _create_current_maps(G, voltages, finitegrounds, cfg; nodemap = Matrix{Float64}(), hbmeta = RasterMeta())

    node_currents = get_node_currents(G, voltages, finitegrounds)

    if cfg["data_type"] == "network"

        branch_currents = _get_branch_currents(G, voltages, true)
        branch_currents = abs.(branch_currents)
        return node_currents, branch_currents

    else

        current_map = zeros(eltype(G), hbmeta.nrows, hbmeta.ncols)
        for j = 1:size(nodemap, 2)
            for i = 1:size(nodemap, 1)
                idx = nodemap[i,j]
                if idx == 0
                    continue
                else
                    current_map[i,j] = node_currents[idx]
                end
            end
        end

        return current_map, spzeros(0,0)
    end
end

function get_node_currents(G, voltages, finitegrounds)

    node_currents_pos = _get_node_currents_posneg(G, voltages, finitegrounds, true)
    node_currents_neg = _get_node_currents_posneg(G, voltages, finitegrounds, false)
    node_currents = map((x,y) -> x > y ? x : y, node_currents_pos, node_currents_neg)

end

function _get_node_currents_posneg(G::SparseMatrixCSC{T,V},
                            voltages, finitegrounds, pos) where {T,V}

    branch_currents = _get_branch_currents(G, voltages, pos)
    branch_currents = branch_currents - branch_currents'
    dropnonzeros!(branch_currents)

	if finitegrounds[1]!= -9999
        finiteground_currents = finitegrounds .* voltages
        if pos
            map!(x -> x < 0 ? -x : 0, finiteground_currents, finiteground_currents)
        else
            map!(x -> x > 0 ? x : 0, finiteground_currents, finiteground_currents)
        end
        n = size(G, 1)
        branch_currents = branch_currents + spdiagm(0 => finiteground_currents)
    end

    s = vec(sum(branch_currents, dims=1))

    s
end

function dropnonzeros!(G)
    for i = 1:size(G, 1)
        for j in nzrange(G, i)
            row = G.rowval[j]
            val = G.nzval[j]
            if val < 0
                G.nzval[j] = 0
            end
        end
    end
    dropzeros!(G)
end

function _get_branch_currents(G::SparseMatrixCSC{T,V}, voltages, pos) where {T,V}

    branch_currents = _get_branch_currents_posneg(G, voltages, pos)

    # Make sparse matrix with branch_currents as right upper triangle
    N = size(G, 1)
    n = size(branch_currents, 1)
    I = zeros(V, n)
    J = zeros(V, n)
    k = 1
    for i = 1:N
        for j in nzrange(G, i)
            row = G.rowval[j]
            if i > row
                I[k] = row
                J[k] = i
                k += 1
            end
        end
    end
    @assert n + 1 == k

    B = sparse(I, J, branch_currents, N, N)

    B
	# branch_currents
end

function _get_branch_currents_posneg(G::SparseMatrixCSC{T,V},
                                v::Vector{T}, pos) where {T,V}

    #=I,J,V = findnz(G)
    mask = I .< J
    @show sum(mask)
    vdiff = zeros(T, sum(mask))
    if pos
        vdiff = v[I[mask]] - v[J[mask]]
        #for (i,v) in enumerate(find(mask))
        #    vdiff[i] = v[I[v]] - v[J[v]]
        #end
    else
        vdiff = v[J[mask]] - v[I[mask]]
        #for (i,v) in enumerate(find(mask))
        #    vdiff[i] = v[J[v]] - v[I[v]]
        #end
    end
    map!(x -> x < 0 ? -x : 0, V, V)

    branch_currents = vdiff .* V[mask]
    maxcur = maximum(branch_currents)
    map!(x -> abs(x / maxcur) < 1e-8 ? 0 : x, branch_currents, branch_currents)=#

    n = count_upper(G)
    b = zeros(T, n)
    k = 1
    if pos
        for i = 1:size(G, 1)
            for j in nzrange(G, i)
                row = G.rowval[j]
                val = G.nzval[j]
                if i > row
                    b[k] = abs(val) * (v[row] - v[i])
                    k += 1
                end
            end
        end
    else
        for i = 1:size(G, 1)
            for j in nzrange(G, i)
                row = G.rowval[j]
                val = G.nzval[j]
                if i > row
                    b[k] = abs(val) * (v[i] - v[row])
                    k += 1
                end
            end
        end
    end
    # @show n, k
    @assert n + 1 == k
    maxcur = maximum(b)
    # map!(x -> abs(x / maxcur) < 1e-8 ? 0 : x, b, b)
    for i = 1:size(b, 1)
        if abs(b[i]/ maxcur) < 1e-8
            b[i] = 0
        end
    end

    # @show sum(abs2, b - branch_currents)
    # branch_currents
    b
end

function count_upper(G)
    n = 0
    for i = 1:size(G, 1)
        for j in nzrange(G,i)
            row = G.rowval[j]
            if i > row
                n += 1
            end
        end
    end
    n
end

function process_grid!(cmap, cellmap, hbmeta; log_transform = false,
                                set_null_to_nodata = false)
    if log_transform
        map!(x -> x > 0 ? log10(x) : float(hbmeta.nodata), cmap, cmap)
    end

    if set_null_to_nodata
        for i in eachindex(cmap)
            if cellmap[i] == 0
                cmap[i] = hbmeta.nodata
            end
        end
    end

end

function write_grid(cmap, name, cfg, hbmeta;
                        voltage = false, cum = false,
                        max = false)

    str = "curmap"
    if cum
        str = "cum_$(str)"
    elseif max
        str = "max_$(str)"
    elseif voltage
        str = "voltmap"
    end

    pref = split(cfg["output_file"], ".out")[1]
	if hbmeta.file_type == FILE_TYPE_GEOTIFF
		filename = "$(pref)_$(str)$(name).tif"
		write_geotiff(cmap,filename,hbmeta)
	else
		filename = "$(pref)_$(str)$(name).asc"
		write_aagrid(cmap,filename,hbmeta)
	end
end

function write_grid(cmap, name, cfg, hbmeta, cellmap;
                        voltage = false, cum = false, max = false,
                        log_transform = false, set_null_to_nodata = false)

    pref = split(cfg["output_file"], ".out")[1]

    if log_transform
        map!(x -> x > 0 ? log10(x) : float(hbmeta.nodata), cmap, cmap)
    end

    if set_null_to_nodata
        for i in eachindex(cmap)
            if cellmap[i] == 0
                cmap[i] = hbmeta.nodata
            end
        end
    end

    str = "curmap"
    if cum
        str = "cum_$(str)"
    elseif max
        str = "max_$(str)"
    elseif voltage
        str = "voltmap"
    end

    filename = "$(pref)_$(str)$(name).asc"
	if hbmeta.file_type == FILE_TYPE_GEOTIFF
		filename = "$(pref)_$(str)$(name).tif"
		write_geotiff(cmap,filename,hbmeta)
	else
		filename = "$(pref)_$(str)$(name).asc"
		write_aagrid(cmap,filename,hbmeta)
	end
end

function write_aagrid(cmap,filename,hbmeta)
	f = open(filename, "w")

	write(f, "ncols         $(hbmeta.ncols)\n")
	write(f, "nrows         $(hbmeta.nrows)\n")
	write(f, "xllcorner     $(hbmeta.xllcorner)\n")
	write(f, "yllcorner     $(hbmeta.yllcorner)\n")
	write(f, "cellsize      $(hbmeta.cellsize)\n")
	write(f, "NODATA_value  $(hbmeta.nodata)\n")

	writedlm(f, round.(cmap, digits=8), ' ')
	close(f)
end

function write_geotiff(cmap, filename, hbmeta)
	ga = GeoArray(permutedims(Array{Float32}(cmap)))
	ga.crs = hbmeta.crs
	ga.f = hbmeta.affine_map
	GeoArrays.write!(filename, ga)
end

function write_volt_maps(name, output, component_data, flags, cfg)

    voltages = flags.is_advanced ? output : output.voltages

    if !flags.is_raster

        cc = component_data.cc
        write_voltages(cfg["output_file"], name, voltages, cc)

    else

        # Desired data
        # voltages = output.voltages
        cc = component_data.cc
        hbmeta = component_data.hbmeta
        nodemap = component_data.local_nodemap
        set_null_voltages_to_nodata = flags.outputflags.set_null_voltages_to_nodata

        vm = _create_voltage_map(voltages, nodemap, hbmeta)
        write_grid(vm, name, cfg, hbmeta, component_data.cellmap, voltage = true,
                        set_null_to_nodata = set_null_voltages_to_nodata)
    end
end

function write_voltages(output, name, voltages::Vector{T}, cc) where {T}

    volt_arr = zeros(T, size(voltages, 1), 2)
    volt_arr[:,1] = cc
    volt_arr[:,2] = voltages

    pref = split(output, ".out")[1]
    writedlm("$(pref)_voltages$(name).txt", volt_arr)

end

function _create_voltage_map(voltages::Vector{T}, nodemap, hbmeta) where {T}
    voltmap = zeros(T, hbmeta.nrows, hbmeta.ncols)
    for j = 1:size(nodemap, 2)
        for i = 1:size(nodemap, 1)
            idx = nodemap[i,j]
            if idx == 0
                continue
            else
                voltmap[i,j] = voltages[idx]
            end
        end
    end
    voltmap
end

alloc_map(hbmeta) = zeros(hbmeta.nrows, hbmeta.ncols)

function accum_voltages!(base, newvolt, nodemap, hbmeta)
    voltmap = _create_voltage_map(newvolt, nodemap, hbmeta)
    for i in eachindex(base)
        base[i] += voltmap[i]
    end
end

function accum_currents!(base, newcurr, cfg, G, voltages, finitegrounds, nodemap, hbmeta)
    node_currents, branch_currents = _create_current_maps(G, voltages, finitegrounds, cfg,
                                                            nodemap = nodemap, hbmeta = hbmeta)

    for i in eachindex(base)
        base[i] += node_currents[i]
    end
end

function save_resistances(r, cfg)
    pref = split(cfg["output_file"], ".out")[1]
    filename = "$(pref)_resistances.out"
    filename_3col = "$(pref)_resistances_3columns.out"
    rcol = compute_3col(r)
    open(filename, "w") do f
        writedlm(f, r, ' ')
    end
    open(filename_3col, "w") do f
        writedlm(f, rcol, ' ')
    end
end

function write_cum_maps(cum, cellmap::Matrix{T}, cfg, hbmeta, write_max, write_cum) where T

    if write_cum || cfg["write_cur_maps"] in TRUELIST
        cum_curr = zeros(T, size(cellmap)...)
        for i = 1:nprocs()
            cum_curr .+= cum.cum_curr[i]
        end
        postprocess_cum_curmap!(cum_curr)
        write_grid(cum_curr, "", cfg, hbmeta, cum = true)
    end

    if write_max
        max_curr = fill(T(-9999), size(cellmap)...)
        for i = 1:nprocs()
            max_curr .= max.(cum.max_curr[i], max_curr)
        end
        postprocess_cum_curmap!(max_curr)
        write_grid(max_curr, "", cfg, hbmeta, max = true)
    end

end


text:
export iauXy06
"""
X,Y coordinates of celestial intermediate pole from series based
on IAU 2006 precession and IAU 2000A nutation.

This function is part of the International Astronomical Union's
SOFA (Standards Of Fundamental Astronomy) software collection.

Status:  canonical model.

Given:
   date1,date2  double     TT as a 2-part Julian Date (Note 1)

Returned:
   x,y          double     CIP X,Y coordinates (Note 2)

Notes:

   1. The TT date date1+date2 is a Julian Date, apportioned in any
   convenient way between the two arguments.  For example,
   JD(TT)=2450123.7 could be expressed in any of these ways,
   among others:

         date1          date2

      2450123.7           0.0       (JD method)
      2451545.0       -1421.3       (J2000 method)
      2400000.5       50123.2       (MJD method)
      2450123.5           0.2       (date & time method)

   The JD method is the most natural and convenient to use in
   cases where the loss of several decimal digits of resolution
   is acceptable.  The J2000 method is best matched to the way
   the argument is handled internally and will deliver the
   optimum resolution.  The MJD method and the date & time methods
   are both good compromises between resolution and convenience.

   2. The X,Y coordinates are those of the unit vector towards the
   celestial intermediate pole.  They represent the combined effects
   of frame bias, precession and nutation.

   3. The fundamental arguments used are as adopted in IERS Conventions
   (2003) and are from Simon et al. (1994) and Souchay et al.
   (1999).

   4. This is an alternative to the angles-based method, via the SOFA
   function iauFw2xy and as used in iauXys06a for example.  The two
   methods agree at the 1 microarcsecond level (at present), a
   negligible amount compared with the intrinsic accuracy of the
   models.  However, it would be unwise to mix the two methods
   (angles-based and series-based) in a single application.

Called:
   iauFal03     mean anomaly of the Moon
   iauFalp03    mean anomaly of the Sun
   iauFaf03     mean argument of the latitude of the Moon
   iauFad03     mean elongation of the Moon from the Sun
   iauFaom03    mean longitude of the Moon's ascending node
   iauFame03    mean longitude of Mercury
   iauFave03    mean longitude of Venus
   iauFae03     mean longitude of Earth
   iauFama03    mean longitude of Mars
   iauFaju03    mean longitude of Jupiter
   iauFasa03    mean longitude of Saturn
   iauFaur03    mean longitude of Uranus
   iauFane03    mean longitude of Neptune
   iauFapa03    general accumulated precession in longitude

References:

   Capitaine, N., Wallace, P.T. & Chapront, J., 2003,
   Astron.Astrophys., 412, 567

   Capitaine, N. & Wallace, P.T., 2006, Astron.Astrophys. 450, 855

   McCarthy, D. D., Petit, G. (eds.), 2004, IERS Conventions (2003),
   IERS Technical Note No. 32, BKG

   Simon, J.L., Bretagnon, P., Chapront, J., Chapront-Touze, M.,
   Francou, G. & Laskar, J., Astron.Astrophys., 1994, 282, 663

   Souchay, J., Loysel, B., Kinoshita, H., Folgueira, M., 1999,
   Astron.Astrophys.Supp.Ser. 135, 111

   Wallace, P.T. & Capitaine, N., 2006, Astron.Astrophys. 459, 981

This revision:  2013 August 21

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# void iauXy06(double date1, double date2, double *x, double *y)

function iauXy06(date1::Real, date2::Real)
   ref_x = Ref{Float64}(0.0)
   ref_y = Ref{Float64}(0.0)

   ccall((:iauXy06, libsofa_c), Cvoid,
         (Cdouble, Cdouble, Ref{Cdouble}, Ref{Cdouble}),
         convert(Float64, date1), convert(Float64, date2),
         ref_x, ref_y)

   return ref_x[], ref_y[]
end


text:
module TestGMRES

using IterativeSolvers
using Test
using LinearMaps
using LinearAlgebra
using Random
using SparseArrays

#GMRES
@testset "GMRES" begin

Random.seed!(1234321)
n = 10

@testset "Matrix{$T}" for T in (Float32, Float64, ComplexF32, ComplexF64)
    A = rand(T, n, n)
    b = rand(T, n)
    F = lu(A)
    reltol = √eps(real(T))

    # Test optimality condition: residual should be non-increasing
    x, history = gmres(A, b, log=true, restart=3, maxiter=10, reltol=reltol);
    @test isa(history, ConvergenceHistory)
    @test all(diff(history[:resnorm]) .<= 0.0)

    # Left exact preconditioner
    x, history = gmres(A, b, Pl=F, maxiter=1, restart=1, reltol=reltol, log=true)
    @test history.isconverged
    @test norm(F \ (A * x - b)) / norm(b) ≤ reltol

    # Right exact preconditioner
    x, history = gmres(A, b, Pl=Identity(), Pr=F, maxiter=1, restart=1, reltol=reltol, log=true)
    @test history.isconverged
    @test norm(A * x - b) / norm(b) ≤ reltol
end

@testset "SparseMatrixCSC{$T, $Ti}" for T in (Float64, ComplexF64), Ti in (Int64, Int32)
    A = sprand(T, n, n, 0.5) + I
    b = rand(T, n)
    F = lu(A)
    reltol = √eps(real(T))

    # Test optimality condition: residual should be non-increasing
    x, history = gmres(A, b, log = true, restart = 3, maxiter = 10);
    @test all(diff(history[:resnorm]) .<= 0.0)

    # Left exact preconditioner
    x, history = gmres(A, b, Pl=F, maxiter=1, restart=1, log=true)
    @test history.isconverged
    @test norm(F \ (A * x - b)) / norm(b) ≤ reltol

    # Right exact preconditioner
    x, history = gmres(A, b, Pl = Identity(), Pr=F, maxiter=1, restart=1, log=true)
    @test history.isconverged
    @test norm(A * x - b) / norm(b) ≤ reltol
end

@testset "Linear operator defined as a function" begin
    A = LinearMap(cumsum!, 100; ismutating=true)
    b = rand(100)
    reltol = 1e-5

    x = gmres(A, b; reltol=reltol, maxiter=2000)
    @test norm(A * x - b) / norm(b) ≤ reltol
end

@testset "Off-diagonal in hessenberg matrix exactly zero" begin
    A = Matrix(1.0I, 2, 2)
    b = [1.0, 2.2]
    x = gmres(A, b)
    @test all(x .== b)
end

@testset "Termination criterion" begin
    for T in (Float32, Float64, ComplexF32, ComplexF64)
        A = T[ 2 -1  0
              -1  2 -1
               0 -1  2]
        n = size(A, 2)
        b = ones(T, n)
        x0 = A \ b
        perturbation = 10 * sqrt(eps(real(T))) * T[(-1)^i for i in 1:n]

        # If the initial residual is small and a small relative tolerance is used,
        # many iterations are necessary
        x = x0 + perturbation
        initial_residual = norm(A * x - b)
        x, ch = gmres!(x, A, b, log=true)
        @test 2 ≤ niters(ch) ≤ n

        # If the initial residual is small and a large absolute tolerance is used,
        # no iterations are necessary
        x = x0 + perturbation
        initial_residual = norm(A * x - b)
        x, ch = gmres!(x, A, b, abstol=2*initial_residual, reltol=zero(real(T)), log=true)
        @test niters(ch) == 0
    end
end
end

end # module


text:
#=

# Compile SPIR-V shaders

Vulkan requires shaders to be in SPIR-V form. We will showcase how to compile shaders into SPIR-V using `glslang` from Julia.

Let's first define a GLSL shader:

=#

shader_code = """
#version 450

layout(location = 0) in vec2 pos;
layout(location = 0) out vec4 color;

void main() {
    color = vec4(pos, 0.0, 1.0);
}
""";

# We will use `glslang` to compile this code into a SPIR-V binary file. First, get the path to the binary.

using glslang_jll: glslangValidator
glslang = glslangValidator(identity)

# Fill in an `IOBuffer` with the shader code (to be used as standard input).

_stdin = IOBuffer()
write(_stdin, shader_code)
seekstart(_stdin)
flags = ["--stdin"]

# Specify which kind of shader we are compiling. Here, we have Vulkan-flavored GLSL, and we will indicate that it is a vertex shader. If your shader is written in HLSL instead of GLSL, you should also provide an additional `"-D"` flag.

push!(flags, "-V", "-S", "vert");

# Use a temporary file as output, as `glslang` does not yet support outputting code to stdout.

path = tempname()
push!(flags, "-o", path)

# Run `glslang`.

run(pipeline(`$glslang $flags`, stdin = _stdin))

# Read the code and clean the temporary file.

spirv_code = reinterpret(UInt32, read(path))
rm(path)
@assert first(spirv_code) == 0x07230203 # SPIR-V magic number
spirv_code

max_stars_repo_name: JuliaGPU/VulkanAbstraction

max_issues_repo_name: JuliaGPU/VulkanAbstraction

max_forks_repo_name: JuliaGPU/VulkanAbstraction


text:
using ZMQ

context = ZMQ.Context()
server = ZMQ.Socket(context, ZMQ.REP)
ZMQ.bind(server, "tcp://127.0.0.1:5555")

cycles = 0
while true
    global cycles

    @info "Waiting for request... (cycles: $(cycles))"
    request = ZMQ.recv(server, Int)
    if cycles > 3 && rand(0:3) == 0
        @info "Simulating a crash"
        break
    elseif cycles > 3 && rand(0:3) == 0
        @info "Simulating CPU overload"
        sleep(2)
    end

    @info "Normal request $(request)"
    sleep(1)
    ZMQ.send(server, request)

    cycles += 1
end

max_stars_repo_name: jwdevantier/julia-zmqtest

max_issues_repo_name: jwdevantier/julia-zmqtest

max_forks_repo_name: jwdevantier/julia-zmqtest


text:
# Guard aginst github #62
using JDF
using Dates
col = [Date(1999,1,1), missing]
df = (d = col, a = [1, missing])
JDF.save("plsdel-date-w-missing.jdf", df)
DataFrame(JDF.load("plsdel-date-w-missing.jdf"))
rm("plsdel-date-w-missing.jdf", recursive=true)


text:
using VoronoiDelaunay
import VoronoiDelaunay: _pushunfixed!, _flipa!, _flipb!, _flipc!
import GeometricalPredicates
import GeometricalPredicates: incircle, intriangle
using Test

@testset "VoronoiDelaunay tests" begin
    @testset begin
        tess = DelaunayTessellation2D(100)
        @test findindex(tess, Point2D(1.1, 1.1)) == 2
        @test findindex(tess, Point2D(1.9, 1.9)) == 3
        @test findindex(tess, Point2D(1.9, 1.9)) == 3
    end

    pa = Point2D(GeometricalPredicates.min_coord, GeometricalPredicates.min_coord)
    pb = Point2D(GeometricalPredicates.min_coord, GeometricalPredicates.max_coord)
    pc = Point2D(GeometricalPredicates.max_coord, GeometricalPredicates.min_coord)
    pd = Point2D(GeometricalPredicates.max_coord, GeometricalPredicates.max_coord)
    pp = Point2D(1.1,1.1)

    @testset begin
        tess = DelaunayTessellation2D(100)
        @test findindex(tess, Point2D(1.9, 1.9)) == 3
        @test findindex(tess, Point2D(1.1, 1.1)) == 2
        @test findindex(tess, Point2D(1.1, 1.1)) == 2

        @test tess._trigs[2]._neighbour_a == 3

        _pushunfixed!(tess, Point2D(1.1, 1.1))

        @test tess._trigs[2]._neighbour_a == 3
        @test tess._trigs[2]._neighbour_b == 4
        @test tess._trigs[2]._neighbour_c == 5

        @test tess._trigs[3]._neighbour_a == 2
        @test tess._trigs[3]._neighbour_b == 1
        @test tess._trigs[3]._neighbour_c == 1

        @test tess._trigs[4]._neighbour_a == 2
        @test tess._trigs[4]._neighbour_b == 1
        @test tess._trigs[4]._neighbour_c == 5

        @test tess._trigs[5]._neighbour_a == 2
        @test tess._trigs[5]._neighbour_b == 4
        @test tess._trigs[5]._neighbour_c == 1

        @test getc(tess._trigs[5]) == pp
        @test geta(tess._trigs[2]) == pp
        @test getb(tess._trigs[4]) == pp

        @test getb(tess._trigs[5]) == pb
        @test getb(tess._trigs[2]) == pb
        @test getc(tess._trigs[3]) == pb

        @test getc(tess._trigs[4]) == pc
        @test getc(tess._trigs[2]) == pc
        @test getb(tess._trigs[3]) == pc

        @test geta(tess._trigs[5]) == pa
        @test geta(tess._trigs[4]) == pa

        @test geta(tess._trigs[3]) == pd

        @test findindex(tess, Point2D(1.01, 1.1)) == 5
        @test findindex(tess, Point2D(1.1, 1.01)) == 4

        @test findindex(tess, Point2D(1.11, 1.11)) == 2
        @test findindex(tess, Point2D(1.6, 1.6)) == 3
        @test findindex(tess, Point2D(1.11, 1.1101)) == 2
        @test findindex(tess, Point2D(1.6, 1.601)) == 3
        @test findindex(tess, Point2D(1.11, 1.11)) == 2
        @test findindex(tess, Point2D(1.6, 1.6)) == 3

        p2 = Point2D(1.9,1.9)
        _pushunfixed!(tess, p2)

        @test geta(tess._trigs[7]) == pd
        @test getb(tess._trigs[7]) == pc
        @test getc(tess._trigs[7]) == p2

        @test tess._trigs[7]._neighbour_a == 3
        @test tess._trigs[7]._neighbour_b == 6
        @test tess._trigs[7]._neighbour_c == 1

        @test geta(tess._trigs[6]) == pd
        @test getb(tess._trigs[6]) == p2
        @test getc(tess._trigs[6]) == pb

        @test tess._trigs[6]._neighbour_a == 3
        @test tess._trigs[6]._neighbour_b == 1
        @test tess._trigs[6]._neighbour_c == 7

        @test geta(tess._trigs[3]) == p2
        @test getb(tess._trigs[3]) == pc
        @test getc(tess._trigs[3]) == pb

        @test tess._trigs[3]._neighbour_a == 2
        @test tess._trigs[3]._neighbour_b == 6
        @test tess._trigs[3]._neighbour_c == 7

        @test geta(tess._trigs[5]) == pa
        @test getb(tess._trigs[5]) == pb
        @test getc(tess._trigs[5]) == pp

        @test tess._trigs[5]._neighbour_a == 2
        @test tess._trigs[5]._neighbour_b == 4
        @test tess._trigs[5]._neighbour_c == 1

        @test geta(tess._trigs[4]) == pa
        @test getb(tess._trigs[4]) == pp
        @test getc(tess._trigs[4]) == pc

        @test tess._trigs[4]._neighbour_a == 2
        @test tess._trigs[4]._neighbour_b == 1
        @test tess._trigs[4]._neighbour_c == 5

        @test geta(tess._trigs[2]) == pp
        @test getb(tess._trigs[2]) == pb
        @test getc(tess._trigs[2]) == pc

        @test tess._trigs[2]._neighbour_a == 3
        @test tess._trigs[2]._neighbour_b == 4
        @test tess._trigs[2]._neighbour_c == 5

        @test tess._last_trig_index == 7

        _flipa!(tess, 2, 3)

        @test geta(tess._trigs[2]) == pp
        @test getb(tess._trigs[2]) == p2
        @test getc(tess._trigs[2]) == pc

        @test tess._trigs[2]._neighbour_a == 7
        @test tess._trigs[2]._neighbour_b == 4
        @test tess._trigs[2]._neighbour_c == 3

        @test geta(tess._trigs[3]) == pp
        @test getb(tess._trigs[3]) == pb
        @test getc(tess._trigs[3]) == p2

        @test tess._trigs[3]._neighbour_a == 6
        @test tess._trigs[3]._neighbour_b == 2
        @test tess._trigs[3]._neighbour_c == 5

        @test geta(tess._trigs[5]) == pa
        @test getb(tess._trigs[5]) == pb
        @test getc(tess._trigs[5]) == pp

        @test geta(tess._trigs[6]) == pd
        @test getb(tess._trigs[6]) == p2
        @test getc(tess._trigs[6]) == pb

        @test geta(tess._trigs[4]) == pa
        @test getb(tess._trigs[4]) == pp
        @test getc(tess._trigs[4]) == pc

        @test geta(tess._trigs[7]) == pd
        @test getb(tess._trigs[7]) == pc
        @test getc(tess._trigs[7]) == p2

        @test tess._trigs[4]._neighbour_a == 2
        @test tess._trigs[4]._neighbour_b == 1
        @test tess._trigs[4]._neighbour_c == 5

        @test tess._trigs[5]._neighbour_a == 3
        @test tess._trigs[5]._neighbour_b == 4
        @test tess._trigs[5]._neighbour_c == 1

        @test tess._trigs[6]._neighbour_a == 3
        @test tess._trigs[6]._neighbour_b == 1
        @test tess._trigs[6]._neighbour_c == 7

        @test tess._trigs[7]._neighbour_a == 2
        @test tess._trigs[7]._neighbour_b == 6
        @test tess._trigs[7]._neighbour_c == 1

        _flipb!(tess, 3, 2)

        @test geta(tess._trigs[3]) == pp
        @test getb(tess._trigs[3]) == pb
        @test getc(tess._trigs[3]) == pc

        @test geta(tess._trigs[2]) == pc
        @test getb(tess._trigs[2]) == pb
        @test getc(tess._trigs[2]) == p2

        @test tess._trigs[3]._neighbour_a == 2
        @test tess._trigs[3]._neighbour_b == 4
        @test tess._trigs[3]._neighbour_c == 5

        @test tess._trigs[2]._neighbour_a == 6
        @test tess._trigs[2]._neighbour_b == 7
        @test tess._trigs[2]._neighbour_c == 3

        @test tess._trigs[6]._neighbour_a == 2
        @test tess._trigs[7]._neighbour_a == 2
        @test tess._trigs[5]._neighbour_a == 3
        @test tess._trigs[4]._neighbour_a == 3

        _flipc!(tess,2,3)

        @test geta(tess._trigs[2]) == pp
        @test getb(tess._trigs[2]) == pb
        @test getc(tess._trigs[2]) == p2

        @test geta(tess._trigs[3]) == pc
        @test getb(tess._trigs[3]) == pp
        @test getc(tess._trigs[3]) == p2

        @test tess._trigs[2]._neighbour_a == 6
        @test tess._trigs[2]._neighbour_b == 3
        @test tess._trigs[2]._neighbour_c == 5

        @test tess._trigs[3]._neighbour_a == 2
        @test tess._trigs[3]._neighbour_b == 7
        @test tess._trigs[3]._neighbour_c == 4

        @test tess._trigs[6]._neighbour_a == 2
        @test tess._trigs[7]._neighbour_a == 3
        @test tess._trigs[5]._neighbour_a == 2
        @test tess._trigs[4]._neighbour_a == 3
    end

    @testset begin
        tess = DelaunayTessellation2D(100)
        pp = Point2D(1.45, 1.49)
        _pushunfixed!(tess, pp)
        _flipa!(tess,2,3)

        @test geta(tess._trigs[2]) == pp
        @test getb(tess._trigs[2]) == pd
        @test getc(tess._trigs[2]) == pc

        @test geta(tess._trigs[3]) == pp
        @test getb(tess._trigs[3]) == pb
        @test getc(tess._trigs[3]) == pd
    end

    @testset begin
        point_arr = Point2D[]
        n=1000
        tess = DelaunayTessellation2D(n*10)
        for i in 1:n
            push!(point_arr, Point2D(rand()+1.0, rand()+1.0))
        end
        push!(tess, point_arr)
        @test tess._last_trig_index == n*2+3

        for p in point_arr
            for t in tess._trigs[2:tess._last_trig_index]
                i = incircle(t, p)
                if i > 0 && ((p == geta(t)) || (p == getb(t)) || (p == getc(t)))
                    i = 0
                end
                @test i <= 0
            end
        end
    end

    @testset begin
        point_arr = Point2D[]
        n=10
        tess = DelaunayTessellation2D(n*n*10)
        for x in range(1.001,stop=1.999,length=n)
            for y in range(1.001,stop=1.999,length=n)
                push!(point_arr, Point2D(x,y))
            end
        end
        push!(tess, point_arr)
        @test tess._last_trig_index == n*n*2+3

        for p in point_arr
            for t in tess._trigs[2:tess._last_trig_index]
                i = incircle(t, p)
                if i > 0 && ((p == geta(t)) || (p == getb(t)) || (p == getc(t)))
                    i = 0
                end
                @test i <= 0
            end
        end
    end

    # Iterator test
    @testset begin
    point_arr = Point2D[]
    n=1000
    tess = DelaunayTessellation2D(n*10)
    for i in 1:n
        push!(point_arr, Point2D(rand()+1.0, rand()+1.0))
    end
    push!(tess, point_arr)
    p = Point2D(rand()+1.0, rand()+1.0)
    counter = 0
    for t in tess
        if intriangle(t, p) == 1
            counter += 1
        end
    end
    @test counter == 1 # p can be contained only in one triangle
end

end
# that's it for today!


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

using Markdown
using InteractiveUtils

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

# ╔═╡ 073afcbe-4cbf-11eb-33a8-a327232257eb
begin
	import Pkg
	Pkg.activate(mktempdir())
	Pkg.add([
			"CSV",
			"DataFrames",
			"Plots",
			"PlutoUI",
	])
	using CSV
	using DataFrames
	using Dates
	using InteractiveUtils
	using Markdown
	using Statistics
	using Plots
	using PlutoUI
end

# ╔═╡ 1b720d94-4f8f-11eb-0baf-bbcd93231972
md"""
# United States COVID-19 Visualization
This is a visualization of weekly reported COVID-19 cases against total reported cases. The plot is done on a log-scale causing exponential growth to be represented linearly. Since time is not explicitly represented on the graph, direct comparisons can be drawn between states who experienced similar case trends at different times. The visualization also provides an option to normalize for population data to give a better comparison for states with vastly different populations.

This project was created by Jeremiah DeGreeff with inspiration from Problem Set 1 for [MIT Course 6.S083 Spring 2020](https://www.eecs.mit.edu/academics-admissions/academic-information/subject-updates-spring-2020/6s083) and [this visualization](https://aatishb.com/covidtrends/) created by Aatish Bhatia.
"""

# ╔═╡ 7ec0992a-4f8e-11eb-3500-53545a706222
md"""
## Data Preparation
"""

# ╔═╡ 2bd8a378-4f8f-11eb-3e67-ed3ed64e820c
begin
	data_url = "https://raw.githubusercontent.com/nytimes/covid-19-data/master/us-states.csv"
	md"""
	US State COVID-19 Data Source: [New York Times](https://github.com/nytimes/covid-19-data)
	"""
end

# ╔═╡ fd95c222-4cbf-11eb-0690-57cc206d122a
begin
	covid_data_filename = "data/covid_state_data.csv"
	download(data_url, covid_data_filename)
	covid_data = CSV.read(covid_data_filename, DataFrame)
	indexed_states = collect(covid_data[:, :state])
	unique_states = sort(unique(indexed_states))
	indexed_dates = collect(covid_data[:, :date])
	unique_dates = sort(unique(indexed_dates))
	covid_data
end

# ╔═╡ 4f0a03d2-4f8f-11eb-3f36-edbec0b6a9c3
md"""
US State Population Data Source: [United States Census Bureau](https://www.census.gov/programs-surveys/popest/technical-documentation/research/evaluation-estimates.html)
"""

# ╔═╡ aadc3bac-4cc5-11eb-298c-1182973f7011
begin
	population_data_filename = "data/state_population_data_2020.csv"
	population_data = CSV.read(population_data_filename, DataFrame)
	population_states = collect(population_data[:, :state])
	normalized_covid_data = covid_data[[s ∈ population_states for s in indexed_states], :]
	transform!(normalized_covid_data, [:state, :cases] => ((states, cases) -> 1e5 * cases ./ (s -> population_data[population_states.==s, :population][1]).(states)) => :cases)
	normalized_indexed_states = collect(normalized_covid_data[:, :state])
	normalized_indexed_dates = collect(normalized_covid_data[:, :date])
	population_data
end

# ╔═╡ 9e089168-4cc1-11eb-023a-8d7cfd681e88
# Function to extract case count for a state on a day
function get_case_count(state, date, normalize=false)
	if normalize
		@assert state in population_states
    	cases = normalized_covid_data[(normalized_indexed_states.==state).*(normalized_indexed_dates.==date), :cases]
	else
		@assert state in unique_states
    	cases = covid_data[(indexed_states.==state).*(indexed_dates.==date), :cases]
	end
	@assert length(cases) < 2
    return length(cases) == 1 ? cases[1] : 0
end;

# ╔═╡ 70496912-4f8e-11eb-28ca-c50ae0af115e
md"""
## Visualization
"""

# ╔═╡ b9c78b8a-4edf-11eb-0a36-c72d60e4680a
begin
	normalize_checkbox = @bind normalize CheckBox(default=true)
	md"""
	Normalize by population $(normalize_checkbox)
	"""
end

# ╔═╡ dd4ae0b4-4cc2-11eb-0dad-7de8cb6c055b
begin
	date_slider = @bind date_index Slider(1:length(unique_dates), default=1, show_value=true)
	md"""
	Days since $(unique_dates[1]): $(date_slider)
	"""
end

# ╔═╡ e7871fa2-4f9c-11eb-0d97-f7040e7c285c
begin
	double_rate_slider = @bind double_rate Slider(1:100, default=7, show_value=true)
	md"""
	Exponential growth doubling rate (in days): $(double_rate_slider)
	"""
end

# ╔═╡ 68593592-4f9e-11eb-2fac-bb30829cb8f0
begin
	legend_checkbox = @bind legend CheckBox(default=true)
	md"""
	Show legend $(legend_checkbox)
	"""
end

# ╔═╡ 3fbcd01c-4cc2-11eb-2613-190f71f8b3c7
md"""
## Plotting Library
"""

# ╔═╡ 513019c0-4f9a-11eb-1086-b795d60a5b50
# Function to apply a function over a dictionary of data sequences
apply_over_dict(func, dict, args=()) = Dict(k => isempty(args) ? func(v) : func(v, args) for (k,v) ∈ dict)

# ╔═╡ 54673912-4f99-11eb-36ca-9b0124ed444b
# Function to convert aggregate data to weekly data
aggreggate_to_weekly_change(aggregate_data) = [aggregate_data[i] -  aggregate_data[max(1, i - 6)] for i ∈ 1:length(aggregate_data)]

# ╔═╡ 76fb1292-4f95-11eb-1824-4d5437b91463
# Function to prepare version of data for a log range
truncate(data, min_value = 0) = map(x -> x > min_value ? x : NaN, data)

# ╔═╡ 6bb8b566-4cc2-11eb-330b-13eae1cc5858
# Function to determine upper bound for a log range
max_value(data) = 10^ceil(Int64, log(10, maximum([d for (_, s) ∈ data for d ∈ s if !isnan(d)])))

# ╔═╡ 7a669466-4cc2-11eb-1359-f3f2d060452b
# Function to plot weekly change vs aggregate on a log scale
function plot_covid(aggregate_data, dates, date_index; min_value = 10, double_rate = missing, normalize = false)
	p = plot(xscale=:log10, yscale=:log10, leg=:bottomright)
	title!("Trajectory of COVID-19 Confirmed Cases ($(dates[date_index]))")
	xlabel!(normalize ? "Total Confirmed Cases per 100,000 Residents" : "Total Confirmed Cases")
	ylabel!(normalize ? "Recent Confirmed Cases per 100,000 Residents" : "Recent Confirmed Cases")

    # Weekly change data is the number of new cases in the past week
    weekly_change_data = apply_over_dict(aggreggate_to_weekly_change, aggregate_data)

	aggregate_data = apply_over_dict(truncate, aggregate_data, min_value)
    weekly_change_data = apply_over_dict(truncate, weekly_change_data)

	min_x = 10^floor(Int64, log(10, min_value))
	if all(s -> all(isnan, s), values(aggregate_data))
		x_range = (min_x, min_x * 10)
		y_range = (1, 10)
	else
		x_range = (min_x, max_value(aggregate_data))
		y_range = (1, max_value(weekly_change_data))
	end
	xlims!(x_range)
	ylims!(y_range)

	for (i, c) in enumerate(keys(aggregate_data))
		plot!(aggregate_data[c][1:date_index], weekly_change_data[c][1:date_index], color=i, label=false)
		if !isnan(aggregate_data[c][date_index]) && !isnan(weekly_change_data[c][date_index])
			scatter!(aggregate_data[c][date_index:date_index], weekly_change_data[c][date_index:date_index], color=i, label=c, markerstrokewidth=0)
		end
	end

	if !ismissing(double_rate)
		# Exponential data provides a reference for what pure exponential growth with a particular growth rate would look like
		exponential_aggregate = [2^(t / double_rate) for t ∈ 0:ceil(log2(x_range[2]) * double_rate)]
		exponential_weekly_change = aggreggate_to_weekly_change(exponential_aggregate)

		exponential_aggregate = truncate(exponential_aggregate)
		exponential_weekly_change = truncate(exponential_weekly_change)

		plot!(exponential_aggregate, exponential_weekly_change, color=:gray, lw=2, label="Exponential Growth (d=$(double_rate))")
	end
	p
end

# ╔═╡ 2180a1f4-4f8b-11eb-073b-4f671634dc37
md"""
## PlutoUI Library Overrides
"""

# ╔═╡ 318885e6-4f8b-11eb-2b78-49f08e7074fc
begin
	import Markdown: htmlesc, withtag

    struct MultiSelect
        options::Array{Pair{<:AbstractString,<:Any},1}
        default::Union{Missing, AbstractVector{AbstractString}}
        size::Int
    end

	MultiSelect(options::Array{<:AbstractString,1}; default=missing, size=5) = MultiSelect([o => o for o in options], default, size)
    MultiSelect(options::Array{<:Pair{<:AbstractString,<:Any},1}; default=missing, size=5) = MultiSelect(options, default, size)

	function Base.show(io::IO, ::MIME"text/html", select::MultiSelect)
        withtag(io, Symbol("select multiple"), :size => string(select.size)) do
            for o in select.options
                print(io, """<option value="$(htmlesc(o.first))"$(!ismissing(select.default) && o.first ∈ select.default ? " selected" : "")>""")
                if showable(MIME"text/html"(), o.second)
                    show(io, MIME"text/html"(), o.second)
                else
                    print(io, o.second)
                end
                print(io, "</option>")
            end
        end
    end

	Base.get(select::MultiSelect) = ismissing(select.default) ? Any[] : select.default
end

# ╔═╡ b7aa216e-4f8d-11eb-339c-394d2c23ee4f
begin
	state_options = normalize ? population_states : unique_states
	select = @bind states MultiSelect(state_options, size=length(state_options))
    md"""
    Select the states to display.

    $(select)
    """
end

# ╔═╡ c60b7e00-4cc1-11eb-2e78-c1a55563de15
begin
	if !isempty(states)
		aggregate_data = Dict(s => [get_case_count(s, d, normalize) for d ∈ unique_dates] for s ∈ states)
		p = plot_covid(aggregate_data, unique_dates, date_index; double_rate=double_rate, normalize=normalize)
		if !legend plot!(p, leg=false) end
		p |> as_svg
	end
end

# ╔═╡ Cell order:
# ╟─073afcbe-4cbf-11eb-33a8-a327232257eb
# ╟─1b720d94-4f8f-11eb-0baf-bbcd93231972
# ╟─7ec0992a-4f8e-11eb-3500-53545a706222
# ╟─2bd8a378-4f8f-11eb-3e67-ed3ed64e820c
# ╟─fd95c222-4cbf-11eb-0690-57cc206d122a
# ╟─4f0a03d2-4f8f-11eb-3f36-edbec0b6a9c3
# ╟─aadc3bac-4cc5-11eb-298c-1182973f7011
# ╟─9e089168-4cc1-11eb-023a-8d7cfd681e88
# ╟─70496912-4f8e-11eb-28ca-c50ae0af115e
# ╟─b9c78b8a-4edf-11eb-0a36-c72d60e4680a
# ╟─b7aa216e-4f8d-11eb-339c-394d2c23ee4f
# ╟─dd4ae0b4-4cc2-11eb-0dad-7de8cb6c055b
# ╟─e7871fa2-4f9c-11eb-0d97-f7040e7c285c
# ╟─68593592-4f9e-11eb-2fac-bb30829cb8f0
# ╟─c60b7e00-4cc1-11eb-2e78-c1a55563de15
# ╟─3fbcd01c-4cc2-11eb-2613-190f71f8b3c7
# ╟─513019c0-4f9a-11eb-1086-b795d60a5b50
# ╟─54673912-4f99-11eb-36ca-9b0124ed444b
# ╟─76fb1292-4f95-11eb-1824-4d5437b91463
# ╟─6bb8b566-4cc2-11eb-330b-13eae1cc5858
# ╟─7a669466-4cc2-11eb-1359-f3f2d060452b
# ╟─2180a1f4-4f8b-11eb-073b-4f671634dc37
# ╟─318885e6-4f8b-11eb-2b78-49f08e7074fc

max_stars_repo_name: jrdegreeff/COVID-Visualization

max_issues_repo_name: jrdegreeff/COVID-Visualization

max_forks_repo_name: jrdegreeff/COVID-Visualization


text:
import Base: show, dump

function show(io::IO, intrng::Array{Intrng,1})
    write(io, """A 1.0 #regions in PPM
# low field   high field  bias        slope\n""")
    entries = sort(intrng, by=(Δ->Δ[1]), rev=true)
    for (i,(δ1, δ2)) in enumerate(entries)
        write(io, " $δ1  $δ2  0.0  0.0  # for region $i\n")
    end
end

function dump(path, s::Spectrum)
    if !isdir(path)
        mkdir(path)
    end

    # Save each child proc
    for (n,p) in s.procs
        dump(joinpath(path, "pdata", string(n)), p)
    end

    # Save fid
    write(open(joinpath(path,"fid"), "w"), Int32.(round.(s.fid)))

    # TODO: Save acqu/acqus
end

function dump(path, s::Spectrum, templatepath)
    cp(templatepath, path)
    dump(path, s)
end

function dump(path, p::ProcessedSpectrum)
    if !isdir(path)
        mkdir(path)
    end

    # Save 1r and 1i
    open(joinpath(path, "1r"), "w") do f
        write(f, Int32.(round.(p.re_ft)))
    end
    open(joinpath(path, "1i"), "w") do f
        write(f, Int32.(round.(p.im_ft)))
    end

    # Save intrng
    open(joinpath(path, "intrng"), "w") do f
        show(f, p.intrng)
    end

    # Save title
    open(joinpath(path, "title"), "w") do f
        write(f, p.title)
    end

    # TODO: Save proc/procs
end


text:
using Test
using HyperbolicDrawSimpleGraphs

@test 1==1


text:
module CombinatorTraceTypes

include("../src/Jaynes.jl")
using .Jaynes
using Gen

model1 = @jaynes function foo(i::Int, x::Float64)::Float64
    y ~ Normal(x, 1.0)
    z ~ Normal(y, 3.0)
    z
end (DefaultPipeline)
display(model1)
uc = Gen.Unfold(model1)

model2 = @jaynes (i::Int, x::Float64) -> begin
    y ~ uc(i, x)
end (DefaultPipeline)
display(model2)

tr = simulate(model2, (10, 5.0))

end # module


text:
"""
    TrajectoryOptimization
Primary module for setting up and solving trajectory optimization problems.
"""
module TrajectoryOptimization

using ForwardDiff
using DocStringExtensions
using Interpolations
using RecipesBase
using LinearAlgebra
using Statistics
using Random
using SparseArrays
using StaticArrays
using Logging
using Formatting
using BenchmarkTools
using Parameters
using Rotations
using MathOptInterface
using Quaternions
using UnsafeArrays

const MOI = MathOptInterface
const MAX_ELEM = 170

import Base.copy

export
    Dynamics,
    Problems,
    Controllers

export
    state_dim,
    control_dim

# Primary types
export
    Model,
    QuadraticCost,
    LQRCost,
    LQRObjective,
    GenericCost,
    Trajectory

export
    Problem,
    iLQRSolver,
    iLQRSolverOptions,
    AugmentedLagrangianSolver,
    AugmentedLagrangianSolverOptions,
    AugmentedLagrangianProblem,
    ALTROSolverOptions,
    DIRCOLSolver,
    DIRCOLSolverOptions,
    ProjectedNewtonSolver,
    ProjectedNewtonSolverOptions,
    SequentialNewtonSolver,
    Discrete,
    Continuous,
    Constraint,
    BoundConstraint,
    Equality,
    Inequality,
    Objective,
    Constraints

export
    rk3,
    rk4,
    midpoint,
    add_constraints!,
    goal_constraint,
    initial_controls!,
    initial_state!,
    circle_constraint,
    sphere_constraint

# Primary methods
export
    solve,
    solve!,
    rollout!,
    rollout,
    forwardpass!,
    backwardpass!,
    cost,
    max_violation,
    max_violation_direct,
    infeasible_control,
    line_trajectory,
    evaluate!,
    jacobian!

export
    get_sizes,
    num_constraints,
    get_num_constraints,
    get_num_controls,
    init_results,
    to_array,
    get_N,
    to_dvecs,
    quat2rot,
    sphere_constraint,
    circle_constraint,
    plot_trajectory!,
    plot_vertical_lines!,
    convergence_rate,
    plot_obstacles,
    evals,
    reset,
    reset_evals,
    final_time,
    total_time,
    count_constraints,
    inequalities,
    equalities,
    bounds,
    labels,
    terminal,
    stage,
    interp_rows

# Static methods
export
    convertProblem

include("utils.jl")
include("rotations.jl")
include("logger.jl")
include("knotpoint.jl")
include("expansions.jl")
include("model.jl")
include("costfunctions.jl")
include("objective.jl")
include("solver_opts.jl")
include("solvers.jl")
include("abstract_constraint.jl")
include("constraints.jl")
include("dynamics_constraints.jl")
include("integration.jl")
include("dynamics.jl")

include("cost.jl")
include("static_methods.jl")
include("constraint_vals.jl")
include("constraint_sets.jl")
include("problem.jl")
# include("solvers/silqr/silqr_solver.jl")
# include("solvers/silqr/silqr_methods.jl")
include("solvers/ilqr/ilqr.jl")
include("solvers/ilqr/ilqr_solve.jl")
include("solvers/ilqr/backwardpass.jl")
include("solvers/ilqr/rollout.jl")
include("solvers/augmented_lagrangian/sal_solver.jl")
include("solvers/augmented_lagrangian/sal_methods.jl")
include("solvers/direct/primals.jl")
include("solvers/direct/pn.jl")
include("solvers/direct/pn_methods.jl")
include("solvers/altro/altro_solver.jl")

include("solvers/direct/moi.jl")
include("solvers/direct/copy_blocks.jl")
include("solvers/direct/direct_constraints.jl")

include("problems.jl")
include("controllers.jl")

write_ipopt_options()
Logging.global_logger(default_logger(true))
end


text:
using Pkg.Artifacts
using vSmartMOM.Absorption
using HTTP

deploy = true

isos = Absorption.iso_info

hitran_url_base    = "https://hitran.org/lbl/api?numin=0&numax=150000&iso_ids_list="
local_download_url = "/net/fluo/data1/ftp/XYZT_hitran/"
target_url = "http://web.gps.caltech.edu/~cfranken/hitran_2016/"
# Loop over all current HITRAN molecules:
for molec=1:49
    if haskey(isos,(molec,1))
        locs = [k for (k,v) in isos if k[1]==molec]

        global_ids = ""
        for i in locs
            global_ids *= ","*string(isos[i][1])
            molec_name = isos[i][5]
        end
        println(global_ids[2:end])
        url = hitran_url_base*global_ids[2:end]

        filename = "hitran_molec_id_" * string(molec)*"_" * isos[molec,1][5] * ".par"
        name = splitext(filename)[1]
        @info("Generating artifact for $(filename)")
        HTTP.open(:GET, url) do http
            open(local_download_url*filename, "w") do file
                write(file, http)
            end
        end
        hash = create_artifact() do artifact_dir
            # Copy in weights
            cp(joinpath(local_download_url, filename), joinpath(artifact_dir, filename))
        end

        # Spit tarballs to be hosted out to local temporary directory:
        if deploy
            tarball_hash = archive_artifact(hash, joinpath(local_download_url, "$(name).tar.gz"))

            # Calculate tarball url
            tarball_url = target_url*name*".tar.gz"
            # Bind this to an Artifacts.toml file
            @info("Binding $(name) in Artifacts.toml...")
            bind_artifact!(joinpath(@__DIR__, "Artifacts.toml"), name, hash; download_info=[(tarball_url, tarball_hash)], lazy=true, force=true)
        end
    end
end


text:
function __init__()
    @require TaylorModels = "314ce334-5f6e-57ae-acf6-00b6e903104a" load_taylormodels()
    @require IntervalMatrices = "5c1f47dc-42dd-5697-8aaa-4d102d140ba9" include("init_IntervalMatrices.jl")
end

function load_taylormodels()
    eval(load_taylormodels_overapproximation())
end