text:
module RobustStatsTutorial
#--------------------------------------------
# SOURCE CODE FOR ROBUST STATISTICS TUTORIAL
#--------------------------------------------

#Draft order
#Plot of data from t-distribution with DoF = 2. Set specific seed. Show how much better we can do by eliminating a couple of observations.
#Simulate empirical densities of sample means, sample medians, and sample trimmed means, for t-distribution with DoF = 2
#Discuss fat-tails of t-distribution with DoF = 2. Explain why kurtosis, although the common characterisation, is not an ideal measure of how fat the tails are. Possibly plot sample kurtosis as a function of sample size.
#Introduce alternative robust measure of fat-tails.
#Provide measures of fat-tails for t-distribution with DoF = 2.
#Provide measures of fat-tails for unconditional risky asset returns.
#Show how fat-tails in unconditional risky asset returns can be duplicated using a conditionally Normal model, i.e. heteroskedasticity in Normal sequence
#Something to do with estimating first moment of risky asset returns (?)
#Move onto OLS model for risky asset returns
#Compare least squares estimated forecast model for risky asset returns with robust loss estimation method
#Also compare to case where returns are standardised for volatility first
#Might need to compare to other robust methods of estimation too, ie ridge regression, lasso (? not robust...)

#Load relevant modules
using Base.Dates, Compose, Gadfly, StatsBase, KernelDensity, Distributions, DependentBootstrap, NLopt
#WARNING: DependentBootstrap is not currently an official package (although I am the author). To install, use:
#pkg.clone("https://github.com/colintbowers/DependentBootstrap.jl.git")

#Fixed output directory based on Linux OS. Will need to be adjusted for Windows or Mac users.
const outputDir = "/home/"*ENV["USER"]*"/robust_stats_tutorial_output/"::ASCIIString
const dataDir = "/home/"*ENV["USER"]*"/.julia/v"*string(VERSION.major)*"."*string(VERSION.minor)*"/RobustStatsTutorial/data/"::ASCIIString

#Fixed theme override for all plots and colour order for multivariate plots
const defaultThemeOverride = Theme(key_title_font_size=14pt, key_label_font_size=14pt, minor_label_font_size=14pt, major_label_font_size=17pt, key_title_font_size=17pt, default_color=color("blue"), default_point_size=4pt)::Theme
const colourVec = ["blue", "green", "red", "black", "purple", "dark blue", "darkgreen", "gray", "brown", "cyan",
				   "violetred", "blue2", "orange", "green2", "darkred", "gray20", "chocolate", "brown2", "darkorange", "gray40",
				   "cadetblue", "violet", "green4", "brown4", "blue4"]::Vector{ASCIIString} #A vector of 25 colours to use in plots

#Security list used throughout the tutorial
const secList = ["AMP", "ANZ", "BHP", "CBA", "CCL", "JBH", "LLC", "NAB", "RIO", "SUN", "TLS", "TOL", "WBC", "WES", "WOW"]::Vector{ASCIIString}

#Create output directory
!isdir(outputDir) && mkdir(outputDir)

#--------------------------------------------
# SIMPLE EXAMPLE WITH t-DISTRIBUTION
#--------------------------------------------
#Plot of data from t-distribution with DoF = 2. Set specific seed. Show how much better we can do by eliminating a couple of observations.
function t_dist_simple_example()
    numObs = 20 #Number of observations to use in this subsection
    tDist = TDist(2) #Initiate t-distribution with 2 degrees of freedom
    srand(78) #Specify a range that results in a good dataset for demonstrating the point (yes I mined this from the first 100 integers)
    tData = rand(tDist, numObs) #Simulate iid from t-distribution
    tDataSort = sort(tData)
    tDataMean = mean(tData)
    p = 0.2 #Set a trimming proportion
    tDataTrimMean = tmean(tDataSort, p, sorted=true) #Get trimmed mean
	println("tDist data sample mean = " * string(tDataMean))
	println("tDist data trim mean = " * string(tDataTrimMean))
    (numLower, numUpper) = tmean_num_cut(tDataSort, p/2, p/2) #Get number of observations cut below and above
    yLower = mean(tDataSort[numLower:numLower+1]) #Get the midpoint of the last lower cut observation and the first lower kept observation
    yUpper = mean(tDataSort[end-numUpper:end-numUpper+1]) #Get the midpoint of the last upper kept observation and the first upper cut observation
    println("Plotting data")
    dataPlot1 = plot(x=collect(1:numObs), y=tData, Geom.point, defaultThemeOverride)
    dataPlot2 = plot(x=collect(1:numObs), y=tData, yintercept=[yLower, yUpper], Geom.point, Geom.hline, defaultThemeOverride)
    draw_local(dataPlot1, "t_Dist_Data", dirPath=outputDir, fileType=:svg)
    draw_local(dataPlot2, "t_Dist_Data_With_Trim_Cutoff", dirPath=outputDir, fileType=:svg)
    println("Routine complete")
end

#--------------------------------------------
# ESTIMATOR DENSITIES WITH t-DISTRIBUTION
#--------------------------------------------
#Simulate empirical densities of sample means, sample medians, and sample trimmed means, for t-distribution with DoF = 2
function t_dist_estimator_densities( ; numIter::Int=5000, numObs::Int=50, tMeanProp1::Float64=0.1, tMeanProp2::Float64=0.5)
    tDist = TDist(2) #Initiate t-distribution with 2 degrees of freedom
    estVecVec = Vector{Float64}[ Array(Float64, numIter) for k = 1:4 ] #Pre-allocate structure to hold sample estimators
    println("Simulating estimators")
    for n = 1:numIter
        #Simulate t-distribution data, and then get mean, two trimmed means, and the median for each dataset
        tData = rand(tDist, numObs)
        tDataSort = sort(tData)
        estVecVec[1][n] = mean(tData)
        estVecVec[2][n] = tmean(tDataSort, tMeanProp1, sorted=true)
        estVecVec[3][n] = tmean(tDataSort, tMeanProp2, sorted=true)
        estVecVec[4][n] = median_sorted(tDataSort)
    end
    println("Building kernel densities")
    kDVec = KernelDensity.UnivariateKDE{FloatRange{Float64}}[ kde(estVecVec[k], boundary=(-1.0, 1.0)) for k = 1:4 ] #Get kernel density for each estimator
    layerVec = Vector{Gadfly.Layer}[ layer(x=collect(kDVec[k].x), y=kDVec[k].density, Geom.line, adjust_default_theme_color(defaultThemeOverride, colourVec[k])) for k = 1:4 ] #Construct plot layers using kernel densities
    println("Plotting kernel densities")
    kernelPlot1 = plot(layerVec..., Guide.xlabel("Estimator value"), Guide.ylabel("Density"), Guide.manual_color_key(default_legend(["mean", string(tMeanProp1)*" trimmed mean", string(tMeanProp2)*" trimmed mean", "median"])...), defaultThemeOverride)
    draw_local(kernelPlot1, "t_Dist_Estimator_Densities", dirPath=outputDir, fileType=:svg)
    println("Routine complete")
end

#--------------------------------------------
# MEASURING TAIL FATNESS
#--------------------------------------------
#Discuss fat-tails of t-distribution with DoF = 2. Explain why kurtosis, although the common characterisation, is not an ideal measure of how fat the tails are. Possibly plot sample kurtosis as a function of sample size.
#Introduce alternative robust measure of fat-tails and show it is stable relative to sample kurtosis
#Specifically, show that sample kurtosis diverges as sample size increases if population kurtosis is undefined
function tail_fatness_and_kurtosis( ; numIter::Int=100, numObsVec::Vector{Int}=collect(20:20:200))
    #Simulate the kurtosis and robust kurtosis as a function of sample size for t-distribution with 2 degrees of freedom
	tDist = TDist(2) #Initiate t-distribution with 2 degrees of freedom
    simKurtosisMat = Array(Float64, numIter, length(numObsVec))
	simRobustKurtosisMat = Array(Float64, numIter, length(numObsVec))
    for k = 1:length(numObsVec)
		for m = 1:numIter
			tData = rand(tDist, numObsVec[k])
        	simKurtosisMat[m, k] = kurtosis(tData)
			simRobustKurtosisMat[m, k] = hogg_robust_kurt!(tData, sorted=false, numerTail=0.05, denomTail=0.5)
		end
    end
	simKurtosis = vec(mean(simKurtosisMat, 1))
	simRobustKurtosis = vec(mean(simRobustKurtosisMat, 1))
    estVec = Vector{Float64}[simKurtosis, simRobustKurtosis]
    layerVec = Vector{Gadfly.Layer}[ layer(x=numObsVec, y=estVec[k], Geom.line, adjust_default_theme_color(defaultThemeOverride, colourVec[k])) for k = 1:2 ] #Construct plot layers for the two estimators
    estPlot1 = plot(layerVec..., Guide.xlabel("Number of observations"), Guide.ylabel("Simulated average estimator value"), Guide.manual_color_key(default_legend(["Sample kurtosis", "Robust measure"])...), defaultThemeOverride)
    println("Plotting estimators...")
    draw_local(estPlot1, "t_Dist_Kurtosis_Versus_Robust_Kurtosis", dirPath=outputDir, fileType=:svg)
	println("Routine complete")
end

#--------------------------------------------
# TAIL FATNESS OF FINANCIAL DATA VERSUS NORMAL
#--------------------------------------------
#Compare robust measure of fat-tails for lots of different stocks to value under Normal distribution and value under t-distribution with 2 DoF
#Discuss how we can generate unconditional fat-tails using a conditional Normal model with time-varying variance. Plot robust measure of fat-tails for returns standardised by realised variance type estimator. Show how it is close to Normal.
#Conclude by comparing mean of financial returns to trimmed mean and median.
function tail_fatness_financial_data( ; scaleMethod::Symbol=:historicvariance, blockLength::Float64=4.0)
	#Get robust measure of return data and compare it to Normal and t-Dist with 2 DoF
	secRet = read_local(secList, :return)
    hoggEst = Float64[ hogg_robust_kurt(secRet[j]) for j = 1:length(secList) ]
    tDist = TDist(2)
    tDistHoggEst = mean(Float64[ hogg_robust_kurt(rand(tDist, 100)) for k = 1:1000 ])
    normalHoggEst = mean(Float64[ hogg_robust_kurt(randn(100)) for k = 1:1000 ])
	println("Drawing plot 1")
	estPlot1 = plot(x=secList, y=hoggEst, yintercept=[normalHoggEst, tDistHoggEst], Geom.point, Geom.hline, Guide.xlabel("Ticker code"), Guide.ylabel("Robust kurtosis"), defaultThemeOverride)
    draw_local(estPlot1, "Robust_Kurtosis_of_Daily_Financial_Returns", dirPath=outputDir, fileType=:svg)
	#Standardise returns by a variance estimate and then repeat above exercise
	if scaleMethod == :historicvariance
		secVar = Vector{Float64}[ historical_variance(secRet[k]) for k = 1:length(secRet) ]
		secRetStd = Vector{Float64}[ secRet[j][end-length(secVar[j])+1:end] ./ sqrt(secVar[j]) for j = 1:length(secList) ]
	else
		error("Invalid scaleMethod")
	end
	hoggEstStd = Float64[ hogg_robust_kurt(secRetStd[j]) for j = 1:length(secList) ]
	println("Drawing plot 2")
	estPlot2 = plot(x=secList, y=hoggEstStd, yintercept=[normalHoggEst, tDistHoggEst], Geom.point, Geom.hline, Guide.xlabel("Ticker code"), Guide.ylabel("Robust kurtosis"), defaultThemeOverride)
    draw_local(estPlot2, "Robust_Kurtosis_of_Standardised_Daily_Financial_Returns", dirPath=outputDir, fileType=:svg)
	#Compare mean and trimmed mean on resampled financial returns (use NAB (fattest tails) and CCL (thinnest tails))
	println("Drawing mean comparison plots (full period)")
	pNABFull = kernel_plot_robust_mean_versus_trimmed_mean(secRet, "NAB", blockLength=blockLength)
	draw_local(pNABFull, "NAB_Bootstrapped_Location_Estimator_Density_(full_period)", dirPath=outputDir, fileType=:svg)
	pCCLFull = kernel_plot_robust_mean_versus_trimmed_mean(secRet, "CCL", blockLength=blockLength)
	draw_local(pCCLFull, "CCL_Bootstrapped_Location_Estimator_Density_(full_period)", dirPath=outputDir, fileType=:svg)
	#Compare mean and trimmed mean on short horizon financial returns for large robust kurtosis versus small robust kurtosis
	println("Estimating rolling robust kurtosis")
	rNAB = secRet[find(secList .== "NAB")[1]]
	rollingHoggNAB = Float64[ hogg_robust_kurt(rNAB[j:j+99]) for j = 1:length(rNAB)-99 ]
	iSmlHogg = indmin(rollingHoggNAB)
	iLrgHogg = indmax(rollingHoggNAB)
	rNABSml = rNAB[iSmlHogg:iSmlHogg+99]
	rNABLrg = rNAB[iLrgHogg:iLrgHogg+99]
	println("Drawing mean comparison plots (short horizon)")
	pNABSml = kernel_plot_robust_mean_versus_trimmed_mean(rNABSml, secStr="NAB Short Horizon (Robust kurt = " * string(rollingHoggNAB[iSmlHogg])[1:3] * ")", blockLength=blockLength)
	pNABLrg = kernel_plot_robust_mean_versus_trimmed_mean(rNABLrg, secStr="NAB Short Horizon (Robust kurt = " * string(rollingHoggNAB[iLrgHogg])[1:3] * ")", blockLength=blockLength)
	draw_local(pNABSml, "NAB_Bootstrapped_Location_Estimator_Density_(short_horizon_thin_tail)", dirPath=outputDir, fileType=:svg)
	draw_local(pNABLrg, "NAB_Bootstrapped_Location_Estimator_Density_(short_horizon_fat_tail)", dirPath=outputDir, fileType=:svg)
	println("Routine complete")
end
function kernel_plot_robust_mean_versus_trimmed_mean(secRet::Vector{Vector{Float64}}, secStr::String ; numResample::Int=1000, blockLength::Float64=4.0)
	iSec = find(secList .== secStr)
	length(iSec) != 1 && error("Unable to find " * secStr * " data")
	r = secRet[iSec[1]]
	return(kernel_plot_robust_mean_versus_trimmed_mean(r, secStr=secStr, numResample=numResample, blockLength=blockLength))
end
function kernel_plot_robust_mean_versus_trimmed_mean(r::Vector{Float64} ; numResample::Int=1000, blockLength::Float64=4.0, secStr::String="", bpMult::Int=10000)
	rCent = r - mean(r) #Centred so re-sampled data has true mean of zero
	rBoot = dbootstrapdata(rCent, blockLength=blockLength, numResample=numResample)
	bootEst = Vector{Float64}[Float64[ bpMult*mean(rBoot[:, m]) for m = 1:numResample ], Float64[ bpMult*tmean(rBoot[:, m], 0.4) for m = 1:numResample ]]
	kDVec = KernelDensity.UnivariateKDE{FloatRange{Float64}}[ kde(bootEst[k]) for k = 1:2 ]
	layerVec = Vector{Gadfly.Layer}[ layer(x=collect(kDVec[k].x), y=kDVec[k].density, Geom.line, adjust_default_theme_color(defaultThemeOverride, colourVec[k])) for k = 1:2 ]
	secStr == "" ? (localTitle = "Location estimator densities") : (localTitle = "Location estimator densities for " * secStr * " return data")
	kernelPlot = plot(layerVec..., Guide.xlabel("Estimator value (basis points)"), Guide.ylabel("Density"), Guide.title(localTitle), Guide.manual_color_key(default_legend(["Mean", "Trimmed mean (0.4)"])...), defaultThemeOverride)
end

#-------------------------------------------
# ROBUST PREDICTION OF FINANCIAL RETURNS
#-------------------------------------------
function return_prediction_financial_data( ; blockLength::Float64=4.0, numResample::Int=1000, varProxyMethod::Symbol=:historicvariance, numObs::Int=2000, rSquared::Float64=0.05, a::Float64=0.0, lag::Int=17, tCost::Float64=0.0, trimSimTail::Float64=0.0, numAsset::Int=20, tradingSimRSquared::Float64=0.01, includeConstantInSim::Bool=true)
	#Get variance proxy
	println("Reading data")
	if varProxyMethod == :realisedvariance
		v = read_local(["NAB"], :realisedvariance)[1]
	elseif varProxyMethod == :historicvariance
		rNAB = read_local(["NAB"], :return)[1]
		v = historical_variance(rNAB, lag)
	else
		error("Invalid varProxyMethod")
	end
	#Simuate returns, signal, error term, constant and coefficient
	println("Simulating and bootstrapping")
	(r, s, e, a, b) = simulate_signal_and_returns(v, rSquared=rSquared, a=a)
	#Resample r, x, and e
	inds = dbootstrapindex(length(r), blockLength, numResample=numResample)
	rBoot = r[inds]
	sBoot = s[inds]
	eBoot = e[inds]
	hoggVec = Float64[ hogg_robust_kurt(rBoot[:, m]) for m = 1:numResample ]
	qProb = Float64[0.05, 0.25, 0.5, 0.75, 0.95]
	hoggVecQ = quantile(hoggVec, qProb)
	#Run regressions on the resampled data using multiple methods and get plot of estimated coefficients
	println("Estimating coefficients on bootstrapped data")
	coefLS = Array(Vector{Float64}, numResample)
	coefLAD = Array(Vector{Float64}, numResample)
	rSqLS = Array(Float64, numResample)
	rSqLAD = Array(Float64, numResample)
	rHatLS = Array(Float64, length(r), numResample)
	rHatLAD = Array(Float64, length(r), numResample)
	for m = 1:numResample
		y = rBoot[:, m]
		x = [ones(Float64, length(y)) sBoot[:, m]]
		coefLS[m] = x \ y
		(coefLAD[m], _, _) = least_absolute_deviation(y, x)
		rHatLS[:, m] = x * coefLS[m]
		rHatLAD[:, m] = x * coefLAD[m]
		eLS = y - rHatLS[:, m]
		eLAD = y - rHatLAD[:, m]
		tSS = sumabs2(y)
		rSqLS[m] = 1 - sumabs2(eLS) / tSS
		rSqLAD[m] = 1 - sumabs2(eLAD) / tSS
	end
	#Get kernel density plot of estimated coefficients (and include true value)
	println("Building coefficient kernel densities")
	aArr = Vector{Float64}[Float64[ coefLS[m][1] for m = 1:numResample ], Float64[ coefLAD[m][1] for m = 1:numResample ]]
	bArr = Vector{Float64}[Float64[ coefLS[m][2] for m = 1:numResample ], Float64[ coefLAD[m][2] for m = 1:numResample ]]
	aKDVec = KernelDensity.UnivariateKDE{FloatRange{Float64}}[ kde(aArr[k]) for k = 1:2 ]
	bKDVec = KernelDensity.UnivariateKDE{FloatRange{Float64}}[ kde(bArr[k]) for k = 1:2 ]
	aLayerVec = Vector{Gadfly.Layer}[ layer(x=collect(aKDVec[k].x), y=aKDVec[k].density, Geom.line, adjust_default_theme_color(defaultThemeOverride, colourVec[k])) for k = 1:2 ]
	bLayerVec = Vector{Gadfly.Layer}[ layer(x=collect(bKDVec[k].x), y=bKDVec[k].density, Geom.line, adjust_default_theme_color(defaultThemeOverride, colourVec[k])) for k = 1:2 ]
	aKernelPlot = plot(aLayerVec..., Guide.xlabel("Regression constant"), Guide.ylabel("Density"), Guide.title("Density of regression constant (true value = " * string(a) * ")"), Guide.manual_color_key(default_legend(["Least squares", "Least absolute deviations"])...), defaultThemeOverride)
	bKernelPlot = plot(bLayerVec..., Guide.xlabel("Regression coefficient"), Guide.ylabel("Density"), Guide.title("Density of regression coefficient (true value = " * string(b) * ")"), Guide.manual_color_key(default_legend(["Least squares", "Least absolute deviations"])...), defaultThemeOverride)
	println("Plotting coefficient kernel densities")
	draw_local(aKernelPlot, "Regression_constant_density", dirPath=outputDir, fileType=:svg)
	draw_local(bKernelPlot, "Regression_coefficient_density", dirPath=outputDir, fileType=:svg)
	#Get kernel density plot of r-squares
	println("Building rSquared kernel density")
	rSqArr = Vector{Float64}[rSqLS, rSqLAD]
	rSqKDVec = KernelDensity.UnivariateKDE{FloatRange{Float64}}[ kde(rSqArr[k]) for k = 1:2 ]
	rSqLayerVec = Vector{Gadfly.Layer}[ layer(x=collect(rSqKDVec[k].x), y=rSqKDVec[k].density, Geom.line, adjust_default_theme_color(defaultThemeOverride, colourVec[k])) for k = 1:2 ]
	rSqKernelPlot = plot(rSqLayerVec..., Guide.xlabel("R-squared"), Guide.ylabel("Density"), Guide.title("Density of regression R-squared (true value = " * string(rSquared) * ")"), Guide.manual_color_key(default_legend(["Least squares", "Least absolute deviations"])...), defaultThemeOverride)
	println("Plotting rSquared kernel density")
	draw_local(rSqKernelPlot, "Regression_rSquared_density", dirPath=outputDir, fileType=:svg)
	println("Hogg robust estimator quantiles (" * string(qProb) * ") = " * string(hoggVecQ))
	#Perform simulations using the estimated models and plot density of terminal values
	println("Performing simple trading simulations")
	println("Total number of observations = " * string(length(v)))
	iEstEnd = Int(round(length(v) / 2))
	terminalValLS = Array(Float64, numResample)
	terminalValLAD = Array(Float64, numResample)
	avgAbsErrorLS = Array(Float64, numResample)
	avgAbsErrorLAD = Array(Float64, numResample)
	for m = 1:numResample
		rMat = Array(Float64, length(r), numAsset)
		sMat = Array(Float64, length(r), numAsset)
		eMat = Array(Float64, length(r), numAsset)
		aVec = Array(Float64, numAsset)
		bVec = Array(Float64, numAsset)
		for j = 1:numAsset
			(rMat[:, j], sMat[:, j], eMat[:, j], aVec[j], bVec[j]) = simulate_signal_and_returns(v, rSquared=tradingSimRSquared, a=0.0)
		end
		coefLSVec = Array(Vector{Float64}, numAsset)
		coefLADVec = Array(Vector{Float64}, numAsset)
		rHatLSMat = Array(Float64, length(r), numAsset)
		rHatLADMat = Array(Float64, length(r), numAsset)
		for j = 1:numAsset
			y = rMat[:, j]
			includeConstantInSim ? (x = [ones(Float64, length(y)) sMat[:, j]]) : (x = sMat[:, j]'')
			yEst = y[1:iEstEnd]
			xEst = x[1:iEstEnd, :]
			ySim = y[iEstEnd+1:end]
			xSim = x[iEstEnd+1:end, :]
			coefLSVec[j] = xEst \ yEst
			(coefLADVec[j], _, _) = least_absolute_deviation(yEst, xEst)
			any(isnan(coefLADVec[j])) && (coefLADVec[j] = coefLSVec[j])
			rHatLSMat[:, j] = x * coefLSVec[j]
			rHatLADMat[:, j] = x * coefLADVec[j]
		end
		avgAbsErrorLS[m] = mean(abs(rHatLSMat[iEstEnd+1:end, :] - rMat[iEstEnd+1:end, :]))
		avgAbsErrorLAD[m] = mean(abs(rHatLADMat[iEstEnd+1:end, :] - rMat[iEstEnd+1:end, :]))
		terminalValLS[m] = very_simple_trading_sim(rMat[iEstEnd+1:end, :], rHatLSMat[iEstEnd+1:end, :])[end]
		terminalValLAD[m] = very_simple_trading_sim(rMat[iEstEnd+1:end, :], rHatLADMat[iEstEnd+1:end, :])[end]
	end
	println("Mean absolute forecast error (LS) = " * string(mean(avgAbsErrorLS)))
	println("Mean absolute forecast error (LAD) = " * string(mean(avgAbsErrorLAD)))
	println("Proportion of samples when LAD has smaller absolute forecast error = " * string(sum(avgAbsErrorLAD .< avgAbsErrorLS) / numResample))
	println("Mean terminal value (LS) = " * string(mean(terminalValLS)))
	println("Mean terminal value (LAD) = " * string(mean(terminalValLAD)))
	println("Median terminal value (LS) = " * string(median(terminalValLS)))
	println("Median terminal value (LAD) = " * string(median(terminalValLAD)))
	sort!(terminalValLS)
	sort!(terminalValLAD)
	iLower = Int(floor(trimSimTail * length(terminalValLS)))
	iLower == 0 && (iLower = 1)
	iUpper = length(terminalValLS) - iLower + 1
	terminalValLS = terminalValLS[iLower:iUpper]
	terminalValLAD = terminalValLAD[iLower:iUpper]
	println("Plotting terminal value kernel densities")
	tV = Vector{Float64}[terminalValLS, terminalValLAD]
	tVKDVec = KernelDensity.UnivariateKDE{FloatRange{Float64}}[ kde(tV[k]) for k = 1:2 ]
	tVLayerVec = Vector{Gadfly.Layer}[ layer(x=collect(tVKDVec[k].x), y=tVKDVec[k].density, Geom.line, adjust_default_theme_color(defaultThemeOverride, colourVec[k])) for k = 1:2 ]
	tVKernelPlot = plot(tVLayerVec..., Guide.xlabel("Terminal value"), Guide.ylabel("Density"), Guide.title("Density of terminal values from simple simulations"), Guide.manual_color_key(default_legend(["Least squares", "Least absolute deviations"])...), defaultThemeOverride)
	draw_local(tVKernelPlot, "Simulated_portfolio_terminal_values", dirPath=outputDir, fileType=:svg)
	println("Routine complete")
end
function simulate_signal_and_returns(v::Vector{Float64} ; rSquared::Float64=0.05, a::Float64=0.0)
	!(0.0 < rSquared < 1.0) && error("Invalid rSquared")
	e = sqrt(v) .* randn(length(v))
	x = randn(length(v))
	residSumSq = sumabs2(e)
	totalSumSq = residSumSq / (1 - rSquared)
	explainSumSq = totalSumSq - residSumSq
	b = sqrt(explainSumSq / sumabs2(x)) #This ensures sample rSquared will equal input rSquared (Note, sample rSquared will become population rSquared for re-sampled data)
	#b = sqrt(explainSumSq / length(v)) #This alternative would ensure population rSquared will equal input rSquared
	r = a + b*x + e
	return(r, x, e, a, b)
end
function least_absolute_deviation(y::Vector{Float64}, x::Matrix{Float64})
	#Build Opt object accepted by NLOpt
	fObj = ((param, grad) -> sumabs(y - x*param)) #Get local anonymous objective function
	initVal = zeros(Float64, size(x, 2))
	opt = Opt(:LN_COBYLA, length(initVal)) #Derivative-free convergence method (yes I'm being lazy here)
	xtol_rel!(opt, 1e-6)
	min_objective!(opt, fObj)
	#Perform optimisation
	(objFuncOpt, paramOpt, flag) = optimize(opt, initVal)
	checkNLoptFlag(flag)
	return(paramOpt, objFuncOpt, y - x*paramOpt) #Return optimal parameter, objective function at optimum, and estimated errors
end
function checkNLoptFlag(flag::Symbol)
	flag == :SUCCESS && return(true)
	flag == :XTOL_REACHED && return(true)
	if flag == :ROUNDOFF_LIMITED
		println("WARNING: Optimisation terminated due to floating point limit. Final iteration parameters used.")
		return(true)
	end
	flag == :MAXEVAL_REACHED && error("Maximum number of evaluations reached in NLopt")
	flag == :MAXTIME_REACHED && error("Maximum optimisation time reached in NLopt")
	flag == :FAILURE && error("Generic failure in NLopt")
	flag == :INVALID_ARGS && error("Invalid arguments to NLopt")
	flag == :OUT_OF_MEMORY && error("NLopt ran out of memory")
	flag == :FORCED_STOP && error("NLopt halted due to forced termination")
	error("Unexpected flag from NLopt. Flag = " * string(flag))
end
#Function for performing a very simple trading simulation. Either hold cash, or trade an asset if forecast return > 0. tCost is in basis points
function very_simple_trading_sim(r::Vector{Float64}, f::Vector{Float64} ; tCost::Float64=50.0)
	length(f) != length(r) && error("Length mismatch")
	N = length(f)
	portValue = Array(Float64, N)
	portValue[1] = 1000000.0
	for n = 2:N
		if f[n] > 0.0
			portValue[n] = portValue[n-1] * (r[n] + 1) - (1/10000) * tCost * portValue[n-1] #notional transaction cost
		else
			portValue[n] = portValue[n-1]
		end
	end
	return(portValue)
end
function very_simple_trading_sim(r::Matrix{Float64}, f::Matrix{Float64})
	size(r) != size(f) && error("Size mismatch")
	portValue = Array(Float64, size(f, 1))
	portValue[1] = 1000000.0
	for n = 2:size(r, 1)
		inds = (vec(f[n, :]) .> 0.0)
		if !any(inds)
			portValue[n] = portValue[n-1]
		else
			portValue[n] = port_value_update(portValue[n-1], vec(r[n, :])[inds])
		end
	end
	return(portValue)
end
port_value_update(p::Float64, r::Float64) = p * (r + 1)
port_value_update(p::Vector{Float64}, r::Vector{Float64}) = (length(p) == length(r)) ? sum(Float64[ port_value_update(p[j], r[j]) for j = 1:length(p) ]) : error("Length mismatch")
port_value_update(p::Float64, r::Vector{Float64}) = port_value_update((p/length(r))*ones(Float64, length(r)), r)

#-------------------------------------------
# ASYMMETRIC DISTRIBUTIONS IN FINANCIAL RETURNS
#-------------------------------------------
function asymmetric_distribution_effect( ; blockLength::Float64=5.0, numResample::Int=1000, varProxyMethod::Symbol=:historicvariance, numObs::Int=2000, rSquared::Float64=0.05, a::Float64=0.0, lag::Int=17, numAsset::Int=20, tradingSimRSquared::Float64=0.01, includeConstantInSim::Bool=true)
	#Get robust measure of return data and compare it to Normal and t-Dist with 2 DoF
	secRet = read_local(secList, :return)
    robustSkew = Float64[ robust_skew(secRet[j]) for j = 1:length(secList) ]
	logNormalRobustSkew = mean(Float64[ robust_skew(rand(LogNormal(0, 1), 100)) for k = 1:1000 ])
	println("Drawing plot 1")
	estPlot1 = plot(x=secList, y=robustSkew, yintercept=[0.0, logNormalRobustSkew], Geom.point, Geom.hline, Guide.xlabel("Ticker code"), Guide.ylabel("Robust skewness"), defaultThemeOverride)
    draw_local(estPlot1, "Robust_Skewness_of_Daily_Financial_Returns", dirPath=outputDir, fileType=:svg)
	#Compare mean and trimmed mean on resampled financial returns (use WOW (Largest right skew) and ANZ (Largest left skew))
	println("Drawing mean comparison plots (full period)")
	pWOWFull = kernel_plot_robust_mean_versus_trimmed_mean(secRet, "WOW", blockLength=blockLength)
	draw_local(pWOWFull, "WOW_Bootstrapped_Location_Estimator_Density_(full_period)", dirPath=outputDir, fileType=:svg)
	pANZFull = kernel_plot_robust_mean_versus_trimmed_mean(secRet, "ANZ", blockLength=blockLength)
	draw_local(pANZFull, "ANZ_Bootstrapped_Location_Estimator_Density_(full_period)", dirPath=outputDir, fileType=:svg)
	#Compare mean and trimmed mean on log-normal data
	kPlotLN = kernel_plot_robust_mean_versus_trimmed_mean(rand(LogNormal(0, 1), 100), bpMult=1, blockLength=1.0)
	draw_local(kPlotLN, "LogNormal_Bootstrapped_Location_Estimator_Density", dirPath=outputDir, fileType=:svg)
	println("Routine complete")
end

#--------------------------------------------
# ROBUST FUNCTIONS USED THROUGHOUT
#--------------------------------------------
#Number of observations to cut when trimming
tmean_num_cut{T}(x::AbstractVector{T}, p::Float64) = Int(floor(p*length(x)))
tmean_num_cut{T}(x::AbstractVector{T}, lp::Float64, up::Float64) = (tmean_num_cut(x, lp), tmean_num_cut(x, up))
#Trimmed mean
function tmean!{T}(x::AbstractVector{T}, lp::Float64=0.1, up::Float64=0.1 ; sorted::Bool=false)
    length(x) < 2 && error("Input must contain at least 2 observations")
    !(0.0 <= lp <= up) && error("Invalid lower portion")
	!(lp <= up <= 1.0) && error("Invalid upper portion")
    (numLower, numUpper) = tmean_num_cut(x, lp, up)
	!sorted && sort!(x)
	return(mean(sub(x, numLower+1:length(x)-numUpper)))
end
tmean!{T}(x::AbstractVector{T}, p::Float64=0.2 ; sorted::Bool=false) = tmean!(x, p/2, p/2, sorted=sorted)
tmean{T}(x::AbstractVector{T}, lp::Float64=0.1, up::Float64=0.1 ; sorted::Bool=false) = sorted ? tmean!(x, lp, up, sorted=true) : tmean!(deepcopy(x), lp, up, sorted=false)
tmean{T}(x::AbstractVector{T}, p::Float64=0.2 ; sorted::Bool=false) = tmean(x, p/2, p/2, sorted=sorted)
#Median when input is already sorted
function median_sorted{T}(x::AbstractVector{T})
    length(x) < 2 && error("Input must contain at least 2 observations")
    if iseven(length(x))
        i = Int(length(x)/2)
        return(mean(x[i:i+1]))
    else
        return(x[Int((length(x)+1)/2)])
    end
end
#Hogg's robust measure of tail-fatness.
function hogg_robust_kurt!{T}(x::Vector{T} ; sorted::Bool=false, numerTail::Float64=0.05, denomTail::Float64=0.5)
    length(x) < 10 && error("Input must contain at least 10 observations")
    !(0.0 < numerTail <= 0.2) && error("Non-sensible value for numerTail")
    !(0.3 < denomTail <= 0.5) && error("Non-sensible value for denomTail")
    numObsNumer = Int(ceil(numerTail * length(x)))
    numObsDenom = Int(ceil(denomTail * length(x)))
    !sorted && sort!(x)
    numerLeft = mean(sub(x, 1:numObsNumer))
    numerRight = mean(sub(x, length(x)-numObsNumer+1:length(x)))
    denomLeft = mean(sub(x, 1:numObsDenom))
    denomRight = mean(sub(x, length(x)-numObsDenom+1:length(x)))
    return((numerRight - numerLeft) / (denomRight - denomLeft))
end
hogg_robust_kurt{T}(x::Vector{T} ; sorted::Bool=false, numerTail::Float64=0.05, denomTail::Float64=0.5) = sorted ? hogg_robust_kurt!(x, sorted=true, numerTail=numerTail, denomTail=denomTail) : hogg_robust_kurt!(deepcopy(x), sorted=false, numerTail=numerTail, denomTail=denomTail)
#Robust skewness
function robust_skew{T<:Number}(x::Vector{T}; qProb::Float64=0.75)
	length(x) < 3 && error("Input must contain at least 10 observations")
	!(0.5 < qProb < 1.0) && error("Invalid quantile for robust skew procedure")
	qEst = quantile(x, Float64[1 - qProb, 0.5, qProb])
	return((qEst[3] + qEst[1] - 2*qEst[2]) / (qEst[3] - qEst[1]))
end

#--------------------------------------------
# ROLLING WINDOW HISTORICAL VARIANCE
#--------------------------------------------
historical_variance{T<:Number}(r::AbstractVector{T}, lagWindow::Int=100) = (0 < lagWindow < length(r)) ? Float64[ var(sub(r, n-lagWindow+1:n)) for n = lagWindow:length(r) ] : error("Invalid lagWindow")

#--------------------------------------------
# DATA READING FUNCTIONS USED THROUGHOUT
#--------------------------------------------
function read_local(secList::Vector{ASCIIString}, dataType::Symbol ; checkLength::Int=0)
    if dataType == :return
        x = Vector{Float64}[ vec(readcsv(dataDir*"ASX_"*secList[j]*"_Return-Trade_uMarketCloseAuction_Bow.csv", Float64)) for j = 1:length(secList) ]
    elseif dataType == :realisedvariance
        x = Vector{Float64}[ vec(readcsv(dataDir*"ASX_"*secList[j]*"_RealisedVariance-5Minute-BidAskMidpoint_uFirstToLast_Bow.csv", Float64)) for j = 1:length(secList) ]
    else
        error("Invalid data type")
    end
    if checkLength != 0
        any(Int[ length(x[j]) for j = 1:length(x) ] .!= checkLength) && error("Invalid data length")
    end
    return(x)
end
function read_local(dataType::Symbol)
    if dataType == :calendar
        x = vec(readcsv(dataDir*"Calendar.csv", DateTime))
    else
        error("Invalid data type")
    end
    return(x)
end

#--------------------------------------------
# PLOTTING FUNCTIONS USED THROUGHOUT
#--------------------------------------------
#Function for drawing plot to saved image file
function draw_local(p::Plot, fileName::ASCIIString ; dirPath::ASCIIString="", fileType::Symbol=:png, width::Measure=40cm, height::Measure=20cm)
	dirPath[end] != '/' && (dirPath = dirPath * "/")
    !isdir(dirPath) && error("Input directory does not exist")
    filePath = dirPath * fileName * "." * string(fileType)
    if fileType == :svg     ; draw(SVG(filePath, width, height), p)
	elseif fileType == :png ; draw(PNG(filePath, width, height), p)
	elseif fileType == :pdf ; draw(PDF(filePath, width, height), p)
	else                    ; error("Invalid fileType symbol")
	end
	return(filePath)
end
function adjust_default_theme_color(x::Theme, colourString::ASCIIString)
    xC = deepcopy(x)
    xC.default_color = parse(Compose.Colorant, colourString)
    return(xC)
end
function default_legend(legendLabel::Vector{ASCIIString})
	legendTitle = "Legend:"
    length(legendLabel) > length(colourVec) && error("Default legend function cannot handle more than " * string(length(colourVec)) * " series")
    return(legendTitle, legendLabel, deepcopy(colourVec[1:length(legendLabel)]))
end

#--------------------------------------------
# OTHER COMMON FUNCTIONS
#--------------------------------------------
arr_to_mat{T}(x::Vector{Vector{T}}) = T[ x[j][k] for j = 1:length(x), k = 1:length(x[1])]
mat_to_arr{T}(x::Matrix{T}) = Vector{T}[ x[:, j] for j = 1:size(x, 2) ]
julia_version_dir() = "v"*string(VERSION.major)*"."*string(VERSION.minor)

end # module


text:
using MsgPack2
using Documenter

makedocs(modules=[MsgPack2],
         sitename="MsgPack2",
         authors="Jarrett Revels and other contributors",
         pages=["API Documentation" => "index.md"])

deploydocs(repo="github.com/beacon-biosignals/MsgPack2.jl.git")

ext: jl


text:
@testset "variables.jl" begin
    @testset "MOI.add_variable()" begin
        opt = SeaPearl.Optimizer()
        push!(opt.moimodel.variables, SeaPearl.MOIVariable("", 1, 2, MOI.VariableIndex(1)))

        @test MOI.add_variable(opt) == MOI.VariableIndex(2)
        @test length(opt.moimodel.variables) == 2
        @test isnothing(opt.moimodel.variables[2].min)
        @test isnothing(opt.moimodel.variables[2].max)
        @test opt.moimodel.variables[2].vi == MOI.VariableIndex(2)
        @test opt.moimodel.variables[2].cp_identifier == ""
    end
end


text:
dims(H::Hermitian{<:Number,<:Tensor}) = dims(parent(H))

dim(H::Hermitian{<:Number,<:Tensor},
    i::Int) = dim(parent(H),i)

matrix(H::Hermitian{<:Number,<:Tensor}) = Hermitian(matrix(parent(H)))

inds(H::Hermitian{<:Number,<:Tensor}) = inds(parent(H))

ind(H::Hermitian{<:Number,<:Tensor},
    i::Int) = ind(parent(H),i)

nnzblocks(H::Hermitian{<:Number,<:Tensor}) = nnzblocks(parent(H))

nzblocks(H::Hermitian{<:Number,<:Tensor}) = nzblocks(parent(H))

nblocks(H::Hermitian{<:Number,<:Tensor}) = nblocks(parent(H))

nzblock(H::Hermitian{<:Number,<:Tensor},
        i::Int) = nzblock(parent(H), i)

blockview(H::Hermitian{<:Number,<:Tensor},
          i::Int) = Hermitian(blockview(parent(H), i))


text:
"""NumPy array protocol type string (typestr) format

A string providing the basic type of the homogenous array. The basic string format
consists of 3 parts: a character describing the byteorder of the data
(<: little-endian, >: big-endian, |: not-relevant), a character code giving the basic
type of the array, and an integer providing the number of bytes the type uses.

https://zarr.readthedocs.io/en/stable/spec/v2.html#data-type-encoding
"""

include("MaxLengthStrings.jl")
using .MaxLengthStrings: MaxLengthString

primitive type ASCIIChar <: AbstractChar 8 end
ASCIIChar(x::UInt8) = reinterpret(ASCIIChar, x)
ASCIIChar(x::Integer) = ASCIIChar(UInt8(x))
UInt8(x::ASCIIChar) = reinterpret(UInt8, x)
Base.codepoint(x::ASCIIChar) = UInt8(x)
Base.show(io::IO, x::ASCIIChar) = print(io, Char(x))
Base.zero(::Union{ASCIIChar,Type{ASCIIChar}}) = ASCIIChar(Base.zero(UInt8))

typestr(t::Type) = string('<', 'V', sizeof(t))
typestr(t::Type{>:Missing}) = typestr(Base.nonmissingtype(t))
typestr(t::Type{Bool}) = string('<', 'b', sizeof(t))
typestr(t::Type{<:Signed}) = string('<', 'i', sizeof(t))
typestr(t::Type{<:Unsigned}) = string('<', 'u', sizeof(t))
typestr(t::Type{Complex{T}} where T<:AbstractFloat) = string('<', 'c', sizeof(t))
typestr(t::Type{<:AbstractFloat}) = string('<', 'f', sizeof(t))
typestr(t::Type{MaxLengthString{N,UInt32}}) where {N,T} = string('<', 'U', N)
typestr(t::Type{MaxLengthString{N,UInt8}}) where {N,T} = string('<', 'S', N)

const typestr_regex = r"^([<|>])([tbiufcmMOSUV])(\d+)$"
const typemap = Dict{Tuple{Char, Int}, DataType}(
    ('b', 1) => Bool,
    ('S', 1) => ASCIIChar,
    ('U', 1) => Char,
)
sizemapf(x::Type{<:Number}) = sizeof(x)
sizemapf(x::Type{<:Complex{T}}) where T = sizeof(T)
typecharf(::Type{<:Signed}) = 'i'
typecharf(::Type{<:Unsigned}) = 'u'
typecharf(::Type{<:AbstractFloat}) = 'f'
typecharf(::Type{<:Complex}) = 'c'
foreach([Float16,Float32,Float64,Int8,Int16,Int32,Int64,Int128,
  UInt8,UInt16,UInt32,UInt64,UInt128,
  Complex{Float16},Complex{Float32},Complex{Float64}]) do t
    typemap[(typecharf(t),sizemapf(t))] = t
end

function typestr(s::AbstractString)
    m = match(typestr_regex, s)
    if m === nothing
        throw(ArgumentError("$s is not a valid numpy typestr"))
    else
        byteorder, typecode, typesize = m.captures
        if byteorder == ">"
            throw(ArgumentError("Big-endian data not yet supported"))
        end
        tc, ts = first(typecode), parse(Int, typesize)
        if (tc in ('U','S')) && ts > 1
          return MaxLengthString{ts,tc=='U' ? UInt32 : UInt8}
        end
        # convert typecode to Char and typesize to Int
        typemap[(tc,ts)]
    end
end

"""Metadata configuration of the stored array

Each array requires essential configuration metadata to be stored, enabling correct
interpretation of the stored data. This metadata is encoded using JSON and stored as the
value of the “.zarray” key within an array store.

https://zarr.readthedocs.io/en/stable/spec/v2.html#metadata
"""
struct Metadata{T, N, C}
    zarr_format::Int
    shape::Ref{NTuple{N, Int}}
    chunks::NTuple{N, Int}
    dtype::String  # structured data types not yet supported
    compressor::C
    fill_value::Union{T, Nothing}
    order::Char
    filters::Nothing  # not yet supported
end

#To make unit tests pass with ref shape
import Base.==
function ==(m1::Metadata{T,N,C}, m2::Metadata{T,N,C}) where {T,N,C}
  m1.zarr_format == m2.zarr_format &&
  m1.shape[] == m2.shape[] &&
  m1.chunks == m2.chunks &&
  m1.dtype == m2.dtype &&
  m1.compressor == m2.compressor &&
  m1.fill_value == m2.fill_value &&
  m1.order == m2.order &&
  m1.filters == m2.filters
end

"Construct Metadata based on your data"
function Metadata(A::AbstractArray{T, N}, chunks::NTuple{N, Int};
        zarr_format::Integer=2,
        compressor::C=BloscCompressor(),
        fill_value::Union{T, Nothing}=nothing,
        order::Char='C',
        filters::Nothing=nothing
    ) where {T, N, C}
    T2 = fill_value === nothing ? T : Union{T,Missing}
    Metadata{T2, N, C}(
        zarr_format,
        size(A),
        chunks,
        typestr(eltype(A)),
        compressor,
        fill_value,
        order,
        filters
    )
end

Metadata(s::Union{AbstractString, IO}) = Metadata(JSON.parse(s))

"Construct Metadata from Dict"
function Metadata(d::AbstractDict)
    # create a Metadata struct from it

    compdict = d["compressor"]
    compressor = getCompressor(compdict)

    T = typestr(d["dtype"])
    N = length(d["shape"])
    C = typeof(compressor)

    fv = fill_value_decoding(d["fill_value"], T)

    TU = fv === nothing ? T : Union{T,Missing}

    Metadata{TU, N, C}(
        d["zarr_format"],
        NTuple{N, Int}(d["shape"]) |> reverse,
        NTuple{N, Int}(d["chunks"]) |> reverse,
        d["dtype"],
        compressor,
        fv,
        first(d["order"]),
        d["filters"]
    )
end

"Describes how to lower Metadata to JSON, used in json(::Metadata)"
function JSON.lower(md::Metadata)
    Dict{String, Any}(
        "zarr_format" => md.zarr_format,
        "shape" => md.shape[] |> reverse,
        "chunks" => md.chunks |> reverse,
        "dtype" => md.dtype,
        "compressor" => JSON.lower(md.compressor),
        "fill_value" => fill_value_encoding(md.fill_value),
        "order" => md.order,
        "filters" => md.filters
    )
end

# Fill value encoding and decoding as described in
# https://zarr.readthedocs.io/en/stable/spec/v2.html#fill-value-encoding

fill_value_encoding(v) = v
fill_value_encoding(::Nothing)=nothing
function fill_value_encoding(v::AbstractFloat)
    if isnan(v)
        "NaN"
    elseif isinf(v)
        v>0 ? "Infinity" : "-Infinity"
    else
        v
    end
end

Base.eltype(::Metadata{T}) where T = T

# this correctly parses "NaN" and "Infinity"
fill_value_decoding(v::AbstractString, T::Type{<:Number}) = parse(T, v)
fill_value_decoding(v::Nothing, T) = v
fill_value_decoding(v, T) = T(v)
fill_value_decoding(v, T::Type{ASCIIChar}) = v == "" ? nothing : v


text:
@testset "Durou" begin
    p, q = synthetic_gradient(SynthSphere(), radius = 50)
    Z = convert_gradient(Durou(), p, q)
    @test Z[38,38] ≈ 1.5742826661339246
    @test Z[1,1] ≈ 0.44970037692631165
    @test Z[76,76] ≈ 47.228212708897075
end


text:
#3D primitives
typealias Primitives3D  Union{AbstractGeometry{3}, AbstractMesh}
#2D primitives AKA sprites, since they are shapes mapped onto a 2D rectangle
typealias Sprites Union{AbstractGeometry{2}, Shape, Char, Type}
typealias AllPrimitives Union{AbstractGeometry, Shape, Char}

function assemble(data, shaders, primitive, instanced=nothing)
    if instanced == nothing
    end

end

function convert2gl!(buffer_preference, dict)
    view = Dict{String, String}()
    map!(dict) do key_value
        k,v = key_value
        new_val = if isa(v, AbstractArray)
            gl_convert(buffer_preference, v)
        elseif isa(v, AbstractMesh)
            NativeMesh(v)
        else
            gl_convert(v)
        end
        return k => new_val
    end

end


text:
macro F_YR()#(V, U, Ax, r, SA) # rear tire force
	code =quote
		if length(Ax) != 1
			Fyr = zeros(length(Ax),1)
			for ii in eachindex(Ax)
				Fyr[ii] = 259.95925391136984 + ((PD2*(FzR0 + (Ax[ii] - V[ii]*R[ii])*KZX)^2 + PD1*(FzR0 + (Ax[ii] - V[ii]*R[ii])*KZX))*sin(PC1*atan((((PK1*sin(2*atan(PK2*(FzR0 + (Ax[ii] - V[ii]*R[ii])*KZX))))/(((PD2*PC1*(FzR0 + (Ax[ii] - V[ii]*R[ii])*KZX)^2 + PD1*PC1*(FzR0 + (Ax[ii] - V[ii]*R[ii])*KZX)) + ((PD2*PC1*(FzR0 + (Ax[ii] - V[ii]*R[ii])*KZX)^2 + PD1*PC1*(FzR0 + (Ax[ii] - V[ii]*R[ii])*KZX)))/(((PD2*PC1*(FzR0 + (Ax[ii] - V[ii]*R[ii])*KZX)^2 + PD1*PC1*(FzR0 + (Ax[ii] - V[ii]*R[ii])*KZX))^2)^(0.5)))+EP))*((atan((V[ii] - lb*R[ii])/(U[ii]+EP))) + PH2*(FzR0 + (Ax[ii] - V[ii]*R[ii])*KZX) + PH1)) - ((PE2*(FzR0 + (Ax[ii] - V[ii]*R[ii])*KZX) + PE1)*(1 - PE3*(((atan((V[ii] - lb*R[ii])/(U[ii]+EP))) + PH2*(FzR0 + (Ax[ii] - V[ii]*R[ii])*KZX) + PH1))/((((atan((V[ii] - lb*R[ii])/(U[ii]+EP))) + PH2*(FzR0 + (Ax[ii] - V[ii]*R[ii])*KZX) + PH1)^2 )^(0.5))))*((((PK1*sin(2*atan(PK2*(FzR0 + (Ax[ii] - V[ii]*R[ii])*KZX))))/(((PD2*PC1*(FzR0 + (Ax[ii] - V[ii]*R[ii])*KZX)^2 + PD1*PC1*(FzR0 + (Ax[ii] - V[ii]*R[ii])*KZX)) + ((PD2*PC1*(FzR0 + (Ax[ii] - V[ii]*R[ii])*KZX)^2 + PD1*PC1*(FzR0 + (Ax[ii] - V[ii]*R[ii])*KZX)))/(((PD2*PC1*(FzR0 + (Ax[ii] - V[ii]*R[ii])*KZX)^2 + PD1*PC1*(FzR0 + (Ax[ii] - V[ii]*R[ii])*KZX))^2)^(0.5)))+EP))*((atan((V[ii] - lb*R[ii])/(U[ii]+EP))) + PH2*(FzR0 + (Ax[ii] - V[ii]*R[ii])*KZX) + PH1)) - atan((((PK1*sin(2*atan(PK2*(FzR0 + (Ax[ii] - V[ii]*R[ii])*KZX))))/(((PD2*PC1*(FzR0 + (Ax[ii] - V[ii]*R[ii])*KZX)^2 + PD1*PC1*(FzR0 + (Ax[ii] - V[ii]*R[ii])*KZX)) + ((PD2*PC1*(FzR0 + (Ax[ii] - V[ii]*R[ii])*KZX)^2 + PD1*PC1*(FzR0 + (Ax[ii] - V[ii]*R[ii])*KZX)))/(((PD2*PC1*(FzR0 + (Ax[ii] - V[ii]*R[ii])*KZX)^2 + PD1*PC1*(FzR0 + (Ax[ii] - V[ii]*R[ii])*KZX))^2)^(0.5)))+EP))*((atan((V[ii] - lb*R[ii])/(U[ii]+EP))) + PH2*(FzR0 + (Ax[ii] - V[ii]*R[ii])*KZX) + PH1)))))) + (PV2*(FzR0 + (Ax[ii] - V[ii]*R[ii])*KZX)^2 + PV1*(FzR0 + (Ax[ii] - V[ii]*R[ii])*KZX)));
			end
		else
			Fyr = Float64
			Fyr = 259.95925391136984 + ((PD2*(FzR0 + (Ax - V*R)*KZX)^2 + PD1*(FzR0 + (Ax - V*R)*KZX))*sin(PC1*atan((((PK1*sin(2*atan(PK2*(FzR0 + (Ax - V*R)*KZX))))/(((PD2*PC1*(FzR0 + (Ax - V*R)*KZX)^2 + PD1*PC1*(FzR0 + (Ax - V*R)*KZX)) + ((PD2*PC1*(FzR0 + (Ax - V*R)*KZX)^2 + PD1*PC1*(FzR0 + (Ax - V*R)*KZX)))/(((PD2*PC1*(FzR0 + (Ax - V*R)*KZX)^2 + PD1*PC1*(FzR0 + (Ax - V*R)*KZX))^2)^(0.5)))+EP))*((atan((V - lb*R)/(U+EP))) + PH2*(FzR0 + (Ax - V*R)*KZX) + PH1)) - ((PE2*(FzR0 + (Ax - V*R)*KZX) + PE1)*(1 - PE3*(((atan((V - lb*R)/(U+EP))) + PH2*(FzR0 + (Ax - V*R)*KZX) + PH1))/((((atan((V - lb*R)/(U+EP))) + PH2*(FzR0 + (Ax - V*R)*KZX) + PH1)^2 )^(0.5))))*((((PK1*sin(2*atan(PK2*(FzR0 + (Ax - V*R)*KZX))))/(((PD2*PC1*(FzR0 + (Ax - V*R)*KZX)^2 + PD1*PC1*(FzR0 + (Ax - V*R)*KZX)) + ((PD2*PC1*(FzR0 + (Ax - V*R)*KZX)^2 + PD1*PC1*(FzR0 + (Ax - V*R)*KZX)))/(((PD2*PC1*(FzR0 + (Ax - V*R)*KZX)^2 + PD1*PC1*(FzR0 + (Ax - V*R)*KZX))^2)^(0.5)))+EP))*((atan((V - lb*R)/(U+EP))) + PH2*(FzR0 + (Ax - V*R)*KZX) + PH1)) - atan((((PK1*sin(2*atan(PK2*(FzR0 + (Ax - V*R)*KZX))))/(((PD2*PC1*(FzR0 + (Ax - V*R)*KZX)^2 + PD1*PC1*(FzR0 + (Ax - V*R)*KZX)) + ((PD2*PC1*(FzR0 + (Ax - V*R)*KZX)^2 + PD1*PC1*(FzR0 + (Ax - V*R)*KZX)))/(((PD2*PC1*(FzR0 + (Ax - V*R)*KZX)^2 + PD1*PC1*(FzR0 + (Ax - V*R)*KZX))^2)^(0.5)))+EP))*((atan((V - lb*R)/(U+EP))) + PH2*(FzR0 + (Ax - V*R)*KZX) + PH1)))))) + (PV2*(FzR0 + (Ax - V*R)*KZX)^2 + PV1*(FzR0 + (Ax - V*R)*KZX)));
		end
		Fyr
	end
	return esc(code)
end


text:
#= simple.jl
Philip Barrett, pobarrett@gmail.com
26may2016, Chicago

Defines basic Polygon operations:
    * Vector Addition
    * SetsetSum (of both Polygon pairs and arrays)
    * Scalar multiplication
    * Weighted setSum
=#

"""
    add( poly::Polygon, pt::Vector )
Adds a vector to all points in a Polygon
"""
function (+)( poly::Polygon, u::Vector{Float64} )
  return Polygon( pts = poly.pts + ones(size(poly.pts)[1]) * u' )
end

function (+)( u::Vector{Float64}, poly::Polygon )
  return poly + u
end

"""
    setSum( poly1::Polygon, poly2::Polygon, dirs, outer=true )
Addition functionality.  Provides either an inner or out approximation
"""
# TODO: Improve the inner approximation when the serach direction is exactly
# a face normal
function setSum( poly1::Polygon, poly2::Polygon,
                    dirs::Matrix{Float64}, outer::Bool=true )
  dists1 = poly1.pts * dirs'
  dists2 = poly2.pts * dirs'
      # The distances in each direction
  maxdists1 = maximum( dists1, 1 )
  maxdists2 = maximum( dists2, 1 )
      # The maximum distances
  maxdists = vec( maxdists1 + maxdists2 )

  if( outer )
    return Polygon( dirs=dirs, dists=maxdists )
  end

  Z1 = zeros( dirs )
  Z2 = zeros( dirs )
        # Initiate output points

  for( i in 1:size(dirs)[1])
      Z1[i,:] = poly1.pts[ dists1[:,i] .== maxdists1[i], : ][1,:]
      Z2[i,:] = poly2.pts[ dists2[:,i] .== maxdists2[i], : ][1,:]
  end
  return Polygon( pts = Z1 + Z2 )
end

"""
    setSum( poly1::Array{Polygon,1}, poly2::Array{Polygon,1}, dirs, outer=true )
Adds an array of Polygons
"""
function setSum( polys::Array{Polygon,1}, dirs::Matrix{Float64},
                    outer::Bool=true )

  N = length(polys)
      # Number of Polygons to add
  dists = zeros( size(dirs)[1] )
      # Initate the distances
  innerpts = zeros(dirs)
      # Initialize the points for the inner approximation
  for( i in 1:N )
    thisdists = polys[i].pts * dirs'
    maxdists = maximum( thisdists, 1 )'
    dists += maximum( polys[i].pts * dirs', 1 )'
        # Sum the distances
    if( !outer )
      for( j in 1:size(dirs)[1])
        innerpts[j,:] += polys[i].pts[ thisdists[:,j] .== maxdists[j], : ][1,:]
      end
    end
  end

  if( outer )
    return Polygon( dirs=dirs, dists=vec( dists ) )
  end
  return Polygon( pts=innerpts )
end

"""
    times( poly1::Array{Polygon,1}, poly2::Array{Polygon,1}, dirs, outer=true )
Scalar multiplication
"""
function (*)( k::Float64, poly::Polygon )
  return Polygon( pts = k * poly.pts )
end

function (.*)( v::Vector{Float64}, poly::Polygon )
# println( "ones(size(poly.pts)[1]):\n", ones(size(poly.pts)[1]) )
# println( "v':\n", v' )
  return Polygon( pts = ( poly.pts .* ( ones(size(poly.pts)[1]) * v' ) ) )
end

"""
    weightedSum( polys::Array{Polygon,1}, wts::Vector dirs, outer=true )
Computes a weighted setSum of Polygons
"""
function weightedSum( polys::Array{Polygon,1}, wts::Vector{Float64},
                        dirs::Matrix{Float64}, outer::Bool=true )
  N = length( polys )
      # Number of Polygons
  polys2 = [ (wts[i] * polys[i])::Polygon for i in 1:N ]
      # Scale polygons according to weights
  return( setSum( polys2, dirs, outer ) )
end


text:
# Sample from arbitrary arrays

################################################################
#
#  A variety of algorithms for sampling without replacement
#
#  They are suited for different cases.
#
#  Particularly,
#  - Fisher-Yates sampler is suited for general cases
#    where n is not overly large
#
#  - Self avoiding sampler is suited for cases where k << n
#
################################################################

function pick2!(a::AbstractArray, x::AbstractArray)
    # Pick a pair of values without replacement

    n0 = length(a)
    i1 = randi(n0)
    i2 = randi(n0 - 1)
    if i2 == i1
        i2 = n0
    end

    x[1] = a[i1]
    x[2] = a[i2]
end

## A sampler that implements without-replacement sampling
## via Fisher-Yates shuffling
##
immutable FisherYatesSampler
    n::Int
    seq::Vector{Int}   # Internal sequence for shuffling

    FisherYatesSampler(n::Int) = new(n, [1:n])
end

function rand!(s::FisherYatesSampler, a::AbstractArray, x::AbstractArray)
    # draw samples without-replacement to x

    n::Int = s.n
    k::Int = length(x)
    if k > n
        throw(ArgumentError("Cannot draw more than n samples without replacement."))
    end

    seq::Vector{Int} = s.seq
    for i = 1:k
        j = randi(i, n)
        sj = seq[j]
        x[i] = a[sj]
        seq[j] = seq[i]
        seq[i] = sj
    end
    x
end

fisher_yates_sample!(a::AbstractArray, x::AbstractArray) = rand!(FisherYatesSampler(length(a)), a, x)

function self_avoid_sample!{T}(a::AbstractArray{T}, x::AbstractArray)
    # This algorithm is suitable when length(x) << length(a)

    s = Set{T}()
    # sizehint(s, length(x))
    rgen = RandIntSampler(length(a))

    # first one
    idx = rand(rgen)
    x[1] = a[idx]
    add!(s, idx)

    # remaining
    for i = 2:length(x)
        idx = rand(rgen)
        while contains(s, idx)
            idx = rand(rgen)
        end
        x[i] = a[idx]
        add!(s, idx)
    end
    x
end

###########################################################
#
#   Interface functions
#
###########################################################

sample(a::AbstractArray) = a[randi(length(a))]

function sample!(a::AbstractArray, x::AbstractArray; replace=true)
    n = length(a)
    k = length(x)

    if !isempty(x)
        if replace   # with replacement
            s = RandIntSampler(n)
            for i = 1:k
                x[i] = a[rand(s)]
            end

        else  # without replacement
            if k > n
                throw(ArgumentError("n exceeds the length of x"))
            end

            if k == 1
                x[1] = sample(a)

            elseif k == 2
                pick2!(a, x)

            elseif n < k * max(k, 100)
                fisher_yates_sample!(a, x)

            else
                self_avoid_sample!(a, x)
            end
        end
    end
    x
end

function sample{T}(a::AbstractArray{T}, n::Integer; replace=true)
    sample!(a, Array(T, n); replace=replace)
end

function sample{T}(a::AbstractArray{T}, dims::Dims; replace=true)
    sample!(a, Array(T, dims); replace=replace)
end

################################################################
#
#  Weighted sampling
#
################################################################

function wsample(a::AbstractArray, w::AbstractArray{Float64}, wsum::Float64)
    n = length(w)
    t = rand() * wsum

    i = 1
    s = w[1]

    while i < n && s < t
        i += 1
        s += w[i]
    end
    a[i]
end

wsample(a::AbstractArray, w::AbstractArray{Float64}) = wsample(a, w, sum(w))

function wsample!(a::AbstractArray, w::AbstractArray{Float64}, x::AbstractArray;
    wsum::Float64=NaN)

    n = length(a)
    if length(w) != n
        throw(ArgumentError("Inconsistent argument dimensions."))
    end

    _wsum::Float64 = isnan(wsum) ? sum(w) : wsum
    for i = 1:length(x)
        x[i] = wsample(a, w, _wsum)
    end
    x
end

function wsample{T}(a::AbstractArray{T}, w::AbstractArray{Float64}, n::Integer; wsum::Float64=NaN)
    wsample!(a, w, Array(T, n); wsum=wsum)
end

function wsample{T}(a::AbstractArray{T}, w::AbstractArray{Float64}, dims::Dims; wsum::Float64=NaN)
    wsample!(a, w, Array(T, dims); wsum=wsum)
end


text:
"""
    restrict(img[, dims]) -> imgr

Reduce the size of `img` by approximately two-fold along the dimensions listed in
`dims`, or all spatial coordinates if `dims` is not specified.

# Output

The type of output array `imgr` depends on the input type:

- If `img` is an `OffsetArray`, then output array `imgr` will also be an `OffsetArray`.
- If `img` is not an `OffsetArray`, then output array `imgr` will be an `Array` type
  even if it has offset indices.

The size of `imgr` is approximately `1/2` of the original size. More specifically:

- if `Nₖ = size(img, k)` is odd, then `size(imgr, k) == (Nₖ+1) ÷ 2`.
- if `Nₖ = size(img, k)` is even, then `size(imgr, k) == (Nₖ÷2) + 1`.

# Examples

The optional argument `dims` can be a `Tuple` or `Integer`:

```julia
A = rand(5, 5) # size: (5, 5)

restrict(A) # size: (3, 3)

restrict(A, 1) # size: (3, 5)
restrict(A, 2) # size: (5, 3)

restrict(A, (1, )) # size: (3, 5)
restrict(A, (1, 2)) # size: (3, 3)
```

Unless the input array is 1-based, the origin will be halfed:

```julia
julia> using ImageBase, OffsetArrays

julia> Ao = OffsetArray(rand(5, 4), 5, 6);

julia> Ar = restrict(Ao);

julia> axes(Ao)
(OffsetArrays.IdOffsetRange(values=6:10, indices=6:10), OffsetArrays.IdOffsetRange(values=7:10, indices=7:10))

julia> axes(Ar)
(OffsetArrays.IdOffsetRange(values=3:5, indices=3:5), OffsetArrays.IdOffsetRange(values=4:6, indices=4:6))
```

# Extended help

The term `restrict` is taken from the coarsening operation of algebraic multigrid
methods; it is the adjoint of "prolongation" (which is essentially interpolation).
`restrict` anti-aliases the image as it goes, so is better than a naive summation
over 2x2 blocks.
The implementation of `restrict` has been tuned for performance, and should be
a fast method for constructing [pyramids](https://en.wikipedia.org/wiki/Pyramid_(image_processing)).

If `l` is the size of `img` along a particular dimension, `restrict` produces an
array of size `(l+1)÷2` for odd `l`,
and `l÷2 + 1` for even `l`. See the example below for an explanation.

See also `ImageTransformations.imresize`.

# Example

```julia
a_course = [0, 1, 0.3]
```
If we were to interpolate this at the halfway points, we'd get
```julia
a_fine = [0, 0.5, 1, 0.65, 0.3]
```
Note that `a_fine` is obtained from `a_course` via the *prolongation* operator `P` as
`P*a_course`, where
```julia
P = [1   0   0;      # this line "copies over" the first point
     0.5 0.5 0;      # this line takes the mean of the first and second point
     0   1   0;      # copy the second point
     0   0.5 0.5;    # take the mean of the second and third
     0   0   1]      # copy the third
```
`restrict` is the adjoint of prolongation. Consequently,
```julia
julia> restrict(a_fine)
3-element Array{Float64,1}:
 0.125
 0.7875
 0.3125

julia> (P'*a_fine)/2
3-element Array{Float64,1}:
 0.125
 0.7875
 0.3125
```
where the division by 2 approximately preserves the mean intensity of the input.

As we see here, for odd-length `a_fine`, restriction is the adjoint of interpolation at half-grid points.
When `length(a_fine)` is even, restriction is the adjoint of interpolation at 1/4 and 3/4-grid points.
This turns out to be the origin of the `l->l÷2 + 1` behavior.

One consequence of this definition is that the edges move towards zero:
```julia
julia> restrict(ones(11))
6-element Array{Float64,1}:
 0.75
 1.0
 1.0
 1.0
 1.0
 0.75
```
In some applications (e.g., image registration), you may find it useful to trim the edges.
"""
restrict(img::AbstractArray, ::Tuple{}) = img

restrict(A::AbstractArray) = restrict(A, coords_spatial(A))
function restrict(A::AbstractArray, dims::Dims)
    restrict(restrict(A, dims[1]), Base.tail(dims))
end

function restrict(A::AbstractArray{T,N}, dim::Integer) where {T,N}
    if Base.has_offset_axes(A)
        # For type stability, we cannot return `OffsetArray` in this method
        A = OffsetArrays.no_offset_view(A)
    end

    indsA = axes(A)
    newinds = ntuple(i->i==dim ? restrict_indices(indsA[dim]) : indsA[i], Val(N))
    out = Array{restrict_eltype(first(A)), N}(undef, last.(newinds))
    restrict!(out, A, dim)
    out
end
function restrict(A::OffsetArray{T,N}, dim::Integer) where {T,N}
    indsA = axes(A)
    newinds = map(UnitRange, ntuple(i->i==dim ? restrict_indices(indsA[dim]) : indsA[i], Val(N)))
    # By shifting it back to normal array, the internal for loop becomes faster because
    # it requires less indexing computation
    OffsetArray(restrict(A.parent, dim), newinds)
end

function restrict_eltype(A::AbstractArray)
    # infer the restrict_eltype on `eltype(A)` while preserving the container type
    # TODO: maybe there's more efficient way than the `similar` here..
    typeof(similar(A, _restrict_eltype(eltype(A)), ntuple(_->1, ndims(A))))
end
restrict_eltype(x) = _restrict_eltype(typeof(x))

for CT in (:AbstractGray, :AbstractRGB, :TransparentGray, :TransparentRGB)
    @eval _restrict_eltype(::Type{C}) where C<:$CT = __restrict_eltype(C)
end
_restrict_eltype(::Type{T}) where T = typeof(one(T)/4 + one(T)/2)
_restrict_eltype(::Type{C}) where C<:Color = __restrict_eltype(RGB{eltype(C)})
_restrict_eltype(::Type{C}) where C<:Colorant = __restrict_eltype(ARGB{eltype(C)})
function __restrict_eltype(::Type{C}) where C
    BT = base_colorant_type(C)
    isconcretetype(BT) && return floattype(BT)
    BT{_restrict_eltype(eltype(C))}
end

function restrict!(out::AbstractArray{T,N}, A::AbstractArray, dim) where {T,N}
    # a no-op for singleton dimension
    if size(A, dim) == 1 # this includes `dim > N` cases
        return copyto!(out, A)
    end
    indsout, indsA = axes(out), axes(A)
    ndims(out) == ndims(A) || throw(DimensionMismatch("input and output must have the same number of dimensions"))
    for d = 1:length(indsA)
        target = d==dim ? restrict_indices(indsA[d]) : indsA[d]
        indsout[d] == target || error("input and output must have corresponding indices; to be consistent with the input indices,\ndimension $d should be $target, got $(indsout[d])")
    end
    indspre, indspost = indsA[1:dim-1], indsA[dim+1:end]
    _restrict!(out, indsout[dim], A, indspre, indsA[dim], indspost)
end

@generated function _restrict!(out, indout, A,
                               indspre::NTuple{Npre,AbstractUnitRange},
                               indA,
                               indspost::NTuple{Npost,AbstractUnitRange}) where {Npre,Npost}
    Ipre = [Symbol(:ipre_, d) for d = 1:Npre]
    Ipost = [Symbol(:ipost_, d) for d = 1:Npost]
    quote
        $(Expr(:meta, :noinline))
        T = eltype(out)
        if isodd(length(indA))
            half = convert(eltype(T), 0.5)
            quarter = convert(eltype(T), 0.25)
            @nloops $Npost ipost d->indspost[d] begin
                iout = first(indout)
                @nloops $Npre ipre d->indspre[d] begin
                    out[$(Ipre...), iout, $(Ipost...)] = zero(T)
                end
                ispeak = true
                for iA in indA
                    if ispeak
                        @inbounds @nloops $Npre ipre d->indspre[d] begin
                            out[$(Ipre...), iout, $(Ipost...)] +=
                                half*convert(T, A[$(Ipre...), iA, $(Ipost...)])
                        end
                    else
                        @inbounds @nloops $Npre ipre d->indspre[d] begin
                            tmp = quarter*convert(T, A[$(Ipre...), iA, $(Ipost...)])
                            out[$(Ipre...), iout, $(Ipost...)]   += tmp
                            out[$(Ipre...), iout+1, $(Ipost...)] = tmp
                        end
                    end
                    ispeak = !ispeak
                    iout += ispeak
                end
            end
        else
            threeeighths = convert(eltype(T), 0.375)
            oneeighth = convert(eltype(T), 0.125)
            z = zero(T)
            fill!(out, zero(T))
            @nloops $Npost ipost d->indspost[d] begin
                peakfirst = true
                iout = first(indout)
                for iA in indA
                    if peakfirst
                        @inbounds @nloops $Npre ipre d->indspre[d] begin
                            tmp = convert(T, A[$(Ipre...), iA, $(Ipost...)])
                            out[$(Ipre...), iout, $(Ipost...)] += threeeighths*tmp
                            out[$(Ipre...), iout+1, $(Ipost...)] += oneeighth*tmp
                        end
                    else
                        @inbounds @nloops $Npre ipre d->indspre[d] begin
                            tmp = convert(T, A[$(Ipre...), iA, $(Ipost...)])
                            out[$(Ipre...), iout, $(Ipost...)]   += oneeighth*tmp
                            out[$(Ipre...), iout+1, $(Ipost...)] += threeeighths*tmp
                        end
                    end
                    peakfirst = !peakfirst
                    iout += peakfirst
                end
            end
        end
        out
    end
end

# If we're restricting along dimension 1, there are some additional efficiencies possible
@generated function _restrict!(out, indout, A, ::NTuple{0,AbstractUnitRange},
                               indA, indspost::NTuple{Npost,AbstractUnitRange}) where Npost
    Ipost = [Symbol(:ipost_, d) for d = 1:Npost]
    quote
        $(Expr(:meta, :noinline))
        T = eltype(out)
        if isodd(length(indA))
            half = convert(eltype(T), 0.5)
            quarter = convert(eltype(T), 0.25)
            @inbounds @nloops $Npost ipost d->indspost[d] begin
                iout, iA = first(indout), first(indA)
                nxt = convert(T, A[iA+1, $(Ipost...)])
                out[iout, $(Ipost...)] = half*convert(T, A[iA, $(Ipost...)]) + quarter*nxt
                for iA in first(indA)+2:2:last(indA)-2
                    prv = nxt
                    nxt = convert(T, A[iA+1, $(Ipost...)])
                    out[iout+=1, $(Ipost...)] = quarter*(prv+nxt) + half*convert(T, A[iA, $(Ipost...)])
                end
                out[iout+1, $(Ipost...)] = quarter*nxt + half*convert(T, A[last(indA), $(Ipost...)])
            end
        else
            threeeighths = convert(eltype(T), 0.375)
            oneeighth = convert(eltype(T), 0.125)
            z = zero(T)
            @inbounds @nloops $Npost ipost d->indspost[d] begin
                c = d = z
                iA = first(indA)
                for iout = first(indout):last(indout)-1
                    a, b = c, d
                    c, d = convert(T, A[iA, $(Ipost...)]), convert(T, A[iA+1, $(Ipost...)])
                    iA += 2
                    out[iout, $(Ipost...)] = oneeighth*(a+d) + threeeighths*(b+c)
                end
                out[last(indout), $(Ipost...)] = oneeighth*c + threeeighths*d
            end
        end
        out
    end
end

restrict_size(len::Integer) = isodd(len) ? (len+1)>>1 : (len>>1)+1
function restrict_indices(r::AbstractUnitRange)
    f, l = first(r), last(r)
    isodd(f) && return (f+1)>>1:restrict_size(l)
    f>>1 : (isodd(l) ? (l+1)>>1 : l>>1)
end
restrict_indices(r::Base.OneTo) = Base.OneTo(restrict_size(length(r)))
function restrict_indices(r::UnitRange)
    f, l = first(r), last(r)
    isodd(f) && return (f+1)>>1:restrict_size(l)
    f>>1 : (isodd(l) ? (l+1)>>1 : l>>1)
end


text:
@testset "fit moments" begin
    D = LogNormal(1,0.6)
    M = Moments(mean(D), var(D))
    Dfit = fit(LogNormal, M)
    @test D ≈ Dfit
    # handle not giving variance
    @test_throws Exception fit(LogNormal, Moments(3.2))
end;
@testset "fit two quantiles" begin
    qpl = @qp_m(3)
    qpu = @qp_u(5)
    d = fit(LogNormal, qpl, qpu);
    @test quantile.(d, [qpl.p, qpu.p]) ≈ [qpl.q, qpu.q]
    d = fit(LogNormal, qpu, qpl) # sort
    @test quantile.(d, [qpl.p, qpu.p]) ≈ [qpl.q, qpu.q]
end;
@testset "fit to quantilepoint and mean" begin
    d = LogNormal(1,1)
    m = log(mean(d))
    qp = @qp(quantile(d,0.95),0.95)
    dfit = fit_mean_quantile(LogNormal, mean(d), qp)
    @test dfit ≈ d
    dfit = fit(LogNormal, mean(d), qp, Val(:mean))
    @test dfit ≈ d
    # with lower quantile
    qp = @qp(quantile(d,0.05),0.05)
    dfit = fit_mean_quantile(LogNormal, mean(d), qp)
    @test dfit ≈ d
    # very close to mean can give very different results:
    qp = @qp(mean(d)-1e-4,0.95)
    dfit = fit_mean_quantile(LogNormal, mean(d), qp)
    @test mean(dfit) ≈ mean(d) && quantile(dfit, qp.p) ≈ qp.q
end;
@testset "fit to quantilepoint and mode" begin
    d = LogNormal(1,1)
    m = log(mode(d))
    qp = @qp(quantile(d,0.95),0.95)
    dfit = fit_mode_quantile(LogNormal, mode(d), qp)
    @test dfit ≈ d
    dfit = fit(LogNormal, mode(d), qp, Val(:mode))
    @test dfit ≈ d
    # with lower quantile
    qp = @qp(quantile(d,0.025),0.025)
    dfit = fit_mode_quantile(LogNormal, mode(d), qp)
    @test mode(dfit) ≈ mode(d) && quantile(dfit, qp.p) ≈ qp.q
end;
@testset "fit to quantilepoint and median" begin
    d = LogNormal(1,1)
    qp = @qp(quantile(d,0.95),0.95)
    dfit = fit(LogNormal, median(d), qp, Val(:median))
    @test dfit ≈ d
end;

@testset "Σstar" begin
    ss = Σstar(4.5)
    @test ss() == 4.5
end;
@testset "fit to mean and Σstar" begin
    d = LogNormal(1, log(1.2))
    dfit = fit(LogNormal, mean(d), Σstar(1.2))
    @test d == dfit
    @test σstar(d) ≈ 1.2
end;

@testset "fit_mean_relerror" begin
    d = fit_mean_relerror(LogNormal, 10.0, 0.03);
    @test all((mean(d), std(d)/mean(d)) .≈ (10.0, 0.03))
end;


text:
using SMACOF
using DelimitedFiles
using Distances
using Plots
using Shapefile
gr()

script_location = @__FILE__
script_dir = dirname(script_location)

dms2deg(d,m,s) = d + m / 60 + s / 3600

clld_path = joinpath(script_dir, "sgb128_lat_lon.txt")
cities_lat_lon_dms = readdlm(clld_path, skipstart = 7)
cities_lat_lon = zeros(size(cities_lat_lon_dms, 1), 2)
for i in 1:size(cities_lat_lon_dms, 1)
    lat = dms2deg(cities_lat_lon_dms[i, 1:3]...)
    lon = dms2deg(cities_lat_lon_dms[i, 5:7]...)
    cities_lat_lon[i, 1] = lat
    cities_lat_lon[i, 2] = -lon
end
cn_path = joinpath(script_dir, "sgb128_names.txt")
cities_names = readdlm(cn_path, ',', skipstart = 2)[:, 2]
cities_names = map(x -> x[2:end], cities_names)

Δ = dist(cities_lat_lon)
anchors = (idx = 1:128, pos = cities_lat_lon)

Xhist = smacof(Δ, Xinit = rand(size(Δ, 1), 2),  anchors = anchors, verbose = true, return_hist = true).hist
# Xhist = wda_smacof(Δ, η = 0.8, verbose = true, anchors = anchors, ε = 1e-6, return_history = true)[2]

# Make mp4 and then convert to gif
anim_dir = joinpath(dirname(script_dir), "tmp")
animation = Animation(anim_dir, String[])
# shp = Shapefile.shapes(Shapefile.Table("examples/cities/cb_2018_us_state_500k.shp"))
for i ∈ 1:length(Xhist)
    println(i)
    Xi = align(Xhist[i], anchors)
    p = scatter(Xi[:, 2], Xi[:, 1], alpha = 0.7, markersize = 3, legend = nothing);
    # plot!(shp)
    annotate!.(Xi[:,2], Xi[:,1], text.(cities_names, :black, :left, 7))
    xlims!((-130, -60))
    ylims!((20, 55))
    title!("$i")
    frame(animation, p)
end
gif_path = joinpath(script_dir, "cities")
run(`ffmpeg -y -r 15 -i $anim_dir/"%06d.png" -vcodec libx264 -crf 25 "$(gif_path).mp4"`)
run(`ffmpeg -y -i $(gif_path).mp4 -vf "fps=10,scale=1000:-1:flags=lanczos,split[s0][s1];[s0]palettegen[p];[s1][p]paletteuse" -loop 0 $(gif_path).gif`)
run(`rm $(gif_path).mp4`)

# Clear image folder
for f in readdir(anim_dir)
    run(`rm $(joinpath(anim_dir, f))`)
end


text:
module TestDeprecated

using Test
using InitialValues

@testset "Init" begin
    @test_deprecated Init(+)
end

end  # module


text:
using InteractModels, DataFrames, Interact, Test

# Test stoled from Interact.jl docs
@testset "interactive model" begin
    color(i) = colors[i%length(colors)+1]
    colors = [
        "black", "gray", "silver", "maroon", "red", "olive", "yellow",
        "green", "lime", "teal", "aqua", "navy", "blue", "purple", "fuchsia"
    ]
    width, height = 700, 300
    nsamples = 256
    pars = (;
        sample_step=Param(val=0.05, range=0.01:0.001:0.1, label="Sample step", description="The step size between samples"),
        phase=Param(val=0.0, range=0:0.1:2pi, label="Phase", description="Phase of the starting point"),
        radii=Param(val=20,range=0:0.1:60, label="Radus", description="Radius of the circles")
    )

    interface = InteractModel(pars; ncolumns=2, title="slinky") do m
        m = stripparams(m)
        println(m)
        cxs_unscaled = [i * m.sample_step + m.phase for i in 1:nsamples]
        cys = sin.(cxs_unscaled) .* height/3 .+ height/2
        cxs = cxs_unscaled .* width/4pi
        # dom"div"()
        dom"svg:svg[width=$width, height=$height]"(
             (dom"svg:circle[cx=$(cxs[i]), cy=$(cys[i]), r=$(m.radii), fill=$(color(i))]"()
              for i in 1:nsamples)...
        )
    end

    ui_ = ui(interface)
    # Title node matches.
    @test first(ui_.children).children == dom"h1"("slinky").children
    # TODO how to test this more?

    submodel_interface = InteractModel(pars; ncolumns=2, submodel=NamedTuple, title="slinky") do m
        m = stripparams(m)
        println(m)
        cxs_unscaled = [i * m.sample_step + m.phase for i in 1:nsamples]
        cys = sin.(cxs_unscaled) .* height/3 .+ height/2
        cxs = cxs_unscaled .* width/4pi
        # dom"div"()
        dom"svg:svg[width=$width, height=$height]"(
             (dom"svg:circle[cx=$(cxs[i]), cy=$(cys[i]), r=$(m.radii), fill=$(color(i))]"()
              for i in 1:nsamples)...
        )
    end

    ui_ = ui(submodel_interface)
    @test first(ui_.children).children == dom"h1"("slinky").children
    # We can find the group name
    @test ui_.children[3].children[1].children[1].children[1].children == dom"h2"("NamedTuple").children

    # To test manually
    # using Blink
    # w = Blink.Window()
    # body!(w, interface)
    # body!(w, submodel_interface)

    @testset "Test the tables interface and getproperty work on InteractModel too" begin
        df = DataFrame(interface)
        @test Tuple(df.val) == interface.val == interface[:val] == (0.05, 0.0, 20)
        @test Tuple(df.range) == interface.range == interface[:range] == (0.01:0.001:0.1, 0.0:0.1:6.2, 0.0:0.1:60.0)
        @test Tuple(df.label) == interface.label == interface[:label] == ("Sample step", "Phase", "Radus")
    end

end

@testset "empty models still work" begin
    InteractModel((named=:tuple, without=:params); ncolumns=2) do m
        hbox()
    end
    InteractModel((p1=Param(1), p2=Param(2)); submodel=Tuple, ncolumns=2) do m
        hbox()
    end
end


text:
# common tag generators
# NOTE: StrOrTag = Union{String, AbstractHtmlTag}

# ----------------------- horizental line & new line
tag_hr() = quickEmptySingleTag("hr")
tag_br() = quickEmptySingleTag("br")
# ----------------------- sup & sub
tag_sup(content::Union{StrOrTag,Vector}) = quickEmptyPairedTag("sup",content)
tag_sub(content::Union{StrOrTag,Vector}) = quickEmptyPairedTag("sub",content)
# ----------------------- bold font, italic, underline, del
tag_b(content::Union{StrOrTag,Vector}) = quickEmptyPairedTag("b",content)
tag_i(content::Union{StrOrTag,Vector}) = quickEmptyPairedTag("i",content)
tag_u(content::Union{StrOrTag,Vector}) = quickEmptyPairedTag("u",content)
tag_del(content::Union{StrOrTag,Vector}) = quickEmptyPairedTag("del",content)
# ----------------------- header
 -

tag_h1(content::Union{StrOrTag,Vector}) = quickEmptyPairedTag("h1",content)
tag_h2(content::Union{StrOrTag,Vector}) = quickEmptyPairedTag("h2",content)
tag_h3(content::Union{StrOrTag,Vector}) = quickEmptyPairedTag("h3",content)
tag_h4(content::Union{StrOrTag,Vector}) = quickEmptyPairedTag("h4",content)
tag_h5(content::Union{StrOrTag,Vector}) = quickEmptyPairedTag("h5",content)
tag_h6(content::Union{StrOrTag,Vector}) = quickEmptyPairedTag("h6",content)

# ----------------------- font block
"""
    tag_font( content::StrOrTag ; color::Union{String,Nothing} = nothing, face::Union{String,Nothing} = nothing, size::Union{String,Nothing} = nothing )

add a paired font tag to your content; returns a PairedHtmlTag instance;
not recommended since HTML 4.x;
StrOrTag = Union{String, AbstractHtmlTag}.
"""
function tag_font( content::StrOrTag ; color::StrOrNo = nothing, face::StrOrNo = nothing, size::StrOrNo = nothing )
    local tmppairedattrs = Dict{String,String}()
    isa(color, Nothing) ? nothing : tmppairedattrs["color"] = add_doublequotations(color)
    isa(face , Nothing) ? nothing : tmppairedattrs["face"] = add_doublequotations(face)
    isa(size , Nothing) ? nothing : tmppairedattrs["size"] = add_doublequotations(size)
    return PairedHtmlTag( "font", nothing, nothing, nothing,
        pairedAttributes = tmppairedattrs,
        content = AbstractHtmlTag[StrOrTag2Tag(content),] )
end # tag_font

# ----------------------- division block

"""
    tag_div( content::StrOrTag ; id::StrOrNo = nothing, class::StrOrNo = nothing, style::StrOrNo = nothing )

creates a
content
 tag PairedHtmlTag.
please organize your contents before packaging them in a div block.
"""
function tag_div( content::StrOrTag ; id::StrOrNo = nothing, class::StrOrNo = nothing, style::StrOrNo = nothing )
    return PairedHtmlTag( "div", id, class, style,
        content = AbstractHtmlTag[StrOrTag2Tag(content),] )
end # tag_div

# ----------------------- Hyperlink
"""
    tag_a( href::String, content::StrOrTag ; id::StrOrNo = nothing, class::StrOrNo = nothing, style::StrOrNo = nothing )

generating paired content tag string, returns a PairedHtmlTag.
"""
function tag_a( href::String, content::StrOrTag ;
    id::StrOrNo = nothing, class::StrOrNo = nothing, style::StrOrNo = nothing )
    # --------
    # preprare paired attributes
    local pairedattrs = Dict(
        "href" => add_doublequotations(href),
    )
    # make a PairedHtmlTag
    return PairedHtmlTag( "a", id, class, style,
        pairedAttributes = pairedattrs,
        content = AbstractHtmlTag[StrOrTag2Tag(content)] )
end # tag_a

# ----------------------- numerate list

"""
    tag_ol( items::Vector{T} where T <: StrOrTag ; id::StrOrNo = nothing, class::StrOrNo = nothing, style::StrOrNo = nothing, type::String = "1", start::Int = 1, flag_reversed::Bool = false )

1. :type can be one of [ "A", "a", "1", "I", "i" ]
2. :items is a list of listing items
3. :start should be an integer to indicate where to start counting
4. return a PairedHtmlTag
"""
function tag_ol( items::Vector{T} where T <: StrOrTag ;
    id::StrOrNo = nothing, class::StrOrNo = nothing, style::StrOrNo = nothing,
    type::String = "1", start::Int = 1, flag_reversed::Bool = false )
    # prepare list items
    local content::Vector = [ quickEmptyPairedTag("li", x) for x in items ]
    # new a tag instance
    return PairedHtmlTag(
        "ol", id, class, style,
            singleAttributes = String[( flag_reversed ? "reversed" : "" ),],
            pairedAttributes = Dict(
                "type" => add_doublequotations(type),
                "start" => add_doublequotations(string(start)),
            ),
            content = content
    )
end # tag_ol

# ----------------------- non-order list

"""
    tag_ul( items::Vector{T} where T <: StrOrTag ; id::StrOrNo = nothing, class::StrOrNo = nothing, style::StrOrNo = nothing, type::StrOrNo = nothing )

1. :type can be one of [ "default" "disc", "square", "circle" ]
2. :items is a list of listing items
4. return a PairedHtmlTag
"""
function tag_ul( items::Vector{T} where T <: StrOrTag ;
    id::StrOrNo = nothing, class::StrOrNo = nothing, style::StrOrNo = nothing,
    type::StrOrNo = nothing )
    # prepare list items
    local content::Vector = [ quickEmptyPairedTag("li", x) for x in items ]
    # new a tag instance
    local tmppairedattrs::Dict = isa(type, Nothing) ? Dict{String,String}() : Dict{String,String}( "type" => add_doublequotations(type), )
    return PairedHtmlTag( "ul", id, class, style,
        pairedAttributes = tmppairedattrs, content = content )
end # tag_ol

#

max_stars_repo_name: Clpr/HtmlConstructor

max_issues_repo_name: Clpr/HtmlConstructor

max_forks_repo_name: Clpr/HtmlConstructor


text:
using RunTests
using Base.Test

@testmodule EmitFunctionTests begin

using MetaTools

module MetaTestModule1
  module MetaTestModule2
    module MetaTestModule3
      function hello()
      end
    end
    function hello()
    end
  end
  function hello()
  end
end

test_functions = [(emit,),
                  (x->emit(ParsedFunction(emit(x))),)]

@parameterize test_functions function test_empty_method_with_no_args(test_function)
  pfunc = ParsedFunction(name=:hello1)
  func = eval(test_function(pfunc))
  @test isgeneric(func) == true
  @test func.env.name == :hello1
  @test func.env.max_args == 0
  @test func() == nothing
end

@parameterize test_functions function test_empty_method_with_one_untyped_arg(test_function)
  pfunc = ParsedFunction(name=:hello2,
                         args=[ParsedArgument(name=:x)])
  func = eval(test_function(pfunc))
  @test isgeneric(func) == true
  @test func.env.name == :hello2
  @test func.env.max_args == 1
  @test func(10) == nothing
end

@parameterize test_functions function test_empty_method_with_two_untyped_args(test_function)
  pfunc = ParsedFunction(name=:hello3,
                         args=[ParsedArgument(name=:x),
                               ParsedArgument(name=:y)])
  func = eval(test_function(pfunc))
  @test isgeneric(func) == true
  @test func.env.name == :hello3
  @test func.env.max_args == 2
  @test func(10, 20) == nothing
end

@parameterize test_functions function test_empty_method_with_one_typed_arg(test_function)
  pfunc = ParsedFunction(name=:hello4,
                         args=[ParsedArgument(name=:x, typ=:Integer)])
  func = eval(test_function(pfunc))
  @test isgeneric(func) == true
  @test func.env.name == :hello4
  @test func.env.max_args == 1
  @test func(10) == nothing
end

@parameterize test_functions function test_empty_method_with_two_typed_args(test_function)
  pfunc = ParsedFunction(name=:hello5,
                         args=[ParsedArgument(name=:x, typ=:Integer),
                               ParsedArgument(name=:y, typ=:Integer)])
  func = eval(test_function(pfunc))
  @test isgeneric(func) == true
  @test func.env.name == :hello5
  @test func.env.max_args == 2
  @test func(10, 20) == nothing
end

@parameterize test_functions function test_empty_method_with_one_union_typed_arg(test_function)
  pfunc = ParsedFunction(name=:hello6,
                         args=[ParsedArgument(:(x::Union(FloatingPoint, Integer)))])
  func = eval(test_function(pfunc))
  @test isgeneric(func) == true
  @test func.env.name == :hello6
  @test func.env.max_args == 1
  @test func(10) == nothing
end

@parameterize test_functions function test_empty_method_with_one_untyped_default_arg(test_function)
  pfunc = ParsedFunction(name=:hello7,
                         args=[ParsedArgument(:(x=100))])
  func = eval(test_function(pfunc))
  @test isgeneric(func) == true
  @test func.env.name == :hello7
  @test func.env.max_args == 1
  @test func() == nothing
  @test func(100) == nothing
end

@parameterize test_functions function test_empty_method_with_two_untyped_default_args(test_function)
  pfunc = ParsedFunction(name=:hello8,
                         args=[ParsedArgument(:(x=100)),
                               ParsedArgument(:(y=100))])
  func = eval(test_function(pfunc))
  @test isgeneric(func) == true
  @test func.env.name == :hello8
  @test func.env.max_args == 2
  @test func() == nothing
  @test func(100) == nothing
  @test func(100, 200) == nothing
end

@parameterize test_functions function test_empty_method_with_one_typed_default_arg(test_function)
  pfunc = ParsedFunction(name=:hello9,
                         args=[ParsedArgument(:(x::Integer=100))])
  func = eval(test_function(pfunc))
  @test isgeneric(func) == true
  @test func.env.name == :hello9
  @test func.env.max_args == 1
  @test func() == nothing
  @test func(100) == nothing
end

@parameterize test_functions function test_empty_method_with_two_typed_default_args(test_function)
  pfunc = ParsedFunction(name=:hello10,
                         args=[ParsedArgument(:(x::Integer=100)),
                               ParsedArgument(:(y::Integer=100))])
  func = eval(test_function(pfunc))
  @test isgeneric(func) == true
  @test func.env.name == :hello10
  @test func.env.max_args == 2
  @test func() == nothing
  @test func(100) == nothing
  @test func(100, 200) == nothing
end

@parameterize test_functions function test_empty_method_with_one_untyped_vararg(test_function)
  pfunc = ParsedFunction(name=:hello11,
                         args=[ParsedArgument(:(x...))])
  func = eval(test_function(pfunc))
  @test isgeneric(func) == true
  @test func.env.name == :hello11
  @test func.env.max_args == 0
  @test func() == nothing
  @test func(100) == nothing
  @test func(100, 200) == nothing
  @test func(100, 200, 300) == nothing
end

@parameterize test_functions function test_empty_method_with_one_typed_vararg(test_function)
  pfunc = ParsedFunction(name=:hello12,
                         args=[ParsedArgument(:(x::Integer...))])
  func = eval(test_function(pfunc))
  @test isgeneric(func) == true
  @test func.env.name == :hello12
  @test func.env.max_args == 0
  @test func() == nothing
  @test func(100) == nothing
  @test func(100, 200) == nothing
  @test func(100, 200, 300) == nothing
end

@parameterize test_functions function test_empty_method_with_one_untyped_default_vararg(test_function)
  pfunc = ParsedFunction(name=:hello13,
                         args=[ParsedArgument(:(x...=1))])
  func = eval(test_function(pfunc))
  @test isgeneric(func) == true
  @test func.env.name == :hello13
  @test func.env.max_args == 0
  @test func() == nothing
  @test func(100) == nothing
  @test func(100, 200) == nothing
  @test func(100, 200, 300) == nothing
end

@parameterize test_functions function test_empty_method_with_one_typed_default_vararg(test_function)
  pfunc = ParsedFunction(name=:hello14,
                         args=[ParsedArgument(:(x::Integer...=1))])
  func = eval(test_function(pfunc))
  @test isgeneric(func) == true
  @test func.env.name == :hello14
  @test func.env.max_args == 0
  @test func() == nothing
  @test func(100) == nothing
  @test func(100, 200) == nothing
  @test func(100, 200, 300) == nothing
end

@parameterize test_functions function test_empty_method_with_one_untyped_kwarg(test_function)
  pfunc = ParsedFunction(name=:hello15,
                        kwargs=[ParsedArgument(:(x=10))])
  func = eval(test_function(pfunc))
  @test isgeneric(func) == true
  @test func.env.name == :hello15
  @test func.env.max_args == 0
  @test func() == nothing
  @test func(x=100) == nothing
end

@parameterize test_functions function test_empty_method_with_two_untyped_kwargs(test_function)
  pfunc = ParsedFunction(name=:hello16,
                        kwargs=[ParsedArgument(:(x=10)),
                                 ParsedArgument(:(y=20))])
  func = eval(test_function(pfunc))
  @test isgeneric(func) == true
  @test func.env.name == :hello16
  @test func.env.max_args == 0
  @test func() == nothing
  @test func(x=100) == nothing
  @test func(y=200) == nothing
  @test func(x=100, y=200) == nothing
end

@parameterize test_functions function test_empty_method_with_one_typed_kwarg(test_function)
  pfunc = ParsedFunction(name=:hello17,
                        kwargs=[ParsedArgument(:(x::Integer=10))])
  func = eval(test_function(pfunc))
  @test isgeneric(func) == true
  @test func.env.name == :hello17
  @test func.env.max_args == 0
  @test func() == nothing
  @test func(x=100) == nothing
end

@parameterize test_functions function test_empty_method_with_two_typed_kwargs(test_function)
  pfunc = ParsedFunction(name=:hello18,
                        kwargs=[ParsedArgument(:(x::Integer=10)),
                                 ParsedArgument(:(y::Integer=20))])
  func = eval(test_function(pfunc))
  @test isgeneric(func) == true
  @test func.env.name == :hello18
  @test func.env.max_args == 0
  @test func() == nothing
  @test func(x=100) == nothing
  @test func(y=200) == nothing
  @test func(x=100, y=200) == nothing
end

@parameterize test_functions function test_empty_method_with_one_union_typed_kwarg(test_function)
  pfunc = ParsedFunction(name=:hello19,
                        kwargs=[ParsedArgument(:(x::Union(FloatingPoint, Integer)=10))])
  func = eval(test_function(pfunc))
  @test isgeneric(func) == true
  @test func.env.name == :hello19
  @test func.env.max_args == 0
  @test func() == nothing
  @test func(x=100) == nothing
end

@parameterize test_functions function test_empty_method_with_one_untyped_varkwarg(test_function)
  pfunc = ParsedFunction(name=:hello20,
                        kwargs=[ParsedArgument(:(x...))])
  func = eval(test_function(pfunc))
  @test isgeneric(func) == true
  @test func.env.name == :hello20
  @test func.env.max_args == 0
  @test func() == nothing
  @test func(x=100) == nothing
  @test func(y=200) == nothing
  @test func(x=100, y=200) == nothing
end

@parameterize test_functions function test_empty_method_with_one_typed_varkwarg(test_function)
  pfunc = ParsedFunction(name=:hello21,
                        kwargs=[ParsedArgument(:(x::Array...))])
  func = eval(test_function(pfunc))
  @test isgeneric(func) == true
  @test func.env.name == :hello21
  @test func.env.max_args == 0
  @test func() == nothing
  @test func(x=100) == nothing
  @test func(y=200) == nothing
  @test func(x=100, y=200) == nothing
end

@parameterize test_functions function test_empty_method_with_a_mix_of_args_and_kwargs(test_function)
  pfunc = ParsedFunction(name=:hello22,
                         args=[ParsedArgument(:x),
                               ParsedArgument(:(y::Integer)),
                               ParsedArgument(:(z=99))],
                        kwargs=[ParsedArgument(:(a::Integer=100)),
                                 ParsedArgument(:(b=200)),
                                 ParsedArgument(:(c...))])
  func = eval(test_function(pfunc))
  @test isgeneric(func) == true
  @test func.env.name == :hello22
  @test func.env.max_args == 3
  @test func(10, 20) == nothing
  @test func(10, 20, 30) == nothing
  @test func(10, 20, a=99) == nothing
  @test func(10, 20, 30, b=99, foobar=100) == nothing
end

@parameterize test_functions function test_empty_method_with_namespace_1_deep(test_function)
  pfunc = ParsedFunction(name=:hello,
                         namespace=[:MetaTestModule1])
  func = eval(test_function(pfunc))
  @test (func === MetaTestModule1.hello) == true
  @test isgeneric(func) == true
  @test func.env.name == :hello
  @test func.env.max_args == 0
  @test func() == nothing
end

@parameterize test_functions function test_empty_method_with_namespace_2_deep(test_function)
  pfunc = ParsedFunction(name=:hello,
                         namespace=[:MetaTestModule1, :MetaTestModule2])
  func = eval(test_function(pfunc))
  @test (func === MetaTestModule1.MetaTestModule2.hello) == true
  @test isgeneric(func) == true
  @test func.env.name == :hello
  @test func.env.max_args == 0
  @test func() == nothing
end

@parameterize test_functions function test_empty_method_with_namespace_3_deep(test_function)
  pfunc = ParsedFunction(name=:hello,
                         namespace=[:MetaTestModule1, :MetaTestModule2, :MetaTestModule3])
  func = eval(test_function(pfunc))
  @test (func === MetaTestModule1.MetaTestModule2.MetaTestModule3.hello) == true
  @test isgeneric(func) == true
  @test func.env.name == :hello
  @test func.env.max_args == 0
  @test func() == nothing
end

@parameterize test_functions function test_empty_method_with_no_type_parameter(test_function)
  pfunc = ParsedFunction(name=:hello23,
                         types=Symbol[])
  func = eval(test_function(pfunc))
  @test isgeneric(func) == true
  @test func.env.name == :hello23
  @test func.env.max_args == 0
  @test func() == nothing
end

@parameterize test_functions function test_empty_method_with_one_type_parameter(test_function)
  pfunc = ParsedFunction(name=:hello24,
                         types=[:T],
                         args=[ParsedArgument(:(x::T))])
  func = eval(test_function(pfunc))
  @test isgeneric(func) == true
  @test func.env.name == :hello24
  @test func.env.max_args == 1
  @test func(10) == nothing
end

@parameterize test_functions function test_empty_method_with_two_type_parameters(test_function)
  pfunc = ParsedFunction(name=:hello25,
                         types=[:T1, :T2],
                         args=[ParsedArgument(:(x::T1)),
                               ParsedArgument(:(y::T2))])
  func = eval(test_function(pfunc))
  @test isgeneric(func) == true
  @test func.env.name == :hello25
  @test func.env.max_args == 2
  @test func(10, 20) == nothing
end

@parameterize test_functions function test_empty_method_with_namespace_3_deep_and_one_type_parameter(test_function)
  pfunc = ParsedFunction(name=:hello,
                         types=[:T],
                         namespace=[:MetaTestModule1, :MetaTestModule2, :MetaTestModule3])
  func = eval(test_function(pfunc))
  @test (func === MetaTestModule1.MetaTestModule2.MetaTestModule3.hello) == true
  @test isgeneric(func) == true
  @test func.env.name == :hello
  @test func.env.max_args == 0
  @test func() == nothing
end
#############################################################################

@parameterize test_functions function test_empty_anonymous_function_with_no_args(test_function)
  pfunc = ParsedFunction()
  func = eval(test_function(pfunc))
  @test func() == nothing
end

@parameterize test_functions function test_empty_anonymous_function_with_one_untyped_arg(test_function)
  pfunc = ParsedFunction(args=[ParsedArgument(:x)])
  func = eval(test_function(pfunc))
  @test func(10) == nothing
end

@parameterize test_functions function test_empty_anonymous_function_with_two_untyped_args(test_function)
  pfunc = ParsedFunction(args=[ParsedArgument(:x),
                               ParsedArgument(:y)])
  func = eval(test_function(pfunc))
  @test isgeneric(func) == false
  @test func(10, 20) == nothing
end

@parameterize test_functions function test_empty_anonymous_function_with_one_typed_arg(test_function)
  pfunc = ParsedFunction(args=[ParsedArgument(:(x::Integer))])
  func = eval(test_function(pfunc))
  @test isgeneric(func) == false
  @test func(10) == nothing
end

@parameterize test_functions function test_empty_anonymous_function_with_two_typed_args(test_function)
  pfunc = ParsedFunction(args=[ParsedArgument(:(x::Integer)),
                               ParsedArgument(:(y::Integer))])
  func = eval(test_function(pfunc))
  @test isgeneric(func) == false
  @test func(10, 20) == nothing
end

@parameterize test_functions function test_empty_anonymous_function_with_one_union_typed_arg(test_function)
  pfunc = ParsedFunction(args=[ParsedArgument(:(x::Union(FloatingPoint, Integer)))])
  func = eval(test_function(pfunc))
  @test isgeneric(func) == false
  @test func(10) == nothing
end

@parameterize test_functions function test_empty_anonymous_function_with_one_untyped_vararg(test_function)
  pfunc = ParsedFunction(args=[ParsedArgument(:(x...))])
  func = eval(test_function(pfunc))
  @test isgeneric(func) == false
  @test func() == nothing
  @test func(100) == nothing
  @test func(100, 200) == nothing
  @test func(100, 200, 300) == nothing
end

@parameterize test_functions function test_empty_anonymous_function_with_one_typed_vararg(test_function)
  pfunc = ParsedFunction(args=[ParsedArgument(:(x::Integer...))])
  func = eval(test_function(pfunc))
  @test isgeneric(func) == false
  @test func() == nothing
  @test func(100) == nothing
  @test func(100, 200) == nothing
  @test func(100, 200, 300) == nothing
end

@parameterize test_functions function test_empty_anonymous_function_with_a_mix_of_args(test_function)
  pfunc = ParsedFunction(args=[ParsedArgument(:x),
                               ParsedArgument(:(y::Integer)),
                               ParsedArgument(:(z...))])
  func = eval(test_function(pfunc))
  @test isgeneric(func) == false
  @test func(10, 20) == nothing
  @test func(10, 20, 30) == nothing
  @test func(10, 20, 30, 40) == nothing
end

#############################################################################

@parameterize test_functions function test_method_with_a_numeric_literal_body(test_function)
  pfunc = ParsedFunction(name=:hello100,
                         body=Expr(:block, 10))
  func = eval(test_function(pfunc))
  @test isgeneric(func) == true
  @test func.env.name == :hello100
  @test func.env.max_args == 0
  @test func() == 10
end

@parameterize test_functions function test_anonymous_function_with_a_numeric_literal_body(test_function)
  pfunc = ParsedFunction(body=Expr(:block, 10))
  func = eval(test_function(pfunc))
  @test isgeneric(func) == false
  @test func() == 10
end
#----------------------------------------------------------------------------

@parameterize test_functions function test_method_with_a_symbol_body(test_function)
  pfunc = ParsedFunction(name=:hello101,
                         args=[ParsedArgument(:x)],
                         body=Expr(:block, :x))
  func = eval(test_function(pfunc))
  @test isgeneric(func) == true
  @test func.env.name == :hello101
  @test func.env.max_args == 1
  @test func(10) == 10
  @test func(100) == 100
end

@parameterize test_functions function test_anonymous_function_with_a_symbol_body(test_function)
  pfunc = ParsedFunction(args=[ParsedArgument(:x)],
                         body=Expr(:block, :x))
  func = eval(test_function(pfunc))
  @test isgeneric(func) == false
  @test func(10) == 10
  @test func(100) == 100
end
#----------------------------------------------------------------------------

@parameterize test_functions function test_method_with_a_expression_body(test_function)
  pfunc = ParsedFunction(name=:hello102,
                         args=[ParsedArgument(:x)],
                         body=:(x + 10))
  func = eval(test_function(pfunc))
  @test isgeneric(func) == true
  @test func.env.name == :hello102
  @test func.env.max_args == 1
  @test func(10) == 20
  @test func(100) == 110
end

@parameterize test_functions function test_anonymous_function_with_a_expression_body(test_function)
  pfunc = ParsedFunction(args=[ParsedArgument(:x)],
                         body=:(x + 10))
  func = eval(test_function(pfunc))
  @test isgeneric(func) == false
  @test func(10) == 20
  @test func(100) == 110
end
#----------------------------------------------------------------------------

@parameterize test_functions function test_method_with_a_multiline_body(test_function)
  pfunc = ParsedFunction(name=:hello103,
                         args=[ParsedArgument(:x)],
                         body=quote
                           y = x + 10
                           z = x - 10
                           return z//y
                         end)
  func = eval(test_function(pfunc))
  @test isgeneric(func) == true
  @test func.env.name == :hello103
  @test func.env.max_args == 1
  @test func(100) == 9//11
end

@parameterize test_functions function test_anonymous_function_with_a_multiline_body(test_function)
  pfunc = ParsedFunction(args=[ParsedArgument(:x)],
                         body=quote
                           y = x + 10
                           z = x - 10
                           return z//y
                         end)
  func = eval(test_function(pfunc))
  @test isgeneric(func) == false
  @test func(100) == 9//11
end

end