text:
# A novel disease mechanism leading to the expression of a disallowed gene in the pancreatic beta-cell identified by non-coding, regulatory variants controlling HK1
## https://doi.org/10.1101/2021.12.03.21267240
## Notebook to generate all genomics figure panels
# Note repository contains source data for all gene expression and quantifications of ChIP-seq and ATAC-seq data at particular locus.
# Repository does not contain alignment files, but code is supplied to regenerate quantifications if suitable alignment files are present.

include(joinpath("..", "src", "project.jl"));

ENV["GKSwstype"] = "100" ### this needed for GR when running on server without a display terminal, comment out should you want GR to display figures rather than saving

mutations, region, distinctmutations, refseq, altseqs, hk1coords, genecoords, datasets, datagroup, FM, hicdata, gtex, betadiffstage, betadiffsc, islet_scrnaseq, betadifflate, isletcooler = loaddata(false, false);

# Motif Analysis

motifmeta, motifs = loadtransfac();
motscan = motifscannall(refseq, altseqs, distinctmutations, motifs)
motscan.MotifFamily = motiffamily.(motscan.MotifName);

motmatch_tier1 = filtermatches(motscan, minscore=0.6)
motmatch_tier2 = filtermatches(motscan, minscore=0.45, maxscore=0.6);

motsel_tier1 = motif_rank1(motmatch_tier1) ## table of rank 1 matches for plotting
motsel_tier2 = motif_rank1(motmatch_tier2) ## table of rank 1 matches for plotting

motifmatches = motifmatchtable(motmatch_tier1, motmatch_tier2, region); ## final table

# Figure 3

fig_3a = gp_hk1_gck(betadiffsc)
fig_3b = calcisletpanel("genelocus", datagroup, hk1coords, FM, genecoords)
fig_3c = calcatacdiff_subpanel("criticalregion_4k", datagroup, hk1coords, FM, genecoords);
fig_3d = coolerpanel(hk1coords["widecontactview"]..., isletcooler, hicdata, hk1coords)

fig_3e_left = disruptmotifscatter(@subset(motmatch_tier1, :ID .!= "71108660_G_A", :Rank .== 1))
fig_3e_right = disruptedmotif_seqlogo(region, @subset(mutations, :start .!= 71108660), @subset(motsel_tier1, :ID .!= "71108660_G_A"))
fig_3e = plot(fig_3e_left, fig_3e_right, layout=@layout([a b{0.75w}]), size=(800, 200), left_margin=5mm, top_margin=5mm)

savefigure("fig_3a", p=fig_3a)
savefigure("fig_3b", p=fig_3b)
savefigure("fig_3c", p=fig_3c)
savefigure("fig_3d", p=fig_3d)
savefigure("fig_3e", p=fig_3e)

# Supplementary Figure 3

fig_s3_a = gtexboxplot(gtex.genetpm, rev=true, orderfield=:MedianOrder)
fig_s3_b = gtexisoformheatmap(gtex.tpm, num_iso=3, c=cgrad(:blues, rev=false))
fig_s3_c = plot_hk1_isoforms(gtex.tpm, genecoords; num_iso=3)

fig_s3_d_l = barbetadiff("HK1", betadifflate)
fig_s3_d_r = barbetadiff("GCK", betadifflate)
fig_s3_e   = scboxplot("HK1", islet_scrnaseq)

fig_s3_de = plot(fig_s3_d_l, fig_s3_d_r, fig_s3_e, layout=(1, 3), size=(1050, 300), fontfamily="helvetica", ylims=(-0.0065, 0.225))

savefigure("fig_s3_a", p=fig_s3_a)
savefigure("fig_s3_b", p=fig_s3_b)
savefigure("fig_s3_c", p=fig_s3_c)
savefigure("fig_s3_de", p=fig_s3_de)

# Supplemental Figure 5

fig_s5_a = hk1_loop_summary(hicdata, hk1coords)
fig_s5_b = coolerpanel(hk1coords["widecontactview"]..., isletcooler, hicdata, hk1coords);
fig_s5_c = coolerpanel(hk1coords["ins_contactview"]..., isletcooler, hicdata, hk1coords, plotroi=false);
fig_s5_d = coolerpanel(hk1coords["sst_contactview"]..., isletcooler, hicdata, hk1coords, plotroi=false);
fig_s5_e = coolerpanel(hk1coords["gcg_contactview"]..., isletcooler, hicdata, hk1coords, plotroi=false);
savefigure("fig_s5_a", p=fig_s5_a)
savefigure("fig_s5_b", p=fig_s5_b)
savefigure("fig_s5_c", p=fig_s5_c)
savefigure("fig_s5_d", p=fig_s5_d)
savefigure("fig_s5_e", p=fig_s5_e)

# Supplemental Figure 6

fig_s6_a = calc_snatachpanel("genelocus", datagroup, hk1coords, FM, genecoords)
fig_s6_b = calc_atacdiff_supp("genelocus", datagroup, hk1coords, FM, genecoords)
fig_s6_c = calc_progentior_panel("genelocus", datagroup, hk1coords, FM, genecoords)
fig_s6_d = calc_histonepanel("roi_200k", datagroup, hk1coords, FM, genecoords)
savefigure("fig_s6_a", p=fig_s6_a)
savefigure("fig_s6_b", p=fig_s6_b)
savefigure("fig_s6_c", p=fig_s6_c)
savefigure("fig_s6_d", p=fig_s6_d)

# Supplemental Figure 7

motifgenes = @subset(islet_scrnaseq.normcounts, occursin.(r"NFAT|NKX2-2|FOXA2", :Gene), :Gene .∉ Ref(["NFATC2IP"])).Gene
sort!(motifgenes, by= x -> ifelse(occursin("NFAT", x), (1, last(x)), (2, x)))

fig_s7_a = scboxplotgenes(motifgenes, islet_scrnaseq)

motifgenesgp = @subset(betadiffsc.data, occursin.(r"NFAT[C]*[0-9]$|NKX2.2|FOXA2", :Gene)).Gene
sort!(motifgenesgp, by= x -> ifelse(occursin("NFAT", x), (1, last(x)), (2, x)))
fig_s7_b_top    = plotgprgenes(filter(g ->  occursin("NFAT", g), motifgenesgp), betadiffsc)
fig_s7_b_bottom = plotgprgenes(filter(g -> !occursin("NFAT", g), motifgenesgp), betadiffsc)

phs_top    = [barbetadiff(g, betadifflate, title=g) for g in filter(g ->  occursin("NFAT", g), motifgenes)]
phs_bottom = [barbetadiff(g, betadifflate, title=g) for g in filter(g -> !occursin("NFAT", g), motifgenes)]
plot!.(phs_top[2:end], leg=false)
fig_s7_c_top    = plot(phs_top..., layout=(1, length(phs_top)), size=(1200, 200), link=:y)
fig_s7_c_bottom = plot(phs_bottom..., layout=(1, length(phs_bottom)), size=(550, 200), link=:y, leg=:outertopright)

savefigure("fig_s7_a", p=fig_s7_a)
savefigure("fig_s7_b_top", p=fig_s7_b_top)
savefigure("fig_s7_b_bottom", p=fig_s7_b_bottom)
savefigure("fig_s7_c_top", p=fig_s7_c_top)
savefigure("fig_s7_c_bottom", p=fig_s7_c_bottom)

# Supplemental Figure 8

pat11_mots = @subset(innerjoin(motmatch_tier1, mutations[!, [:ID, :PatientGroup]], on=:ID), :PatientGroup .== 11, :Rank .== 1)
fig_s8a_left = disruptmotifscatter(pat11_mots)
fig_s8a_right = disruptedmotif_seqlogo(region, @subset(mutations, :PatientGroup .== 11), pat11_mots)
plot!(first.(@subset(mutations, :PatientGroup .== 11).RelReg) .+ [1, -1], [1, 1].*-0.2, c=:red, ls=:dash, lab="")
fig_s8a = plot(fig_s8a_left, fig_s8a_right, layout=@layout([a b{0.75w}]), size=(800, 200), left_margin=5mm, top_margin=5mm)
fig_s8b = scboxplotgenes(filter(g -> occursin(Regex(string("^", first(@subset(pat11_mots, :ID .== "71108660_G_A").MotifFamily))), g), islet_scrnaseq.normcounts.Gene), islet_scrnaseq, size=(500, 315))

matchtier2 = @subset(motmatch_tier2, :MotifFamily .∉ Ref(Set(motsel_tier1.MotifFamily)), :Rank .== 1)
motsel2 = @subset(motsel_tier2, :MotifFamily .∉ Ref(Set(motsel_tier1.MotifFamily)))

fig_s8c_left  = disruptmotifscatter([@subset(motmatch_tier1, :ID .!= "71108660_G_A", :Rank .== 1) ; matchtier2], xlims=(0.2, 0.8), ylims=(-.2, 0.8))
fig_s8c_right = disruptedmotif_seqlogo(region, @subset(mutations, :start .!= 71108660), motsel2)
fig_s8c = plot(fig_s8c_left, fig_s8c_right, layout=@layout([a b{0.75w}]), size=(800, 200), left_margin=5mm, top_margin=5mm)

fig_s8d = scboxplotgenes(filter(g -> occursin(motsel2.MotifFamily[1], g), islet_scrnaseq.normcounts.Gene)[1:2], islet_scrnaseq, size=(500, 315))
fig_s8e = scboxplotgenes(filter(g -> occursin(Regex(string("^", motsel2.MotifFamily[2], "[0-9]\$")), g), islet_scrnaseq.normcounts.Gene), islet_scrnaseq, size=(1250, 315))

savefigure("fig_s8a", p=fig_s8a)
savefigure("fig_s8b", p=fig_s8b)
savefigure("fig_s8c", p=fig_s8c)
savefigure("fig_s8d", p=fig_s8d)
savefigure("fig_s8e", p=fig_s8e)


text:
# data_handling.jl
# New way to load data into memory to plug into Julia
# Also includes training/testing set separation and evaluation of method results

MANHATTAN_NODES = 6134

function createReducedJLDs()
	"""
	Converts reduced CSV files into JLD format with a DataFrame
	Takes in no arguments and returns nothing (designed to be run once)
	"""
	for j = 11:12
		println(j)
		f = open("../../Full_data/reduced_trip_data_$j.csv")
		df2 = DataFrame(pTime = DateTime[], dTime=DateTime[], pX=Float64[], pY=Float64[], dX=Float64[], dY=Float64[],)
		for (i, ln) in enumerate(eachline(f))
			if i == 1
				continue
			end
			if i%1000000 == 0
				println(i)
			end
			s = split(ln,",")
			p = DateTime(s[1][2:end-1], "y-m-d H:M:S")
			d = DateTime(s[2][2:end-1], "y-m-d H:M:S")
			push!(df2, [p,d, float(s[3]), float(s[4]), float(s[5]), float(strip(s[6], '\n'))])
		end
	    close(f)
	    save("../Ride_data/new_trip_data_$j.jld","df",df2)
	end
	return nothing
end

function loadRides(;year::Int=2013, startTime::Int=12, endTime::Int=14, weekdays::Bool = true, startMonth::Int=1, endMonth::Int=1, loadFromCache = true)
	"""
	Load appropriate set of rides, between startTime and endTime, startMonth and endMonth, on weekdays or weekends. Supports caching to save time.
	"""
	# Initialize output file name
	fileName = "../Ride_data/$(year)/rides_$(startTime)$(endTime)_m$(startMonth)$(endMonth)"
	if weekdays
		fileName = string(fileName, "_wd.jld")
	else
		fileName = string(fileName, "_we.jld")
	end
	# Load from file if available
	if isfile(fileName)
		df2 = load(fileName, "df")
		return df2
	end
	# Create dataframe
	df2 = DataFrame()
	# Load data
	for i = startMonth:endMonth
		df = load("../Ride_data/$(year)/reduced_jld/trip_data_$i.jld", "df")
		if weekdays
			df2 = vcat(df2,df[((dayofweek(df[:,:pTime]) .<= 5) & (startTime .<= hour(df[:,:pTime]) .<= endTime) & (dayofweek(df[:,:dTime]) .<= 5) & (startTime .<= hour(df[:,:dTime]) .<= endTime)),:])
		else
			df2 = vcat(df2,df[((dayofweek(df[:,:pTime]) .> 5) & (startTime .<= hour(df[:,:pTime]) .<= endTime) & (dayofweek(df[:,:dTime]) .> 5) & (startTime .<= hour(df[:,:dTime]) .<= endTime)),:])
		end
	end
	# Filter out bad rides
	function goodRide(i::Int)
	    return (df2[i,:dTime] - df2[i,:pTime]).value > 60*1000 && (abs(df2[i,:pX] - df2[i,:dX]) + abs(df2[i,:pY] - df2[i,:dY])) >= 140
	end
	df2 = df2[Bool[goodRide(i) for i = 1:nrow(df2)],:]
	save(fileName, "df", df2)
	return df2
end

function splitTrainingAndTestingSets(df::DataFrame, inverseTestingSetFraction::Int; split_type::AbstractString = "random")
	"""
	Takes in DataFrame with rides. Returns two DataFrames, one for training and the other for testing
	split_type can be random or localized
	"""
	if split_type == "random"
		shuf = shuffle(collect(1:nrow(df)))
		testDf = df[shuf[1:round(Int,length(shuf)/inverseTestingSetFraction)],:]
		trainDf= df[shuf[(round(Int,length(shuf)/inverseTestingSetFraction) + 1):end],:]
	else
		p = (1000, -1200, -1000, -3200)
		rad = 2000
		trainDf = df[(df[:,:pX]- p[1]).^2 + (df[:,:pY] - p[2]).^2 + (df[:,:dX] - p[3]).^2 + (df[:,:dY] - p[4]).^2 .> rad^2, :]
		testDf = df[(df[:,:pX] -p[1]).^2 + (df[:,:pY] - p[2]).^2 + (df[:,:dX] - p[3]).^2 + (df[:,:dY] - p[4]).^2 .<= rad^2, :]
	end
	return trainDf, testDf
end

function loadTrainingAndTestingSets(;year::Int=2013, startTime::Int=12, endTime::Int=14, weekdays::Bool = true, startMonth::Int=1, endMonth::Int=1, loadFromCache::Bool = true, split_training_type::AbstractString = "random")
	"""
	Loads dataframe for rides, splits off 20% for permanent testing set (NEVER to be touched), then splits training set into 80% tr-learning and 20% tr-testing sets
	"""
	fileName = "../Ride_data/$(year)/rides_$(startTime)$(endTime)_m$(startMonth)$(endMonth)"
	if weekdays
		fileName = string(fileName, "_wd")
	else
		fileName = string(fileName, "_we")
	end
	if split_training_type == "random"
		testingFileName = string(fileName, "_testing.jld")
		trainLearnFileName = string(fileName, "_tr_learn.jld")
		trainTestFileName = string(fileName, "_tr_test.jld")
	elseif split_training_type == "localized"
		testingFileName = string(fileName, "_loc_testing.jld")
		trainLearnFileName = string(fileName, "_loc_tr_learn.jld")
		trainTestFileName = string(fileName, "_loc_tr_test.jld")
	end
	if isfile(testingFileName) && isfile(trainLearnFileName) && isfile(trainTestFileName) && loadFromCache
		testDf = load(testingFileName, "df")
		trainLearnDf = load(trainLearnFileName, "df")
		trainTestDf = load(trainTestFileName, "df")
	else
		df = loadRides(startTime=startTime, endTime=endTime, weekdays=weekdays,startMonth=startMonth, endMonth=endMonth, loadFromCache=true)
		trainDf, testDf = splitTrainingAndTestingSets(df, 5, split_type = "random")
		trainLearnDf, trainTestDf = splitTrainingAndTestingSets(trainDf, 0, split_type = split_training_type)
		save(testingFileName, "df", testDf)
		save(trainLearnFileName, "df", trainLearnDf)
		save(trainTestFileName, "df", trainTestDf)
	end
	return trainLearnDf, trainTestDf, testDf
end

function mapRidesToNodes(trainDf::DataFrame, nodePositions::Array{Coordinates, 1}, method::AbstractString; radius::Int = 140)
	"""
	Given the learning set of rides, an array of the node locations and the appropriate method, converts dataframe full of rides into travelTime matrix
	"""
	highwayNodes = load("Cities/Saved/highwayNodes.jld", "highwayNodes")
	travelTimes = zeros(length(nodePositions), length(nodePositions))
	numRides = zeros(Int, (length(nodePositions), length(nodePositions)))
	if method == "nearest"
		# Put nodes in KDTree
		nodeTree, nodePairs = buildNodeKDTree(nodePositions)
		# Map rides to nodes, without overlap
		# Loop through rides
		for i = 1:nrow(trainDf)
			# Choose method
			if i % 50000 == 0
				println(i)
			end
			# Nearest neighbor search
			index, dist = knn(nodeTree, [trainDf[i,:pX], trainDf[i,:pY], trainDf[i,:dX], trainDf[i,:dY]], 1)
			index = index[1]
			newTime = (trainDf[i,:dTime] - trainDf[i,:pTime]).value/1000
			if !(nodePairs[index,1] in highwayNodes) && !(nodePairs[index,2] in highwayNodes)
				travelTimes[nodePairs[index,1], nodePairs[index,2]] = (travelTimes[nodePairs[index,1], nodePairs[index,2]] * numRides[nodePairs[index,1], nodePairs[index,2]] + newTime)/(numRides[nodePairs[index,1],nodePairs[index,2]]+1)
				numRides[nodePairs[index,1], nodePairs[index,2]] += 1
			end
		end
	elseif method == "average"
		# Put rides in KDTree
		treeL = zeros(4,nrow(trainDf))
		for i in 1:4
			treeL[i,:] = trainDf[:,(i+2)]
		end
		println("--Constructing KDTree")
		@time treeL = KDTree(treeL)
		for i = 1:length(nodePositions), j=1:length(nodePositions)
			if !(i in highwayNodes) && !(j in highwayNodes) && i != j
				vec = [nodePositions[i].x, nodePositions[i].y, nodePositions[j].x, nodePositions[j].y]
				index = inball(treeL, vec, float(radius))
				for idx in index
					newTime = (trainDf[idx,:dTime] - trainDf[idx,:pTime]).value/1000
					travelTimes[i,j] = (travelTimes[i,j] * numRides[i,j] + newTime)/(numRides[i,j] + 1)
					numRides[i,j] += 1
				end
			end
		end
	end
	return travelTimes, numRides
end

function loadInputTravelTimes(nodePositions::Array{Coordinates}, method::AbstractString; year::Int=2013, startTime::Int=12, endTime::Int=14, weekdays::Bool = true, startMonth::Int=1, endMonth::Int=1, loadFromCache = true, radius::Int=140, split_training_type::AbstractString = "random")
	"""
	Top level function, called before fast_LP to load travel times for data.
	"""
	println("**** Loading training and testing sets ****")
	# Select file name
	fileName = "../Ride_data/$(year)/Input_travel_times/travel_times_$(startTime)$(endTime)_m$(startMonth)$(endMonth)_r$(radius)_$(method)"
	if weekdays
		fileName = string(fileName, "_wd")
	else
		fileName = string(fileName, "_we")
	end
	if split_training_type == "random"
		trainLearnFileName = string(fileName, "_tr_learn.jld")
		trainTestFileName = string(fileName, "_tr_test.jld")
	else
		trainLearnFileName = string(fileName, "_loc_tr_learn.jld")
		trainTestFileName = string(fileName, "_loc_tr_test.jld")
	end
	trainLearnDf, trainTestDf, testDf = loadTrainingAndTestingSets(startTime=startTime, endTime=endTime, weekdays=weekdays, startMonth=startMonth, endMonth=endMonth, split_training_type=split_training_type)
	if isfile(trainLearnFileName) && loadFromCache
		travelTimes = load(trainLearnFileName, "travelTimes")
		numRides = load(trainLearnFileName, "numRides")
	else
		travelTimes, numRides = mapRidesToNodes(trainLearnDf, nodePositions, method, radius=radius)
		save(trainLearnFileName, "travelTimes", travelTimes, "numRides", numRides)
	end
	if isfile(trainTestFileName) && loadFromCache
		testingTravelTimes = load(trainTestFileName, "travelTimes")
		testingNumRides = load(trainTestFileName, "numRides")
	else
		testingTravelTimes, testingNumRides = mapRidesToNodes(trainTestDf, nodePositions, method, radius=radius)
		save(trainTestFileName, "travelTimes", testingTravelTimes, "numRides", testingNumRides)
	end
	return travelTimes, numRides, trainTestDf, testingTravelTimes, testingNumRides, testDf
end

function buildNodeKDTree(nodePositions::Array{Coordinates,1})
	"""
	Given an array of node positions, returns a KDTree containing the nodePairs in 4D space. Also returns array mapping index in KDTree to pair of nodes.
	"""
	nodeTree = zeros(4, length(nodePositions)^2 - length(nodePositions))
	index = 1
	nodePairs = zeros(Int, (length(nodePositions)^2 - length(nodePositions), 2))
	for i = 1:length(nodePositions), j=1:length(nodePositions)
		if i != j
			nodeTree[1, index] = nodePositions[i].x
			nodeTree[2, index] = nodePositions[i].y
			nodeTree[3, index] = nodePositions[j].x
			nodeTree[4, index] = nodePositions[j].y
			nodePairs[index, 1] = i
			nodePairs[index, 2] = j
			index += 1
		end
	end
	@time nodeTree = KDTree(nodeTree)
	return nodeTree, nodePairs
end

function computeAverageTestingError(testingTravelTimes::Array{Float64,2}, algorithmOutput::Array{Float64, 2}; num_nodes::Int = MANHATTAN_NODES)
	"""
	Given a matrix of TABs from testing Data and calculated by our algorithm, returns some error measures
	"""
	num_results = 0
	average_squared_error = 0
	average_relative_error = 0
	average_bias = 0
	for i = 1:num_nodes, j=1:num_nodes
		if testingTravelTimes[i,j] > 0
			average_squared_error = (num_results * average_squared_error + (testingTravelTimes[i,j] - algorithmOutput[i,j]) ^ 2) / (num_results + 1)
			average_relative_error = (num_results * average_relative_error + abs(testingTravelTimes[i,j] - algorithmOutput[i,j])/testingTravelTimes[i,j])/(num_results + 1)
			average_bias = (num_results * average_bias + (algorithmOutput[i,j] - testingTravelTimes[i,j]))/(num_results + 1)
			num_results += 1
		end
	end
	println("------------------------------------------------")
	println("Average squared error: \t\t\t\t\t", average_squared_error)
	println("Average relative error: \t\t\t\t", average_relative_error)
	println("Average bias: \t\t\t\t\t\t", average_bias)
	println("------------------------------------------------")
	return average_squared_error, average_relative_error, average_bias
end

function computeTestingError(rideTestingSet::DataFrame, nodeTree::KDTree{Float64}, nodePairs::Array{Int,2}, algorithmOutput::Array{Float64,2}, method::AbstractString; radius::Int = 140)
	"""
	Uses testing set to compute error on our method.
	"""
	# Initialize vars
	num_results = 0
	average_squared_error = 0
	average_relative_error = 0
	average_bias = 0
	# Loop through testing set and choose method
	for i = 1:nrow(testDf)
		if method == "nearest"
			# Find nearest node pair to this ride
			index, dist = knn(nodeTree, [testDf[i,:pX], testDf[i,:pY], testDf[i,:dX], testDf[i,:dY]], 1)
			idx = index[1]
			# Calculate time of ride in seconds
			testingTime = (testDf[i,:dTime] - testDf[i,:pTime]).value/1000 # don't forget to convert
			# Calculate predicted time
			predictedTime = algorithmOutput[nodePairs[idx,1],nodePairs[idx,2]]
			# Update errors
			if predictedTime != 0
				average_squared_error = (num_results * average_squared_error + (testingTime - predictedTime) ^ 2) / (num_results + 1)
				average_relative_error = (num_results * average_relative_error + abs(testingTime - predictedTime)/testingTime) / (num_results + 1)
				average_bias = (num_results * average_bias + (predictedTime - testingTime))/(num_results + 1)
				num_results += 1
			end
		elseif method == "average"
			index = inball(nodeTree, [testDf[i,:pX], testDf[i,:pY], testDf[i,:dX], testDf[i,:dY]], float(radius))
			# Compute testing time in seconds
			testingTime = (testDf[i,:dTime] - testDf[i,:pTime]).value/1000
			# Compute predicted time by averaging neighboring node pairs (within radius)
			predictedTime = 0.
			numPoints = 0
			for idx in index
				predictedTime = (predictedTime * numPoints + algorithmOutput[nodePairs[idx,1],nodePairs[idx,2]])/(numPoints + 1)
				numPoints += 1
			end
			# Update errors (if any neighbors were found within 140 meters)
			if predictedTime != 0
				average_squared_error = (num_results * average_squared_error + (testingTime - predictedTime) ^ 2) / (num_results + 1)
				average_relative_error = (num_results * average_relative_error + abs(testingTime - predictedTime)/testingTime) / (num_results + 1)
				average_bias = (num_results * average_bias + (predictedTime - testingTime))/(num_results + 1)
				num_results += 1
			end
		end
	end
	println("------------------------------------------------")
	println("Average squared error: \t\t\t\t\t", average_squared_error)
	println("Average relative error: \t\t\t\t", average_relative_error)
	println("Average bias: \t\t\t\t\t\t", average_bias)
	println("------------------------------------------------")
	return average_squared_error, average_relative_error, average_bias
end

ext: jl

max_stars_repo_name: sebmart/NYC-timing

max_issues_repo_name: sebmart/NYC-timing

max_forks_repo_name: sebmart/NYC-timing


text:
using BlockArrays, Test

@testset "broadcast" begin
    @testset "BlockArray" begin
        A = BlockArray(randn(6), 1:3)

        @test BlockArrays.BroadcastStyle(typeof(A)) == BlockArrays.BlockStyle{1}()

        @test exp.(A) == exp.(Vector(A))
        @test blocksizes(A) == blocksizes(exp.(A))

        @test A+A isa BlockArray
        @test blocksizes(A + A) == blocksizes(A .+ A) == blocksizes(A)
        @test blocksizes(A .+ 1) == blocksizes(A)

        A = BlockArray(randn(6,6), 1:3,1:3)

        @test BlockArrays.BroadcastStyle(typeof(A)) == BlockArrays.BlockStyle{2}()

        @test exp.(A) == exp.(Matrix(A))
        @test blocksizes(A) == blocksizes(exp.(A))

        @test blocksizes(A + A) == blocksizes(A .+ A) == blocksizes(A)
        @test blocksizes(A .+ 1) == blocksizes(A)
    end

    @testset "PseudoBlockArray" begin
        A = PseudoBlockArray(randn(6), 1:3)

        @test BlockArrays.BroadcastStyle(typeof(A)) == BlockArrays.PseudoBlockStyle{1}()

        @test exp.(A) == exp.(Vector(A))
        @test blocksizes(A) == blocksizes(exp.(A))

        @test A+A isa PseudoBlockArray
        @test blocksizes(A + A) == blocksizes(A .+ A) == blocksizes(A)
        @test blocksizes(A .+ 1) == blocksizes(A)

        B = PseudoBlockArray(randn(6,6), 1:3,1:3)

        @test BlockArrays.BroadcastStyle(typeof(B)) == BlockArrays.PseudoBlockStyle{2}()

        @test exp.(B) == exp.(Matrix(B))
        @test blocksizes(B) == blocksizes(exp.(B))

        @test blocksizes(B + B) == blocksizes(B .+ B) == blocksizes(B)
        @test blocksizes(B .+ 1) == blocksizes(B)
        @test blocksizes(A .+ 1 .+ B) == blocksizes(B)
        @test A .+ 1 .+ B == Vector(A) .+ 1 .+ B == Vector(A) .+ 1 .+ Matrix(B)
    end

    @testset "Mixed" begin
        A = BlockArray(randn(6), 1:3)
        B = PseudoBlockArray(randn(6), 1:3)

        @test A + B isa BlockArray
        @test B + A isa BlockArray

        @test blocksizes(A + B) == blocksizes(A)

        C = randn(6)

        @test A + C isa BlockVector{Float64}
        @test C + A isa BlockVector{Float64}
        @test B + C isa PseudoBlockVector{Float64}
        @test C + B isa PseudoBlockVector{Float64}

        @test blocksizes(A+C) == blocksizes(C+A) == blocksizes(A)
        @test blocksizes(B+C) == blocksizes(C+B) == blocksizes(B)

        A = BlockArray(randn(6,6), 1:3, 1:3)
        D = Diagonal(ones(6))
        @test blocksizes(A + D) == blocksizes(A)
        @test blocksizes(B .+ D) == BlockArrays.BlockSizes([1,2,3],[6])
    end

    @testset "Mixed block sizes" begin
        A = BlockArray(randn(6), 1:3)
        B = BlockArray(randn(6), fill(2,3))
        @test blocksizes(A+B) == BlockArrays.BlockSizes([1,1,1,1,2])

        A = BlockArray(randn(6,6), 1:3, 1:3)
        B = BlockArray(randn(6,6), fill(2,3), fill(3,2))

        @test blocksizes(A+B) == BlockArrays.BlockSizes([1,1,1,1,2], 1:3)
    end
end


text:
using DiffEqFlux, StaticArrays, Test

fd = FastDense(2,25,tanh)
pd = initial_params(fd)
fd(ones(2),pd)

f1 = FastDense(2,25,tanh)
f2 = FastDense(25,2,tanh)
p1 = initial_params(f1)
p2 = initial_params(f2)
@test FastChain(f1,f2)(ones(2),[p1;p2]) == f2(f1(ones(2),p1),p2)

f = FastChain(FastDense(2,25,tanh),FastDense(25,2,tanh))
p = initial_params(f)
@test f(ones(2),p) == f2(f1(ones(2),p[1:length(p1)]),p[length(p1)+1:end])

fd1 = FastDense(2,25,tanh,bias=false)
pd1 = initial_params(fd1)
fd1(ones(2),pd1)

f3 = FastDense(2,25,tanh,bias=false)
f4 = FastDense(25,2,tanh,bias=false)
p3 = initial_params(f3)
p4 = initial_params(f4)

@test FastChain(f3,f4)(ones(2),[p3;p4]) == f4(f3(ones(2),p3),p4)

f5 = FastChain(FastDense(2,25,tanh,bias=false),FastDense(25,2,tanh,bias=false))
p5 = initial_params(f5)
@test f5(ones(2),p5) == f4(f3(ones(2),p5[1:length(p3)]),p5[length(p3)+1:end])

fs = StaticDense(2,25,tanh)
x = rand(2)

@test fs(x,pd) ≈ fd(x,pd)
fdgrad = Flux.Zygote.gradient((x,p)->sum(fd(x,p)),x,pd)
fsgrad = Flux.Zygote.gradient((x,p)->sum(fs(x,p)),x,pd)
@test fdgrad[1] ≈ fsgrad[1]
@test fdgrad[2] ≈ fsgrad[2] rtol=1e-5

# Now test vs Zygote
struct TestDense{F,F2} <: DiffEqFlux.FastLayer
  out::Int
  in::Int
  σ::F
  initial_params::F2
  function TestDense(in::Integer, out::Integer, σ = identity;
                 initW = Flux.glorot_uniform, initb = Flux.zeros)
    initial_params() = vcat(vec(initW(out, in)),initb(out))
    new{typeof(σ),typeof(initial_params)}(out,in,σ,initial_params)
  end
end
(f::TestDense)(x,p) = f.σ.(reshape(p[1:(f.out*f.in)],f.out,f.in)*x .+ p[(f.out*f.in+1):end])
ts = StaticDense(2,25,tanh)
testgrad = Flux.Zygote.gradient((x,p)->sum(ts(x,p)),x,pd)
@test fdgrad[1] ≈ testgrad[1]
@test fdgrad[2] ≈ testgrad[2] rtol=1e-5

fsgrad = Flux.Zygote.gradient((x,p)->sum(fs(x,p)),SVector{2}(x),SVector{75}(pd))
@test fsgrad[1] isa SArray
@test fsgrad[2] isa SArray

layer = FastDense(3, 4, bias=false)
p = initial_params(layer)

rand_loss(p) = begin
    y = sum(abs2, layer(randn(3, 5), p))
end

y, back = Zygote.pullback(rand_loss, p)
grad = back(1.0)[1]
@test length(grad) == length(p)

layer = StaticDense(3, 4, bias=false)
p = initial_params(layer)

rand_loss(p) = begin
    y = sum(abs2, layer(randn(3), p))
end

y, back = Zygote.pullback(rand_loss, p)
grad = back(1.0)[1]
@test length(grad) == length(p)


text:
const PP = PhysicalParam(4)

# const NORMAL_BEHAVIOR = IDMMOBILBehavior("normal", PP.v_med, PP.l_car, 1)
# const THREE_BEHAVIORS = IDMMOBILBehavior[
#                     NORMAL_BEHAVIOR,
#                     IDMMOBILBehavior("cautious", PP.v_slow+0.5, PP.l_car, 2),
#                     IDMMOBILBehavior("aggressive", PP.v_fast, PP.l_car, 3)
#                     ]
#
# const NINE_BEHAVIORS = IDMMOBILBehavior[IDMMOBILBehavior(x[1],x[2],x[3],idx) for (idx,x) in
#                                                  enumerate(Iterators.product(["cautious","normal","aggressive"],
#                                                         [PP.v_slow+0.5;PP.v_med;PP.v_fast],
#                                                         [PP.l_car]))]

const UNIFORM_MEAN = IDMMOBILBehavior(IDMParam(1.4, 2.0, 1.5, 33.3, 2.0, 4.0), MOBILParam(0.5, 2.0, 0.1), 1)

const NORMAL_IDM = IDMParam(1.4, 2.0, 1.5, 120/3.6, 2.0, 4.0) # 33.3
const TIMID_IDM = IDMParam(1.0, 1.0, 1.8, 100/3.6, 4.0, 4.0) # 27.7
const AGGRESSIVE_IDM = IDMParam(2.0, 3.0, 1.0, 140/3.6, 1.0, 4.0) # 38.9
agents_behavior(idm, idx) = IDMMOBILBehavior(idm, MOBILParam(0.6, idm.b, 0.1), idx)

const NORMAL = agents_behavior(NORMAL_IDM, 1)
const TIMID = agents_behavior(TIMID_IDM, 2)
const AGGRESSIVE = agents_behavior(AGGRESSIVE_IDM, 3)

function generate_nine_behaviors()
    bs = IDMMOBILBehavior[]
    velocities = [120/3.6, 100/3.6, 140/3.6]
    for (i,v) in enumerate(velocities)
        push!(bs, agents_behavior(IDMParam(1.4, 2.0, 1.5, v, 2.0, 4.0), 1+(i-1)*3))
        push!(bs, agents_behavior(IDMParam(1.0, 1.0, 1.8, v, 4.0, 4.0), 2+(i-1)*3))
        push!(bs, agents_behavior(IDMParam(2.0, 3.0, 1.0, v, 1.0, 4.0), 3+(i-1)*3))
    end
    return bs
end

const NINE_BEHAVIORS = generate_nine_behaviors()

const DEFAULT_BEHAVIORS = Dict{String, Any}(
    # "3_even" => DiscreteBehaviorSet(THREE_BEHAVIORS, WeightVec(ones(3))),
    # "9_even" => DiscreteBehaviorSet(NINE_BEHAVIORS, WeightVec(ones(9))),
    "agents" => DiscreteBehaviorSet([NORMAL, TIMID, AGGRESSIVE], Weights(ones(3))),
    "9_agents" => DiscreteBehaviorSet(NINE_BEHAVIORS, Weights(ones(3)))
)

const DEFAULT_PROBLEM_PARAMS = Dict{Symbol, Any}( #NOTE VALUES ARE NOT VECTORS like in linked
    :behaviors => "agents",
    :p_appear => 1.0,
    :vel_sigma => 0.5,
    :brake_terminate_thresh => Inf,
    :max_dist => Inf,
    :lane_terminate => false,

    :rmodel_type => :NoCrashRewardModel,
    :lambda => 1.0,
    :brake_threshold => 2.5,
    :target_lane => 4,
)

const INITIAL_RELEVANT = [:behaviors]

type TestSet
    solver_key::String
    problem_params::Dict{Symbol, Any}
    solver_problem_params::Dict{Symbol, Any}
    linked_problem_params::Dict{Symbol, AbstractVector}
    N::Int # number of initial conditions
    nb_problems::Int
    rng_seed::UInt32
    key::String
end

"""
Generate a test set.

The keyword arguments should be either
    1) fields of the TestSet type, or
    2) problem parameters.

If the the argument in is a problem parameter, if it begins with "solver_", it
will only be applied to the problem given to the solver, otherwise it will be
applied to both problems. The arguments will be applied in order, so solver_
parameter values should come later.
"""
function TestSet(ts::TestSet; kwargs...)
    ts = deepcopy(ts)
    fn = fieldnames(TestSet)
    for (k,v) in kwargs
        if k in fn
            setfield!(ts, k, v)
        elseif isa(v, AbstractVector)
            if ts.nb_problems == 1
                ts.nb_problems = length(v)
            end
            if length(v) != ts.nb_problems
                error("""
                    All linked parameters in a TestSet must have the same number of values.
                    Argument $k had had $(length(v)) entries while others had $(ts.nb_problems) entries.
                    """)
            end
            ts.linked_problem_params[k] = v
            delete!(ts.problem_params, k)
            delete!(ts.solver_problem_params, k)
        elseif haskey(ts.problem_params, k)
            ts.problem_params[k] = v
            ts.solver_problem_params[k] = v
        else
            solver_match = match(r"solver_(.*)", string(k))
            if solver_match != nothing
                sk = Symbol(solver_match[1])
                if haskey(ts.solver_problem_params, sk)
                    ts.solver_problem_params[sk] = v
                end
            else
                warn("Unrecognized TestSet argument $k")
            end
        end
    end
    return ts::TestSet
end

function TestSet(;kwargs...)
    return TestSet(TestSet(randstring()); kwargs...)
end

function TestSet(key::AbstractString)
    return TestSet("",
            deepcopy(DEFAULT_PROBLEM_PARAMS),
            deepcopy(DEFAULT_PROBLEM_PARAMS),
            Dict{String, AbstractVector}(),
            500, 1, 1, key)
end

# function TestSet(;key::AbstractString=randstring(), kwargs...)
#=
function TestSet(;kwargs...)
    ts = TestSet("",
            deepcopy(default_problem_params),
            deepcopy(default_problem_params),
            deepcopy(linked_problem_params),
            500,
            1,
            randstring())
            # key)
    return TestSet(ts, kwargs...)
end
=#

function gen_initials(tests::AbstractVector, initials::Dict=Dict{String,Any}();
                      behaviors::Dict{String,Any}=get(initials, "behaviors", DEFAULT_BEHAVIORS),
                      #generate_physical=false,
                      generate_physical=true,
                      rng::AbstractRNG=MersenneTwister(rand(UInt32)))
    initials=Dict{String,Any}(k=>v for (k,v) in initials)
    for t in tests
        add_initials!(initials, t, behaviors=behaviors, rng=rng, generate_physical=generate_physical)
    end
    return initials
end

function gen_initials_cz(tests::TestSet, initials::Dict=Dict{String,Any}();
                      behaviors::Dict{String,Any}=get(initials, "behaviors", DEFAULT_BEHAVIORS),
                      generate_physical=true,
                      #generate_physical=false,
                      rng::AbstractRNG=MersenneTwister(rand(UInt32)))
    initials=Dict{String,Any}(k=>v for (k,v) in initials)
    #for t in tests
    print("begin ini")
    add_initials!(initials, tests, behaviors=behaviors, rng=rng, generate_physical=generate_physical)
    #end
    print("end ini")
    return initials
end

function gen_dmodel(row, behaviors::BehaviorGenerator)
    nb_lanes = 4
    nb_cars = 10
    pp = PhysicalParam(nb_lanes,lane_length=100.)
    dmodel = NoCrashIDMMOBILModel(nb_cars, pp,
                              behaviors=DEFAULT_BEHAVIORS["agents"],
                              lane_terminate=false)

    dmodel.behaviors = behaviors
    # p_appear
    dmodel.p_appear = row[:p_appear]
    # vel_sigma
    dmodel.vel_sigma = row[:vel_sigma]
    # max_dist
    dmodel.max_dist = row[:max_dist]
    # lane_terminate
    dmodel.lane_terminate = row[:lane_terminate]
    # brake_terminate_thresh
    dmodel.brake_terminate_thresh = row[:brake_terminate_thresh]

    return dmodel
end

function gen_rmodel(row)
#=
    if Symbol(row[:rmodel_type]) == :TargetLaneReward
        return TargetLaneReward(row[:target_lane])
    elseif Symbol(row[:rmodel_type]) == :NoCrashRewardModel
        target_lane_reward = 10.
        rmodel = NoCrashRewardModel(target_lane_reward*10., target_lane_reward,2.5,row[:target_lane])
        rmodel.cost_dangerous_brake = row[:lambda]*rmodel.reward_in_target_lane
        rmodel.brake_penalty_thresh = row[:brake_threshold]
    end
=#
########I changed the Rmodel
    target_lane_reward = 10.
    rmodel = NoCrashRewardModel(target_lane_reward*10., target_lane_reward,2.5,row[:target_lane])
    rmodel.cost_dangerous_brake = row[:lambda]*rmodel.reward_in_target_lane
    rmodel.brake_penalty_thresh = row[:brake_threshold]
    return rmodel
end

function gen_problem(row, behaviors::Dict{String,Any})
    _discount = 1.0
    dmodel = gen_dmodel(row, behaviors[row[:behaviors]])
    rmodel = gen_rmodel(row)
    return MLMDP{MLState, MLAction, typeof(dmodel), typeof(rmodel)}(dmodel, rmodel, _discount)
end

function gen_initial_physical(base_problem, N; rng::AbstractRNG=MersenneTwister(rand(UInt32)))
    # return [MLPhysicalState(initial_state(base_problem, rng)) for i in 1:N]
    return [MLPhysicalState(relaxed_initial_state(base_problem, 200, rng)) for i in 1:N]
end

function add_initials!(objects::Dict{String, Any},
                       ts::TestSet;
                       behaviors::Dict{String,Any}=get(objects, "behaviors"),
                       rng::AbstractRNG=MersenneTwister(rand(UInt32)),
                       generate_physical=false)
                       #generate_physical=true)

    new_table = DataFrame()

    for p in ts.linked_problem_params
        if isempty(new_table)
            new_table = DataFrame(Dict(p))
        else
            new_table = join(new_table, DataFrame(Dict(p)), kind=:cross)
        end
    end

    # run through and check to see if all are the same
    different = false
    for (k,v) in ts.problem_params
        if ts.solver_problem_params[k] != v
            different = true
            break
        end
    end

    if different
        dict = Dict(k=>typeof(v)[v,ts.solver_problem_params[k]] for (k,v) in ts.problem_params)
        specific_table = DataFrame(dict)
    else
        vectors = Dict(k=>DataArray(typeof(v)[v]) for (k,v) in ts.problem_params)
        specific_table = DataFrame(vectors)
    end

    if nrow(new_table) == 0 # there were no linked params
        new_table = specific_table
    else
        new_table = join(new_table, specific_table, kind=:cross)
    end
    param_list = names(new_table)

    problems = get(objects, "problems", Dict{String,Any}())
    if haskey(objects, "param_table")
        param_table = objects["param_table"]
        for p in param_list
            if !(p in names(param_table))
                default = DEFAULT_PROBLEM_PARAMS[p]
                warn("adding $p to param table with default value $default.")
                param_table[p] = default
            end
        end
        param_table = join(param_table, new_table, on=param_list, kind=:outer)
    else
        param_table = new_table
        param_table[:problem_key] = DataArray(String, nrow(param_table))
        param_table[:state_list_key] = DataArray(String, nrow(param_table))
    end

    # go through and make sure each row has a problem
    for row in eachrow(param_table)
        if first(isna(row[:problem_key]))
            key = randstring(rng)
            problems[key] = gen_problem(row, behaviors)
            row[:problem_key] = key
        end
    end

    initial_states = get(objects, "initial_states", Dict{String,Any}())
    state_lists = get(objects, "state_lists", Dict{String,Any}())

    if haskey(objects, "initial_physical_states")
        initial_physical_states = objects["initial_physical_states"]
    elseif generate_physical
        initial_physical_states = gen_initial_physical(first(values(problems)), ts.N, rng=rng)
    elseif isinteractive()
        println("\n\nNo initial physical states found! Press Enter to generate them, Ctrl-C to cancel.")
        readline(STDIN)
        initial_physical_states = gen_initial_physical(gen_base_problem(), ts.N, rng=rng)
    else
        error("""No initial physical states found!""")
    end
     print("check st")
    for g in groupby(param_table, INITIAL_RELEVANT)
        if isna(first(g[:state_list_key]))
            key = randstring(rng)
            pk = g[1,:problem_key]
            p = problems[pk]
            these_states = assign_keys([initial_state(p, initial_physical_states[i], rng) for i in 1:ts.N])
            merge!(initial_states, these_states)
            g[:,:state_list_key] = key
            state_lists[key] = collect(keys(these_states))
        end
    end
    print("check en")
    objects["param_table"] = param_table
    objects["problems"] = problems
    objects["state_lists"] = state_lists
    objects["initial_states"] = initial_states
    objects["initial_physical_states"] = initial_physical_states
    objects["behaviors"] = behaviors

    return objects
end

function find_row(table::DataFrame, vals::Dict{Symbol,Any})
    for i in 1:nrow(table)
        found = true
        for (k,v) in vals
            if table[i,k] != v
                found = false
                break
            end
        end
        if found
            return table[i,:]
        end
    end
end

function setup_stats(tests::AbstractVector, objects::Dict{String,Any})
    param_table = objects["param_table"]
    state_lists = objects["state_lists"]

    nb_sims = sum([t.N*t.nb_problems for t in tests])

    stats = DataFrame(
        id=1:nb_sims,
        uuid=UInt128[Base.Random.uuid4() for i in 1:nb_sims],
        solver_key=DataArray(String,nb_sims),
        problem_key=DataArray(String,nb_sims),
        solver_problem_key=DataArray(String,nb_sims),
        test_key=DataArray(String,nb_sims),
        initial_key=DataArray(String,nb_sims),
        rng_seed=DataArray(Int,nb_sims),
        time=ones(nb_sims).*time(),
    )

    id = 0
    for t in tests
        test_key = t.key
        solver_key = t.solver_key
        for i in 1:t.nb_problems
            # find the row
            d = Dict{Symbol, Any}()
            solver_d = Dict{Symbol, Any}()
            for (k,v) in t.problem_params
                d[k] = v
                solver_d[k] = t.solver_problem_params[k]
            end
            for (k,v) in t.linked_problem_params
                d[k] = v[i]
                solver_d[k] = v[i]
            end
            row = find_row(param_table, d)
            ep_key = first(row[:problem_key])
            solver_row = find_row(param_table, solver_d)
            sp_key = first(solver_row[:problem_key])
            is_keys = state_lists[first(row[:state_list_key])]
            for (is_i, is_key) in enumerate(take(is_keys, t.N))
                id += 1
                stats[:solver_key][id] = solver_key
                stats[:problem_key][id] = ep_key
                stats[:solver_problem_key][id] = sp_key
                stats[:initial_key][id] = is_key
                stats[:rng_seed][id] = is_i+t.rng_seed
                stats[:test_key][id] = test_key
            end
        end
    end
    return stats
end

"""
Run the simulations in tests.
"""
function evaluate(tests::AbstractVector, objects::Dict{String,Any};
                  parallel=true,
                  desc="Progress: ",
                  metrics::AbstractVector=[],
                  max_steps=10000)

    stats = setup_stats(tests, objects)
    sort!(stats, cols=:uuid) # so that they are randomly distributed so that the time estimate will be better

    sims = run_simulations(stats, objects,
                           parallel=parallel,
                           desc=desc, max_steps=max_steps)
    fill_stats!(stats, objects, sims, metrics=metrics)

    results = deepcopy(objects)
    results["stats"] = stats
    results["histories"] = sims
    results["tests"] = Dict{String,Any}([(t.key, t) for t in tests])
    return results

end

function merge_tests!(t1::Dict, t2::Dict)
    for (k,t) in t2
        if haskey(t1, k)
            this = t1[k]
            @assert this.solver_key == t.solver_key
            careful_merge!(this.problem_params, t.problem_params)
            careful_merge!(this.solver_problem_params, t.solver_problem_params)
            @assert this.N == t.N
            @assert this.rng_seed == t.rng_seed
            @assert this.key == t.key
            for p in keys(t.linked_problem_params)
                this.linked_problem_params[p] = cat(1,this.linked_problem_params[p],t.linked_problem_params[p])
            end
            this.nb_problems += t.nb_problems
        else
            t1[k] = t
        end
    end
    return t1
end

max_stars_repo_name: zhcao/HighwayExit

max_issues_repo_name: zhcao/HighwayExit

max_forks_repo_name: zhcao/HighwayExit


text:
import Base: show

type AR{T<:Float64}
    target::Vector{T}
    features::Matrix{T}   # design matrix
    params::Vector{T}
    errors::Vector{T}
    se::Vector{T}         # standard errors
end

function show(io::IO, arm::AR)
    tstats = arm.params ./ arm.se
    df     = size(arm.features,1) - size(arm.features,2) # rows (or n) - columns (or p+1)
    pt     = cdf(TDist(df), tstats)
    cpt    = 100 * ccdf(TDist(df), tstats)
#     println(io, @sprintf("coefficients:             %.4f  %.4f  %.4f", arm.params[1], arm.params[2], arm.params[3]))
#     println(io, @sprintf("standard errors:          %.4f  %.4f  %.4f", arm.se[1], arm.se[2], arm.se[3]))
#     println(io, @sprintf("t values:                 %.4f  %.4f  %.4f", tstats[1], tstats[2], tstats[3]))
#     println(io, @sprintf("P(>|t|):                  %.4f  %.4f  %.4f", pt[1], pt[2], pt[3]))
    println(io, "Meaningful coefficients:")
    cpt[1]  > 95 ?  println(io, @sprintf("Intercept:  %.1f", cpt[1]))  : nothing
    cpt[2]  > 95 ?  println(io, @sprintf("Phi 1:      %.1f", cpt[2]))  : nothing
    cpt[3]  > 95 ?  println(io, @sprintf("Phi 2:      %.1f", cpt[3]))  : nothing
    cpt[4]  > 95 ?  println(io, @sprintf("Phi 3:      %.1f", cpt[4]))  : nothing
    cpt[5]  > 95 ?  println(io, @sprintf("Phi 4:      %.1f", cpt[5]))  : nothing
    cpt[6]  > 95 ?  println(io, @sprintf("Phi 5:      %.1f", cpt[6]))  : nothing
    cpt[7]  > 95 ?  println(io, @sprintf("Phi 6:      %.1f", cpt[7]))  : nothing
    cpt[8]  > 95 ?  println(io, @sprintf("Phi 7:      %.1f", cpt[8]))  : nothing
    cpt[9]  > 95 ?  println(io, @sprintf("Phi 8:      %.1f", cpt[9]))  : nothing
    cpt[10] > 95 ? println(io, @sprintf("Phi 9:      %.1f", cpt[10]))  : nothing
    cpt[11] > 95 ? println(io, @sprintf("Phi 10:     %.1f", cpt[11]))  : nothing

    println(io, "")
end

#function ar(v::Vector{Float64}, p::Int)
function ar(ta::TimeArray{Float64,1}, p::Int)

# p = 10
# ta = c["Settle"]

    pc = percentchange(ta, method="log")
    v  = flipud(pc.values)
    Xs = ones(length(v)-p, p+1)
    for i in 1:length(v)-p
        Xs[i,2:end] = flipud(v[i:i+p-1])'
    end
    y = v[p+1:end]

    # sigs = findall(abs(StatsBase.autocor(v)) .> .05) # find rows that pass .05 confidence
    # sig  = sigs[sigs.<11] # don't care about lags more than 10

    sig = 1:p+1 #canonical method to include all features up to and including the last one that shows serial correlation

    X = Xs[:,sig] # the first column remains 1s, some put this as the last

    # Normal Equations
    # ß = inv(X'X)*X'*y

    # QR Decomposition
    Q,R = qr(X)

    beta    = inv(R)Q'y
    epsilon = y-X*beta
    mse     = mean(epsilon.^2)
    ∑       = mse * inv(X'X)
    stderr  = sqrt(diag(∑)) # diag(sqrtm(∑))

    AR(y,X,beta,epsilon,stderr)
end

ar(ta::TimeArray{Float64,1}) = ar(ta,10)


text:
using Colors, Zygote

target = RGB(1, 0, 0)
colour = RGB(1, 1, 1)

function update_color(c, Δ, η = 0.001)
    return RGB(
        c.r - η*Δ.r,
        c.g - η*Δ.g,
        c.b - η*Δ.b,
    )
end

for idx in 1:51
    global colour, target

    # Calculate gradients
    grads = Zygote.gradient(colour) do y
        colordiff(target, y)
    end
    # Update colour
    colour = update_color(colour, grads[1])
    if idx % 5 == 1
        @info idx, colour
    end
end

max_stars_repo_name: MikeInnes/zygote-paper

max_issues_repo_name: MikeInnes/zygote-paper

max_forks_repo_name: MikeInnes/zygote-paper


text:
# Adders

function HalfAdder(; name)
    @named x1 = DigitalPin()
    @named x2 = DigitalPin()
    @named y1 = DigitalPin()
    @named y2 = DigitalPin()
    @variables sum(t), carry(t)

    eqs = [
        y1.val ~ _xor(x1.val, x2.val)
        y2.val ~ _and(x1.val, x2.val)
        sum ~ y1.val
        carry ~ y2.val
    ]
    ODESystem(eqs, t, [sum, carry], [], systems=[x1, x2, y1, y2], name=name)
end

function FullAdder(; name)
    @named x1 = DigitalPin()
    @named x2 = DigitalPin()
    @named x3 = DigitalPin()
    @named y1 = DigitalPin()
    @named y2 = DigitalPin()
    @variables sum(t), carry(t)

    eqs = [
        y1.val ~ _xor(x1.val, x2.val, x3.val)
        y2.val ~ _or(_and(x3.val, _xor(x1.val, x2.val)), _and(x1.val, x2.val))
        sum ~ y1.val
        carry ~ y2.val
    ]
    ODESystem(eqs, t, [sum, carry], [], systems=[x1, x2, x3, y1, y2], name=name)
end

# Multiplexers

# This selects data from the `N` input ports (`d₀` to `dₙ₋₁`)
# using values of `n` select lines, where `N = 2^n`
function MUX(; name, N=4)
    n = log2(N)
    try n = Int(n) catch(e) throw("`N` must be a power of 2") end
    s = map(0:n-1) do i
        DigitalPin(; name=Symbol(:s, i))
    end
    d = map(0:N-1) do i
        DigitalPin(; name=Symbol(:d, i))
    end
    @named y = DigitalPin()

    nodes = Num[]
    for i in 1:N
        bin = digits!(zeros(Int64,n), i-1, base=2)
        statelist = Term{Real, Nothing}[]
        for j in 1:n
            varstate = bin[j] == 0 ? _not(s[j].val) : s[j].val
            push!(statelist, varstate)
        end
        push!(nodes, _and(statelist..., d[i].val))
    end

    eqs = Equation[
        y.val ~ _or(nodes...)
    ]

    ODESystem(eqs, t, [], [], systems=[d..., s..., y], name=name)
end

# This selects one of the `N` output ports (`y₀` to `yₙ₋₁`)
# to transmit data `d` using values of `n` select lines, where `N = 2^n`
function DEMUX(; name, N=4)
    n = log2(N)
    try n = Int(n) catch(e) throw("`N` must be a power of 2") end
    @named d = DigitalPin()
    s = map(0:n-1) do i
        DigitalPin(; name=Symbol(:s, i))
    end
    y = map(0:N-1) do i
        DigitalPin(; name=Symbol(:y, i))
    end

    eqs = Equation[]
    for i in 1:N
        bin = digits!(zeros(Int64,n), i-1, base=2)
        statelist = Term{Real, Nothing}[]
        for j in 1:n
            varstate = bin[j] == 0 ? _not(s[j].val) : s[j].val
            push!(statelist, varstate)
        end
        push!(eqs, y[i].val ~ _and(statelist..., d.val))
    end

    ODESystem(eqs, t, [], [], systems=[d, s..., y...], name=name)
end

# Encoder-Decoder

# Encodes `N` inputs to `n` outputs, where `N = 2^n`
function Encoder(; name, N=4)
    n = log2(N)
    try n = Int(n) catch(e) throw("`N` must be a power of 2") end
    d = map(0:N-1) do i
        DigitalPin(; name=Symbol(:d, i))
    end
    y = map(0:n-1) do i
        DigitalPin(; name=Symbol(:y, i))
    end

    nodes = Vector{Term{Real, Nothing}}[]
    i = 0
    for j in 1:n
        counter = 1
        statelist = Term{Real, Nothing}[]
        while i < N
            while counter <= 2^j
                counter > 2^(j-1) && push!(statelist, d[i+1].val)
                counter += 1
                i = i+1
            end
            counter = 1
        end
        i = 0
        push!(nodes, statelist)
    end

    eqs = Equation[]
    for i in n:-1:1
        push!(eqs, y[i].val ~ _or(nodes[i]...))
    end

    ODESystem(eqs, t, [], [], systems=[d..., y...], name=name)
end

# Decodes `n` inputs to `N` outputs, where `N = 2^n`
function Decoder(; name, n=2)
    N = 2^n
    d = map(0:n-1) do i
        DigitalPin(; name=Symbol(:d, i))
    end
    y = map(0:N-1) do i
        DigitalPin(; name=Symbol(:y, i))
    end

    nodes = Vector{Term{Real, Nothing}}[]
    for i in 1:N
        bin = digits!(zeros(Int64, n), i-1, base=2)
        statelist = Term{Real, Nothing}[]
        for j in 1:n
            varst = bin[j] == 0 ? _not(d[j].val) : d[j].val
            push!(statelist, varst)
        end
        push!(nodes, statelist)
    end

    eqs = Equation[]
    for i in N:-1:1
        push!(eqs, y[i].val ~ _and(nodes[i]...))
    end

    ODESystem(eqs, t, [], [], systems=[d..., y...], name=name)
end


text: include("truncatedconjugategradient.jl")

max_stars_repo_name: JuliaNLSolvers/NLSolvers

max_issues_repo_name: JuliaNLSolvers/NLSolvers

max_forks_repo_name: JuliaNLSolvers/NLSolvers


text:
test_ex(x, y) = @test prettify(x) == prettify(y)
test_ex(x::Dict, y) = @test prettify(to_expr(x)) == prettify(y)

test_wrap(f, value, ex; kwargs...) = test_ex(wrap(f(value); kwargs...), ex)
test_wrap_handle(name, ex) = test_wrap(handle_by_name, name, ex)
test_wrap_struct(name, ex) = test_wrap(struct_by_name, name, ex)
test_wrap_func(name, ex; kwargs...) = test_wrap(func_by_name, name, ex; kwargs...)
test_wrap_bitmask(name, ex) = test_wrap(bitmask_by_name, name, ex)
test_add_constructor(f, name, ex; kwargs...) = test_ex(add_constructor(f(name); kwargs...), ex)
test_add_constructors(f, name, exs; kwargs...) = test_ex.(add_constructors(f(name); kwargs...), exs)
test_struct_add_constructor(args...) = test_add_constructor(struct_by_name, args...)
test_handle_add_constructor(name, ex; kwargs...) = test_ex(first(add_constructors(handle_by_name(name); kwargs...)), ex)
test_handle_add_constructors(args...; kwargs...) = test_add_constructors(handle_by_name, args...; kwargs...)
test_extend_from_vk(name, ex) = test_ex(extend_from_vk(struct_by_name(name)), :(from_vk(T::Type{$(VulkanGen.remove_vk_prefix(name))}, x::$name) = $ex))
test_extend_handle_constructor(name, ex; kwargs...) = test_ex(extend_handle_constructor(create_func_no_batch(handle_by_name(name)); kwargs...), ex)

@testset "Wrapping" begin
    @testset "Handles" begin
        test_wrap_handle(:VkInstance, :(
            mutable struct Instance <: Handle
                vks::VkInstance
                refcount::RefCounter
                destructor
                Instance(vks::VkInstance, refcount::RefCounter) = new(vks, refcount, undef)
            end
        ))

        test_wrap_handle(:VkPhysicalDevice, :(
            mutable struct PhysicalDevice <: Handle
                vks::VkPhysicalDevice
                instance::Instance
                refcount::RefCounter
                destructor
                PhysicalDevice(vks::VkPhysicalDevice, instance::Instance, refcount::RefCounter) = new(vks, instance, refcount, undef)
            end
        ))

        test_wrap_handle(:VkDevice, :(
            mutable struct Device <: Handle
                vks::VkDevice
                physical_device::PhysicalDevice
                refcount::RefCounter
                destructor
                Device(vks::VkDevice, physical_device::PhysicalDevice, refcount::RefCounter) = new(vks, physical_device, refcount, undef)
            end
        ))
    end

    @testset "Structs" begin

        test_wrap_struct(:VkPhysicalDeviceProperties, :(
            struct PhysicalDeviceProperties <: ReturnedOnly
                api_version::VersionNumber
                driver_version::VersionNumber
                vendor_id::UInt32
                device_id::UInt32
                device_type::VkPhysicalDeviceType
                device_name::String
                pipeline_cache_uuid::String
                limits::PhysicalDeviceLimits
                sparse_properties::PhysicalDeviceSparseProperties
            end))

        test_wrap_struct(:VkApplicationInfo, :(
            struct ApplicationInfo <: VulkanStruct{true}
                vks::VkApplicationInfo
                deps::Vector{Any}
            end))

        test_wrap_struct(:VkExtent2D, :(
            struct Extent2D <: VulkanStruct{false}
                vks::VkExtent2D
            end))

        test_wrap_struct(:VkExternalBufferProperties, :(
            struct ExternalBufferProperties <: ReturnedOnly
                s_type::VkStructureType
                p_next::Ptr{Cvoid}
                external_memory_properties::ExternalMemoryProperties
            end
        ))

        test_wrap_struct(:VkPipelineExecutableInternalRepresentationKHR, :(
            struct PipelineExecutableInternalRepresentationKHR <: ReturnedOnly
                s_type::VkStructureType
                p_next::Ptr{Cvoid}
                name::String
                description::String
                is_text::Bool
                data_size::UInt
                p_data::Ptr{Cvoid}
            end
        ))

        test_wrap_struct(:VkDescriptorSetAllocateInfo, :(
            struct DescriptorSetAllocateInfo <: VulkanStruct{true}
                vks::VkDescriptorSetAllocateInfo
                deps::Vector{Any}
                descriptor_pool::DescriptorPool
            end
        ))

        test_wrap_struct(:VkAccelerationStructureBuildGeometryInfoKHR, :(
            struct AccelerationStructureBuildGeometryInfoKHR <: VulkanStruct{true}
                vks::VkAccelerationStructureBuildGeometryInfoKHR
                deps::Vector{Any}
                src_acceleration_structure::OptionalPtr{AccelerationStructureKHR}
                dst_acceleration_structure::OptionalPtr{AccelerationStructureKHR}
            end
        ))
    end

    @testset "API functions" begin
        test_wrap_func(:vkEnumeratePhysicalDevices, :(
            function enumerate_physical_devices(instance::Instance)::Result{Vector{PhysicalDevice},VulkanError}
                pPhysicalDeviceCount = Ref{UInt32}()
                @repeat_while_incomplete begin
                    @check vkEnumeratePhysicalDevices(instance, pPhysicalDeviceCount, C_NULL)
                    pPhysicalDevices = Vector{VkPhysicalDevice}(undef, pPhysicalDeviceCount[])
                    @check vkEnumeratePhysicalDevices(instance, pPhysicalDeviceCount, pPhysicalDevices)
                end
                PhysicalDevice.(pPhysicalDevices, identity, instance)
            end
        ))

        test_wrap_func(:vkGetPhysicalDeviceProperties, :(
            function get_physical_device_properties(physical_device::PhysicalDevice)::PhysicalDeviceProperties
                pProperties = Ref{VkPhysicalDeviceProperties}()
                vkGetPhysicalDeviceProperties(physical_device, pProperties)
                from_vk(PhysicalDeviceProperties, pProperties[])
            end
        ))

        test_wrap_func(:vkEnumerateInstanceVersion, :(
            function enumerate_instance_version()::Result{VersionNumber,VulkanError}
                pApiVersion = Ref{UInt32}()
                @check vkEnumerateInstanceVersion(pApiVersion)
                from_vk(VersionNumber, pApiVersion[])
            end
        ))

        test_wrap_func(:vkEnumerateInstanceExtensionProperties, :(
            function enumerate_instance_extension_properties(; layer_name = C_NULL)::Result{Vector{ExtensionProperties},VulkanError}
                pPropertyCount = Ref{UInt32}()
                @repeat_while_incomplete begin
                    @check vkEnumerateInstanceExtensionProperties(layer_name, pPropertyCount, C_NULL)
                    pProperties = Vector{VkExtensionProperties}(undef, pPropertyCount[])
                    @check vkEnumerateInstanceExtensionProperties(layer_name, pPropertyCount, pProperties)
                end
                from_vk.(ExtensionProperties, pProperties)
            end
        ))

        test_wrap_func(:vkGetGeneratedCommandsMemoryRequirementsNV, :(
            function get_generated_commands_memory_requirements_nv(device::Device, info::GeneratedCommandsMemoryRequirementsInfoNV)::MemoryRequirements2
                pMemoryRequirements = Ref{VkMemoryRequirements2}()
                vkGetGeneratedCommandsMemoryRequirementsNV(device, info, pMemoryRequirements)
                from_vk(MemoryRequirements2, pMemoryRequirements[])
            end
        ))

        test_wrap_func(:vkGetInstanceProcAddr, :(get_instance_proc_addr(name::AbstractString; instance = C_NULL)::FunctionPtr = vkGetInstanceProcAddr(instance, name)))
        test_wrap_func(:vkGetInstanceProcAddr, :(get_instance_proc_addr(name::AbstractString, fun_ptr::FunctionPtr; instance = C_NULL)::FunctionPtr = vkGetInstanceProcAddr(instance, name, fun_ptr)); with_func_ptr=true)

        test_wrap_func(:vkGetPhysicalDeviceSurfacePresentModesKHR, :(
            function get_physical_device_surface_present_modes_khr(physical_device::PhysicalDevice, surface::SurfaceKHR)::Result{Vector{VkPresentModeKHR},VulkanError}
                pPresentModeCount = Ref{UInt32}()
                @repeat_while_incomplete begin
                    @check vkGetPhysicalDeviceSurfacePresentModesKHR(physical_device, surface, pPresentModeCount, C_NULL)
                    pPresentModes = Vector{VkPresentModeKHR}(undef, pPresentModeCount[])
                    @check vkGetPhysicalDeviceSurfacePresentModesKHR(physical_device, surface, pPresentModeCount, pPresentModes)
                end
                pPresentModes
            end
        ))

        test_wrap_func(:vkGetRandROutputDisplayEXT, :(
            function get_rand_r_output_display_ext(physical_device::PhysicalDevice, dpy::vk.Display, rr_output::vk.RROutput)::Result{DisplayKHR,VulkanError}
                pDisplay = Ref{VkDisplayKHR}()
                @check vkGetRandROutputDisplayEXT(physical_device, Ref(dpy), rr_output, pDisplay)
                DisplayKHR(pDisplay[], identity, physical_device)
            end
        ))

        test_wrap_func(:vkRegisterDeviceEventEXT, :(
            function register_device_event_ext(device::Device, device_event_info::DeviceEventInfoEXT; allocator = C_NULL)::Result{Fence,VulkanError}
                pFence = Ref{VkFence}()
                @check vkRegisterDeviceEventEXT(device, device_event_info, allocator, pFence)
                Fence(pFence[], (x->destroy_fence(device, x; allocator)), device)
            end
        ))

        test_wrap_func(:vkCreateInstance, :(
            function create_instance(create_info::InstanceCreateInfo; allocator = C_NULL)::Result{Instance,VulkanError}
                pInstance = Ref{VkInstance}()
                @check vkCreateInstance(create_info, allocator, pInstance)
                Instance(pInstance[], x -> destroy_instance(x; allocator))
            end
        ))

        test_wrap_func(:vkCreateDebugReportCallbackEXT, :(
            function create_debug_report_callback_ext(instance::Instance, create_info::DebugReportCallbackCreateInfoEXT, fun_ptr_create::FunctionPtr, fun_ptr_destroy::FunctionPtr; allocator = C_NULL)::Result{DebugReportCallbackEXT,VulkanError}
                pCallback = Ref{VkDebugReportCallbackEXT}()
                @check vkCreateDebugReportCallbackEXT(instance, create_info, allocator, pCallback, fun_ptr_create)
                DebugReportCallbackEXT(pCallback[], (x->destroy_debug_report_callback_ext(instance, x, fun_ptr_destroy; allocator)), instance)
            end
        ); with_func_ptr=true)

        test_wrap_func(:vkCreateGraphicsPipelines, :(
            function create_graphics_pipelines(device::Device, create_infos::AbstractArray; pipeline_cache = C_NULL, allocator = C_NULL)::Result{Tuple{Vector{Pipeline}, VkResult}, VulkanError}
                pPipelines = Vector{VkPipeline}(undef, pointer_length(create_infos))
                @check vkCreateGraphicsPipelines(device, pipeline_cache, pointer_length(create_infos), create_infos, allocator, pPipelines)
                (Pipeline.(pPipelines, x -> destroy_pipeline(device, x; allocator), device), _return_code)
            end
        ))

        test_wrap_func(:vkAllocateDescriptorSets, :(
            function allocate_descriptor_sets(device::Device, allocate_info::DescriptorSetAllocateInfo)::Result{Vector{DescriptorSet},VulkanError}
                pDescriptorSets = Vector{VkDescriptorSet}(undef, allocate_info.vks.descriptorSetCount)
                @check vkAllocateDescriptorSets(device, allocate_info, pDescriptorSets)
                DescriptorSet.(pDescriptorSets, identity, getproperty(allocate_info, :descriptor_pool))
            end
        ))

        test_wrap_func(:vkMergePipelineCaches, :(
            merge_pipeline_caches(device::Device, dst_cache::PipelineCache, src_caches::AbstractArray)::Result{VkResult,VulkanError} = @check(vkMergePipelineCaches(device, dst_cache, pointer_length(src_caches), src_caches))
        ))

        test_wrap_func(:vkGetFenceFdKHR, :(
            function get_fence_fd_khr(device::Device, get_fd_info::FenceGetFdInfoKHR)::Result{Int,VulkanError}
                pFd = Ref{Int}()
                @check vkGetFenceFdKHR(device, get_fd_info, pFd)
                pFd[]
            end
        ))

        test_wrap_func(:vkGetDeviceGroupPeerMemoryFeatures, :(
            function get_device_group_peer_memory_features(device::Device, heap_index::Integer, local_device_index::Integer, remote_device_index::Integer)::PeerMemoryFeatureFlag
                pPeerMemoryFeatures = Ref{VkPeerMemoryFeatureFlags}()
                vkGetDeviceGroupPeerMemoryFeatures(device, heap_index, local_device_index, remote_device_index, pPeerMemoryFeatures)
                pPeerMemoryFeatures[]
            end
        ))

        test_wrap_func(:vkUpdateDescriptorSets, :(
            update_descriptor_sets(device::Device, descriptor_writes::AbstractArray, descriptor_copies::AbstractArray)::Cvoid = vkUpdateDescriptorSets(device, pointer_length(descriptor_writes), descriptor_writes, pointer_length(descriptor_copies), descriptor_copies)
        ))

        test_wrap_func(:vkCmdSetViewport, :(
            cmd_set_viewport(command_buffer::CommandBuffer, viewports::AbstractArray)::Cvoid = vkCmdSetViewport(command_buffer, 0, pointer_length(viewports), viewports)
        ))

        test_wrap_func(:vkCmdSetLineWidth, :(
            cmd_set_line_width(command_buffer::CommandBuffer, line_width::Real)::Cvoid = vkCmdSetLineWidth(command_buffer, line_width)
        ))

        test_wrap_func(:vkDestroyInstance, :(
            destroy_instance(instance::Instance; allocator = C_NULL)::Cvoid = vkDestroyInstance(instance, allocator)
        ))

        test_wrap_func(:vkMapMemory, :(
            function map_memory(device::Device, memory::DeviceMemory, offset::Integer, size::Integer; flags = 0)::Result{Ptr{Cvoid},VulkanError}
                ppData = Ref{Ptr{Cvoid}}()
                @check vkMapMemory(device, memory, offset, size, flags, ppData)
                ppData[]
            end
        ))

        test_wrap_func(:vkEnumeratePhysicalDeviceQueueFamilyPerformanceQueryCountersKHR, :(
            function enumerate_physical_device_queue_family_performance_query_counters_khr(physical_device::PhysicalDevice, queue_family_index::Integer)::Result{Tuple{Vector{PerformanceCounterKHR}, Vector{PerformanceCounterDescriptionKHR}},VulkanError}
                pCounterCount = Ref{UInt32}()
                @repeat_while_incomplete begin
                    @check vkEnumeratePhysicalDeviceQueueFamilyPerformanceQueryCountersKHR(physical_device, queue_family_index, pCounterCount, C_NULL, C_NULL)
                    pCounters = Vector{VkPerformanceCounterKHR}(undef, pCounterCount[])
                    pCounterDescriptions = Vector{VkPerformanceCounterDescriptionKHR}(undef, pCounterCount[])
                    @check vkEnumeratePhysicalDeviceQueueFamilyPerformanceQueryCountersKHR(physical_device, queue_family_index, pCounterCount, pCounters, pCounterDescriptions)
                end
                (from_vk.(PerformanceCounterKHR, pCounters), from_vk.(PerformanceCounterDescriptionKHR, pCounterDescriptions))
            end
        ))

        test_wrap_func(:vkGetPipelineCacheData, :(
            function get_pipeline_cache_data(device::Device, pipeline_cache::PipelineCache)::Result{Tuple{UInt,Ptr{Cvoid}},VulkanError}
                pDataSize = Ref{UInt}()
                @repeat_while_incomplete begin
                    @check vkGetPipelineCacheData(device, pipeline_cache, pDataSize, C_NULL)
                    pData = Libc.malloc(pDataSize[])
                    @check vkGetPipelineCacheData(device, pipeline_cache, pDataSize, pData)
                    if _return_code == VK_INCOMPLETE
                        Libc.free(pData)
                    end
                end
                (pDataSize[], pData)
            end
        ))

        test_wrap_func(:vkWriteAccelerationStructuresPropertiesKHR, :(
            write_acceleration_structures_properties_khr(device::Device, acceleration_structures::AbstractArray, query_type::VkQueryType, data_size::Integer, data::Ptr{Cvoid}, stride::Integer)::Result{VkResult, VulkanError} = @check vkWriteAccelerationStructuresPropertiesKHR(device, pointer_length(acceleration_structures), acceleration_structures, query_type, data_size, data, stride)
        ))

        test_wrap_func(:vkGetQueryPoolResults, :(
            get_query_pool_results(device::Device, query_pool::QueryPool, first_query::Integer, query_count::Integer, data_size::Integer, data::Ptr{Cvoid}, stride::Integer; flags = 0)::Result{VkResult, VulkanError} = @check vkGetQueryPoolResults(device, query_pool, first_query, query_count, data_size, data, stride, flags)
        ))

        test_wrap_func(:vkGetFenceStatus, :(
            get_fence_status(device::Device, fence::Fence)::Result{VkResult, VulkanError} = @check(vkGetFenceStatus(device, fence))
        ))

        test_wrap_func(:vkGetSwapchainImagesKHR, :(
            function get_swapchain_images_khr(device::Device, swapchain::SwapchainKHR)::Result{Vector{Image}, VulkanError}
                pSwapchainImageCount = Ref{UInt32}()
                @repeat_while_incomplete begin
                        @check vkGetSwapchainImagesKHR(device, swapchain, pSwapchainImageCount, C_NULL)
                        pSwapchainImages = Vector{VkImage}(undef, pSwapchainImageCount[])
                        @check vkGetSwapchainImagesKHR(device, swapchain, pSwapchainImageCount, pSwapchainImages)
                    end
                Image.(pSwapchainImages, identity, device)
            end
        ))
    end

    @testset "Additional constructors" begin
        @testset "Struct constructors" begin
            test_struct_add_constructor(:VkInstanceCreateInfo, :(
                function InstanceCreateInfo(enabled_layer_names::AbstractArray, enabled_extension_names::AbstractArray; next=C_NULL, flags=0, application_info=C_NULL)
                    next = cconvert(Ptr{Cvoid}, next)
                    application_info = cconvert(Ptr{VkApplicationInfo}, application_info)
                    enabled_layer_names = cconvert(Ptr{Cstring}, enabled_layer_names)
                    enabled_extension_names = cconvert(Ptr{Cstring}, enabled_extension_names)
                    deps = [
                        next,
                        application_info,
                        enabled_layer_names,
                        enabled_extension_names
                    ]
                    vks = VkInstanceCreateInfo(VK_STRUCTURE_TYPE_INSTANCE_CREATE_INFO,
                                            unsafe_convert(Ptr{Cvoid}, next),
                                            flags,
                                            unsafe_convert(Ptr{VkApplicationInfo}, application_info),
                                            pointer_length(enabled_layer_names),
                                            unsafe_convert(Ptr{Cstring}, enabled_layer_names),
                                            pointer_length(enabled_extension_names),
                                            unsafe_convert(Ptr{Cstring}, enabled_extension_names),
                                            )
                    InstanceCreateInfo(vks, deps)
                end
            ))

        test_struct_add_constructor(:VkSubpassSampleLocationsEXT, :(
            function SubpassSampleLocationsEXT(subpass_index::Integer, sample_locations_info::SampleLocationsInfoEXT)
                SubpassSampleLocationsEXT(VkSubpassSampleLocationsEXT(subpass_index, sample_locations_info.vks))
            end
        ))

        test_struct_add_constructor(:VkDebugUtilsMessengerCreateInfoEXT, :(
            function DebugUtilsMessengerCreateInfoEXT(message_severity::DebugUtilsMessageSeverityFlagEXT, message_type::DebugUtilsMessageTypeFlagEXT, pfn_user_callback::FunctionPtr; next = C_NULL, flags = 0, user_data = C_NULL)
                next = cconvert(Ptr{Cvoid}, next)
                user_data = cconvert(Ptr{Cvoid}, user_data)
                deps = [next, user_data]
                vks = VkDebugUtilsMessengerCreateInfoEXT(VK_STRUCTURE_TYPE_DEBUG_UTILS_MESSENGER_CREATE_INFO_EXT, unsafe_convert(Ptr{Cvoid}, next), flags, message_severity, message_type, pfn_user_callback, unsafe_convert(Ptr{Cvoid}, user_data))
                DebugUtilsMessengerCreateInfoEXT(vks, deps)
            end
        ))

        test_struct_add_constructor(:VkDescriptorSetAllocateInfo, :(
            function DescriptorSetAllocateInfo(descriptor_pool::DescriptorPool, set_layouts::AbstractArray; next = C_NULL)
                next = cconvert(Ptr{Cvoid}, next)
                set_layouts = cconvert(Ptr{VkDescriptorSetLayout}, set_layouts)
                deps = [next, set_layouts]
                vks = VkDescriptorSetAllocateInfo(VK_STRUCTURE_TYPE_DESCRIPTOR_SET_ALLOCATE_INFO, unsafe_convert(Ptr{Cvoid}, next), descriptor_pool, pointer_length(set_layouts), unsafe_convert(Ptr{VkDescriptorSetLayout}, set_layouts))
                DescriptorSetAllocateInfo(vks, deps, descriptor_pool)
            end
        ))

        test_struct_add_constructor(:VkPipelineShaderStageCreateInfo, :(
            function PipelineShaderStageCreateInfo(stage::ShaderStageFlag, _module::ShaderModule, name::AbstractString; next = C_NULL, flags = 0, specialization_info = C_NULL)
                next = cconvert(Ptr{Cvoid}, next)
                name = cconvert(Cstring, name)
                specialization_info = cconvert(Ptr{VkSpecializationInfo}, specialization_info)
                deps = [next, name, specialization_info]
                vks = VkPipelineShaderStageCreateInfo(VK_STRUCTURE_TYPE_PIPELINE_SHADER_STAGE_CREATE_INFO, unsafe_convert(Ptr{Cvoid}, next), flags, VkShaderStageFlagBits(stage.val), _module, unsafe_convert(Cstring, name), unsafe_convert(Ptr{VkSpecializationInfo}, specialization_info))
                PipelineShaderStageCreateInfo(vks, deps, _module)
            end
        ))

        test_struct_add_constructor(:VkDescriptorImageInfo, :(
            function DescriptorImageInfo(sampler::Sampler, image_view::ImageView, image_layout::VkImageLayout)
                DescriptorImageInfo(VkDescriptorImageInfo(sampler, image_view, image_layout), sampler, image_view)
            end
        ))

        test_struct_add_constructor(:VkDescriptorSetLayoutBinding, :(
            function DescriptorSetLayoutBinding(binding::Integer, descriptor_type::VkDescriptorType, stage_flags::ShaderStageFlag; descriptor_count = 0, immutable_samplers = C_NULL)
                immutable_samplers = cconvert(Ptr{VkSampler}, immutable_samplers)
                deps = [immutable_samplers]
                vks = VkDescriptorSetLayoutBinding(binding, descriptor_type, descriptor_count, stage_flags, unsafe_convert(Ptr{VkSampler}, immutable_samplers))
                DescriptorSetLayoutBinding(vks, deps)
            end
        ))

        end

        @testset "Handle constructors" begin
            test_handle_add_constructor(:VkInstance, :(
                Instance(enabled_layer_names::AbstractArray, enabled_extension_names::AbstractArray; allocator = C_NULL, next=C_NULL, flags=0, application_info=C_NULL) = unwrap(create_instance(enabled_layer_names, enabled_extension_names; allocator, next, flags, application_info))
            ))

            test_handle_add_constructor(:VkDebugReportCallbackEXT, :(
                    DebugReportCallbackEXT(instance::Instance, pfn_callback::FunctionPtr, fun_ptr_create::FunctionPtr, fun_ptr_destroy::FunctionPtr; allocator = C_NULL, next = C_NULL, flags = 0, user_data = C_NULL) =    unwrap(create_debug_report_callback_ext(instance, pfn_callback, fun_ptr_create, fun_ptr_destroy; allocator, next, flags, user_data))
                ),
                with_func_ptr=true,
            )

            test_handle_add_constructor(:VkDeferredOperationKHR, :(
                DeferredOperationKHR(device::Device; allocator = C_NULL) = create_deferred_operation_khr(device; allocator)
            ))

            test_handle_add_constructors(:VkRenderPass, [
                :(
                    RenderPass(device::Device, attachments::AbstractArray, subpasses::AbstractArray, dependencies::AbstractArray; allocator = C_NULL, next = C_NULL, flags = 0) = unwrap(create_render_pass(device, attachments, subpasses, dependencies; allocator, next, flags))
                ),
                :(
                    RenderPass(device::Device, attachments::AbstractArray, subpasses::AbstractArray, dependencies::AbstractArray, correlated_view_masks::AbstractArray; allocator = C_NULL, next = C_NULL, flags = 0) = unwrap(create_render_pass_2(device, attachments, subpasses, dependencies, correlated_view_masks; allocator, next, flags))
                ),
            ])
        end

        @testset "Extended handle constructors" begin
            test_extend_handle_constructor(:VkInstance, :(
                create_instance(enabled_layer_names::AbstractArray, enabled_extension_names::AbstractArray; allocator = C_NULL, next=C_NULL, flags=0, application_info=C_NULL) = create_instance(InstanceCreateInfo(enabled_layer_names, enabled_extension_names; next, flags, application_info); allocator)
            ))

            test_extend_handle_constructor(:VkDebugReportCallbackEXT, :(
                create_debug_report_callback_ext(instance::Instance, pfn_callback::FunctionPtr, fun_ptr_create::FunctionPtr, fun_ptr_destroy::FunctionPtr; allocator = C_NULL, next = C_NULL, flags = 0, user_data = C_NULL) = create_debug_report_callback_ext(instance, DebugReportCallbackCreateInfoEXT(pfn_callback; next, flags, user_data), fun_ptr_create, fun_ptr_destroy; allocator)
                ),
                with_func_ptr=true,
            )
        end
    end

    @testset "New methods for `from_vk`" begin
        test_extend_from_vk(:VkLayerProperties, :(
            T(from_vk(String, x.layerName), from_vk(VersionNumber, x.specVersion), from_vk(VersionNumber, x.implementationVersion), from_vk(String, x.description))
        ))

        test_extend_from_vk(:VkQueueFamilyProperties, :(
            T(x.queueFlags, x.queueCount, x.timestampValidBits, from_vk(Extent3D, x.minImageTransferGranularity))
        ))

        test_extend_from_vk(:VkPhysicalDeviceMemoryProperties, :(
            T(x.memoryTypeCount, from_vk.(MemoryType, x.memoryTypes), x.memoryHeapCount, from_vk.(MemoryHeap, x.memoryHeaps))
        ))

        test_extend_from_vk(:VkDisplayPlaneCapabilities2KHR, :(
            T(x.sType, x.pNext, from_vk(DisplayPlaneCapabilitiesKHR, x.capabilities))
        ))

        test_extend_from_vk(:VkDrmFormatModifierPropertiesListEXT, :(
            T(x.sType, x.pNext, unsafe_wrap(Vector{DrmFormatModifierPropertiesEXT}, x.pDrmFormatModifierProperties, x.drmFormatModifierCount; own = true))
        ))

        test_extend_from_vk(:VkPhysicalDeviceGroupProperties, :(
            T(x.sType, x.pNext, x.physicalDeviceCount, PhysicalDevice.(x.physicalDevices), from_vk(Bool, x.subsetAllocation))
        ))
    end

    @testset "Bitmask flags" begin
        test_wrap_bitmask(:VkQueryPipelineStatisticFlagBits, :(
            @bitmask_flag QueryPipelineStatisticFlag::UInt32 begin
                QUERY_PIPELINE_STATISTIC_INPUT_ASSEMBLY_VERTICES_BIT = 1
                QUERY_PIPELINE_STATISTIC_INPUT_ASSEMBLY_PRIMITIVES_BIT = 2
                QUERY_PIPELINE_STATISTIC_VERTEX_SHADER_INVOCATIONS_BIT = 4
                QUERY_PIPELINE_STATISTIC_GEOMETRY_SHADER_INVOCATIONS_BIT = 8
                QUERY_PIPELINE_STATISTIC_GEOMETRY_SHADER_PRIMITIVES_BIT = 16
                QUERY_PIPELINE_STATISTIC_CLIPPING_INVOCATIONS_BIT = 32
                QUERY_PIPELINE_STATISTIC_CLIPPING_PRIMITIVES_BIT = 64
                QUERY_PIPELINE_STATISTIC_FRAGMENT_SHADER_INVOCATIONS_BIT = 128
                QUERY_PIPELINE_STATISTIC_TESSELLATION_CONTROL_SHADER_PATCHES_BIT = 256
                QUERY_PIPELINE_STATISTIC_TESSELLATION_EVALUATION_SHADER_INVOCATIONS_BIT = 512
                QUERY_PIPELINE_STATISTIC_COMPUTE_SHADER_INVOCATIONS_BIT = 1024
            end
        ))

        test_wrap_bitmask(:VkSparseMemoryBindFlagBits, :(
            @bitmask_flag SparseMemoryBindFlag::UInt32 begin
                SPARSE_MEMORY_BIND_METADATA_BIT = 1
            end
        ))
    end

    @testset "Manual tweaks" begin
        test_struct_add_constructor(:VkWriteDescriptorSet, :(
            function WriteDescriptorSet(dst_set::DescriptorSet, dst_binding::Integer, dst_array_element::Integer, descriptor_type::VkDescriptorType, image_info::AbstractArray, buffer_info::AbstractArray, texel_buffer_view::AbstractArray; next = C_NULL, descriptor_count = max(pointer_length(image_info), pointer_length(buffer_info), pointer_length(texel_buffer_view)))
                next = cconvert(Ptr{Cvoid}, next)
                image_info = cconvert(Ptr{VkDescriptorImageInfo}, image_info)
                buffer_info = cconvert(Ptr{VkDescriptorBufferInfo}, buffer_info)
                texel_buffer_view = cconvert(Ptr{VkBufferView}, texel_buffer_view)
                deps = [next, image_info, buffer_info, texel_buffer_view]
                vks = VkWriteDescriptorSet(VK_STRUCTURE_TYPE_WRITE_DESCRIPTOR_SET, unsafe_convert(Ptr{Cvoid}, next), dst_set, dst_binding, dst_array_element, descriptor_count, descriptor_type, unsafe_convert(Ptr{VkDescriptorImageInfo}, image_info), unsafe_convert(Ptr{VkDescriptorBufferInfo}, buffer_info), unsafe_convert(Ptr{VkBufferView}, texel_buffer_view))
                WriteDescriptorSet(vks, deps, dst_set)
            end
        ))
    end
end


text:
"""
Solar position and radiation functions.

* [`eccentricity`](@ref)
* [`atmos_refrac`](@ref)
* [`solar_angle`](@ref)
"""
module Solar

include("../docstring_style.jl")

import Dates

export eccentricity,
    atmos_refrac,
    solar_angle

"""
Convert `Dates.Millisecond` to number of days (`Float64`).

# Examples

```jldoctest
julia> ms2day(Dates.Millisecond(86400_000))
1.0
```
"""
ms2day(ms::Dates.Millisecond) = ms.value * 1e-3 / 86400.0

"""
Compute the earth's orbital eccentricity. `t` can be either a number or a
`Dates.DateTime` type. If `t` is a number, it is treated as in thousands of
Julian years from J2000.

# References

* [SBC94] Simon, J. L., Bretagnon, P., Chapront, J., Chapront-Touzé, M.,
  Francou, G., and Laskar, J. (1994). *Astron. Astrophys.*, 282, 663--683.
  Retrieved from <http://adsabs.harvard.edu/abs/1994A%26A...282..663S> on
  2018-06-12.

# Examples

```jldoctest
julia> println("eccentricity in 2018: ", eccentricity((2018 - 2000) * 1e-3))
eccentricity in 2018: 0.01670106351746052

julia> println("eccentricity in 2018: ", eccentricity(Dates.DateTime(2018)))
eccentricity in 2018: 0.016701062941387256
```

# Known issues

!!! note

    `eccentricity(dt::Dates.DateTime)` gives more accurate result for a known
    date because the number of days in a leap year differs from that in a
    normal year.
"""
eccentricity(t) = abs(t) <= 1 ?
    0.016_708_634_2 - 4.203_654e-4 * t - 1.267_34e-5 * t * t +
    1.444e-7 * t^3 - 2e-10 * t^4 + 3e-10 * t^5 :
    error("Applicable range: year 1000 -- 3000 AD.")

eccentricity(dt::Dates.DateTime) =
    eccentricity(ms2day(dt - Dates.DateTime(2000)) / 365250.)

"""
Convert the time part of `Dates.DateTime` to the fraction of a day (`Float64`).

# Examples

```jldoctest
julia> partial_day(Dates.DateTime(2017, 6, 5, 11, 15))
0.46875
```
"""
partial_day(dt::Dates.DateTime) = Dates.hour(dt) / 24. +
    Dates.minute(dt) / 1440. +
    (Dates.second(dt) + Dates.millisecond(dt) * 1e-3) / 86400.

"""
Calculate the atmospheric refraction effect from solar elevation angle
(`angle`). Both the argument and the return are in degrees.

# Examples

```jldoctest
julia> atmos_refrac(40.)
0.019200724191321145

julia> atmos_refrac(-1.0)
0.3305949063659434
```
"""
function atmos_refrac(angle)
    if angle > 85.
        return 0.
    elseif 5. < angle <= 85.
        return @evalpoly(
            1. / tand(angle), 0., 58.1, 0., -0.07, 0., 8.6e-5) / 3600.
    elseif -0.575 < angle <= 5.
        return @evalpoly(
            angle, 1735., -518.2, 103.4, -12.79, 0.711) / 3600.
    else
        return -20.774 / (3600. * tand(angle))
    end
end

"""
Calculate solar angles and sunrise/sunset times.

The original program was created by NOAA Global Radiation Group
<http://www.esrl.noaa.gov/gmd/grad/solcalc/calcdetails.html>. Ported to Julia
by Wu Sun on 2017-09-24.

# Parameters

* `dt`: `Dates.DateTime`.
* `lat`: Latitude (-90 -- +90°).
* `lon`: Longitude (-180 -- +180°).
* `timezone` (*optional*): Time zone (-12 -- +12). Default is 0.

# Returns

`solar_angle_results::Dict` --- Unpack the keys to get the results.

* "solar noon": Local solar noon, in fraction of a day.
* "sunrise": Sunrise time, in fraction of a day.
* "sunset": Sunset time, in fraction of a day.
* "sunlight duration": Sunlight duration in fraction of a day.
* "hour angle": Hour angle [deg].
* "solar zenith angle": Solar zenith angle [deg].
* "solar elevation angle": Solar elevation angle [deg].
* "solar azimuth angle": Solar azimuth angle [deg].
* "atmospheric refraction": Atmospheric refraction [deg].
* "corrected solar zenith angle": Solar zenith angle corrected for
  atmospheric refraction [deg].
* "corrected solar elevation angle": solar elevation angle corrected
  for atmospheric refraction [deg].

# Examples

```jldoctest
julia> solar_angle(Dates.DateTime(2017, 9, 24), 34, -118, -7)
Dict{String,Float64} with 11 entries:
  "sunset"                          => 0.782313
  "atmospheric refraction"          => 0.00403875
  "corrected solar elevation angle" => -55.0082
  "solar zenith angle"              => 145.012
  "solar elevation angle"           => -55.0122
  "solar azimuth angle"             => 340.545
  "hour angle"                      => 168.99
  "solar noon"                      => 0.530585
  "sunrise"                         => 0.278857
  "sunlight duration"               => 0.503456
  "corrected solar zenith angle"    => 145.008
```
"""
function solar_angle(dt::Dates.DateTime, lat, lon, timezone=0.)
    # NASA truncated Julian date method: JD = 2440000.5 at 1968-05-24T00:00:00
    timedelta = dt - Dates.DateTime(1968, 5, 24)  # Dates.Millisecond
    julian_day = 2440000.5 + ms2day(timedelta) - timezone / 24.
    julian_century = (julian_day - 2451545.) / 36525.

    # geometric mean longitude of the sun [deg]
    geom_mean_lon_sun = (280.46646 + julian_century *
                         (36000.76983 + julian_century * 0.0003032)) % 360.0
    # geometric mean anomaly of the sun [deg]
    geom_mean_anom_sun = 357.52911 + julian_century *
        (35999.05029 - 0.0001537 * julian_century)
    # eccentricity of the earth orbit
    eccent_earth_orbit = eccentricity(dt - Dates.Millisecond(timezone * 3.6e6))
    # solar equator of center
    sun_eq_ctr = sind(geom_mean_anom_sun) *
        (1.914602 - julian_century * (0.004817 + 1.4e-5 * julian_century)) +
        sind(geom_mean_anom_sun * 2) *
        (0.019993 - 0.000101 * julian_century) +
        sind(geom_mean_anom_sun * 3) * 0.000289  # sind := sin o deg2rad
    # solar true longitude [deg]
    sun_true_lon = geom_mean_lon_sun + sun_eq_ctr
    # solar true anomaly [deg]
    sun_true_anom = geom_mean_anom_sun + sun_eq_ctr
    # solar apparent longitude [deg]
    sun_app_lon = sun_true_lon - 0.00569 - 0.00478 *
        sind(125.04 - 1934.136 * julian_century)
    # mean obliquity of the ecliptic [deg]
    mean_obliq_ecliptic = 23. +
        (26. + (21.448 - julian_century *
                (46.815 + julian_century *
                 (0.00059 - julian_century * 0.001813))) / 60.) / 60.
    # corrected obliquity [deg]
    obliq_corr = mean_obliq_ecliptic + 0.00256 *
        cosd(125.04 - 1934.136 * julian_century)
    # solar declination angle [deg]
    sun_declin = asind(sind(obliq_corr) * sind(sun_app_lon))

    #= DEPRECATED
    # solar rad vector [AU]
    sun_rad_vector = 1.000001018 *
        (1. - eccent_earth_orbit * eccent_earth_orbit) /
        (1. + eccent_earth_orbit * cosd(sun_true_anom))
    # cosd := cos o deg2rad
    sun_rt_ascen = atand(cosd(obliq_corr) * sind(sun_app_lon) /
                         cosd(sun_app_lon))   # atand := rad2deg o atan
    =#

    var_y = let
        y = tand(obliq_corr * 0.5)
        y * y
    end
    # equation of time [minute]
    eq_time = 4. * rad2deg(
        var_y * sind(2. * geom_mean_lon_sun) -
        2. * eccent_earth_orbit * sind(geom_mean_anom_sun) +
        4. * eccent_earth_orbit * var_y * sind(geom_mean_anom_sun) *
        cosd(2. * geom_mean_lon_sun) -
        0.5 * var_y * var_y * sind(4. * geom_mean_lon_sun) -
        1.25 * eccent_earth_orbit * eccent_earth_orbit *
        sind(2. * geom_mean_anom_sun))
    HA_sunrise = acosd(cosd(90.833) / (cosd(lat) * cosd(sun_declin)) -
                       tand(lat) * tand(sun_declin))
    # local solar noon [day]
    solar_noon_local = (720. - 4. * lon - eq_time + timezone * 60.) / 1440.
    # local sunrise [day]
    sunrise_local = solar_noon_local - HA_sunrise * 4. / 1440.
    # local sunset [day]
    sunset_local = solar_noon_local + HA_sunrise * 4. / 1440.
    # sunlight duration [day]
    sunlight_duration = 8. * HA_sunrise / 1440.
    # true solar time [minute]
    true_solar_time_minute = mod(
        partial_day(dt) * 1440. + eq_time + 4. * lon - 60. * timezone, 1440.)
    # hour angle [deg]
    hour_angle = true_solar_time_minute / 4. +
        (true_solar_time_minute < 0. ? 180. : -180.)
    # solar zenith angle [deg]
    solar_zenith_angle = acosd(sind(lat) * sind(sun_declin) +
                               cosd(lat) * cosd(sun_declin) * cosd(hour_angle))
    # solar elevation angle [deg]
    solar_elev_angle = 90. - solar_zenith_angle
    # solar azimuth angle [deg]
    solar_azimuth_angle = let
        x = acosd((sind(lat) * cosd(solar_zenith_angle) - sind(sun_declin)) /
                  (cosd(lat) * sind(solar_zenith_angle)))
        hour_angle > 0 ? mod(x + 180., 360.) : mod(540 - x, 360.)
    end
    # approximate atmospheric refraction effect [deg]
    approx_atmos_refrac = atmos_refrac(solar_elev_angle)
    # corrected zenith angle and elevation angles [deg]
    solar_zenith_angle_corr = solar_zenith_angle - approx_atmos_refrac
    solar_elev_angle_corr = solar_elev_angle + approx_atmos_refrac

    solar_angle_results = Dict(
        "solar noon" => solar_noon_local,
        "sunrise" => sunrise_local,
        "sunset" => sunset_local,
        "sunlight duration" => sunlight_duration,
        "hour angle" => hour_angle,
        "solar zenith angle" => solar_zenith_angle,
        "solar elevation angle" => solar_elev_angle,
        "solar azimuth angle" => solar_azimuth_angle,
        "atmospheric refraction" => approx_atmos_refrac,
        "corrected solar zenith angle" => solar_zenith_angle_corr,
        "corrected solar elevation angle" => solar_elev_angle_corr)
    return solar_angle_results
end

end  # module


text:
facts("CANDECOMP") do
    r = 2
    T = rand_kruskal3(r, (10, 20, 30), false)
    initial_guess = tuple([randn(dim, r) for dim in size(T)]...)

    context("ALS (Alternating least squares)") do
        @time factors = candecomp(T, r, initial_guess, tol=1e-6, compute_error=true, method=:ALS)
        @fact length(factors.factors) --> ndims(T)
        @fact map(size, factors.factors) --> (collect(zip(size(T), (r, r, r)))...)
        @fact rank(factors) --> r
        @fact rel_residue(factors) --> less_than(1e-5)
    end

    context("SGSD (Simultaneous generalized Schur decomposition)") do
        @time factors = candecomp(T, r, initial_guess, compute_error=true, method=:SGSD)
        @fact length(factors.factors) --> ndims(T)
        @fact map(size, factors.factors) --> (collect(zip(size(T), (r, r, r)))...)
        @fact rank(factors) --> r
        @fact length(factors.lambdas) --> r
        @fact rel_residue(factors) --> less_than(1e-5)
    end

end