text:
module Mortchartgen

using CSV, DataFrames, DataStructures, JSON, Loess, Mustache, ODBC, Statistics, PyCall, Missings

const bo = PyNULL()
const bp = PyNULL()
const bml  = PyNULL()
const bpal = PyNULL()

function __init__()
	copy!(bo, pyimport("bokeh"))
	copy!(bp, pyimport("bokeh.plotting"))
	copy!(bml, pyimport("bokeh.models"))
	copy!(bpal, pyimport("bokeh.palettes"))
end

mainpath = normpath(@__DIR__, "..")
datapath = normpath(mainpath, "data")
chartpath = normpath(mainpath, "charts")
tmpoutpath = normpath(tempdir(), "mout.html")
mkpath(chartpath)
conf = JSON.parsefile(normpath(datapath, "chartgen.json"),
	dicttype=DataStructures.OrderedDict)
tables = Dict(:deaths => "Deaths", :pop => "Pop")
wpp = CSV.read(normpath(datapath, "wppconverted.csv"))
dfarrmatch(col, arr) = map((x) -> in(x, arr), Vector(col))
ctrycodes = map((x)->parse(Int, x), collect(keys(conf["countries"])))
perc_round(value) = replace("$(round(value, 4))", ".", ",")
dthalias(language) = uppercasefirst(conf["deaths"]["alias"][language])
obend = "svg"

function caalias(cause, language)
	if cause=="pop"
		return conf["pop"]["alias"][language]
	else
		return conf["causes"][cause]["alias"][language]
	end
end

function cgen_frames(causes = keys(conf["causes"]))
	conn_config = conf["settings"]["conn_config"]
	#conn = MySQL.connect(conn_config["host"], conn_config["user"],
	#	conn_config["password"], db = conn_config["database"],
	#	unix_socket = conn_config["unix_socket"])
	conn = ODBC.DSN(conn_config["dsn"], conn_config["user"], conn_config["password"])
	qstr = *("select Sex,Year,Country,Admin1,",
		join(map((x)->"$(tables[:pop])$x,", 1:25)),
		"$(tables[:pop])26 from $(tables[:pop])")
	popframe = ODBC.query(conn, qstr)
	rename!(popframe,
		Dict(zip(map((x)->Symbol("$(tables[:pop])$(x)"), 1:26),
			 map((x)->Symbol("Age$(x)"), 1:26))))
	dthframes = DataStructures.OrderedDict()
	for cause in causes
		causeexpr = conf["causes"][cause]["causeexpr"]
		selstat = *("Sex,Year,List,Country,Admin1,",
			join(map((x)->"$(tables[:deaths])$x,", 1:25)),
			"$(tables[:deaths])26")
		#Prepared statements give incorrect results.
		qstr = *("select ", selstat,
			" from $(tables[:deaths]) where ",
			"(case when List='07A' then Cause REGEXP '$(causeexpr["07A"])' ",
			"when List='08A' then Cause REGEXP '$(causeexpr["08A"])' ",
			"when List REGEXP '09(B|N)' then Cause REGEXP '$(causeexpr["09B"])' ",
			"when List REGEXP '10(M|[3-4])' then Cause REGEXP '$(causeexpr["10"])' ",
			"when List='101' then Cause REGEXP '$(causeexpr["101"])' end) ",
			"order by Year")
		dthframe = aggregate(ODBC.query(conn, qstr),
				[:Sex, :Year, :List, :Country, :Admin1], sum)
		rename!(dthframe,
			Dict(zip(map((x)->Symbol("$(tables[:deaths])$(x)_sum"), 1:26),
				 map((x)->Symbol("Age$(x)"), 1:26))))
		dthframes[cause] = hcat(DataFrame(Cause = fill(cause, size(dthframe)[1])),
			dthframe)
	end
	ODBC.disconnect!(conn)
	frames = Dict(:deaths => vcat(values(dthframes)...), :pop => popframe)
	for tblkey in keys(tables)
		frame = frames[tblkey]
		for r in eachrow(frame)
			if r[:Country]==4100 && r[:Year]<1990
				r[:Country] = 4085
			elseif r[:Country]==3150 && r[:Year]>1974 && !ismissing(r[:Admin1])
				r[:Country] = 0
			end
		end
		delete!(frame, :Admin1)
		ncols = size(frame)[2]
		frame_trimmed = frame[dfarrmatch(frame[:Country], ctrycodes), :]
		frame_trimmed_long = stack(frame_trimmed, ncols-25:ncols)
		frames[tblkey] = frame_trimmed_long
	end
	for ctrycode in ctrycodes
		ctry = conf["countries"]["$ctrycode"]
		if haskey(ctry, "wppstart")
			wppctry = wpp[((wpp[:Country].==ctrycode).&
				(wpp[:Year].>=ctry["wppstart"])), :]
			frames[:pop] = vcat(frames[:pop], wppctry)
		end
	end
	frames
end

function save_frames(framedict)
	for tblkey in keys(tables)
		writetable(normpath(datapath, *(string(tblkey), ".csv")),
			framedict[tblkey])
	end
end

function load_frames()
	framedict = Dict()
	for tblkey in keys(tables)
		framedict[tblkey] = CSV.read(normpath(datapath, *(string(tblkey), ".csv")))
		framedict[tblkey][:variable] = map(Symbol, framedict[tblkey][:variable])
	end
	framedict
end

subframe_sray(df, sex, countries, agelist, years) = df[((df[:Sex].==sex)
	.& (dfarrmatch(df[:Country], countries)) .& (dfarrmatch(df[:variable], agelist))
	.& (dfarrmatch(df[:Year], years))), :]
dfgrp_agemean(df, grpcol, f = mean) = by(df, grpcol, x -> DataFrame(value = f(x[:value])))
dfgrp_sum(df, grpcol, f = sum) = by(df, grpcol,
x -> DataFrame(value = f(x[:value]), value_1 = f(x[:value_1])))

function listchanges(country, years, framedict)
	countryframe = subframe_sray(framedict[:deaths], 2, country, [:Age1], years)
	listframe = countryframe[countryframe[:Cause].=="all", [:Year; :List]]
	nrows = size(listframe)[1]
	lcomp = DataFrame(Year = listframe[:Year][2:nrows],
		List = listframe[:List][2:nrows], Listprev = listframe[:List][1:nrows-1])
	lcomp[lcomp[:List].!=lcomp[:Listprev], :]
end

function grpprop(numframe_sub, denomframe_sub, grpcol, agemean)
	numdenomframe_sub =
		join(numframe_sub, denomframe_sub, on = grpcol, makeunique = true)
	if agemean
		propfr_agesp = DataFrame()
		propfr_agesp[grpcol] = numdenomframe_sub[grpcol]
		propfr_agesp[:value] = numdenomframe_sub[:value]./numdenomframe_sub[:value_1]
		return dfgrp_agemean(propfr_agesp, grpcol)
	else
		numdenomgrp = dfgrp_sum(numdenomframe_sub, grpcol)
		propfr_agegr = DataFrame()
		propfr_agegr[grpcol] = numdenomgrp[grpcol]
		propfr_agegr[:value] = numdenomgrp[:value]./numdenomgrp[:value_1]
		return propfr_agegr
	end
end

function propgrp(numframe, denomframe, sex, countries, agelist, years, agemean, grpcol)
	numframe_sub = subframe_sray(numframe, sex, countries, agelist, years)
	denomframe_sub = subframe_sray(denomframe, sex, countries, agelist, years)
	grpprop(numframe_sub, denomframe_sub, grpcol, agemean)
end

function ageslice(sage, eage, agemean, language)
	infsymb = "\u03c9"
	ages = map((x)->Symbol("Age$x"), 1:25)
	agest = [0;0:4;5:5:90;95]
	ageend = [infsymb;0:4;9:5:94;infsymb]
	if agemean
		agemeanstr = " $(conf["agemean"]["alias"][language])"
	else
		agemeanstr = ""
	end
	Dict(:agelist => ages[sage:eage],
		:alias => "$(agest[sage])\u2013$(ageend[eage])$agemeanstr",
		:agest => agest[sage], :ageend => ageend[eage],
		:color => bpal.Category20[20][mod1(sage, 20)])
end

function propplotframes(ca1, ca2, framedict, language)
	deaths = framedict[:deaths]
	ca1frame = deaths[deaths[:Cause].==ca1, :]
	ca1alias = caalias(ca1, language)
	ca2alias = caalias(ca2, language)
	if ca2=="pop"
		ca2frame = framedict[:pop]
	else
		ca2frame = deaths[deaths[:Cause].==ca2, :]
	end
	Dict(:ca1frame => ca1frame, :ca1alias => ca1alias,
		:ca2frame => ca2frame, :ca2alias => ca2alias)
end

function listlabels(country, years, minvals, framedict)
	lch = listchanges(country, years, framedict)
	lchdata = bml[:ColumnDataSource](data = Dict("year" => lch[:Year],
		"list" => lch[:List]))
	listlabels = bml[:LabelSet](x = "year", y = minimum(minvals),
		text = "list", text_color = "red", angle = pi/2, render_mode = "canvas",
		source = lchdata)
end

function preplot(outfile)
	bp[:reset_output]()
	bp[:output_file](outfile)
end

function postplot(figure, showplot)
	if showplot
		bp[:show](figure)
	else
		bp[:save](figure)
	end
end

function propplot_sexesyrs(ca1, ca2, sexes, country, sage, eage, years, agemean,
	framedict, language, outfile, showplot, apctype = "pa", acbcoh = 1900)
	ages = ageslice(sage, eage, agemean, language)
	if (apctype == "pa" || apctype == "pc")
		xdict = (sage, eage)
		xcol = :Year
		ylab = ages[:alias]
		ages = [(sage, eage)]
	elseif (apctype == "ap" || apctype == "ac")
		if apctype == "ap"
			xdict = years
			ylab = string(years)
		else
			xdict = acbcoh
			ylab = string(acbcoh)
		end
		xcol = :Age
		ages = collect(map((a)->(a,a), sage:eage))
	end
	meta = Dict(:ca1 => ca1, :ca2 => ca2, :country => country,
		:sage => sage, :eage => eage, :agemean => agemean,
		:apctype => apctype, :years => years, :acbcoh => acbcoh,
		:xcol => xcol, :ylab => ylab, :plclass => "sexesyrs")
	propframes = Dict()
	for sex in sexes
		propframe = propplot_agesyrs(ca1, ca2, sex, country, ages,
			years, agemean, framedict, language, "linear", "",
			false, apctype)[xdict][:frame]
		yrfloatarr = convert(Array{Float64}, propframe[xcol])
		valarr = convert(Array, propframe[:value])
		propsm = Loess.predict(loess(yrfloatarr, valarr), yrfloatarr)
		propframes[sex] = Dict(:propframe => propframe, :propsm => propsm)
	end
	propframes[:meta] = meta
	if outfile != ""
		propframes_render(propframes, language, sexes, outfile, framedict, showplot)
	end
	propframes
end

function propframes_render(propframes, language, propframeds, outfile, framedict, showplot,
			  x_axis_type = "linear", y_axis_type = "linear")
	preplot(outfile)
	meta = propframes[:meta]
	plclass = meta[:plclass]
	ca1 = meta[:ca1]
	ca2 = meta[:ca2]
	years = meta[:years]
	apctype = meta[:apctype]
	ca1alias = caalias(ca1, language)
	ca2alias = caalias(ca2, language)
	legends = []
	minvals = []
	if plclass == "sexesyrs"
		country = meta[:country]
		xcol = meta[:xcol]
		ylab = meta[:ylab]
		ctryalias = conf["countries"][string(country)]["alias"][language]
		figtitle = "$(dthalias(language)) $ca1alias/$ca2alias $ctryalias"
		p = bp[:figure](output_backend = obend, title = figtitle, y_axis_label = ylab,
			toolbar_location = "below", toolbar_sticky = false,
			plot_width = 600, plot_height = 600)
		for sex in propframeds
			sexalias = conf["sexes"][string(sex)]["alias"][language]
			col = conf["sexes"][string(sex)]["color"]
			propframed = propframes[sex]
			propframe = propframed[:propframe]
			propsm = propframed[:propsm]
			sexcirc = p[:circle](propframe[xcol], propframe[:value], color = col)
			sexline = p[:line](propframe[xcol], propsm, color = col)
			legends = vcat(legends,
				[(sexalias, [sexcirc]); ("$sexalias loess", [sexline])]...)
			minvals = vcat(minvals, minimum(propframe[:value]))
		end
		legend = bml[:Legend](items = legends, location = (0, -30))
		p[:add_layout](legend, "right")
		if apctype == "pa"
			p[:add_layout](listlabels(country, years, minvals, framedict))
		end
	elseif plclass == "agesyrs"
		country = meta[:country]
		sex = meta[:sex]
		ctryalias = conf["countries"][string(country)]["alias"][language]
		sexalias = conf["sexes"][string(sex)]["alias"][language]
		figtitle = "$(dthalias(language)) $ca1alias/$ca2alias $sexalias $ctryalias"
		p = bp[:figure](output_backend = obend, title = figtitle, y_axis_type = y_axis_type,
			toolbar_location = "below", toolbar_sticky = false,
			plot_width = 600, plot_height = 600)
		for agetuple in propframeds
			propframed = propframes[agetuple]
			propframe = propframed[:frame]
			ages = propframed[:ages]
			agealias = ages[:alias]
			if (apctype == "pa" || apctype == "pc")
				ageline = p[:line](propframe[:Year], propframe[:value], color = ages[:color])
				legends = vcat(legends, (ages[:alias], [ageline]))
			end
		end
		legend = bml[:Legend](items = legends, location = (0, -30))
		p[:add_layout](legend, "right")
	end
	p[:add_tools](bml[:CrosshairTool]())
	postplot(p, showplot)
end

function sexesrat(sexesyrsplot, numsex, denomsex, xcol)
	xs = sexesyrsplot[numsex][:propframe][xcol]
	values = sexesyrsplot[numsex][:propframe][:value]./sexesyrsplot[denomsex][:propframe][:value]
	df = DataFrame()
	df[xcol] = xs
	df[:value] = values
	df
end

function srataliases(sexesyrsplots, language)
	ca1 = map((s)->caalias(s[:meta][:ca1], language), sexesyrsplots)
	ca2 = map((s)->caalias(s[:meta][:ca2], language), sexesyrsplots)
	country = map((s)->conf["countries"][string(s[:meta][:country])]["alias"][language],
		sexesyrsplots)
	ages = map((s)->ageslice(s[:meta][:sage], s[:meta][:eage],
		s[:meta][:agemean], language)[:alias], sexesyrsplots)
	years = map((s)->string(s[:meta][:years]), sexesyrsplots)
	acbcohs = map((s)->string(s[:meta][:acbcoh]), sexesyrsplots)
	aldict = Dict(:ca1 => ca1, :ca2 => ca2, :country => country, :ages => ages,
		:years => years, :acbcohs => acbcohs)
end

unicount(i) = size(unique(i))[1]

# Takes an array of sexesyrsplots with the same apctype.
function plot_sexesrats(sexesyrsplots, language, outfile, showplot,
	sexesratio = true, numsex = 1, denomsex = 2, y_axis_type = "linear")
	preplot(outfile)
	numsexalias = conf["sexes"][string(numsex)]["alias"][language]
	denomsexalias = conf["sexes"][string(denomsex)]["alias"][language]
	srals = srataliases(sexesyrsplots, language)
	figtitle = "$(dthalias(language))"
	apctype = sexesyrsplots[1][:meta][:apctype]
	if sexesratio
		figtitle = "$figtitle $numsexalias/$denomsexalias"
	else
		figtitle = "$figtitle $numsexalias"
	end
	ylab = ""
	if unicount(srals[:ca1]) == 1
		figtitle = "$figtitle $(srals[:ca1][1])"
	end
	if unicount(srals[:ca2]) == 1
		figtitle = "$figtitle/$(srals[:ca2][1])"
	end
	if unicount(srals[:country]) == 1
		figtitle = "$figtitle, $(srals[:country][1])"
	end
	if (apctype == "pa" || apctype == "pc")
		xcol = :Year
		if unicount(srals[:ages]) == 1
			ylab = srals[:ages][1]
		end
	elseif apctype == "ap"
		xcol = :Age
		if unicount(srals[:years]) == 1
			ylab = srals[:years][1]
		end
	elseif apctype == "ac"
		xcol = :Age
		if unicount(srals[:acbcohs]) == 1
			ylab = srals[:acbcohs][1]
		end
	end
	p = bp[:figure](output_backend = obend, title = figtitle, y_axis_label = ylab,
		toolbar_location = "below", toolbar_sticky = false,
		plot_width = 800, plot_height = 600, y_axis_type = y_axis_type)
	sratlegends = []
	for (i, syplot) in enumerate(sexesyrsplots)
		meta = syplot[:meta]
		col = bpal[:Category20][20][mod1(i, 20)]
		if sexesratio
			sratframe = sexesrat(syplot, numsex, denomsex, xcol)
			yrfloatarr = convert(Array{Float64}, sratframe[xcol])
			valarr = convert(Array, sratframe[:value])
			sratsm = Loess.predict(loess(yrfloatarr, valarr), yrfloatarr)
		else
			sratframe = syplot[numsex][:propframe]
			sratsm = syplot[numsex][:propsm]
		end
		sratcirc = p[:circle](sratframe[xcol], sratframe[:value], color = col)
		sratline = p[:line](sratframe[xcol], sratsm, color = col)
		sratleg = ""
		if unicount(srals[:ca1]) > 1
			sratleg = "$(srals[:ca1][i])"
		end
		if unicount(srals[:ca2]) > 1
			sratleg = "$sratleg/$(srals[:ca2][i])"
		end
		if unicount(srals[:country]) > 1
			sratleg = "$sratleg, $(srals[:country][i])"
		end
		if (apctype == "pa" || apctype == "pc")
			if unicount(srals[:ages]) > 1
				sratleg = "$sratleg, $(srals[:ages][i])"
			end
		elseif apctype == "ap"
			if unicount(srals[:years]) > 1
				sratleg = "$sratleg, $(srals[:years][i])"
			end
		elseif apctype == "ac"
			if unicount(srals[:acbcohs]) > 1
				sratleg = "$sratleg, $(srals[:acbcohs][i])"
			end
		end
		sratlegends = vcat(sratlegends, (sratleg, [sratcirc]))
	end
	legend = bml[:Legend](items = sratlegends, location = (0, -30))
	p[:add_layout](legend, "right")
	p[:add_tools](bml[:CrosshairTool]())
	postplot(p, showplot)
end

function propplot_agesyrs(ca1, ca2, sex, country, agetuples, years, agemean,
	framedict, language, y_axis_type, outfile, showplot, apctype = "pa")
	pframes = propplotframes(ca1, ca2, framedict, language)
	minvals = []
	meta = Dict(:ca1 => ca1, :ca2 => ca2, :sex => sex, :country => country,
		    :apctype => apctype, :years => years,
		    :agemean => agemean, :plclass => "agesyrs")
	propframes = Dict()
	for agetuple in agetuples
		ages = ageslice(agetuple[1], agetuple[2], agemean, language)
		agelist = ages[:agelist]
		propframe = propgrp(pframes[:ca1frame], pframes[:ca2frame],
			sex, country, agelist, years, agemean, :Year)
		if (apctype == "pc" || apctype == "ac")
			propframe[:Year] = propframe[:Year].-ages[:agest]
		end
		minvals = vcat(minvals, minimum(propframe[:value]))
		propframes[agetuple] = Dict(:frame => propframe, :ages => ages)
	end
	if (apctype == "ap" || apctype == "ac")
		framedicts = values(propframes)
		ageframes = map((agdict)->DataFrame(Year = agdict[:frame][:Year],
			value = agdict[:frame][:value],
			Age = agdict[:ages][:agest]),
			framedicts)
		ageframe = vcat(ageframes...)
		propframes = Dict()
		uniyears = sort!(unique(ageframe[:Year]))
		for year in uniyears
			propframe = ageframe[ageframe[:Year].==year, :]
			sort!(propframe, cols=[:Age])
			propframes[year] = Dict(:frame => propframe)
		end
	end
	propframes[:meta] = meta
	if outfile != ""
		propframes_render(propframes, language, agetuples, outfile, framedict, showplot)
	end
	propframes
end

function propscat_yrsctry(ca1, ca2, sex, countries, sage, eage, year1, year2, agemean,
	framedict, language, outfile, showplot)
	preplot(outfile)
	pframes = propplotframes(ca1, ca2, framedict, language)
	ages = ageslice(sage, eage, agemean, language)
	agealias = ages[:alias]
	agelist = ages[:agelist]
	sexalias =  conf["sexes"][string(sex)]["alias"][language]
	figtitle = "$(dthalias(language)) $(pframes[:ca1alias])/$(pframes[:ca2alias])\n$agealias $sexalias"
	meta = Dict(:ca1 => ca1, :ca2 => ca2, :sex => sex, :countries => countries,
		:sage => sage, :eage => eage, :agemean => agemean)
	yr1propframe = propgrp(pframes[:ca1frame], pframes[:ca2frame], sex,
		countries, agelist, year1, agemean, :Country)
	yr2propframe = propgrp(pframes[:ca1frame], pframes[:ca2frame], sex,
		countries, agelist, year2, agemean, :Country)
	propframe = join(yr1propframe, yr2propframe, on=:Country)
	isos = map((c)->conf["countries"][string(c)]["iso3166"], propframe[:Country])
	ctrynames = map((c)->conf["countries"][string(c)]["alias"][language], propframe[:Country])
	scatdata = bml[:ColumnDataSource](data = Dict("year1prop" => propframe[:value],
		"year2prop" => propframe[:value_1], "isos" => isos, "ctrynames" => ctrynames))
	hover = bml[:HoverTool](tooltips =
			[("befolkning", "@ctrynames"),
			("$year1", "@year1prop"),
			("$year2", "@year2prop")])
	p = bp[:figure](output_backend = obend, title = figtitle, x_axis_label = "$year1",
		y_axis_label = "$year2", plot_width = 600, plot_height = 600)
	p[:add_tools](hover)
	p[:circle](x = "year1prop", y = "year2prop", size = 12, source = scatdata)
	isolabels = bml[:LabelSet](x = "year1prop", y = "year2prop", text = "isos",
		level = "glyph", x_offset = 5, y_offset = 5, source = scatdata)
	p[:add_layout](isolabels)
	postplot(p, showplot)
	Dict(:meta => meta, year1 => yr1propframe, year2 => yr2propframe)
end

function propscat_sexesctry(ca1, ca2, countries, sage, eage, year, agemean,
	framedict, language, outfile, showplot)
	pframes = propplotframes(ca1, ca2, framedict, language)
	preplot(outfile)
	ages = ageslice(sage, eage, agemean, language)
	agealias = ages[:alias]
	agelist = ages[:agelist]
	figtitle = "$(dthalias(language)) $(pframes[:ca1alias])/$(pframes[:ca2alias])\n$agealias $year"
	meta = Dict(:ca1 => ca1, :ca2 => ca2, :countries => countries,
		:sage => sage, :eage => eage, :agemean => agemean)
	fempropframe = propgrp(pframes[:ca1frame], pframes[:ca2frame], 2,
		countries, agelist, year, agemean, :Country)
	malepropframe = propgrp(pframes[:ca1frame], pframes[:ca2frame], 1,
		countries, agelist, year, agemean, :Country)
	propframe = join(fempropframe, malepropframe, on=:Country)
	femalias =  conf["sexes"]["2"]["alias"][language]
	malealias = conf["sexes"]["1"]["alias"][language]
	isos = map((c)->conf["countries"][string(c)]["iso3166"], propframe[:Country])
	ctrynames = map((c)->conf["countries"][string(c)]["alias"][language], propframe[:Country])
	scatdata = bml[:ColumnDataSource](data = Dict("femprop" => propframe[:value],
		"maleprop" => propframe[:value_1], "isos" => isos, "ctrynames" => ctrynames))
	hover = bml[:HoverTool](tooltips =
			[("befolkning", "@ctrynames"),
			("$femalias", "@femprop"),
			("$malealias", "@maleprop")])
	p = bp[:figure](output_backend = obend, title = figtitle, x_axis_label = femalias,
		y_axis_label = malealias, plot_width = 600, plot_height = 600)
	p[:add_tools](hover)
	p[:circle](x = "femprop", y = "maleprop", size = 12, source = scatdata)
	isolabels = bml[:LabelSet](x = "femprop", y = "maleprop", text = "isos",
		level = "glyph", x_offset = 5, y_offset = 5, source = scatdata)
	p[:add_layout](isolabels)
	postplot(p, showplot)
	Dict(:meta => meta, 2 => fempropframe, 1 => malepropframe)
end

batchages_caflt(cause) = filter((age)->
	age["ca2"]!=cause &&
	((age["startage"]==1 && age["endage"]==1) ||
	(!(haskey(conf["causes"][cause], "lowerage")) ||
	age["endage"]>=conf["causes"][cause]["lowerage"]) &&
	(!(haskey(conf["causes"][cause], "upperage")) ||
	age["startage"]<=conf["causes"][cause]["upperage"])),
	conf["batchages"])

fname_sexesyrs(ca1, ca2, country, sage, eage, agemean) =
*(ca1, ca2, country, "s", string(sage), "e", string(eage),
"mean", string(agemean), ".html")

agedict_bothsexes(sage, eage, agemean, language, ca2, fname) =
Dict("alias" => *(ageslice(sage, eage, agemean, language)[:alias],
	"/$(caalias(ca2, language))"),
	"fname" => fname)
agedict_sex(sage, eage, sex, agemean, language, ca2, fname) =
Dict("alias" => *(ageslice(sage, eage, agemean, language)[:alias],
	"/$(caalias(ca2, language)) ",
	conf["sexes"][string(sex)]["alias"][language]),
	"fname" => fname)

cadict(cause, children, language) = Dict("alias" => uppercasefirst(caalias(cause, language)),
"children" => children, "name" => cause)

function ctryints_flt(countries, year1, year2)
	countries_flt = filter((ctry)->conf["countries"][ctry]["startyear"]<=year1
				&& conf["countries"][ctry]["endyear"]>=year2, countries)
	map((c)->parse(Int, c), countries_flt)
end

fname_yrsctry(ca1, ca2, sex, sage, eage, agemean, year1, year2) =
 *(ca1, ca2, string(sex), "s", string(sage), "e", string(eage),
"mean", string(agemean), "ctries", string(year1), "vs", string(year2), ".html")

causes_flt_skipyrs(causes, year1, year2) =
filter((ca)-> !(haskey(conf["causes"][ca], "skipyrs")) ||
(!(year1 in conf["causes"][ca]["skipyrs"]) &&
!(year2 in conf["causes"][ca]["skipyrs"])),
causes)

fname_sexesctry(ca1, ca2, sage, eage, agemean, year) =
*(ca1, ca2, "s", string(sage), "e", string(eage), "mean", string(agemean),
"sexesctries", string(year), ".html")

function agebatchplot(framedict, age, plottype, language, ca1, child, countries, yr1, yr2, sexes)
	ca2 = age["ca2"]
	sage = age["startage"]
	eage = age["endage"]
	agemean = age["agemean"]
	if plottype == "sexesyrs"
		fname = fname_sexesyrs(ca1, ca2, child, sage, eage, agemean)
		outfile = normpath(chartpath, fname)
		propplot_sexesyrs(ca1, ca2, sexes, parse(Int, child),
			sage, eage, yr1:yr2,
			agemean, framedict, language, outfile, false)
		return agedict_bothsexes(sage, eage, agemean, language, ca2, fname)
	elseif plottype == "yrsctry"
		if ca1 in causes_flt_skipyrs(keys(conf["causes"]), yr1, yr2)
			agesexdicts = []
			ctryints = ctryints_flt(countries, yr1, yr2)
			for sex in sexes
				fname = fname_yrsctry(ca1, ca2, sex, sage, eage,
					agemean, yr1, yr2)
				outfile = normpath(chartpath, fname)
				propscat_yrsctry(ca1, ca2, sex, ctryints,
					sage, eage, yr1, yr2,
					agemean, framedict, language, outfile, false)
				agesexdict = agedict_sex(sage, eage, sex,
					agemean, language, ca2, fname)
				agesexdicts = vcat(agesexdicts, agesexdict)
			end
			return agesexdicts
		end
	elseif plottype == "sexesctry"
		if (ca1 in causes_flt_skipyrs(keys(conf["causes"]), yr1, yr2) &&
			conf["causes"][ca1]["sex"]==[2;1])
			ctryints = ctryints_flt(countries, yr1, yr2)
			year = yr1
			fname = fname_sexesctry(ca1, ca2, sage, eage, agemean, year)
			outfile = normpath(chartpath, fname)
			propscat_sexesctry(ca1, ca2, ctryints,
				sage, eage, year, agemean,
				framedict, language, outfile, false)
			return agedict_bothsexes(sage, eage, agemean, language, ca2, fname)
		end
	end
end

function batchplotalias(child, plottype, language)
	if plottype == "sexesyrs"
		return conf["countries"][child]["alias"][language]
	elseif plottype == "yrsctry"
		return "$(child[1]) vs $(child[2])"
	elseif plottype == "sexesctry"
		return "$child"
	end
end

function batchplotyrs(child, plottype)
	if plottype == "sexesyrs"
		return [conf["countries"][child]["startyear"];
			conf["countries"][child]["endyear"]]
	elseif plottype == "yrsctry"
		return child
	elseif plottype == "sexesctry"
		return [child; child]
	end
end

function childplot(framedict, language, plottype, ca1, countries, years, yrtups)
	ages = batchages_caflt(ca1)
	sexes = conf["causes"][ca1]["sex"]
	childdicts = []
	if plottype == "sexesyrs"
		children = countries
	elseif plottype == "yrsctry"
		children = yrtups
	elseif plottype == "sexesctry"
		children = years
	end
	for child in children
		print("$child\n")
		agedicts = []
		ages = batchages_caflt(ca1)
		yrs = batchplotyrs(child, plottype)
		for age in ages
			print("$age\n")
			agedict = agebatchplot(framedict, age, plottype, language, ca1, child,
				countries, yrs[1], yrs[2], sexes)
			agedicts = vcat(agedicts, agedict)
		end
		if !(agedicts[1]==nothing)
			childdict = Dict(
				"alias" => batchplotalias(child, plottype, language),
				"ages" => agedicts)
			childdicts = vcat(childdicts, childdict)
		end
	end
	childdicts
end

function batchplot(framedict, language, plottype, causes = collect(keys(conf["causes"])),
	countries = collect(keys(conf["countries"])), years = conf["batchsexesyrs"],
	yrtups = conf["batchyrtups"])
	sort!(causes, by=((c)->
		(conf["causes"][c]["causeclass"],
		!(conf["causes"][c]["classtot"]),
		conf["causes"][c]["alias"][language])))
	sort!(countries, by=((c)->conf["countries"][c]["alias"][language]))
	cadicts = []
	if plottype == "sexesctry"
		causes = filter((ca)->conf["causes"][ca]["sex"]==[2;1], causes)
	end
	for ca1 in causes
		print("$plottype, $ca1\n")
		children = childplot(framedict, language, plottype, ca1,
			countries, years, yrtups)
		cadicts = vcat(cadicts, cadict(ca1, children, language))
	end
	Dict(
		"plottype" => plottype,
		"cadicts" => cadicts
	)
end

function writeplotlist(batchplotdict, outfile)
	tpl = read(normpath(datapath, "plotlist.mustache"), String)
	write(outfile, render(tpl, batchplotdict))
end

function ccflt(causeclass, language)
	cavals = filter((c)->c["causeclass"]==causeclass, collect(values(conf["causes"])))
	cadicts = map((c)->Dict("alias" => uppercasefirst(c["alias"][language]),
		"classtot" => c["classtot"], "codedesc" => c["codedesc"][language],
		"note" => c["note"][language]), cavals)
	sort!(cadicts, by=((c)->(!(c["classtot"]), c["alias"])))
end

function writetempl(language, fname, sitepath)
	haktemplpath = normpath(sitepath, "templates")
	mkpath(haktemplpath)
	tpl = read(normpath(datapath, "$fname.mustache"), String)
	maintempl = Dict()
	if (fname == "default" || fname == "site")
		maintempl["maintempldicts"] = map((p)->
			Dict("plottype" => p, "alias" => conf["plottypes"][p]["alias"][language]),
			keys(conf["plottypes"]))
		if fname == "default"
			outpath = haktemplpath
			ext = "html"
			maintempl["language"] = language
			for siteelem in ["dochead"; "homehead"; "sitehead";
				"sitenavhead"; "blognavhead"; "othlanghead"]
				maintempl[siteelem] = conf[siteelem][language]
			end
			maintempl["othlangs"] = map((od)->Dict(od),
				conf["othlangs"][language])
			maintempl["footer"] =
				read(normpath(datapath, "footer-$language.$ext"), String)
		elseif fname == "site"
			outpath = sitepath
			maintempl["indexpagetitle"] = conf["indexpagetitle"][language]
			ext = "hs"
		end
	elseif (fname == "index" || fname == "mortchartdoc")
		if fname == "index"
			ext = "html"
		elseif fname == "mortchartdoc"
			ext = "md"
			maintempl["maintempldicts"] = map((cc)->
				Dict("class" => cc,
				"alias" => conf["causeclasses"][cc]["alias"][language],
				"causes" => ccflt(parse(Int, cc), language)),
				keys(conf["causeclasses"]))
			maintempl["refhead"] = conf["refhead"][language]
		end
		outpath = sitepath
		maintempl["body"] = read(normpath(datapath, "$fname-$language.$ext"), String)
	end
	write(normpath(outpath, "$fname.$ext"), render(tpl, maintempl = maintempl))
end

function writeplotsite(framedict, language,
	sitepath = normpath(mainpath, "mortchart-site"))
	sitesubpaths = DataStructures.OrderedDict(
		"siteroot" => Dict("path" => "", "files" =>
			["default.csl"; "mortchartdoc_biber.bib"]),
		"css" => Dict("path" => "css", "files" => ["default.css"]),
		"images" => Dict("path" => "images", "files" => ["mortchartico.png"]),
		"charts" => Dict("path" => "charts", "files" => []))
	sitefullsubpaths = DataStructures.OrderedDict()
	for subpath in keys(sitesubpaths)
		sitefullsubpath = (normpath(sitepath, sitesubpaths[subpath]["path"]))
		sitefullsubpaths[subpath] = sitefullsubpath
		mkpath(sitefullsubpath)
		for file in sitesubpaths[subpath]["files"]
			cp(normpath(datapath, file), normpath(sitefullsubpath, file),
				remove_destination = true)
		end
	end
	for plottype in keys(conf["plottypes"])
		batchplotdict = batchplot(framedict, language, plottype)
		writeplotlist(batchplotdict, normpath(sitepath, "$plottype.html"))
	end
	chartdestpath = sitefullsubpaths["charts"]
	for chartfile in readdir(chartpath)
		mv(normpath(chartpath, chartfile), normpath(chartdestpath, chartfile),
		remove_destination = true)
	end
	for fname in ["default"; "index"; "site"; "mortchartdoc"]
		writetempl(language, fname, sitepath)
	end
end

end # module


text:
module FixYourWorkaround

using Pkg
using Pkg.Types: VersionSpec, semver_spec
using Test

export package_compatible
export CompatNotFound, PackageNotInCompat, VersionNotCompatible

struct CompatNotFound <: Exception
    message::String
end
Base.showerror(io::IO, e::CompatNotFound) = println(io, e.message)

struct PackageNotInCompat <: Exception
    message::String
end
Base.showerror(io::IO, e::PackageNotInCompat) = println(io, e.message)

struct VersionNotCompatible <: Exception
    message::String
end
Base.showerror(io::IO, e::VersionNotCompatible) = println(io, e.message)

"""
    package_compatible(package_name::String, version::String)

Check to see if package_name@version is inbounds with the compat section in your Project.toml

# Arguments
- `package_name::String`: Name of the package
- `version::String`: Package version to check being inbounds

# Keywords
- `toml_path::String`: Path to the Project.toml file

# Returns
- `Bool`: If the version is in the compat versions

# Throws
- `PackageNotInCompat`: package_name was not found in the compat section
- `CompatNotFound`: Compat section not found in the Project.toml
- `VersionOutsideSpec`: Version is no longer inbounds of the compat versions
"""
function package_compatible(package_name::String, version::String; toml_path=joinpath(@__DIR__, "..", "Project.toml"))
    version = VersionSpec(version)
    toml = Pkg.TOML.parsefile(toml_path)

    if haskey(toml, "compat")
        if haskey(toml["compat"], package_name)
            compat_version = semver_spec(toml["compat"][package_name])

            if !isempty(intersect(version, compat_version))
                return true
            else
                throw(VersionNotCompatible("$package_name@$version is not support"))
            end
        else
            throw(PackageNotInCompat("$package_name not found in compat section"))
        end
    else
        throw(CompatNotFound("Compat section not found in Project.toml"))
    end

    return false
end

package_compatible(package_name::String, version::Int64; kwargs...) = package_compatible(package_name, string(version); kwargs...)
package_compatible(package_name::String, version::Float64; kwargs...) = package_compatible(package_name, string(version); kwargs...)
package_compatible(package_name::String, version::VersionNumber; kwargs...) = package_compatible(package_name, string(version); kwargs...)

end  # module

ext: jl


text:
@testset "Clibrary" begin
	lib1 = Clibrary()
	lib2 = Clibrary()
	@test lib1.handle != C_NULL
	@test lib1 == lib2

	lib3 = Clibrary(Base.libm_name)
	lib4 = Clibrary(Base.libm_name)
	@test lib3.handle != C_NULL
	@test lib3 == lib4

	@test lib1 != lib3
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "ArcadeLearningEnvironmentRoms"
version = v"0.0.1"

# Collection of sources required to complete build
sources = [
    FileSource("http://www.atarimania.com/roms/Roms.rar", "4e35879fbd3da7d008f80f8d3a48360b9513859aa6c694164e67d5a82daca498"),
    FileSource("https://raw.githubusercontent.com/mgbellemare/Arcade-Learning-Environment/v0.6.1/md5.txt", "218673cbeba56f3a7066293c259ae7b31ebece686bf3ff4ae4fe746e7d58a51e"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/
unrar e Roms.rar
unzip ROMS.zip
mkdir $prefix/roms
for f in ROMS/*; do
    md5=$(md5sum "$f" | awk '{print $1}')
    newname=$(grep $md5 md5.txt | awk '{print $2}')
    if [[ $newname != "" ]]; then
        cp "$f" $prefix/roms/$newname
    fi
done
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [AnyPlatform()]

# The products that we will ensure are always built
products = [
    FileProduct("roms", :roms_dir)
]

# Dependencies that must be installed before this package can be built
dependencies = []

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: jbrea/Yggdrasil

max_issues_repo_name: jbrea/Yggdrasil

max_forks_repo_name: jbrea/Yggdrasil


text:
import Knet.Ops20: rnnforw
using Knet.Ops20: RNN
using Knet.KnetArrays: DevArray, KnetArray, Cptr
using CUDA: CuArray, CUDNN, CU_NULL
using AutoGrad: AutoGrad, @primitive1, value, recording, Param, Value

"RNN descriptor"
mutable struct RD; ptr; end

"Dropout descriptor"
mutable struct DD; ptr; states; end

"Keeps an array of 3D tensor descriptors"
mutable struct TDs; pvec::Vector{Cptr}; xDesc::Vector{TD}; end     # Keep xDesc in TDs so it does not get gc'ed

Base.unsafe_convert(::Type{Cptr}, dd::DD)=dd.ptr
Base.unsafe_convert(::Type{Cptr}, rd::RD)=rd.ptr
Base.unsafe_convert(::Type{Ptr{Cptr}}, tds::TDs)=pointer(tds.pvec)

function DD(; atype, handle=CUDNN.handle(), dropout=0.0, seed=0, o...)
    if seed==0; seed=floor(Culonglong,time()); end
    d = Cptr[0]; s = Csize_t[0] # TODO: Can multiple RNNs share dropout descriptors? Can dropout probability be changed?
    CUDNN.cudnnCreateDropoutDescriptor(d)
    CUDNN.cudnnDropoutGetStatesSize(handle,s)
    states = rnnworkspace(s[1], atype)
    @cudnn_retry CUDNN.unsafe_cudnnSetDropoutDescriptor(d[1],handle,dropout,states,bytes(states),seed)
    dd = DD(d[1],states)
    finalizer(x->CUDNN.cudnnDestroyDropoutDescriptor(x.ptr),dd)
    return dd
end

function RD()
    d = Cptr[0]
    @cudnn_retry CUDNN.unsafe_cudnnCreateRNNDescriptor(d)
    rd = RD(d[1])
    finalizer(x->CUDNN.cudnnDestroyRNNDescriptor(x.ptr),rd)
    return rd
end

function RD(hiddenSize,numLayers,dropoutDesc,inputMode,direction,mode,algo,dataType; handle=CUDNN.handle())
    rnnDesc = RD()
    inputMode = CUDNN.cudnnRNNInputMode_t(inputMode)
    direction = CUDNN.cudnnDirectionMode_t(direction)
    mode = CUDNN.cudnnRNNMode_t(mode)
    algo = CUDNN.cudnnRNNAlgo_t(algo)
    dt = CUDNN.cudnnDataType_t(DT(dataType))
    if CUDNN.version() < v"8"
        CUDNN.cudnnSetRNNDescriptor(handle,rnnDesc,hiddenSize,numLayers,dropoutDesc,inputMode,direction,mode,algo,dt)
    else
        CUDNN.cudnnSetRNNDescriptor_v6(handle,rnnDesc,hiddenSize,numLayers,dropoutDesc,inputMode,direction,mode,algo,dt)
    end
    return rnnDesc
end

Base.length(tds::TDs)=length(tds.pvec)

function TDs(x::DevArray{A},::Nothing) where {A} # Treat x: (X,B?,T?) as a 4D array: (1,X,B,T)
    xDesc = TD(A,1,size(x,1),size(x,2)) # we can use a single xDesc
    pvec = Vector{Cptr}(undef, size(x,3))
    pvec[:] .= xDesc.ptr
    return TDs(pvec, [xDesc])
end

function TDs(x::DevArray{A},batchSizes) where {A} # x: (X,B*), batchSizes gives us Bt sizes
    @assert sum(batchSizes) == div(length(x),size(x,1))
    X = size(x,1)
    xs = [ TD(A,1,X,B) for B in batchSizes ]
    ps = [ xd.ptr for xd in xs ]
    return TDs(ps,xs)
end

function TD3(a::DevArray) # Treat a as a 3D array, pad from right
    n = ndims(a)
    if n==3; TD(a)
    elseif n==2; TD(reshape(a, size(a,1), size(a,2), 1))
    elseif n==1; TD(reshape(a, size(a,1), 1, 1))
    else; throw(DimensionMismatch())
    end
end

function FD3(a::DevArray) # Treat a as a 3D array, pad from left
    n = ndims(a)
    if n==3; FD(a)
    elseif n==2; FD(reshape(a, 1, size(a,1), size(a,2)))
    elseif n==1; FD(reshape(a, 1, 1, size(a,1)))
    else; throw(DimensionMismatch())
    end
end

function rnnforw(r::RNN, w, x::Union{DevArray{T},Value{<:DevArray{T}}}, hx=nothing, cx=nothing;
                 handle=CUDNN.handle(), batchSizes=nothing, hy = (hx != nothing), cy = (cx != nothing && r.mode == 2)) where T
    @assert value(w) === value(r.w)
    @assert size(x,1) == r.inputSize
    x3 = reshape(value(x), size(x,1), size(x,2), size(x,3))
    @assert typeof(x3) == typeof(value(w)) "$(typeof(value(w))) weights do not match $(typeof(x)) input. Please use RNN(;atype) option."
    if r.rnnDesc === nothing    # initialize rnn for gpu with first input
        r.dataType = eltype(x3)
        r.dropoutDesc = DD(handle=CUDNN.handle(),dropout=r.dropout,seed=r.seed,atype=typeof(x3))
        r.rnnDesc = RD(r.hiddenSize,r.numLayers,r.dropoutDesc,r.inputMode,r.direction,r.mode,r.algo,r.dataType)
    end
    _rnnforw(w,x,hx,cx; rnn=r,handle=handle,batchSizes=batchSizes,hy=hy,cy=cy)
end

function _rnnforw(w, x, hx, cx; rnn, handle, batchSizes, hy, cy)
    # Input descriptors
    seqLength = batchSizes==nothing ? size(x,3) : length(batchSizes) # (X,B,T) or (X,B+) with batchSizes
    wDesc = FD3(w)              # (1,1,W)
    xtds = TDs(x,batchSizes)    # (1,X,Bt) x T
    isnothing(a) = a === nothing || a === C_NULL || a === CU_NULL
    if hx==nothing; hx=CU_NULL; hxDesc=C_NULL; else; hxDesc=TD3(hx); end # (H,B,L/2L)
    if cx==nothing || rnn.mode != 2; cx=CU_NULL; cxDesc=C_NULL; else; cxDesc=TD3(cx); end

    # Output arrays and descriptors
    ysize = collect(size(x))
    ysize[1] = rnn.hiddenSize * (rnn.direction == 1 ? 2 : 1)
    y = similar(x, ysize...)    # (H/2H,B,T) or (H/2H,B+) -- y mirrors x except for the first dimension
    ytds = TDs(y,batchSizes)    # (1,H/2H,Bt) x T

    # Optionally output hidden and cell of last step
    hyout = cyout = CU_NULL
    hyDesc = cyDesc = C_NULL
    if hy || cy
        firstBatchSize = batchSizes==nothing ? size(x,2) : batchSizes[1]
        hsize = (Int(rnn.hiddenSize), Int(firstBatchSize), Int(rnn.numLayers * (rnn.direction == 1 ? 2 : 1))) # (H,B,L/2L)
        if hy; hyout=similar(y,hsize); hyDesc=TD3(hyout); end
        if cy && rnn.mode==2; cyout=similar(y,hsize); cyDesc=TD3(cyout); end
        if !isnothing(hx) && any(size(hx,i)!=hsize[i] for i=1:3) # compare one by one in case hx is 1-D or 2-D
            throw(DimensionMismatch("size(hx)=$(size(hx)) does not match hsize=$(hsize)"))
        end
        if !isnothing(cx) && rnn.mode == 2 && any(size(cx,i)!=hsize[i] for i=1:3)
            throw(DimensionMismatch("size(cx)=$(size(cx)) does not match hsize=$(hsize)"))
        end
    end

    # workSpace and reserveSpace
    wss = cudnnGetRNNWorkspaceSize(rnn.rnnDesc, xtds; handle=handle)
    ws = rnnworkspace(wss, typeof(value(w)))

    if AutoGrad.recording()
        rss = cudnnGetRNNTrainingReserveSize(rnn.rnnDesc, xtds; handle=handle)
        rs = rnnworkspace(rss, typeof(value(w)))
        @cudnn_retry CUDNN.unsafe_cudnnRNNForwardTraining(handle, rnn.rnnDesc, seqLength, xtds, x, hxDesc, hx, cxDesc, cx, wDesc, w, ytds, y, hyDesc, hyout, cyDesc, cyout, ws, wss, rs, rss)
    else
        rs = nothing
        @cudnn_retry CUDNN.unsafe_cudnnRNNForwardInference(handle, rnn.rnnDesc, seqLength, xtds, x, hxDesc, hx, cxDesc, cx, wDesc, w, ytds, y, hyDesc, hyout, cyDesc, cyout, ws, wss)
    end
    if hyout === CU_NULL; hyout = nothing; end
    if cyout === CU_NULL; cyout = nothing; end
    return y, hyout, cyout, rs, ws
end

function _rnnback(dt, t, w, x, hx, cx; rnn, o...)
    @assert value(rnn.w) === value(w)
    y,hy,cy,rs,ws = value(t)
    dy,dhy,dcy,drs,dws = value(dt)
    rnn=value(rnn); w=value(w); x=value(x); hx=value(hx); cx=value(cx)
    # To prevent dependencies to next iteration we need to clear the Result type from rnn.h,rnn.c
    # We can't do this during forward, because another forward may be run within the same iteration.
    # Doing it here is safe, means the iteration is done and we are taking gradients.
    # Note that this does not work on the cpu and these have to be cleaned by hand.
    # The cpu version is not a primitive and has no back function. (TODO: find better solution)
    rnn.h = value(rnn.h); rnn.c = value(rnn.c)
    _rnnback2(rnn, w, x, y, dy, hx, cx, dhy, dcy, rs, ws; o...)
end

function _rnnback2(r, w, x, y, dy, hx, cx, dhy, dcy, rs, ws;
                   handle=CUDNN.handle(), batchSizes=nothing, o...)
    @assert value(r.w) === value(w)
    # Input descriptors:
    seqLength = batchSizes==nothing ? size(x,3) : length(batchSizes) # (X,B,T) or (X,B+) with batchSizes
    wDesc = FD3(w)              # (1,1,W)
    xtds = TDs(x,batchSizes)    # (X,B,T) -> (1,X,B) x T
    ytds = TDs(y,batchSizes)    # (H/2H,B,T) -> (1,H/2H,B) x T
    # dytds = TDs(dy,batchSizes)  # we use ytds for dytds
    if dy == nothing; dy=zero(y); end
    if hx == nothing; hx=CU_NULL; hxDesc=C_NULL; else; hxDesc=TD3(hx); end
    if cx == nothing || r.mode != 2; cx=CU_NULL; cxDesc=C_NULL; else; cxDesc=TD3(cx); end
    if dhy == nothing; dhy=CU_NULL; dhyDesc=C_NULL; else; dhyDesc=TD3(dhy); end
    if dcy == nothing || r.mode != 2; dcy=CU_NULL; dcyDesc=C_NULL; else; dcyDesc=TD3(dcy); end

    # Output arrays and descriptors:
    dx = similar(x)             # (X,B,T) or (X,B+) with batchSizes
    # dxtds = TDs(dx,batchSizes)  # we use xtds here
    dw = zero(w)               # dw is used additively, so we need zero
    dwDesc = FD3(dw)
    if hx === CU_NULL; dhx=CU_NULL; dhxDesc=C_NULL; else; dhx=similar(hx); dhxDesc=TD3(dhx); end
    if cx === CU_NULL; dcx=CU_NULL; dcxDesc=C_NULL; else; dcx=similar(cx); dcxDesc=TD3(dcx); end

    # workSpace and reserveSpace
    # ws = cudnnWorkSpace()
    wss = bytes(ws)
    rss = bytes(rs)
    @cudnn_retry CUDNN.unsafe_cudnnRNNBackwardData(handle, r.rnnDesc, seqLength, ytds, y, ytds, dy, dhyDesc, dhy, dcyDesc, dcy, wDesc, w, hxDesc, hx, cxDesc, cx, xtds, dx, dhxDesc, dhx, dcxDesc, dcx, ws, wss, rs, rss)
    @cudnn_retry CUDNN.unsafe_cudnnRNNBackwardWeights(handle, r.rnnDesc, seqLength, xtds, x, hxDesc, hx, ytds, y, ws, wss, dwDesc, dw, rs, rss)
    # Update the cache
    if dhx===CU_NULL; dhx=nothing; end
    if dcx===CU_NULL; dcx=nothing; end
    r.dx, r.dhx, r.dcx = dx, dhx, dcx
    return dw
end

@primitive1 _rnnforw(w,x,hx,cx; rnn, o...),dy,y  _rnnback(dy,y,w,x,hx,cx; rnn=rnn, o...) value(rnn).dx value(rnn).dhx value(rnn).dcx

#506: Because r.dx,dhx,dcx may be freed by gcnode, their C_NULL pointers cause trouble in deepcopy.
import Base: deepcopy_internal
function deepcopy_internal(x::RNN, s::IdDict)
    if !haskey(s,x)
        s[x] = RNN(deepcopy_internal(x.w,s), deepcopy_internal(x.h,s), deepcopy_internal(x.c,s), x.inputSize, x.hiddenSize, x.numLayers, x.dropout, x.seed, x.inputMode, x.direction, x.mode, x.algo, x.dataType, deepcopy_internal(x.rnnDesc,s), deepcopy_internal(x.dropoutDesc,s), nothing, nothing, nothing)
    end
    return s[x]
end

function rnnworkspace(n, type)
    n8 = (n-1)÷sizeof(Int)+1
    if type <: KnetArray
        buf = KnetArray{Int}(undef, n8)
    elseif type <: CuArray
        buf = CuArray{Int}(undef, n8)
    else
        error("$type not a known GPU array type.")
    end
    return buf
end

function cudnnGetRNNParamsSize(r::RNN)
    res = Csize_t[0]
    xDesc = TD(r.dataType, 1, r.inputSize, 1)    # xDesc: (1,X,B) where X = inputSize, B is ignored, so assume 1
    dt = CUDNN.cudnnDataType_t(DT(r.dataType))
    CUDNN.cudnnGetRNNParamsSize(CUDNN.handle(), r.rnnDesc, xDesc, res, dt)
    div(res[1], sizeof(r.dataType))
end

# This is buggy, why?
# X,H,L,I = r.inputSize, r.hiddenSize, r.numLayers, rnnids(r)
# biases = L*I
# inputMatrices = (r.inputMode == 1 ? 0 : r.direction == 1 ? I : div(I,2))
# hiddenMatrices = (r.direction == 1 ? (L-1)*I : (L-1)*I + div(I,2))
# biases * H + inputMatrices * X * H + hiddenMatrices * H * H

function cudnnGetRNNWorkspaceSize(rd::RD, tds::TDs; handle=CUDNN.handle())
    res = Csize_t[1]
    CUDNN.cudnnGetRNNWorkspaceSize(handle, rd, length(tds), tds, res)
    return Int(res[1])
end

function cudnnGetRNNTrainingReserveSize(rd::RD, tds::TDs; handle=CUDNN.handle())
    res = Csize_t[1]
    CUDNN.cudnnGetRNNTrainingReserveSize(handle, rd, length(tds), tds, res)
    return Int(res[1])
end

# Return eltype,size
function cudnnGetFilterNdDescriptor(wDesc::FD; nbDimsRequested = 8)
    dataType = Cint[0]
    format = Cint[0]
    nbDims = Cint[0]
    filterDimA = Vector{Cint}(undef,nbDimsRequested)
    CUDNN.cudnnGetFilterNdDescriptor(wDesc, nbDimsRequested, dataType, format, nbDims, filterDimA)
    if nbDims[1] > nbDimsRequested
        cudnnGetFilterNdDescriptor(wDesc::FD; nbDimsRequested = nbDims[1])
    else
        (Float32,Float64,Float16)[1+dataType[1]],
        (filterDimA[nbDims[1]:-1:1]...,)
    end
end

function cudnnGetRNNParam(r::RNN, layer::Integer, id::Integer, par::Integer; useview=false)
    params_are_good =
    ((1 <= par <= 2) &&
     ((r.direction == 0 && 1 <= layer <= r.numLayers) ||
      (r.direction == 1 && 1 <= layer <= 2*r.numLayers)) &&
     ((r.mode == 0 && 1 <= id <= 2) ||
      (r.mode == 1 && 1 <= id <= 2) ||
      (r.mode == 2 && 1 <= id <= 8) ||
      (r.mode == 3 && 1 <= id <= 6)))
    params_are_good || throw(ArgumentError("Bad arguments for rnnparam, please see doc."))
    should_return_nothing =
        ((r.inputMode == 1) &&
         (par == 1) &&
         ((r.mode == 0 && id == 1) ||
          (r.mode == 1 && id == 1) ||
          (r.mode == 2 && 1 <= id <= 4) ||
          (r.mode == 3 && 1 <= id <= 3)) &&
         ((layer == 1) ||
          (r.direction == 1 && layer == 2)))

    i1 = i2 = len = 0
    w = value(r.w)
    @assert isa(w, DevArray)
    T = eltype(w)
    xDesc = TD(T,1,r.inputSize,1)
    wDesc = FD(T,1,1,length(w))
    paramDesc = FD(T,1,1,1,1)
    param = Cptr[0]
    if par == 1 # matrix
        CUDNN.cudnnGetRNNLinLayerMatrixParams(handle, r.rnnDesc, layer-1, xDesc, wDesc, w, id-1, paramDesc, param)
    else # bias
        CUDNN.cudnnGetRNNLinLayerBiasParams(handle, r.rnnDesc, layer-1, xDesc, wDesc, w, id-1, paramDesc, param)
    end
    dt,sz = cudnnGetFilterNdDescriptor(paramDesc)
    if should_return_nothing
        @assert param[1] === C_NULL
        @assert sz == ()
        return nothing
    end
    len = prod(sz)
    i1 = 1 + div(Int(param[1] - pointer(w)), sizeof(T))
    i2 = i1 + len - 1
    if i1 > i2
        @assert should_return_nothing
        nothing
    elseif par == 1 # matrix; weights are transposed
        h = Int(r.hiddenSize)
        reshape(view(r.w, i1:i2),:,h)
    else # bias
        view(r.w, i1:i2)
    end
end

# CuArray specific support: should move this to cuarrays
TD(a::CuArray{T}) where {T} = TD(T, size(a))
FD(a::CuArray{T}) where {T} = FD(T, size(a))
bytes(x::CuArray{T}) where T = length(x)*sizeof(T)

# KnetArray getindex contiguous indices already returns a view.
# We need the following for rnnparam/rnntest to work:
Base.view(A::KnetArray, I::AbstractUnitRange{Int}) = getindex(A, I)

# This supports cpucopy/gpucopy:
import Knet.KnetArrays: _ser
function _ser(x::RNN, s::IdDict, m::Val)
    if !haskey(s,x)
        # we need rd,dd only if there is a gpu, we are not in cpumode,
        # and if we are in jldmode we are loading, not saving
        # if (CUDA.functional() && m != CPUMODE && !(m == JLDMODE && x.rnnDesc != nothing))
        #     dd = DD(dropout=x.dropout,seed=x.seed)
        #     rd = RD(x.hiddenSize,x.numLayers,dd,x.inputMode,x.direction,x.mode,x.algo,x.dataType)
        # else
        #     rd = dd = nothing
        # end

        # 20200806: We no longer need to load/save rd/dd: rnnforw will construct as needed.
        rd = dd = nothing

        # dx, dhx, dcx are temporary fields used by rnnback, they do not need to be copied
        # gcnode sets dx.ptr to C_NULL which breaks serialize, best not to try
        s[x] = RNN(_ser(x.w,s,m), _ser(x.h,s,m), _ser(x.c,s,m), x.inputSize, x.hiddenSize, x.numLayers, x.dropout, x.seed, x.inputMode, x.direction, x.mode, x.algo, x.dataType, rd, dd, nothing, nothing, nothing)
    end
    return s[x]
end

import JLD2

struct JLD2RNN; w; h; c; inputSize; hiddenSize; numLayers; dropout; seed; inputMode; direction; mode; algo; dataType; end
JLD2RNN(x::RNN) = JLD2RNN(x.w, x.h, x.c, x.inputSize, x.hiddenSize, x.numLayers, x.dropout, x.seed, x.inputMode, x.direction, x.mode, x.algo, x.dataType)
RNN(x::JLD2RNN) = RNN(x.w, x.h, x.c, x.inputSize, x.hiddenSize, x.numLayers, x.dropout, x.seed, x.inputMode, x.direction, x.mode, x.algo, x.dataType, nothing, nothing, nothing, nothing, nothing)

JLD2.writeas(::Type{RNN}) = JLD2RNN
JLD2.wconvert(::Type{JLD2RNN}, x::RNN) = JLD2RNN(x)
JLD2.rconvert(::Type{RNN}, x::JLD2RNN) = RNN(x)


text:
using JSON
function dump(net::Network,file_str::String)::Nothing
    dict = Dict("layers" => net.layers,
                "neurons" => net.neurons,
                "weights" => net.weights,
                "biases" => net.biases,
                "z" => net.z,
                "batch" => net.batch,
                "epochs" => net.epochs,
                "η" => net.η)

    open(file_str, "w") do file
        write(file, JSON.json(dict))
    end

    println("Saved to $(file_str)")
    return nothing
end

max_stars_repo_name: ian-double-u/NNfromScratch

max_issues_repo_name: ian-double-u/NNfromScratch

max_forks_repo_name: ian-double-u/NNfromScratch


text:
# Autogenerated wrapper script for ICU_jll for i686-w64-mingw32-cxx11
export libicudata, libicui18n, libicuio, libicutest, libicutu, libicuuc

JLLWrappers.@generate_wrapper_header("ICU")
JLLWrappers.@declare_library_product(libicudata, "icudt69.dll")
JLLWrappers.@declare_library_product(libicui18n, "icuin69.dll")
JLLWrappers.@declare_library_product(libicuio, "icuio69.dll")
JLLWrappers.@declare_library_product(libicutest, "icutest69.dll")
JLLWrappers.@declare_library_product(libicutu, "icutu69.dll")
JLLWrappers.@declare_library_product(libicuuc, "icuuc69.dll")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libicudata,
        "bin\\icudt69.dll",
        nothing,
    )

    JLLWrappers.@init_library_product(
        libicui18n,
        "bin\\icuin69.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libicuio,
        "bin\\icuio69.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libicutest,
        "bin\\icutest69.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libicutu,
        "bin\\icutu69.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libicuuc,
        "bin\\icuuc69.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
export getperm

# computes permutation to cartesian index for elements
# needs lin = reshape(collect(1:length(xC)), (ex, ey, ez))
# but then xC[p[lin[i,j,k]]] works where i,j,k are cartesian
function getperm(xC, yC, zC, ex, ey, ez)
    pz = sortperm(zC)
    tmpY = reshape(yC[pz], (ex*ey, ez))
    tmp_py = [sortperm(tmpY[:,i]) for i in 1:ez]
    py = zeros(Int64,length(pz))
    for i in eachindex(tmp_py)
        n = length(tmp_py[i])
        ii = (i-1) * n + 1
        py[ii : ii + n-1] .= tmp_py[i] .+ ii .- 1
    end
    tmpX = reshape(xC[pz][py], (ex, ey*ez))
    tmp_px = [sortperm(tmpX[:,i]) for i in 1:ey*ez]
    px = zeros(Int64,length(pz))
    for i in eachindex(tmp_px)
        n = length(tmp_px[i])
        ii = (i-1) * n + 1
        px[ii : ii + n-1] .= tmp_px[i] .+ ii .- 1
    end
    p = [pz[py[px[i]]] for i in eachindex(px)]
    return p
end

max_stars_repo_name: CliMA/ImperoHooks

max_issues_repo_name: CliMA/ImperoHooks

max_forks_repo_name: CliMA/ImperoHooks


text:
# io = readstring(STDIN)
io = "what,is,the;meaning,of:life."
i  = 0

readbyte!() = io[global i += 1]
writebyte(c) = print(Char(c))

function odd(prev::Function = () -> false)
    a = readbyte!()
    if !isalpha(a)
        prev()
        writebyte(a)
        return a != '.'
    end

    # delay action until later, in the shape of a closure
    clos() = (writebyte(a); prev())

    return odd(clos)
end

function even()
    while true
        c = readbyte!()
        writebyte(c)
        if !isalpha(c) return c != '.' end
    end
end

evn = false
while evn ? odd() : even()
    evn = !evn
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using SafeTestsets

@safetestset "Likelihood" begin
    cd(@__DIR__)
    using Pkg
    Pkg.activate("../")
    using Test, ACTRModels, Distributions, FFTDists
    include("../models/Semantic_FFT/Semantic_FFT_Model.jl")
    include("../models/Semantic_FFT/model_functions.jl")
    Random.seed!(44154)
    #######################################################################################
    #                                   Generate Data
    #######################################################################################
    blc = 1.0
    δ = 1.0
    parms = (noise = true, τ = 0.0, s = .3, mmp = true)
    stimuli = get_stimuli()
    n_reps = 2_000
    data = map(s -> simulate(parms, s, n_reps; blc, δ), stimuli)

    x = range(.8*blc, 1.2*blc, length=100)
    y = map(x->loglike(data, parms, x, δ), x)
    mxv,mxi = findmax(y)
    blc′ = x[mxi]
    @test blc′≈ blc rtol=.02

    x = range(.8*δ, 1.2*δ, length=100)
    y = map(x->loglike(data, parms, blc, x), x)
    mxv,mxi = findmax(y)
    δ′ = x[mxi]
    @test δ′ ≈ δ rtol=.02
end

@safetestset "Semantic FFT Zero Chains" begin
    using Test, ACTRModels, Distributions, FFTDists, KernelDensity
    include("../models/Semantic_FFT/Semantic_FFT_Model.jl")
    include("../models/Semantic_FFT/model_functions.jl")
    Random.seed!(354301)
    #######################################################################################
    #                                   Generate Data
    #######################################################################################
    blc = 1.0
    δ = 1.0
    parms = (noise = true, τ = 0.0, s = .2, mmp = true)
    stimuli = get_stimuli()
    n_reps = 100_000
    stimulus = stimuli[1]
    data = simulate(parms, stimulus, n_reps; blc, δ)
    sim_dist_no = kde(data.no_rts)
    sim_dist_yes = kde(data.yes_rts)
    p_correct = length(data.yes_rts) / n_reps

    chunks = populate_memory()
    memory = Declarative(;memory=chunks)
    actr = ACTR(;declarative=memory, parms..., blc, δ, noise=false)
    x = .8:.01:1.8

    exact_density_yes = map(i -> zero_chain_yes(actr, stimulus, i) |> exp, x)
    sim_density_yes = pdf(sim_dist_yes, x) * p_correct
    @test .15 > maximum(abs.(exact_density_yes .- sim_density_yes))

    exact_density_no = map(i -> zero_chain_no(actr, stimulus, i), x)
    exact_density_no .= exp.(exact_density_no)
    sim_density_no = pdf(sim_dist_no, x) * (1 - p_correct)
    @test .05 > maximum(abs.(exact_density_no .- sim_density_no))
end

@safetestset "Semantic FFT One Chain" begin
    using Test, ACTRModels, Distributions, FFTDists, KernelDensity
    include("../models/Semantic_FFT/Semantic_FFT_Model.jl")
    include("../models/Semantic_FFT/model_functions.jl")
    Random.seed!(354301)
    #######################################################################################
    #                                   Generate Data
    #######################################################################################
    blc = 1.0
    δ = 1.0
    parms = (noise = true, τ = 0.0, s = .2, mmp = true)
    stimuli = get_stimuli()
    n_reps = 100_000
    stimulus = stimuli[2]
    data = simulate(parms, stimulus, n_reps; blc, δ)
    sim_dist_no = kde(data.no_rts)
    sim_dist_yes = kde(data.yes_rts)
    p_correct = length(data.yes_rts) / n_reps

    chunks = populate_memory()
    memory = Declarative(;memory=chunks)
    actr = ACTR(;declarative=memory, parms..., blc, δ, noise=false)
    x = .8:.01:2.5

    exact_density_yes = map(i -> one_chain_yes(actr, stimulus, i), x)
    exact_density_yes .= exp.(exact_density_yes)
    sim_density_yes = pdf(sim_dist_yes, x) * p_correct
    @test .15 > maximum(abs.(exact_density_yes .- sim_density_yes))

    exact_density_no = map(i -> one_chain_no(actr, stimulus, i), x)
    exact_density_no .= exp.(exact_density_no)
    sim_density_no = pdf(sim_dist_no, x) * (1 - p_correct)
    @test .05 > maximum(abs.(exact_density_no .- sim_density_no))
end

@safetestset "Semantic FFT all no one Chain" begin
    using Test, ACTRModels, Distributions, FFTDists, KernelDensity
    include("../models/Semantic_FFT/Semantic_FFT_Model.jl")
    include("../models/Semantic_FFT/model_functions.jl")
    Random.seed!(354301)
    #######################################################################################
    #                                   Generate Data
    #######################################################################################
    blc = 1.0
    δ = 1.0
    parms = (noise = true, τ = 0.0, s = .2, mmp = true)
    stimuli = get_stimuli()
    n_reps = 100_000
    stimulus = stimuli[3]
    data = simulate(parms, stimulus, n_reps; blc, δ)
    sim_dist_no = kde(data.no_rts)
    @test isempty(data.yes_rts)

    chunks = populate_memory()
    memory = Declarative(;memory=chunks)
    actr = ACTR(;declarative=memory, parms..., blc, δ, noise=false)
    x = .8:.01:2.5

    exact_density_no = map(i -> one_chain_no(actr, stimulus, i), x)
    exact_density_no .= exp.(exact_density_no)
    sim_density_no = pdf(sim_dist_no, x)
    @test .10 > maximum(abs.(exact_density_no .- sim_density_no))
end

@safetestset "Semantic FFT all no two Chains" begin
    using Test, ACTRModels, Distributions, FFTDists, KernelDensity
    include("../models/Semantic_FFT/Semantic_FFT_Model.jl")
    include("../models/Semantic_FFT/model_functions.jl")
    Random.seed!(354301)
    #######################################################################################
    #                                   Generate Data
    #######################################################################################
    blc = 1.0
    δ = 1.0
    parms = (noise = true, τ = 0.0, s = .2, mmp = true)
    stimuli = get_stimuli()
    n_reps = 100_000
    stimulus = stimuli[4]
    data = simulate(parms, stimulus, n_reps; blc=blc, δ=δ)
    sim_dist_no = kde(data.no_rts)
    @test isempty(data.yes_rts)

    chunks = populate_memory()
    memory = Declarative(;memory=chunks)
    actr = ACTR(;declarative=memory, parms..., blc, δ, noise=false)
    x = .8:.01:3.0

    exact_density_no = map(i -> two_chains_no(actr, stimulus, i), x)
    exact_density_no .= exp.(exact_density_no)
    sim_density_no = pdf(sim_dist_no, x)
    @test .10 > maximum(abs.(exact_density_no .- sim_density_no))
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Desription
# ==============================================================================
#
#   Tests related to conversion from direction cosine matrices to Euler angle
#   and axis.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

# File: ./src/conversions/dcm_to_angleaxis.jl
# ===========================================

# Functions: dcm_to_angleaxis
# ---------------------------

@testset "DCM => Euler angle and axis (Float64)" begin
    T = Float64

    # Create a random DCM.
    D = angle_to_dcm(_rand_ang(T), :X) *
        angle_to_dcm(_rand_ang(T), :Y) *
        angle_to_dcm(_rand_ang(T), :Z)

    # Convert to Euler angle and axis.
    av = dcm_to_angleaxis(D)
    @test eltype(av) === T

    # Check if the rotation expressed by D is consistent, which is performed in
    # two steps:
    #
    #   1. A vector aligned with `v` does not change.
    #   2. A vector perpendicular to `v` is rotated by `a`.

    v = av.v
    vr = D * v
    @test vr ≈ v
    @test eltype(v) === eltype(vr) === T

    # Auxiliary vector to obtain a vector perpendicular to `v`.
    aux = @SVector randn(3)
    aux = aux / norm(aux)

    vp = av.v × aux
    vp = vp / norm(vp)
    vpr = D * vp
    @test eltype(vp) === eltype(vpr) === T

    # Compute the angle between vp and vpr in [0, 2π].
    a = acos(vp ⋅ vpr)
    @test a ≈ av.a
end

@testset "DCM => Euler angle and axis (Float32)" begin
    T = Float32

    # Create a random DCM.
    D = angle_to_dcm(_rand_ang(T), :X) *
        angle_to_dcm(_rand_ang(T), :Y) *
        angle_to_dcm(_rand_ang(T), :Z)

    # Convert to Euler angle and axis.
    av = dcm_to_angleaxis(D)
    @test eltype(av) === T

    # Check if the rotation expressed by D is consistent, which is performed in
    # two steps:
    #
    #   1. A vector aligned with `v` does not change.
    #   2. A vector perpendicular to `v` is rotated by `a`.

    v = av.v
    vr = D * v
    @test vr ≈ v
    @test eltype(v) === eltype(vr) === T

    # Auxiliary vector to obtain a vector perpendicular to `v`.
    aux = @SVector randn(T, 3)
    aux = aux / norm(aux)

    vp = av.v × aux
    vp = vp / norm(vp)
    vpr = D * vp
    @test eltype(vp) === eltype(vpr) === T

    # Compute the angle between vp and vpr in [0, 2π].
    a = acos(vp ⋅ vpr)
    @test a ≈ av.a
end

@testset "DCM => Euler angle and axis (Special cases)" begin
    # Float64
    # ==========================================================================

    T = Float64

    D = DCM(T(1) * I)
    av = dcm_to_angleaxis(D)
    @test eltype(av) === T
    @test av.a ≈ 0
    @test av.v ≈ [0, 0, 0]

    D = angle_to_dcm(T(π), :X)
    av = dcm_to_angleaxis(D)
    @test eltype(av) === T
    @test av.a ≈ π
    @test av.v ≈ [1, 0, 0]

    D = angle_to_dcm(T(π), :Y)
    av = dcm_to_angleaxis(D)
    @test eltype(av) === T
    @test av.a ≈ π
    @test av.v ≈ [0, 1, 0]

    D = angle_to_dcm(T(π), :Z)
    av = dcm_to_angleaxis(D)
    @test eltype(av) === T
    @test av.a ≈ π
    @test av.v ≈ [0, 0, 1]

    D = angle_to_dcm(T(2π / 3), :X)
    av = dcm_to_angleaxis(D)
    @test eltype(av) === T
    @test av.a ≈ T(2π / 3)
    @test av.v ≈ [1, 0, 0]

    D = angle_to_dcm(T(-2π / 3), :X)
    av = dcm_to_angleaxis(D)
    @test eltype(av) === T
    @test av.a ≈ T(2π / 3)
    @test av.v ≈ [-1, 0, 0]

    D = angle_to_dcm(T(2π / 3), :Y)
    av = dcm_to_angleaxis(D)
    @test eltype(av) === T
    @test av.a ≈ T(2π / 3)
    @test av.v ≈ [0, 1, 0]

    D = angle_to_dcm(T(-2π / 3), :Y)
    av = dcm_to_angleaxis(D)
    @test eltype(av) === T
    @test av.a ≈ T(2π / 3)
    @test av.v ≈ [0, -1, 0]

    D = angle_to_dcm(T(2π / 3), :Z)
    av = dcm_to_angleaxis(D)
    @test eltype(av) === T
    @test av.a ≈ T(2π / 3)
    @test av.v ≈ [0, 0, 1]

    D = angle_to_dcm(T(-2π / 3), :Z)
    av = dcm_to_angleaxis(D)
    @test eltype(av) === T
    @test av.a ≈ T(2π / 3)
    @test av.v ≈ [0, 0, -1]

    D = angle_to_dcm(T(-2π / 3), :Z) * angle_to_dcm(T(-π / 2), :Y)
    av = dcm_to_angleaxis(D)
    @test eltype(av) === T
    @test av.a ≈ 2.4188584057763776
    @test av.v ≈ [0.6546536707079772, -0.37796447300922736, -0.6546536707079772]

    D = DCM(
        T(-1 / 3), T(-2 / 3), T(-1 / 3),
        T(-1 / 3), T(-1 / 3), T(+2 / 3),
        T(-2 / 3), T(+2 / 3), T(-1 / 3)
    )'
    av = dcm_to_angleaxis(D)
    @test eltype(av) === T
    @test av.a ≈ T(π)
    @test av.v ≈ [T(√3 / 3), T(-√3 / 3), T(-√3 / 3)]

    # Float32
    # ==========================================================================

    T = Float32

    D = DCM(T(1) * I)
    av = dcm_to_angleaxis(D)
    @test eltype(av) === T
    @test av.a ≈ 0
    @test av.v ≈ [0, 0, 0]

    D = angle_to_dcm(T(π), :X)
    av = dcm_to_angleaxis(D)
    @test eltype(av) === T
    @test av.a ≈ π
    @test av.v ≈ [1, 0, 0]

    D = angle_to_dcm(T(π), :Y)
    av = dcm_to_angleaxis(D)
    @test eltype(av) === T
    @test av.a ≈ π
    @test av.v ≈ [0, 1, 0]

    D = angle_to_dcm(T(π), :Z)
    av = dcm_to_angleaxis(D)
    @test eltype(av) === T
    @test av.a ≈ π
    @test av.v ≈ [0, 0, 1]

    D = angle_to_dcm(T(2π / 3), :X)
    av = dcm_to_angleaxis(D)
    @test eltype(av) === T
    @test av.a ≈ T(2π / 3)
    @test av.v ≈ [1, 0, 0]

    D = angle_to_dcm(T(-2π / 3), :X)
    av = dcm_to_angleaxis(D)
    @test eltype(av) === T
    @test av.a ≈ T(2π / 3)
    @test av.v ≈ [-1, 0, 0]

    D = angle_to_dcm(T(2π / 3), :Y)
    av = dcm_to_angleaxis(D)
    @test eltype(av) === T
    @test av.a ≈ T(2π / 3)
    @test av.v ≈ [0, 1, 0]

    D = angle_to_dcm(T(-2π / 3), :Y)
    av = dcm_to_angleaxis(D)
    @test eltype(av) === T
    @test av.a ≈ T(2π / 3)
    @test av.v ≈ [0, -1, 0]

    D = angle_to_dcm(T(2π / 3), :Z)
    av = dcm_to_angleaxis(D)
    @test eltype(av) === T
    @test av.a ≈ T(2π / 3)
    @test av.v ≈ [0, 0, 1]

    D = angle_to_dcm(T(-2π / 3), :Z)
    av = dcm_to_angleaxis(D)
    @test eltype(av) === T
    @test av.a ≈ T(2π / 3)
    @test av.v ≈ [0, 0, -1]

    D = angle_to_dcm(T(-2π / 3), :Z) * angle_to_dcm(T(-π / 2), :Y)
    av = dcm_to_angleaxis(D)
    @test eltype(av) === T
    @test av.a ≈ 2.4188584057763776
    @test av.v ≈ [0.6546536707079772, -0.37796447300922736, -0.6546536707079772]

    D = DCM(
        T(-1 / 3), T(-2 / 3), T(-1 / 3),
        T(-1 / 3), T(-1 / 3), T(+2 / 3),
        T(-2 / 3), T(+2 / 3), T(-1 / 3)
    )'
    av = dcm_to_angleaxis(D)
    @test eltype(av) === T
    @test av.a ≈ T(π)
    @test av.v ≈ [T(√3 / 3), T(-√3 / 3), T(-√3 / 3)]
end


text:
export ConstantRate,
    ExponentialRate,
    LinearRate,
    Rate

import Base: convert

abstract type Rate end

convert(::Type{Rate}, f::Float64) = ConstantRate(f)

struct ConstantRate <: Rate
    α::Float64
end

(rate::ConstantRate)(::Integer) = rate.α

struct LinearRate <: Rate
    initial_rate::Float64
    final_rate::Float64
    slope::Float64
end

(rate::LinearRate)(generation::Integer) =
    max(rate.initial_rate - generation * rate.slope, rate.final_rate)

struct ExponentialRate <: Rate
    initial_rate::Float64
    final_rate::Float64
    decay::Float64
end

(rate::ExponentialRate)(generation::Integer) =
    max(rate.initial_rate * rate.decay ^ (generation - 1), rate.final_rate)


text:
####### Routines related to Chebyshev interpolation or polynomials that are not about
#   evaluation nor interpolation set up.

# scale x ∈ [a,b] to [-1,1].
function scaletoChebyinterval(x::T, a::T, b::T)::T where T
    @assert a <= x <= b
    @assert b > a # disallow empty interval.

    Δ::T = b-a
    return (x-a)/Δ*convert(T,2) - one(T)
end

# scale y ∈ [-1,1] to [a,b].
function scalefromChebyinterval(y::T, a::T, b::T)::T where T
    @assert -one(T) <= y <= one(T)

    Δ::T = b-a
    return (y+1)/2*Δ + a
end

function Chebyshev2nd(x::T, L::Int)::T where T
    two = one(T) + one(T)

    out = zero(T)
    for l = 0:div(L,2)
        out += (-one(Int))^l*binomial(L-l,l)*(two*x)^(L-2*l)
    end

    return out
end

function Chebyshev1stall(x::T, L::Int)::Vector{T} where T
    @assert L >= 2

    # setup.
    two = one(T)+one(T)

    out = Vector{T}(undef,L+1)
    out[1] = one(T)
    out[2] = x
    for i = 3:length(out)
        out[i] = two*x*out[i-1] - out[i-2]
    end

    return out
end

function Chebyshev1st(x::T, L::Int)::T where T
    @assert L >= 0

    @assert isnumericallyin(x,-one(T),one(T))
    clamp(x,-one(T),one(T))

    return cos(convert(T,L)*acos(x))
end

# fill array with function
# equivalent to: collect( 𝝋[j](y[i]) for j = 1:n+1 for i = 1:m+1 ), but transposed.
function getunivariateChebyVmatrix(𝓧::Vector{T}, L::Int) where T
    p = length(𝓧)
    V = Matrix{T}(undef,L+1,p)

    for j = 1:p
        V[:,j] = Chebyshev1stall(𝓧[j],L)
    end

    return V
end

function Chebyshev1stnodes(Lp1::Int)::Vector{BigFloat}
    L = Lp1 - 1
    half = BigFloat("0.5")
    return collect( cos( (BigFloat(L-i)+half)*BigFloat(π)/BigFloat(L+1) ) for i = 0:L+1 )
end

function Chebyshev2ndnodes(L::Int)::Vector{BigFloat}
    return collect( cos( (L-i)*BigFloat(π)/BigFloat(L) ) for i = 0:L+1 )
end

function Chebyshev2ndnode(i::Int, L::Int)::BigFloat
    return cos( (L-i)*BigFloat(π)/BigFloat(L) )
end

# Colleague matrix for a L-degree Chebyshev polynomial of the mth kind.
function constructcolleaguematrix(c::Vector{T}, m::Int) where T
    @assert m == 1 || m == 2

    L = length(c) - 1

    while isnumericallyclose(c[L+1], zero(T), 1e-9) && L > 2
        L -= 1
    end

    # pre-compute.
    two = convert(T,2)
    half = one(T)/two

    C = zeros(T,L,L)
    C[2,1] = half

    # fill out the band.
    for l = 2:L-1
        C[l-1,l] = half
        C[l+1,l] = half
    end
    C[L-1,L] = half

    if m == 1
        C[1,2] = one(T)
    end

    # fill out bottom row.
    for l = 1:L
        C[L,l] = -c[l]/(two*c[L+1])
    end

    C[L,L-1] += half

    return C
end

# polynomial basis is Chebyshev polynomials of the m-th kind.
# coeffs are a_k for m = 1
# coeffs are c_k-1 = k*a_k for m = 2
function findChebyitproots(a::Vector{T}, m::Int)::Vector{T} where T
    if count(!isnumericallyclose(a[i],zero(T)) for i = 1:length(a)) < 2
        # case, a is the zero vector or a has less than 2 non-zero entries.
        #   The latter case means constant term and linear term. No local extrema.
        return Vector{T}(undef,0)
    end

    C = constructcolleaguematrix(a,m)

    return findChebyitproots(C)
end

function findChebyitproots(C::Matrix{T})::Vector{T} where T
    C_eig = LinearAlgebra.eigen(convert(Matrix{Float64},C))
    𝑧 = convert(Vector{T}, real.(C_eig.values))
    unique!(𝑧)
    filter!(xx->isnumericallyin(xx,-one(T),one(T)), 𝑧)
    clamp!(𝑧,-one(T),one(T))

    return 𝑧
end

function findChebyitpextrema(a::Vector{T}) where T
    L = length(a)-1

    c = collect( k*a[k+1] for k = 1:L)

    extrema_positions = findChebyitproots(c,2)

    return extrema_positions
end

# remove roots that do not have a zero derivative.
function findChebyitpextremawithchk(a::Vector{T}, tol::T = convert(T,0.2)) where T
    X = findChebyitpextrema(a)

    𝑑itp_X = evalunivariateChebyshevpolynomialderivative(a,X)

    remove_list = Vector{Int}(undef,0)
    for i = 1:length(X)
        if !isnumericallyclose(𝑑itp_X[i],zero(T),tol)
            push!(remove_list,i)
        end
    end

    deleteat!(X,remove_list)

    return X
end

# N is the number of samples to take between each reference position.
# 𝓧 is the set of reference positions.
function findChebyitpextremainterval(  𝑒::Function,
                                    L::Int,
                                    𝓟_in::Vector{T},
                                    tol_𝓧_spacing, # was defaulted to convert(T,1e-6/length(𝓟_in)
                                    tol_derivative_zero) where T

    ### get roots.
    sample_position_template = collect( Chebyshev2ndnode(k,L) for k = 0:L )
    out, d_set, X_set = getextremafromroots(𝑒, L, 𝓟_in,
                                            sample_position_template,
                                            tol_𝓧_spacing, tol_derivative_zero)

    # check each point in out to see if its a stationary point.
    chkfunc = xx->numericalchkextrema(𝑒, xx, tol_derivative_zero)
    filter!(chkfunc,out)

    sort!(out)
    unique!(out)

    return out, d_set, X_set
end

function discreteinnerprod(f::Function, g::Function, L::Int)

    ν = collect( Chebyshev2ndnode(k,L) for k = 0:L )

    term1 = 1/2 * f(ν[1]) * g(ν[1])

    term2 = BigFloat("0")
    for l = 1:L-1
        term2 += f(ν[l+1]) * g(ν[l+1])
    end

    term3 = 1/2 * f(ν[L+1]) * g(ν[L+1])

    return term1 + term2 + term3
end

max_stars_repo_name: RoyCCWang/FIRRemez

max_issues_repo_name: RoyCCWang/FIRRemez

max_forks_repo_name: RoyCCWang/FIRRemez


text:
"""
    SteepestDescent

Descent method which minimizes the objective function in the direction
of antigradient at each step.
"""
mutable struct SteepestDescent{T<:AbstractFloat,V<:AbstractVector{T}} <: CoreMethod
    x::V
    xpre::V
    g::V
    dir::V
    y::T
    α::T
end

@inline fnval(M::SteepestDescent) = M.y
@inline gradientvec(M::SteepestDescent) = M.g
@inline argumentvec(M::SteepestDescent) = M.x
@inline step_origin(M::SteepestDescent) = M.xpre

function SteepestDescent(x::AbstractVector, α::Real)
    F = float(eltype(x))
    SteepestDescent(similar(x, F), similar(x, F), similar(x, F), similar(x, F), zero(F), convert(F, α))
end

SteepestDescent(x::AbstractVector) = SteepestDescent(x, 1)

"""
`optfn!` must be the 3-arg closure that computes fdf(x + α*d) and overwrites `M`'s gradient
"""
function init!(::SteepestDescent{T}, optfn!, x0; kw...) where {T}
    optfn!(x0, zero(T), x0)
    return
end

@inline reset!(::SteepestDescent) = nothing

function reset!(M::SteepestDescent, x0, α = M.α)
    if length(M.x) != length(x0)
        foreach((M.x, M.xpre, M.g, M.dir)) do v
            resize!(v, length(x0))
        end
    end
    M.α = α
    return
end

function callfn!(M::SteepestDescent, fdf, x, α, d)
    __update_arg!(M, x, α, d)
    y, g = fdf(M.x, M.g)
    __update_grad!(M, g)
    M.y = y
    return y, g
end

@inline function __descent_dir!(M::SteepestDescent)
    M.dir .= .-M.g
    return M.dir
end

function step!(M::SteepestDescent, optfn!; constrain_step = infstep)
    M.x, M.xpre = M.xpre, M.x
    d = __descent_dir!(M)
    xpre = M.xpre
    αmax = constrain_step(xpre, d)
    α = strong_backtracking!(optfn!, xpre, d, M.y, M.g, α=M.α, αmax=αmax, β=1e-4, σ=0.1)
    return α
end

@inline function __update_arg!(M::SteepestDescent, x, α, d)
    map!(M.x, d, x) do a, b
        muladd(α, a, b)
    end
    return
end

@inline function __update_arg!(M::SteepestDescent, x)
    if x !== M.x
        copy!(M.x, x)
    end
    return
end

@inline function __update_xpre!(M::SteepestDescent, x)
    if x !== M.xpre
        copy!(M.xpre, x)
    end
    return
end

@inline function __update_grad!(M::SteepestDescent, g)
    if M.g !== g
        copy!(M.g, g)
    end
    return
end


text:
################################################################################
# Starting Area
################################################################################

mutable struct StartingArea{T}
    # This is a rectangular starting area on state space
    x_nom::VehicleState{T} # Nominal state
    x_min::VehicleState{T} # Bound on the starting state
    x_max::VehicleState{T} # Bound on the starting state
end

function StartingArea()
    return StartingArea(VehicleState(), VehicleState(), VehicleState())
end

function StartingArea(x_nom::VehicleState{T}) where {T}
    return StartingArea{T}(copy(x_nom), copy(x_nom), copy(x_nom))
end

function StartingArea(x_min::VehicleState{T}, x_max::VehicleState{T}) where {T}
    return StartingArea{T}(0.5*(x_min+x_max), x_min, x_max)
end

function randstate(start::StartingArea{T}) where {T}
    x_min = start.x_min
    x_max = start.x_max
    return VehicleState{T}(
        x_min.x + rand(T)*(x_max.x - x_min.x),
        x_min.y + rand(T)*(x_max.y - x_min.y),
        x_min.θ + rand(T)*(x_max.θ - x_min.θ),
        x_min.v + rand(T)*(x_max.v - x_min.v),
        )
end


text:
using Revise

using Pkg
Pkg.activate(".")
Pkg.instantiate()

using RadMod2D

include("./models2D.jl")
include("./plot2D.jl")

using CairoMakie
CairoMakie.activate!(type = "svg")
filetype = ".svg"
# CairoMakie.activate!(type = "png")
# filetype = ".png"

function calc_view2D_algo_n_study(;elemsize = 0.05)
    # calculating blocking with tiles
    # for different tile numbers from 1 to 200 tiles per dimension
    # choose model
    # m = model_circles_in_circle_rand_full(0.1, 1.8, elemsize)
    # m = model_circles_in_circle_rand_half(0.1, 1.8, elemsize)
    # m = model_circles_in_circle_rand_quarter(0.1, 1.8, elemsize)
    m = model_circles_in_circle_cross(0.1, 1.8, 9, elemsize)
    n_study = [1,2,3,4,5,10,15,20,25,30,35,40,45,50,55,60,65,70,75,80,85,90,95,100,120,140,160,180,200]
    times = zeros(size(n_study,1),2)
    analy_tiles = zeros(size(n_study,1),4)
    for i = 1:size(n_study,1)
        n = n_study[i]
        println("my algorithm running with ", n)
        vfmat = zeros(Float64, m.nelem, m.nelem)
        existing_vf!(m, vfmat)
        dx, dy = create_tiles(m, n)
        stats_tiles = @timed t_occ = check_tile_occupation(m, dx, dy, n)
        stats_calc = @timed blocking_vf_with_tiles!(m, vfmat, dx, dy, n, t_occ)
        # calculating_vf!(m, vfmat, normit = false)
        times[i,1] = stats_tiles.time
        times[i,2] = stats_calc.time
        t_max, t_min, t_mean, ratio = tile_occ_analysis(t_occ, printit = false)
        analy_tiles[i,1] = t_max
        analy_tiles[i,2] = t_min
        analy_tiles[i,3] = t_mean
        analy_tiles[i,4] = ratio
    end
    # save results
    header = ["n" "time_tiles" "time_calc" "t_max" "t_min" "t_mean" "ratio"]
    open("file_algo_n_study.txt", "w") do io
        writedlm(io, header)
        writedlm(io, [n_study times analy_tiles])
    end
end

function fig_n_study_mesh2D_allinone()
    m = Matrix{Model}(undef,2,2)
    m[1,1] = model_circles_in_circle_rand_full(0.1, 1.8, 0.005)
    m[1,2] = model_circles_in_circle_rand_half(0.1, 1.8, 0.005)
    m[2,1] = model_circles_in_circle_rand_quarter(0.1, 1.8, 0.005)
    m[2,2] = model_circles_in_circle_cross(0.1, 1.8, 9, 0.005)
    # 2D plot
    fig = Figure(resolution = (300, 300), font = "Arial", fontsize = 10)
    for i = 1:2
        for j = 1:2
            ax = fig[i, j] = Axis(fig)
            linewidth = 1.0
            setup_axis(ax, linewidth)
            ax.xlabel = "X in m"
            ax.ylabel = "Y in m"
            ax.aspect = DataAspect()
            colors = Vector{Any}(undef,18)
            colors[:] .= :black
            plot_model(fig, ax, m[i,j], shownvec = false, shownodes = false, showcom = false, showleg = false, colors = colors, linewidth = linewidth)
        end
    end
    # display(fig)
    save("fig_cyl_in_cyl_n_study_allinone"*filetype, fig)
end

# calc_view2D_algo_n_study(elemsize = 0.0005)
# fig_n_study_mesh2D_allinone()

# include("./test/run_bf_vs_algo.jl")


text:
const CV_CN_MAX = 512
const CV_CN_SHIFT = 3
const CV_DEPTH_MAX = (1 << CV_CN_SHIFT)

const CV_8U = 0
const CV_8S = 1
const CV_16U = 2
const CV_16S = 3
const CV_32S = 4
const CV_32F = 5
const CV_64F = 6

const CV_MAT_DEPTH_MASK = (CV_DEPTH_MAX - 1)
CV_MAT_DEPTH(flags) = ((flags) & CV_MAT_DEPTH_MASK)

CV_MAKETYPE(depth,cn) = (CV_MAT_DEPTH(depth) + (((cn)-1) << CV_CN_SHIFT))
CV_MAKE_TYPE = CV_MAKETYPE

function cpp_to_julia(mat::CxxMat)
    rets = jlopencv_core_Mat_mutable_data(mat)
    if rets[2] == CV_MAKE_TYPE(CV_8U, rets[3])
        dtype = UInt8
    elseif rets[2]==CV_MAKE_TYPE(CV_8S, rets[3])
        dtype = Int8
    elseif rets[2]==CV_MAKE_TYPE(CV_16U, rets[3])
        dtype = UInt16
    elseif rets[2]==CV_MAKE_TYPE(CV_16S, rets[3])
        dtype = Int16
    elseif rets[2]==CV_MAKE_TYPE(CV_32S, rets[3])
        dtype = Int32
    elseif rets[2]==CV_MAKE_TYPE(CV_32F, rets[3])
        dtype = Float32
    elseif rets[2]==CV_MAKE_TYPE(CV_64F, rets[3])
        dtype = Float64
    else
        error("Bad type returned from OpenCV")
    end
    steps = [rets[6]/sizeof(dtype), rets[7]/sizeof(dtype)]
    # println(steps[1]/rets[3], steps[2]/rets[3]/rets[4])
    #TODO: Implement views when steps do not result in continous memory
    arr = Base.unsafe_wrap(Array{dtype, 3}, Ptr{dtype}(rets[1].cpp_object), (rets[3], rets[4], rets[5]))

    #Preserve Mat so that array allocated by C++ isn't deallocated
    return Mat{dtype}(mat, arr)
end

function julia_to_cpp(img::InputArray)
    if typeof(img) <: CxxMat
        return img
    end
    steps = 0
    try
        steps = strides(img)
    catch
        # Copy array since array is not strided
        img = img[:, :, :]
        steps = strides(img)
    end

    if steps[1] <= steps[2] <= steps[3] && steps[1]==1
        steps_a = Array{size_t, 1}()
        ndims_a = Array{Int32, 1}()
        sz = sizeof(eltype(img))
        push!(steps_a, UInt64(steps[3]*sz))
        push!(steps_a, UInt64(steps[2]*sz))
        push!(steps_a, UInt64(steps[1]*sz))

        push!(ndims_a, Int32(size(img)[3]))
        push!(ndims_a, Int32(size(img)[2]))
        if eltype(img) == UInt8
            return CxxMat(2, pointer(ndims_a), CV_MAKE_TYPE(CV_8U, size(img)[1]), Ptr{Nothing}(pointer(img)), pointer(steps_a))
        elseif eltype(img) == UInt16
            return CxxMat(2, pointer(ndims_a), CV_MAKE_TYPE(CV_16U, size(img)[1]), Ptr{Nothing}(pointer(img)), pointer(steps_a))
        elseif eltype(img) == Int8
            return CxxMat(2, pointer(ndims_a), CV_MAKE_TYPE(CV_8S, size(img)[1]), Ptr{Nothing}(pointer(img)), pointer(steps_a))
        elseif eltype(img) == Int16
            return CxxMat(2, pointer(ndims_a), CV_MAKE_TYPE(CV_16S, size(img)[1]), Ptr{Nothing}(pointer(img)), pointer(steps_a))
        elseif eltype(img) == Int32
            return CxxMat(2, pointer(ndims_a), CV_MAKE_TYPE(CV_32S, size(img)[1]), Ptr{Nothing}(pointer(img)), pointer(steps_a))
        elseif eltype(img) == Float32
            return CxxMat(2, pointer(ndims_a), CV_MAKE_TYPE(CV_32F, size(img)[1]), Ptr{Nothing}(pointer(img)), pointer(steps_a))
        elseif eltype(img) == Float64
            return CxxMat(2, pointer(ndims_a), CV_MAKE_TYPE(CV_64F, size(img)[1]), Ptr{Nothing}(pointer(img)), pointer(steps_a))
        end
    else
        # Copy array, invalid config
        return julia_to_cpp(img[:, :, :])
    end
end

function julia_to_cpp(var::Array{T, 1}) where {T <: InputArray}
    ret = CxxWrap.StdVector{CxxMat}()
    for x in var
        push!(ret, julia_to_cpp(x))
    end
    return ret
end

function cpp_to_julia(var::CxxWrap.StdVector{T}) where {T <: CxxMat}
    ret = Array{Mat, 1}()
    for x in var
        push!(ret, cpp_to_julia(x))
    end
    return ret
end


text:
"""
    LikelihoodApproximation

An HMM likelihood approximation.
"""
abstract type LikelihoodApproximation{T} end

"""
    Truncation(hmm)

Constructs a truncation likelihood approximation object for `hmm`.

See also: [`Projection`](@ref)
"""
struct Truncation{T<:HMM} <: LikelihoodApproximation{T}
  hmm::T
end

function (ll::Truncation)(seq::Vector, obs::Vector, range)
  @unpack μ, σ, Σᵨ, Σₐ, W, D = ll.hmm
  H  = view(W,  range, range) * view(D, range, range)
  Sᵨ = view(Σᵨ, range, range)
  Sₐ = view(Σₐ, range, range)
  d  = view(obs, range)
  loglikelihood(seq, d, μ, σ, Sₐ, H, Sᵨ) #/ length(range)
end

"""
    Projection(hmm)

Constructs a projection likelihood approximation object for `hmm`.

See also: [`Truncation`](@ref)
"""
struct Projection{T<:HMM,V<:AbstractVector,M<:AbstractMatrix} <: LikelihoodApproximation{T}
  hmm::T
  μ₁::V
  μ₂::V
  Σ₁₁::M
  Σ₂₂::M
  Σ₂₁::M
end

function Projection(hmm::T) where {T<:HMM}
  @unpack W, D, Σₐ = hmm
  H = W * D

  # response approximation
  μ₁, Σ₁₁ = project(hmm)

  # observation moments
  μ₂ = H * μ₁
  Σ₂₂ = Symmetric(H * Σ₁₁ * H' + Σₐ)

  # observation-response covariance
  Σ₂₁ = H * Σ₁₁

  Projection{T,typeof(μ₁),typeof(Σ₁₁)}(hmm, μ₁, μ₂, Σ₁₁, Σ₂₂, Σ₂₁)
end

function (ll::Projection)(seq::Vector, obs::Vector, range)
  @unpack μ₁, μ₂, Σ₁₁, Σ₂₂, Σ₂₁ = ll
  @unpack μ, σ, Σᵨ = ll.hmm

  # extract slices for conditioning
  Σ₂₁ₖ = view(Σ₂₁,     :, range)
  Σ₁₁ₖ = view(Σ₁₁, range, range)
  Σᵨₖ  = view(Σᵨ,  range, range)
  μ₁ₖ  = view(μ₁,  range)

  # values to condition on
  μₛ  = μ[seq]
  S  = Diagonal(σ[seq])
  Σₛ = S * Σᵨₖ * S

  # approximate moments
  C = Σ₂₁ₖ * inv(Σ₁₁ₖ)
  μₚ = μ₂ + C * (μₛ - μ₁ₖ)
  Σₚ = Σ₂₂ - C * Σ₂₁ₖ' + C * Σₛ * C'

  logpdf(MvNormal(μₚ, Symmetric(Σₚ)), obs)
end

"""
    project(hmm)

Computes Gaussian mixture approximation to `p(m)`. Returns mean and covariance of an `n`-dimensional marginal Gaussian.
"""
function project(hmm::ConvolvedHMM)
  @unpack μ, σ, logP, Σᵨ, logπ₀ = hmm
  n = size(Σᵨ, 1)
  nk = length(μ)

  # transition and stationary distributions
  π₀ = exp.(logπ₀)
  P = exp.(logP)

  # point-wise projected mean
  μₘ = μ' * π₀
  μₚ = μₘ * ones(n)

  # vector of covariances c(|s-t|)
  c = Vector{eltype(σ)}(undef, n)

  # diagonal elements for lag 0
  c[1] = sum((σ.^2 .+ (μ .- μₘ).^2) .* π₀)

  # initialize n-step transition matrix
  Pⁿ = P

  # lags 1 to n-1
  for h in 1:n-1
    c[h+1] = 0
    for k₁ in 1:nk, k₂ in 1:nk
      c[h+1] += (σ[k₁] * σ[k₂] * Σᵨ[1,h+1] + (μ[k₁] - μₘ) * (μ[k₂] - μₘ)) * π₀[k₁] * Pⁿ[k₁,k₂]
    end

    # increment lag
    Pⁿ *= P
  end

  # build covariance matrix
  Σₚ = Matrix{eltype(σ)}(undef, n, n)
  for i in 1:n, j in 1:n
    Σₚ[i,j] = c[abs(j-i)+1]
  end

  μₚ, Σₚ
end


text:
@testset "durations" begin
    cvrp, _, _ = CVRPLIB.readCVRPLIB("CMT6")
    ap = AlgorithmParameters(timeLimit=3)
    result = solve_cvrp(cvrp, ap, verbose=true, use_dist_mtx=false, round=false)
    @show result.cost
    @test result.cost ≈ 555.43 atol=1e-2

    cvrp, _, _ = CVRPLIB.readCVRPLIB("CMT7")
    service_times = cvrp.service_time .* ones(cvrp.dimension)
    service_times[1] = 0.0
    ap = AlgorithmParameters(timeLimit=3)

    result = solve_cvrp(
        cvrp.coordinates[:, 1],
        cvrp.coordinates[:, 2],
        cvrp.demand,
        cvrp.capacity,
        ap;
        verbose = true,
        round = false,
        duration_limit = cvrp.distance,
        service_times = service_times
    )
    @show result.cost
    @test result.cost ≈ 909.68 atol=1e-2

end

@testset "P-n19-k2 by coordinates" begin
    cvrp, _, _ = CVRPLIB.readCVRPLIB("P-n19-k2")
    ap = AlgorithmParameters(timeLimit=1.8)
    result = solve_cvrp(cvrp, ap; verbose=true)
    @test result.cost <= 212 * 1.01
end

@testset "P-n19-k2 by dist_mtx with coordinates" begin
    cvrp, _, _ = CVRPLIB.readCVRPLIB("P-n19-k2")
    ap = AlgorithmParameters(timeLimit=2.1)
    result = solve_cvrp(cvrp, ap; verbose=true, use_dist_mtx=true)
    @test result.cost <= 212 * 1.01
end

@testset "P-n19-k2 by dist_mtx without coordinates" begin
    cvrp, _, _ = CVRPLIB.readCVRPLIB("P-n19-k2")
    cvrp.coordinates = zeros(size(cvrp.coordinates))
    ap = AlgorithmParameters(timeLimit=1.23)
    result = solve_cvrp(cvrp, ap; verbose=true, use_dist_mtx=true)
    @test result.cost <= 212 * 1.01
end

@testset "X-n101-k25" begin
    cvrp, _, _ = CVRPLIB.readCVRPLIB("X-n101-k25")

    ap = AlgorithmParameters(timeLimit=3)
    result = solve_cvrp(cvrp, ap; verbose=true)

    @show result.cost
    @show result.time
    visited_customers = result.routes

    @show typeof(visited_customers), size(visited_customers)

    report_visited_customers = reporting(visited_customers)
    for i in eachindex(report_visited_customers)
        @test report_visited_customers[i] == visited_customers[i] .- 1
    end
end

@testset "A-n32-k5.vrp" begin
    cvrp = CVRPLIB.readCVRP("A-n32-k5.vrp")
    ap = AlgorithmParameters(timeLimit=2.3)
    result = solve_cvrp(cvrp, ap; verbose=true)
    @test result.cost <= 784 * 1.01
end

@testset "x, y, dist_mtx CVRP" begin
    cvrp, _, _ = CVRPLIB.readCVRPLIB("A-n32-k5")
    x = cvrp.coordinates[:, 1]
    y = cvrp.coordinates[:, 2]
    dist_mtx = cvrp.weights
    demand = cvrp.demand
    capacity = cvrp.capacity
    n_vehicles = 5

    ap = AlgorithmParameters(timeLimit=1.8)

    result1 = solve_cvrp(x, y, demand, cvrp.capacity, ap;
                            n_vehicles=n_vehicles, verbose=true)
    result2 = solve_cvrp(dist_mtx, demand, cvrp.capacity, ap; verbose=true)
    result3 = solve_cvrp(dist_mtx, demand, cvrp.capacity, ap; x_coordinates=x, y_coordinates=y, verbose=true, n_vehicles=n_vehicles)

    @test result1.cost == result2.cost == result3.cost
end


text:
import TPT
using Base.Test

using DataFrames
using Plots; pyplot()

println("--- TM HHTB ---")

# Load elemental parameters
p = readtable(joinpath("data", "parameters", "tm_hhtb.csv"))
N = size(p, 1) # number of elements

#
# Prepare a directory to store the results
#
!isdir("results") && mkdir("results")

resdir = joinpath("results", "tm_hhtb")
!isdir(resdir) && mkdir(resdir)

# Convertion of parameters (Hartree units)
X = p[:X]
Vc = p[:Vc] / (5.292^3 * 1e-3) # a.u
r₀ = (Vc * 3/4π) .^ (1/3) # a.u.
Vl = p[:Vl] / (5.292^3 * 1e-3) # a.u.
ρ = 1 / Vl # a.u.
σ = p[:σ] / 5.291e-1 # Å -> a.u.
Ns = p[:Ns]
Nd = p[:Nd]
Ed = p[:Ed] / 27.21 # eV -> a.u.
Wd = p[:Wd] / 27.21 # eV -> a.u.
Rc = p[:Rc] / 5.291e-1 # Å -> a.u.
rc = p[:rc] # Å -> a.u.
a = p[:a] # Å -> a.u.

# some parameters
rmin, rmax = (2, 12)
umin, umax = (-0.04, 0.04)

#
# Puremetals
#
E_band = zeros(N)
E_bond = zeros(N)
E_rep = zeros(N)
E_nfe = zeros(N)

for i in 2:8
  ahs = TPT.AHS(σ = σ[i], ρ = ρ[i])
  # pp = TPT.Ashcroft(Ns[i], Rc[i])
  pp = TPT.BretonnetSilbert(Ns[i], rc[i], a[i])
  nfe = TPT.NFE(ahs, pp, :IU)
  hhtb = TPT.HHTB(Nd[i], Ed[i], Wd[i], r₀[i])
  nfetb = TPT.NFETB(nfe, hhtb)

  u_nfe = TPT.pairpotential(nfe)[1,1]
  u_tb = TPT.pairpotential(hhtb)[1,1]
  u_tot = TPT.pairpotential(nfetb)[1,1]
  plot([u_nfe, u_tb, u_tot], 2, 12, ylims=(umin, umax), xl="r (a.u.)", yl="u (a.u.)", labels=["NFE" "TB" "total"])
  png(joinpath(resdir, "$(i)_$(X[i])_u"))

  E_band[i] = TPT.bandenergy(hhtb)
  E_bond[i] = TPT.internal(ahs, hhtb, TPT.pairpotential_bond(hhtb))
  E_rep[i] = TPT.internal(ahs, hhtb, TPT.pairpotential_rep(hhtb))
  E_nfe[i] = TPT.internal(ahs, hhtb, TPT.pairpotential(nfe))
end

plot(X, [E_bond, E_rep, E_nfe, E_bond + E_rep + E_nfe], labels=["bond" "rep" "nfe" "total"], m=:o, xl="X", yl="E (a.u. / atom)", ylims=:auto)
png(joinpath(resdir, "E"))

#
# Fe-based binary alloys
#
E_alloy = zeros(N)
ΔE_bond = zeros(N)
ΔE_rep = zeros(N)
ΔE_nfe = zeros(N)
x = [0.5, 0.5]
for i in 6, j in 2:8
  i == j && continue

  ρ_alloy = (ρ[i] + ρ[j]) / 2
  ahs = TPT.AHS(ρ = ρ_alloy, σ = [ σ[i], σ[j] ], c = [0.5, 0.5])
  # pp = TPT.Ashcroft([ Ns[i], Ns[j] ], [ Rc[i], Rc[j] ])
  pp = TPT.BretonnetSilbert([ Ns[i], Ns[j] ], [ rc[i], rc[j] ], [ a[i], a[j] ])
  nfe = TPT.NFE(ahs, pp, :IU)
  hhtb = TPT.HHTB(x, [ Nd[i], Nd[j] ], [ Ed[i], Ed[j] ], [ Wd[i], Wd[j] ], [ r₀[i], r₀[j] ])
  nfetb = TPT.NFETB(nfe, hhtb)

  Emin, Emax = TPT.cutoffenergy(hhtb)
  Ef = TPT.fermienergy(hhtb)
  D = TPT.partialdensityofstate(hhtb)
  D_total = TPT.totaldensityofstate(hhtb)

  plot([ D[1], D[2], D_total ], Emin, Emax, ylims=:auto, labels=[ X[i] X[j] "total" ])
  vline!([Ef], label="E_F")
  xlabel!("E (a.u.)")
  ylabel!("DOS (states / a.u. atom)")
  png(joinpath(resdir, "$(i)-$(j)_$(X[i])-$(X[j])_D"))

  u = TPT.pairpotential(nfetb)
  plot([ u[1,1], u[1,2], u[2,2] ], rmin, rmax, labels=["1-1" "1-2" "2-2"], ylims=(umin, umax), xl="r (a.u.)", yl="u (a.u.)")
  png(joinpath(resdir, "$(i)-$(j)_$(X[i])-$(X[j])_u"))

  E_rep_alloy = TPT.internal(ahs, hhtb, TPT.pairpotential_rep(hhtb))
  E_bond_alloy = TPT.internal(ahs, hhtb, TPT.pairpotential_bond(hhtb))
  E_nfe_alloy = TPT.internal(ahs, hhtb, TPT.pairpotential(nfe))
  E_alloy[j] = E_rep_alloy + E_bond_alloy + E_nfe_alloy

  ΔE_bond[j] = E_bond_alloy - x[1]*E_bond[i] - x[2]*E_bond[j]
  ΔE_rep[j] = E_rep_alloy - x[1]*E_rep[i] - x[2]*E_rep[j]
  ΔE_nfe[j] = E_nfe_alloy - x[1]*E_nfe[i] - x[2]*E_nfe[j]
end

plot(X, E_alloy, label="total", m=:o, xl="Fe-X", yl="E (a.u. / atom)", ylims=:auto)
png(joinpath(resdir, "E_alloy"))

plot(X, [ΔE_bond, ΔE_rep, ΔE_nfe, ΔE_bond + ΔE_rep + ΔE_nfe], labels=["bond" "rep" "nfe" "total"], m=:o, xl="Fe-X", yl="Emix (a.u. / atom)", ylims=:auto)
png(joinpath(resdir, "Emix"))

@testset "TM HHTB" begin
end


text:
"""
PlottingCallback

A callback function that plots a number of convergence metrics using [ConvergencePlots.jl](https://github.com/mohamed82008/ConvergencePlots.jl). The following are the fields of the struct:
- `plot`: an instance of `ConvergencePlots.ConvergencePlot`.
- `n`: the number of history points to keep track of.
- `iter`: the number of times the callback was called.
"""
mutable struct PlottingCallback <: Function
    plot
    n::Int
    iter::Int
end
PlottingCallback(n = 1000) = PlottingCallback(nothing, n, 0)

(callback::PlottingCallback)(solution) = nothing

"""
    LazyPlottingCallback

A callback function that plots a number of convergence metrics using [ConvergencePlots.jl](https://github.com/mohamed82008/ConvergencePlots.jl) only when `show` is called. The following are the fields of the struct:
- `plot`: an instance of `ConvergencePlots.ConvergencePlot`.
- `n`: the number of history points to keep track of.
- `iter`: the number of times the callback was called.
"""
mutable struct LazyPlottingCallback <: Function
    plot
    n::Int
    iter::Int
    show_plot::Bool
    save_plot::Union{Nothing, String}
end
LazyPlottingCallback(n = 10^10; show_plot = true, save_plot = nothing) = LazyPlottingCallback(nothing, n, 0, show_plot, save_plot)

(callback::LazyPlottingCallback)(solution::Solution; update = false) = nothing

struct NoCallback <: Function end
(::NoCallback)(args...; kwargs...) = nothing

@require ConvergencePlots = "e9c76d45-d282-4a42-92e1-2fa9292d7569" @eval begin
    using .ConvergencePlots

    """
    (callback::PlottingCallback)(solution::Solution)

    Updates the convergence plots using the current solution.
    """
    function (callback::PlottingCallback)(solution::Solution)
        @unpack convstate = solution
        if callback.iter == 0
            callback.plot = ConvergencePlot(
                callback.n,
                names = ["KKT residual", "|Δf|", "|Δx|_∞", "Infeasibility"],
                #names = ["KKT residual", "IPOPT residual", "|Δf|", "|Δx|_∞", "Infeasibility"],
                options = Dict(
                    "KKT residual" => (color = "black",),
                    #"IPOPT residual" => (color = "blue",),
                    "|Δf|" => (color = "green",),
                    "|Δx|_∞" => (color = "violet",),
                    "Infeasibility" => (color = "red",),
                ),
            )
        else
            @unpack plot = callback
            addpoint!(
                plot,
                Dict(
                    "KKT residual" => convstate.kkt_residual,
                    #"IPOPT residual" => convstate.ipopt_residual,
                    "|Δf|" => convstate.Δf,
                    "|Δx|_∞" => convstate.Δx,
                    "Infeasibility" => convstate.infeas,
                ),
            )
        end
        callback.iter += 1
    end

    """
    (callback::LazyPlottingCallback)(solution::Solution)

    Updates the convergence plots using the current solution.
    """
    function (callback::LazyPlottingCallback)(solution::Solution; update = false)
        @unpack convstate = solution
        if callback.iter == 0
            callback.plot = ConvergencePlot(
                callback.n,
                names = ["KKT residual", "|Δf|", "|Δx|_∞", "Infeasibility"],
                #names = ["KKT residual", "IPOPT residual", "|Δf|", "|Δx|_∞", "Infeasibility"],
                options = Dict(
                    "KKT residual" => (color = "black",),
                    #"IPOPT residual" => (color = "blue",),
                    "|Δf|" => (color = "green",),
                    "|Δx|_∞" => (color = "violet",),
                    "Infeasibility" => (color = "red",),
                ),
                show = update && callback.show_plot,
            )
        else
            @unpack plot = callback
            addpoint!(
                plot,
                Dict(
                    "KKT residual" => convstate.kkt_residual,
                    #"IPOPT residual" => convstate.ipopt_residual,
                    "|Δf|" => convstate.Δf,
                    "|Δx|_∞" => convstate.Δx,
                    "Infeasibility" => convstate.infeas,
                ),
                show = callback.show_plot,
                update = update,
                filename = callback.save_plot,
            )
        end
        callback.iter += 1
    end
end