text:
using Pkg
 Pkg.activate(".")

verbose = true
 if verbose && !isdefined(Main,:RvSpectML)  println("# Loading RvSpecML")    end
 using RvSpectML
 include("shared/scripts.jl")
 if verbose   println("# Loading other packages")    end
 using DataFrames, Query, Statistics, Dates

include("shared/scripts.jl")
include("read_expres_data_101501.jl")

order_list_timeseries = extract_orders(all_spectra,pipeline_plan)

line_list_df = prepare_line_list_pass1(linelist_for_ccf_filename, all_spectra, pipeline_plan,
         v_center_to_avoid_tellurics=ccf_mid_velocity, Δv_to_avoid_tellurics = 21e3)

#need_to!(pipeline_plan, :template)
if need_to(pipeline_plan, :template)  # Compute order CCF's & measure RVs
   if verbose println("# Making template spectra.")  end
   @assert !need_to(pipeline_plan,:extract_orders)
#   @assert !need_to(pipeline_plan,:rvs_ccf_total)
   GC.gc()   # run garbage collector for deallocated memory
   map(i->order_list_timeseries.metadata[i][:rv_est] = 0.0, 1:length(order_list_timeseries) )
   @time ( spectral_orders_matrix, f_mean, var_mean, deriv, deriv2, order_grids )  = make_template_spectra(order_list_timeseries)
   if save_data(pipeline_plan, :template)
      #using CSV
      #CSV.write(joinpath(output_dir,target_subdir * "_template.csv"),DataFrame("λ"=>spectral_orders_matrix.λ,"flux_template"=>f_mean,"var"=>var_mean, "dfluxdlnλ_template"=>deriv,"d²fluxdlnλ²_template"=>deriv2))
      using JLD2, FileIO
      save(joinpath(output_dir,target_subdir * "_matrix.jld2"), Dict("λ"=>spectral_orders_matrix.λ,"spectra"=>spectral_orders_matrix.flux,"var_spectra"=>spectral_orders_matrix.var,"flux_template"=>f_mean,"var"=>var_mean, "dfluxdlnλ_template"=>deriv,"d²fluxdlnλ²_template"=>deriv2))
   end
   dont_need_to!(pipeline_plan, :template);
end

if make_plot(pipeline_plan,:template)
   using Plots
   chunkid = 10
   idx = spectral_orders_matrix.chunk_map[chunkid]
   plt = scatter(spectral_orders_matrix.λ[idx],(f_mean[idx].-1.0)./maximum(abs.((f_mean[idx].-1.0))),markersize=1.0,label=:none)
   scatter!(plt,spectral_orders_matrix.λ[idx],deriv[idx]./maximum(abs.(deriv[idx])),markersize=1.0,label=:none)
   scatter!(plt,spectral_orders_matrix.λ[idx],deriv2[idx]./maximum(abs.(deriv2[idx])),markersize=1.0,label=:none)
   xlabel!("λ (Å)")
   ylabel!("f(λ), f'(λ), f''(λ), all standardized")
   title!("Template spectrum for chunk " * string(chunkid) )
end

#need_to!(pipeline_plan, :dcpca)
if need_to(pipeline_plan, :dcpca)
   if verbose println("# Performing Doppler constrained PCA analysis.")  end
   #@assert !need_to(pipeline_plan,:rvs_ccf_total)
   @assert !need_to(pipeline_plan,:template)
   using MultivariateStats
   Δrvs_dcpca = calc_rvs_from_taylor_expansion(spectral_orders_matrix; mean=f_mean, deriv=deriv).rv
   #map(i->order_list_timeseries.metadata[i][:rv_est] = Δrvs_dcpca[i], 1:length(order_list_timeseries) )
   dcpca_out, M = DCPCA.doppler_constrained_pca(spectral_orders_matrix.flux, deriv, Δrvs_dcpca)
   frac_var_unexplained = 1.0.-cumsum(principalvars(M))./tvar(M)
   num_basis = length(frac_var_unexplained)
   println("# Fraction of variance unexplained = ", frac_var_unexplained[1:min(5,num_basis)])
   if save_data(pipeline_plan, :dcpca)
      using CSV
      CSV.write(joinpath(output_dir,target_subdir * "_dcpca_basis.csv"),  Tables.table(M.proj) )
      CSV.write(joinpath(output_dir,target_subdir * "_dcpca_scores.csv"), Tables.table(dcpca_out) )
   end
   dont_need_to!(pipeline_plan, :dcpca);
end

if make_plot(pipeline_plan, :dcpca)  # Ploting results from DCPCA
  include("../scripts/plots/dcpca.jl")
  #using Plots
  # Set parameters for plotting analysis
  plt = scatter(frac_var_unexplained, xlabel="Number of PCs", ylabel="Frac Variance Unexplained")
  if save_plot(pipeline_plan,:dcpca)   savefig(plt,joinpath(output_dir,target_subdir * "_dcpca_frac_var.png"))   end
  display(plt)
end

if make_plot(pipeline_plan, :dcpca)
  plt_order = 10
  plt = plot_basis_vectors(order_grids, f_mean, deriv, M.proj, idx_plt = spectral_orders_matrix.chunk_map[plt_order], num_basis=2 )
  #xlims!(5761.5,5766)
  if save_plot(pipeline_plan,:dcpca)   savefig(plt,joinpath(output_dir,target_subdir * "_dcpca_basis.png"))   end
  display(plt)
end

if make_plot(pipeline_plan, :dcpca)
  plt = plot_basis_scores(order_list_timeseries.times, Δrvs_dcpca, dcpca_out, num_basis=min(5,num_basis) )
  if save_plot(pipeline_plan,:dcpca)   savefig(plt,joinpath(output_dir,target_subdir * "_dcpca_scores.png"))   end
  display(plt)
end

if make_plot(pipeline_plan, :dcpca)
  plt = plot_basis_scores_cor( Δrvs_dcpca, dcpca_out, num_basis=min(5,num_basis))
  if save_plot(pipeline_plan,:dcpca)   savefig(plt,joinpath(output_dir,target_subdir * "_dcpca_scores_cor.png"))   end
  display(plt)
end

#need_to!(pipeline_plan,:fit_lines)
if need_to(pipeline_plan,:fit_lines)
   if verbose println("# Performing fresh search for lines in template spectra.")  end
   cl = ChunkList(map(grid->ChunkOfSpectrum(spectral_orders_matrix.λ,f_mean, var_mean, grid), spectral_orders_matrix.chunk_map))
   #= # We're done with the spectral_orders_matrix, so we can release the memory now
   spectral_orders_matrix = nothing
   GC.gc()
   need_to!(pipeline_plan,:template)
   =#
   lines_in_template = LineFinder.find_lines_in_chunklist(cl, plan=RvSpectML.LineFinder.LineFinderPlan(min_deriv2=3))  # TODO: Automate threshold for finding a line

   if verbose println("# Finding above lines in all spectra.")  end
   @time fits_to_lines = LineFinder.fit_all_lines_in_chunklist_timeseries(order_list_timeseries, lines_in_template )

   if save_data(pipeline_plan,:fit_lines)
      using CSV
      CSV.write(joinpath(output_dir,target_subdir * "_linefinder_lines.csv"), lines_in_template )
      CSV.write(joinpath(output_dir,target_subdir * "_linefinder_line_fits.csv"), fits_to_lines )
      #CSV.write(joinpath(output_dir,target_subdir * "_linefinder_line_fits_clean.csv"), lines_to_try )
   end
   dont_need_to!(pipeline_plan,:fit_lines);
end

#need_to!(pipeline_plan,:clean_line_list_blends)
if need_to(pipeline_plan,:clean_line_list_blends)
   # Exploratory data analysis of the distribution of line properties over time to figure out how to select "clean" lines
   @assert !need_to(pipeline_plan,:fit_lines)
   fit_distrib = fits_to_lines |> @groupby(_.line_id) |>
            @map( { median_a=median(_.fit_a), median_b=median(_.fit_b), median_depth=median(_.fit_depth), median_σ²=median(_.fit_σ²), median_λc=median(_.fit_λc),
                   std_a=std(_.fit_a), std_b=std(_.fit_b), std_depth=std(_.fit_depth), std_σ²=std(_.fit_σ²), std_λc=std(_.fit_λc), line_id=_.line_id, frac_converged=mean(_.fit_converged)  } ) |>
            @filter(_.frac_converged == 1.0 ) |> DataFrame
                  # , min_telluric_model_all_obs=minimum(_.min_telluric_model_this_obs),
   good_lines = fit_distrib |>
            #@filter(_.min_telluric_model_all_obs == 1.0 ) |> # Already done
            @filter( _.median_depth > 0.1 ) |>
            @filter( _.median_σ² < 0.015) |>
            @filter( _.std_σ² < 0.004) |>
            @filter(_.std_b < 0.06) |>
            @filter( _.std_depth < 0.006 ) |>
            #@filter(_.std_a < 0.025) |>             # Unncessary?
            #@filter( -0.25 < _.median_b < 0.25) |>  # Unncessary?
            DataFrame
      println("# Found ", size(lines_in_template,1), " lines, including ",size(good_lines,1), " good lines, rejected ", size(lines_in_template,1)-size(good_lines,1), " lines.")
      bad_lines = fit_distrib |> @filter(_.std_λc >= 0.002 ) |> DataFrame
      good_lines_high_scatter = good_lines |> @filter(_.std_λc >= 0.002 ) |> DataFrame
      if size(bad_lines,1) >= 1
         println("# ", size(bad_lines,1), " lines have large λc scatter, including ", size(good_lines_high_scatter,1), " good lines.")
      end
      lines_to_try = lines_in_template[first.(good_lines[!,:line_id]),:]
      dont_need_to!(pipeline_plan,:clean_line_list_blends);
 end
 fit_distrib

# Reset steps of pipeline_plan to rerun with new linelist.
need_to!(pipeline_plan, :ccf_total);
 need_to!(pipeline_plan, :rvs_ccf_total);
 need_to!(pipeline_plan, :ccf_orders);
 need_to!(pipeline_plan, :rvs_ccf_orders);
 need_to!(pipeline_plan, :template);
 need_to!(pipeline_plan, :dcpca);

RvSpectML.dont_save_data!(pipeline_plan,:ccf_total)
idx = sortperm(lines_to_try.fit_λc)
line_list_clean_df = DataFrame(:lambda=>lines_to_try.fit_λc[idx], :weight=>lines_to_try.fit_depth[idx] )
(ccfs, v_grid) = ccf_total(order_list_timeseries, line_list_df, pipeline_plan,  mask_scale_factor=1.0,
                           ccf_mid_velocity=ccf_mid_velocity, recalc=true)
(ccfs2, v_grid2) = ccf_total(order_list_timeseries, line_list_clean_df, pipeline_plan,  mask_scale_factor=1.0,
                           ccf_mid_velocity=0, recalc=true)

if make_plot(pipeline_plan, :ccf_total)
   using Plots
   t_idx = 20
   using Plots
   plt = plot(v_grid.-ccf_mid_velocity,ccfs[:,t_idx]./maximum(ccfs[:,t_idx],dims=1),label=:none)
   scatter!(plt,v_grid2,ccfs2[:,t_idx]./maximum(ccfs2[:,t_idx],dims=1),markersize=1.2,label=:none)
   xlabel!("v (m/s)")
   ylabel!("CCF")
   if save_plot(pipeline_plan,:ccf_total)   savefig(plt,joinpath(output_dir,target_subdir * "_ccf2_sum.png"))   end
   display(plt)
end

#need_to!(pipeline_plan,:rvs_ccf_total)
# TODO: Replace using ccf's here with calc_rvs_from_taylor_expansion using new template?
# Δrvs_dcpca = calc_rvs_from_taylor_expansion(spectral_orders_matrix; mean=f_mean, deriv=deriv).rv

rvs_ccf = calc_rvs_from_ccf_total(ccfs, pipeline_plan, v_grid=v_grid, times=order_list_timeseries.times,recalc=true)
rvs_ccf2 = calc_rvs_from_ccf_total(ccfs2, pipeline_plan, v_grid=v_grid2, times=order_list_timeseries.times,recalc=true)

#need_to!(pipeline_plan, :template)
if need_to(pipeline_plan, :template)
   @assert !need_to(pipeline_plan,:rvs_ccf_total)
   @assert !need_to(pipeline_plan,:clean_line_list_blends)
   # Revert to velocities before cleaning for now
   map(i->order_list_timeseries.metadata[i][:rv_est] = rvs_ccf2[i]-mean(rvs_ccf2), 1:length(order_list_timeseries) )
   #chunk_list_df2 = lines_to_try |> @select(:fit_min_λ,:fit_max_λ) |> @rename(:fit_min_λ=>:lambda_lo, :fit_max_λ=>:lambda_hi) |> DataFrame
   expand_chunk_factor = 4
   chunk_list_df2 = lines_to_try |> @select(:fit_λc,:fit_min_λ,:fit_max_λ) |> @map({:lambda_lo=>_.fit_λc-expand_chunk_factor*(_.fit_λc-_.fit_min_λ), :lambda_hi=>_.fit_λc+expand_chunk_factor*(_.fit_max_λ-_.fit_λc)}) |> DataFrame
   chunk_list_timeseries2 = make_chunk_list_timeseries(all_spectra,chunk_list_df2)

   # Check that no NaN's included
   #(chunk_list_timeseries2, chunk_list_df2) = filter_bad_chunks(chunk_list_timeseries2,chunk_list_df2)
   chunk_list_timeseries2 = filter_bad_chunks(chunk_list_timeseries2)
   println(size(chunk_list_df2), " vs ", num_chunks(chunk_list_timeseries2) )

   GC.gc()   # run garbage collector for deallocated memory
   if verbose println("# Making template spectra.")  end
   @time ( spectral_orders_matrix2, f_mean2, var_mean2, deriv_2, deriv2_2, order_grids2 )  = make_template_spectra(chunk_list_timeseries2)

   if save_data(pipeline_plan, :template)
      #using CSV
      #CSV.write(joinpath(output_dir,target_subdir * "_template2.csv"),DataFrame("λ"=>spectral_orders_matrix2.λ,"flux_template"=>f_mean2,"var"=>var_mean2, "dfluxdlnλ_template"=>deriv_2,"d²fluxdlnλ²_template"=>deriv2_2))
      using JLD2, FileIO
      save(joinpath(output_dir,target_subdir * "_matrix2.jld2"), Dict("λ"=>spectral_orders_matrix2.λ,"spectra"=>spectral_orders_matrix2.flux,"var_spectra"=>spectral_orders_matrix2.var,"flux_template"=>f_mean2,"var"=>var_mean2, "dfluxdlnλ_template"=>deriv_2,"d²fluxdlnλ²_template"=>deriv2_2))
   end
   dont_need_to!(pipeline_plan, :template)
end

need_to!(pipeline_plan, :dcpca)
 if need_to(pipeline_plan, :dcpca)
   @assert !need_to(pipeline_plan,:rvs_ccf_total)
   @assert !need_to(pipeline_plan,:template)
   if verbose println("# Performing Doppler constrained PCA analysis.")  end
   using MultivariateStats
   # Second EXPRES observation of 101501 is very weird. Let's try leaving it out of the DCPCA analysis.
   obs_incl_for_dcpca = vcat([1], collect(3:size(spectral_orders_matrix2.flux,2)) )
   spectral_orders_matrix2 = make_spectral_time_series_common_wavelengths_with_selected_times(spectral_orders_matrix2, obs_incl_for_dcpca)
   Δrvs_dcpca2 = calc_rvs_from_taylor_expansion(spectral_orders_matrix2; mean=f_mean2, deriv=deriv_2).rv
   dcpca2_out, M2 = DCPCA.doppler_constrained_pca(spectral_orders_matrix2.flux, deriv_2, Δrvs_dcpca2)
   frac_var_unexplained2 = 1.0.-cumsum(principalvars(M2))./tvar(M2)
   num_basis = length(frac_var_unexplained2)
   println("# Fraction of variance unexplained (orig) = ", frac_var_unexplained[1:min(5,length(frac_var_unexplained))])
   println("# Fraction of variance unexplained (clean) = ", frac_var_unexplained2[1:min(5,num_basis)])
   if save_data(pipeline_plan, :dcpca)
      using CSV
      CSV.write(joinpath(output_dir,target_subdir * "_dcpca_basis2.csv"),  Tables.table(M2.proj) )
      CSV.write(joinpath(output_dir,target_subdir * "_dcpca_scores2.csv"), Tables.table(dcpca2_out) )
   end
   dont_need_to!(pipeline_plan, :dcpca)
end

if make_plot(pipeline_plan, :dcpca)  # Ploting results from DCPCA
  using Plots
  plt = scatter(frac_var_unexplained2, xlabel="Number of PCs", ylabel="Frac Variance Unexplained")
  if save_plot(pipeline_plan,:dcpca)   savefig(plt,joinpath(output_dir,target_subdir * "_dcpca2_frac_var.png"))   end
  display(plt)
end

if make_plot(pipeline_plan, :dcpca)
   plt_line = 10
   plt = plot_basis_vectors(order_grids2, f_mean2, deriv_2, M2.proj, idx_plt = spectral_orders_matrix2.chunk_map[plt_line], num_basis=min(4,num_basis), label=plt_line)
   if save_plot(pipeline_plan,:dcpca)
     savefig(plt,joinpath(output_dir,target_subdir * "_dcpca2_basis.png"))
   end
   # Making animations is slow.  # TODO: Add animations key to pipeline_plan?
   #= suspicous_lines = Int[]
   anim = @animate for plt_line ∈ 1:length(order_grids2)
      local idx = spectral_orders_matrix2.chunk_map[plt_line]
      argmin_pixel = argmin(f_mean2[idx])
      if 0.45*length(idx) <= argmin_pixel <= 0.55*length(idx)
        plot_basis_vectors(order_grids2, f_mean2, deriv_2, M2.proj, idx_plt = spectral_orders_matrix2.chunk_map[plt_line], num_basis=min(4,num_basis), label=plt_line  )
     else
        push!(suspicous_lines,plt_line)
     end
   end
   gif(anim, joinpath(output_dir,"dcpca2_basis.gif"), fps = 25)
   =#

  display(plt)
end

if make_plot(pipeline_plan, :dcpca)
  plt = plot_basis_scores(order_list_timeseries.times[obs_incl_for_dcpca], rvs_ccf2[obs_incl_for_dcpca], dcpca2_out, num_basis=min(4,num_basis) )
  if save_plot(pipeline_plan,:dcpca)   savefig(plt,joinpath(output_dir, target_subdir * "_dcpca2_scores.png"))   end
  display(plt)

end

if make_plot(pipeline_plan, :dcpca)
  plt = plot_basis_scores_cor( rvs_ccf2[obs_incl_for_dcpca], dcpca2_out, num_basis=min(4,num_basis))
  if save_plot(pipeline_plan,:dcpca)   savefig(plt,joinpath(output_dir,target_subdir * "_dcpca2_score_cor.png"))   end
  display(plt)
end


text:
using SimpleMock
using Documenter

makedocs(;
    modules=[SimpleMock],
    authors="Chris de Graaf <chrisadegraaf@gmail.com>",
    repo="https://github.com/christopher-dG/SimpleMock.jl/blob/{commit}{path}#L{line}",
    sitename="SimpleMock.jl",
    format=Documenter.HTML(;
        canonical="https://christopher-dG.github.io/SimpleMock.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/christopher-dG/SimpleMock.jl",
)

ext: jl


text:
#Code written by Skipper af Rosenborg and Esben Bagerx
"""
findMaxCor(standX, bestSolution)
Returns the maximum correlation form the selected variables
"""
function findMaxCor(standX, bestSolution)
	chosenVar = find(bestSolution)
	#Shrink values
	for i = chosenVar
		if bestSolution[i] < 1e-6 && bestSolution[i] > -1e-6
			bestSolution[i] = 0
		end
	end
	chosenVar = find(bestSolution)
	return maximum(cor(standX[:,chosenVar])-Diagonal(cor(standX[:,chosenVar])))
end

"""
Function that returns the zscore column zscore for each column in the Matrix.
Must have at least two columns
"""
function zScoreByColumn(X)
	standX = Array{Float64}(X)
	for i=1:size(X)[2]
		temp = zscore(X[:,i])
		if any(isnan.(temp))

		else
			standX[:,i] = copy(zscore(X[:,i]))
		end
	end
	return standX
end

"""
Function an array with the kMax largest absolute values
and all other values shrunk to 0
"""
function shrinkValuesH(betaVector, kMax, HCArray)
	#Make aboslute copy of betaVector
	absCopy = copy(abs.(betaVector))
	#Create zeroVector
	zeroVector = zeros(betaVector)

	for i in 1:kMax
		#Find index of maximum value in absolute vector
		ind = indmax(absCopy)

		#Replace index in zeroVector with betaVector value of index
		zeroVector[ind] = betaVector[ind]

		#Replace index in absolute vector with 0
		absCopy[ind] = 0

        #Check if index is in HC array and add all
        HCvar = find(HCArray[:,1] .== ind) #Maybe need to check both sides

        if !isempty(HCvar)
            for j in 1:length(HCvar)
				absCopy[Int64(HCArray[HCvar[j],2])] = 0
            end
        end

		#Find place in column
		nColOrig = size(zeroVector)[1]/4

		if ind%nColOrig == 0
			indNum = nColOrig
		else
			indNum = ind%nColOrig
		end
		placeInColumn = floor(ind/nColOrig)

		for j in 0:3
			if j != placeInColumn
				absCopy[Int64(indNum+j*nColOrig)] = 0
			end
		end
	end
	return zeroVector
end

"""
Function to calculate two norm regression error
"""
function twoNormRegressionError(X, y, betaVector)
	return 1/2*norm(y-X*betaVector)^2
end

"""
Vanilla gradient decent which only keeps the kMax biggest values. All other
values are shrunk to 0.
"""
function gradDecent(X, y, L, epsilon, kMax, HC, bSolved, HCArray)
	#If first iteration intialise bSolve randomly
	if countnz(bSolved) < 1
		bSolved = rand(size(X)[2])
	end

    betaVector = shrinkValuesH(bSolved,kmax, HCArray)

	#Initialise parameters
	iter = 0
	oldBetaVector = copy(betaVector)
	curError = twoNormRegressionError(X, y, betaVector) - twoNormRegressionError(X, y, oldBetaVector) + 10000

	counter = 0

	while (iter < 2000 && curError > epsilon)
		if counter == 100
			println("Broke in iteration ",iter," due to counter")
			break
		end

		oldBetaVector = copy(betaVector)

		#Calculate delta(g(beta))
		gradBeta = X'*(y-X*betaVector)

		#Shrink smalles values
		betaVector = copy(shrinkValuesH(betaVector+1/L*gradBeta, kMax, HCArray))
		if find(oldBetaVector) != find(betaVector)
			#println("New beta")
		end

		if find(oldBetaVector) == find(betaVector)
			counter +=1
		else
			count = 0
		end

		curError = abs.(twoNormRegressionError(X, y, oldBetaVector) - twoNormRegressionError(X, y, betaVector))
		iter += 1

		if iter%1000 == 0
			println("Iteration $iter with error on $curError")
		end
	end
	#println("End error: $curError")
	#println("Number of iterations: $iter")
	return  betaVector
end

"""
Function to that returns the R^2 value
"""
function getRSquared(X,y,betaSolve)
    SSres = sum((y[i] - X[i,:]'*betaSolve)^2 for i=1:length(y))
    SSTO = sum((y[i]-mean(y))^2 for i=1:length(y))
    Rsquared = 1-SSres/SSTO
    return Rsquared
end

function getRMSE(X,y,betaSolve)
    SSres = sum((y[i] - X[i,:]'*betaSolve)^2 for i=1:length(y))
	n=length(y)
	RMSE = sqrt(SSres/n)
	return RMSE
end

"""
Function that prints the non zero solutions values of beta
"""
function printNonZeroValues(bSolved)
    for i=1:length(bSolved)
        if !isequal(bSolved[i], 0)
            println("b[$i] = ", bSolved[i])
        end
    end
end

"""
Transforming X as: Squared, log & square-root
"""
function expandWithTransformations(X)
	println("Transforming as: Squared, log & square-root")
	expandedX = copy(Array{Float64}(X))
	xCols = size(expandedX)[2]
	xRows = size(expandedX)[1]

	# squared transformation column
	for i=1:xCols
		insertArray = []
		for j=1:xRows
			push!(insertArray, expandedX[j,i]^2)
		end
		expandedX = hcat(expandedX, insertArray)
	end

	# Natural log transformation column
	for i=1:xCols
		insertArray = []
		for j=1:xRows
			if count(k->(k<0), expandedX[j,i]) > 0
				push!(insertArray, 0)
			elseif count(k->(k==0), expandedX[j,i]) > 0
				push!(insertArray, log(expandedX[j,i]+0.00001))
			else
				push!(insertArray, log(expandedX[j,i]))
			end
		end
		expandedX = hcat(expandedX, insertArray)
	end

	# sqrt transformation column
	for i=1:xCols
		insertArray = []
		for j=1:xRows
			if count(k->(k<=0), expandedX[j,i]) > 0
				push!(insertArray, 0)
			else
				push!(insertArray, sqrt(expandedX[j,i]))
			end
		end
		expandedX = hcat(expandedX, insertArray)
	end

	#Ensure we return a Float64 array
	expandedX = copy(Array{Float64}(expandedX))
	return expandedX
end

"""
Transforming X with time elements -3, -6 and -12
"""
function expandWithTime3612(X)
	println("Transforming with t-1, t-3, ..., t-12")
	expandedX = copy(Array{Float64}(X))
	xCols = size(expandedX)[2]
	xRows = size(expandedX)[1]

	for k=1:12
		# k = how many month we look back: -1:-12 time series
		for i=1:xCols
			# i loops over all columns in the X matrix
			insertArray = []

			# j loops over rows in X matrix
			for j=1:k
				# adding 0-rows for the first k rows as we do not have any
				# previous observations
				push!(insertArray, 0)
			end

			for j=k+1:xRows
				push!(insertArray, expandedX[j-k,i])
			end
			expandedX = hcat(expandedX, insertArray)
		end
	end

	#Ensure we return a Float64 array
	expandedX = copy(Array{Float64}(expandedX))
	return expandedX
end

"""
Transforming X with moving average and momentum signals
"""
function expandWithMAandMomentum(X, Y, originalColumns)
	println("Transforming with momentum (9, 12) and MA(1,2,3 & 9,12) ")
	#the first row is -1 already, so it will be row +2, +5 and +11
	expandedX = copy(Array{Float64}(X))
	xCols = originalColumns
	xRows = size(expandedX)[1]
	#Already have moving average for s=1, since that our base row.

	#s=2 would be the average of the last 2 observations
	MAs2arr = []
	for j=1:2
		push!(MAs2arr,0)
	end
	for j=1:(xRows-2)
		push!(MAs2arr, sum(Y[s] for s=(j):(j+1))/2)
	end

	#s=3 would be the average of the last 3 observations
	MAs3arr = []
	for j=1:3
		push!(MAs3arr,0)
	end
	for j=1:(xRows-3)
		push!(MAs3arr, sum(Y[s] for s=(j):(j+2))/3)
	end

	#l=9 would be the average of the last 9 observations
	MAl9arr = []
	for j=1:9
		push!(MAl9arr,0)
	end
	for j=1:(xRows-9)
		push!(MAl9arr, sum(Y[s] for s=(j):(j+8))/9)
	end

	#l=12 would be the average of the last 12 observations
	MAl12arr = []
	for j=1:12
		push!(MAl12arr,0)
	end
	for j=1:(xRows-12)
		push!(MAl12arr, sum(Y[s] for s=(j):(j+11))/12)
	end

	#Compare s=1 and l=9 and create a column of buy signals
	#then do the same for s=2, 3 vs l=9
	#then l=12 vs. s=1,2,3
	MAMatrix = zeros(xRows, 6)
	#check s=1 and l=9
	for j=10:xRows
		if Y[j-1] >= MAl9arr[j]
			MAMatrix[j, 1] = 1
		end
	end
	#check s=2 and l=9
	for j=10:xRows
		if MAs2arr[j] >= MAl9arr[j]
			MAMatrix[j, 2] = 1
		end
	end
	#check s=3 and l=9
	for j=10:xRows
		if MAs3arr[j] >= MAl9arr[j]
			MAMatrix[j, 3] = 1
		end
	end

	#check s=1 and l=12
	for j=13:xRows
		if Y[j-1] >= MAl12arr[j]
			MAMatrix[j, 4] = 1
		end
	end

	#check s=2 and l=12
	for j=13:xRows
		if MAs2arr[j] >= MAl12arr[j]
			MAMatrix[j, 5] = 1
		end
	end
	#check s=3 and l=12
	for j=13:xRows
		if MAs3arr[j] >= MAl12arr[j]
			MAMatrix[j, 6] = 1
		end
	end

	expandedX = hcat(expandedX, MAMatrix)

	# MOMENTUM
	# If a stock is higher then it was m periods ago, it has MOMENTUM
	# computing for m=9
	Mom9arr = []
	for j=1:10
		push!(Mom9arr,0)
	end
	for j=11:xRows
		if Y[j-1] > Y[j-10]
			push!(Mom9arr, 1)
		else
			push!(Mom9arr, 0)
		end
	end

	expandedX = hcat(expandedX, Mom9arr)

	# computing for m=12
	Mom12arr = []
	for j=1:13
		push!(Mom12arr,0)
	end
	for j=14:xRows
		if Y[j-1] > Y[j-13]
			push!(Mom12arr, 1)
		else
			push!(Mom12arr, 0)
		end
	end

	expandedX = hcat(expandedX, Mom12arr)

	#Ensure we return a Float64 array
	expandedX = copy(Array{Float64}(expandedX))
	return expandedX
end

function expandColNamesTimeToString(featureNames)
	loopLength = length(featureNames)

	for j = 1:12
		for i=1:loopLength
			#println(Symbol(string(featureNames[i])*" t-"*string(j)))
			featureNames = vcat(featureNames, Symbol(string(featureNames[i])*" t-"*string(j)))
		end
	end

	return featureNames
end

"""
Function that returns all the column names including ^2, log(), sqrt()
"""
function expandedColNamesToString(colNames, ExpTrans, TA)
	masterString = ""

	if ExpTrans == 1
		for i=1:length(colNames)
			masterString = masterString*string(colNames[i])*","
		end

		for i=1:length(colNames)
			masterString = masterString*string(colNames[i])*"^2,"
		end

		for i=1:length(colNames)
			masterString = masterString*"ln("*string(colNames[i])*"),"
		end

		for i=1:length(colNames)
			masterString = masterString*"sqrt("*string(colNames[i])*"),"
		end
	end

	if TA == 1
		masterString = masterString * "MA 1-9,"
		masterString = masterString * "MA 2-9,"
		masterString = masterString * "MA 3-9,"
		masterString = masterString * "MA 1-12,"
		masterString = masterString * "MA 2-12,"
		masterString = masterString * "MA 3-12,"
		masterString = masterString * "Mo 9,"
		masterString = masterString * "MA 12"
	end

	return masterString
end

"""
Function to identify which parameters have been selected (transformed or not)
"""
function identifyParameters(betaSolution, colNames)
	originalColumns = Int64(length(betaSolution)/4)
	# test original parameters
	count = 1
	for i=1:originalColumns
		if !isequal(betaSolution[i],0)
			println("Parameter '", colNames[count],"' has been selected with value ", round(betaSolution[i],3))
		end
		count += 1
	end

	# test for squared
	count = 1
	for i=(originalColumns+1):(2*originalColumns)
		if !isequal(betaSolution[i],0)
			println("Parameter '", colNames[count],"^2' has been selected with value ", round(betaSolution[i],3))
		end
		count += 1
	end

	# test for natural log
	count = 1
	for i=(originalColumns*2+1):(3*originalColumns)
		if !isequal(betaSolution[i],0)
			println("Parameter 'ln(", colNames[count],")' has been selected with value ", round(betaSolution[i],3))
		end
		count += 1
	end
	# test for sqrt
	count = 1
	for i=(originalColumns*3+1):(4*originalColumns)
		if !isequal(betaSolution[i],0)
			println("Parameter 'sqrt(", colNames[count],")' has been selected with value ", round(betaSolution[i],3))
		end
		count += 1
	end
end

"""
Functions to make a split in data
"""
function createSampleX(x, inputRows)
	inputX = copy(x)
	outputX = inputX[inputRows,:]
	return outputX
end

function createSampleY(y, inputRows)
	inputY = copy(y)
	outputY = inputY[inputRows,:]
	return outputY
end

"""
Creates a sample of random rows without replacement.
Not useful for financial data, if we want a timeseries.
"""

function selectSampleRows(rowsWanted, nRows)
	rowsSelected = sample(1:nRows, rowsWanted, replace = false)
	return rowsSelected
end

function selectSampleRowsWR(rowsWanted, nRows)
	rowsSelected = rand(1:nRows, rowsWanted)
	return rowsSelected
end

function splitDataIn2(data, rowsWanted, nRows)
	rowsOne = selectSampleRows(rowsWanted, nRows)
	rowsTwo = []
	for i = 1:nRows
		if !(i in rowsOne)
			push!(rowsTwo, i)
		end
	end
	dataSplitOne = data[rowsOne, :]
	dataSplitTwo = data[rowsTwo, :]

	return dataSplitOne, dataSplitTwo
end

"""
Generating samplesize amount of beta-values for a beta distribution
"""
function createBetaDistribution(bSample, standX, standY, k, sampleSize, rowsPerSample, gamma, allCuts, bestZ)
	println("0% through samples")
	empty1 = 0
	empty2 = 0
	continueBool = false
	for i=1:sampleSize
		continueBool = false
		sampleRows = selectSampleRowsWR(rowsPerSample, nRows)
		sampleX = createSampleX(standX, sampleRows)
		sampleY = createSampleY(standY, sampleRows)
		HC = cor(sampleX) # this is added to support function below
		#println("Maximum pairwise correlation is: ",findMaxCor(standX, bestZ))

		bSample[i,:], empty1, empty2 = buildAndSolveStage2(sampleX, sampleY, k, gamma, false, [], allCuts, bestZ, HC)
		if empty1 == 9999
			newStandX = sampleX
			newStandY = sampleY
			newK = k
			newGamma = gamma
			newAllCuts = allCuts
			newBestZ = bestZ
			newHC = HC
		end

		#bSample[i,:] = solveForBeta(sampleX, sampleY, k, gamma, allCuts, bestZ); #This is replaced by above function
		#bSample[i,:] = solveForBetaClosedForm(sampleX, sampleY, k, bestZ)
		if i == floor(sampleSize/2)
			println("50% through samples, $i missing")
		end
	end
	println("100% through samples")
	return bSample
end

"""
Creates a confidence interval based on confidence interval level and nBoot
bootstrapped samples from the created beta distribution
"""
function createConfidenceIntervalArray(sampleInput, nBoot, confLevel)
	bSample = convert(Array{Float64}, sampleInput)
	bCols = size(bSample)[2]
	confIntArray = Matrix(2, bCols)
	for i=1:bCols
		bs = bootstrap(bSample[:,i], mean, BasicSampling(nBoot))
		cil = confLevel
		CiEst = Bootstrap.ci(bs, BasicConfInt(cil))
		confIntArray[1,i] = CiEst[1][2] #lower
		confIntArray[2,i] = CiEst[1][3]
	end
	return confIntArray
end
"""
Based on the confidence intervals created, significance of found parameters can
be tested with this function
"""
function testSignificance(confIntArray99, confIntArray95, confIntArray90, bResult)
	bCols = size(bResult)[1]
	significance = zeros(bCols)
	for i=1:bCols
		if bResult[i] >= confIntArray99[1,i] && bResult[i] <=confIntArray99[2,i] && confIntArray99[1,i] > 0
			significance[i] = 0.99
		elseif bResult[i] >= confIntArray99[1,i] && bResult[i] <=confIntArray99[2,i] && confIntArray99[2,i] < 0
			significance[i] = 0.99
		elseif bResult[i] >= confIntArray95[1,i] && bResult[i] <=confIntArray95[2,i] && confIntArray95[1,i] > 0
			significance[i] = 0.95
		elseif bResult[i] >= confIntArray95[1,i] && bResult[i] <=confIntArray95[2,i] && confIntArray95[2,i] < 0
			significance[i] = 0.95
		elseif bResult[i] >= confIntArray90[1,i] && bResult[i] <=confIntArray90[2,i] && confIntArray90[1,i] > 0
			significance[i] = 0.90
		elseif bResult[i] >= confIntArray90[1,i] && bResult[i] <=confIntArray90[2,i] && confIntArray90[2,i] < 0
			significance[i] = 0.90
		else
			significance[i] = 0
		end
	end
	return significance
end

"""
SUPPORT FUNCTIONS FOR CEO AND LASSO MEAN VARIANCE
"""
function createDataSplits(XArrays, YArrays, t, trainingSize)
    trainingXArrays = Array{Array{Float64, 2}}(industriesTotal)
    trainingYArrays = Array{Array{Float64, 2}}(industriesTotal)

    validationXRows = Array{Array{Float64, 2}}(industriesTotal)
    validationY     = Array{Array{Float64, 2}}(industriesTotal)

    OOSXArrays      = Array{Array{Float64, 2}}(industriesTotal)
    OOSYArrays      = Array{Array{Float64, 2}}(industriesTotal)

    OOSRow          = Array{Array{Float64, 2}}(industriesTotal)
    OOSY            = Array{Array{Float64, 2}}(industriesTotal)
    for i = 1:industriesTotal
        trainingXArrays[i] = XArrays[i][t:(t-1+trainingSize), :]
        trainingYArrays[i] = YArrays[i][t:(t-1+trainingSize), :]

        validationXRows[i] = XArrays[i][(t+trainingSize):(t+trainingSize), :]
        validationY[i]     = YArrays[i][(t+trainingSize):(t+trainingSize), :]

        OOSXArrays[i]      = XArrays[i][(t+1):(t+trainingSize), :]
        OOSYArrays[i]      = YArrays[i][(t+1):(t+trainingSize), :]

        OOSRow[i]          = XArrays[i][(t+trainingSize+1):(t+trainingSize+1), :]
        OOSY[i]            = YArrays[i][(t+trainingSize+1):(t+trainingSize+1), :]
    end
    return trainingXArrays, trainingYArrays, validationXRows, validationY, OOSXArrays, OOSYArrays, OOSRow, OOSY
end

function calculatePvalue(return1N, returnPerfect, returnCEO)
    PValue = 1 - (returnCEO-returnPerfect)/(return1N-returnPerfect)
    return PValue
end

function findPerfectResults(trainingXArrays, Xrow, Yvalues, gamma)
    indexes = 10
    Sigma =  cov(trainingXArrays[1][:,1:10])

    #A=U^(T)U where U is upper triangular with real positive diagonal entries
    F = lufact(Sigma)
    U = F[:U]  #Cholesky factorization of Sigma

    M = JuMP.Model(solver = GurobiSolver(OutputFlag = 0))
    @variables M begin
            w[1:indexes]
            u[1:indexes]
            z
            y
    end

    @objective(M,Min, gamma*y - Yvalues'*w)
    @constraint(M, 0 .<= w)
    @constraint(M, sum(w[i] for i=1:indexes) == 1)
    @constraint(M, norm([2*U'*w;y-1]) <= y+1)
    solve(M)
    wPerfect = getvalue(w)
    forecastRow = (exp(Xrow')-1)*100
    periodPerfectReturn = forecastRow*wPerfect

    return wPerfect, periodPerfectReturn
end

function findPerfectRFRResults(trainX, Xrow, Yvalues, gamma)
    indexes = 11
    Sigma =  cov(trainX)

    #A=U^(T)U where U is upper triangular with real positive diagonal entries
    F = lufact(Sigma)
    U = F[:U]  #Cholesky factorization of Sigma

    M = JuMP.Model(solver = GurobiSolver(OutputFlag = 0))
    @variables M begin
            w[1:indexes]
            u[1:indexes]
            z
            y
    end

    @objective(M,Min, gamma*y - Yvalues'*w)
    @constraint(M, 0 .<= w)
    @constraint(M, sum(w[i] for i=1:indexes) == 1)
    @constraint(M, norm([2*U'*w;y-1]) <= y+1)
    solve(M)
    wPerfect = getvalue(w)
    forecastRow = (exp(Xrow')-1)*100
    periodPerfectReturn = forecastRow*wPerfect

    return wPerfect, periodPerfectReturn
end

function performMVOptimization(expectedReturns, U, gamma, Xrow, Yvalues)
    indexes = 10
    M = JuMP.Model(solver = GurobiSolver(OutputFlag = 0))
    @variables M begin
            w[1:indexes]
            u[1:indexes]
            z
            y
    end

    @objective(M,Min, gamma*y - expectedReturns'*w)
    @constraint(M, 0 .<= w)
    @constraint(M, sum(w[i] for i=1:indexes) == 1)
    @constraint(M, norm([2*U'*w;y-1]) <= y+1)
    solve(M)
    wStar = getvalue(w)

    forecastRow = (exp.(Xrow)-1)*100

    periodReturn = forecastRow'*wStar
    period1NReturn = forecastRow'*w1N

    return period1NReturn, periodReturn, wStar, forecastRow
end

function performMVOptimizationRISK(expectedReturns, U, gamma, Xrow, Yvalues)
    indexes = 11
    M = JuMP.Model(solver = GurobiSolver(OutputFlag = 0))
    @variables M begin
            w[1:indexes]
            u[1:indexes]
            z
            y
    end

    @objective(M,Min, gamma*y - expectedReturns'*w)
    @constraint(M, 0 .<= w)
    @constraint(M, sum(w[i] for i=1:indexes) == 1)
    @constraint(M, norm([2*U'*w;y-1]) <= y+1)
    solve(M)
    wStar = getvalue(w)

	forecastRow = zeros(11)
	forecastRow[1:10] = (exp.(Xrow[1:10])-1)*100
	forecastRow[11]   = (exp.(Xrow[11])-1)

    periodReturn = forecastRow'*wStar
    period1NReturn = forecastRow'*w1N

    return period1NReturn, periodReturn, wStar, forecastRow
end

function generateExpectedReturns(betaArray, trainingXArrays, trainingYArrays, validationXRows)
    industriesTotal = 10

    expectedReturns = zeros(10)
    bRows = 240
    bCols = floor(Int,zeros(industriesTotal))
    for i = 1:industriesTotal
        bCols[i] = size(trainingXArrays[i])[2]
    end

    for i = 1:10
        errorArray = (trainingXArrays[i]*betaArray[i,1:bCols[i]]-trainingYArrays[i])*(1/bRows)
        errorSum = sum(errorArray)
        expectedReturns[i] = validationXRows[i][:]'*betaArray[i,1:bCols[i]]
    end
    return expectedReturns
end

function runCEO(trainingXArrays, trainingYArrays, modelConfigRow, gamma)
    industriesTotal = 10
    periodMean = zeros(10)
    for i = 1:10
        periodMean[i] = mean(trainingXArrays[1][:,i])
    end
    l1, l2, l3, l4 = modelConfigRow
    env = Gurobi.Env()

    maxPredictors = 1412
    bCols = floor(Int,zeros(industriesTotal))
    Sigma =  cov(trainingXArrays[1][:,1:10])

    #A=U^(T)U where U is upper triangular with real positive diagonal entries
    F = lufact(Sigma)
    U = F[:U]  #Cholesky factorization of Sigma

    for i = 1:industriesTotal
        bCols[i] = size(trainingXArrays[i])[2]
    end
    gamma = gamma
    bRows = size(trainingXArrays[1])[1]
    M = JuMP.Model(solver = GurobiSolver(env, OutputFlag = 0, Threads=(nprocs()-1)))

    @variables M begin
            b[1:industriesTotal, 1:maxPredictors] #Beta in LASSO; industry i and predictors 1:1304
            z[1:industriesTotal, 1:maxPredictors] #auxilliary for 1norm
            q[1:industriesTotal] #auxilliary for norm
            f[1:industriesTotal, 1:bRows] #expected return for industry i at time t
            w[1:industriesTotal, 1:bRows] #weights for industry i at time t
            y[1:bRows] #auxilliary variables for MV optimization
    end

    @objective(M,Min,sum(l1*q[i] + l2*ones(maxPredictors)'*z[i,:] for i=1:industriesTotal)+ l3*(sum(gamma*y[t]-sum(f[i,t] for i=1:10) for t=1:bRows)) - l4*(sum(w[i,t]*trainingYArrays[i][t] for i=1:10, t=1:bRows)))

    for i = 1:industriesTotal
        for t = 1:bRows
            @constraint(M, (periodMean[i]+1000)*w[i,t] >= f[i,t])
        end
    end
    for i = 1:industriesTotal
        @constraint(M, norm( [1-q[i];2*(trainingXArrays[i]*b[i,1:bCols[i]]-trainingYArrays[i])] ) <= 1+q[i]) #second order cone constraint SOC
    end
    @constraint(M,  b .<= z)
    @constraint(M, -z .<= b)

    for i = 1:10
        errorArray = (trainingXArrays[i]*b[i,1:bCols[i]]-trainingYArrays[i])*(1/bRows)
        errorSum = sum(errorArray)
        for t = 1:bRows
            prediction = trainingXArrays[i][t, :]'*b[i,1:bCols[i]]
            @constraint(M, prediction + errorSum >= f[i, t]) ##== originally
        end
    end

    @constraint(M, 0 .<= w)
    for t = 1:bRows
        @constraint(M, sum(w[i,t] for i=1:industriesTotal) == 1)
        @constraint(M, norm([2*U'*w[:,t];y[t]-1]) <= y[t]+1)
    end

    @time(solve(M))
    betaArray = getvalue(b)
    #Insert shrinking here
    return betaArray, U
end

function generateXandYs(industries, modelMatrix)
    industriesTotal = length(industries)
    for i = 1:industriesTotal
        industry = industries[i]
        mainData = loadIndexDataLOGReturn(industry, path)

        dateAndRecession = Array(mainData[:,end-1:end])
        mainDataArr = Array(mainData[:,1:end-2])

        colNames = names(mainData)

        nRows = size(mainDataArr)[1]
        nCols = size(mainDataArr)[2]

        fileName = path*"Results/IndexData/CEOTest/"*industry*"/"*string(trainingSize)*"-1/"*string(trainingSize)*"_"

        mainYarr = mainDataArr[:, nCols:nCols]

        mainXarr = mainDataArr[:,1:nCols-1]

        VIX = modelMatrix[i, 1]
        raw = modelMatrix[i, 2]
        expTrans = modelMatrix[i, 3]
        timeTrans = modelMatrix[i, 4]
        TA = modelMatrix[i,5]

        if raw == 1
            fileName = fileName*"Raw"
            mainXarr = mainXarr[:,1:10]
        else
            if VIX == 1
                fileName = fileName*"VIX_Macro"
                mainXarr = mainXarr[:,1:end-1]
            else
                fileName = fileName*"Macro"
            end
        end

        # Transform with time elements
        if timeTrans == 1
            fileName = fileName*"Time"
            mainXarr = expandWithTime3612(mainXarr)
        end

        # Transform #
        if expTrans == 1
            fileName = fileName*"Exp"
            mainXarr = expandWithTransformations(mainXarr)
        end

        if TA == 1
            fileName = fileName*"TA"
            mainXarr = expandWithMAandMomentum(mainXarr, mainYarr, (nCols-1))
        end

        # Standardize
        standX = zScoreByColumn(mainXarr)
        standY = mainYarr

        XArrays[i] = standX
        YArrays[i] = standY
    end
    return XArrays, YArrays
end

max_stars_repo_path: ML/Code to hand in/SupportFunction.jl

max_stars_repo_name: skipperafrosenborg/Thesis

max_issues_repo_path: ML/Code to hand in/SupportFunction.jl

max_issues_repo_name: skipperafrosenborg/Thesis

max_forks_repo_path: ML/Code to hand in/SupportFunction.jl

max_forks_repo_name: skipperafrosenborg/Thesis


text:
module Define

using Tokenize

include("utilities.jl")

export @cdef
include("Cdef.jl")

end


text:
import Avalon: CPU, GPU

device = GPU()
cpu = CPU()

@testset "cuda: conv" begin
    # 1D
    x = rand(7, 3, 10); w = rand(3, 3, 1)
    g = grad((x, w) -> sum(conv1d(x, w)), x, w)[2]
    d_g = grad((x, w) -> sum(conv1d(x, w)), device(x), device(w))[2]
    @test g[2] ≈ cpu(d_g[2])
    @test g[3] ≈ cpu(d_g[3])

    # 2D
    x = rand(7, 7, 3, 10); w = rand(3, 3, 3, 1)
    g = grad((x, w) -> sum(conv2d(x, w)), x, w)[2]
    d_g = grad((x, w) -> sum(conv2d(x, w)), device(x), device(w))[2]
    @test g[2] ≈ cpu(d_g[2])
    @test g[3] ≈ cpu(d_g[3])

    # 3D
    x = rand(7, 7, 7, 3, 10); w = rand(3, 3, 3, 3, 1)
    g = grad((x, w) -> sum(conv3d(x, w)), x, w)[2]
    d_g = grad((x, w) -> sum(conv3d(x, w)), device(x), device(w))[2]
    @test g[2] ≈ cpu(d_g[2])
    @test g[3] ≈ cpu(d_g[3])
end

@testset "cuda: pooling" begin
    x = rand(7, 7, 3, 10);
    g = grad(x -> sum(maxpool2d(x, 2)), x)[2]
    d_g = grad(x -> sum(maxpool2d(x, 2)), device(x))[2]
    @test g[2] ≈ cpu(d_g[2])
end

@testset "cuda: batchnorm" begin
    Random.seed!(898);

    m = BatchNorm2d(3)
    x = rand(10, 10, 3, 5)
    d_m = device(m)
    d_x = device(x)

    y = m(x)
    d_y = d_m(d_x)
    @test isapprox(cpu(d_y), y; rtol=1e-2)

    dy = randn(size(x)...)
    d_dy = device(dy)

    dgamma, dbeta, dx = ∇batchnorm2d(dy, m, m.gamma, m.beta, x)
    d_dgamma, d_dbeta, d_dx = ∇batchnorm2d(d_dy, d_m, d_m.gamma, d_m.beta, d_x)
    @test isapprox(cpu(d_dgamma), dgamma; rtol=1e-2)
    @test isapprox(cpu(d_dbeta), dbeta; rtol=1e-2)
    @test isapprox(cpu(d_dx), dx; rtol=1e-2)

    _, g = grad((m, x) -> sum(m(x)), m, x)
    _, d_g = grad((m, x) -> sum(m(x)), d_m, d_x)

    @test isapprox(g[2].gamma, cpu(d_g[2].gamma); rtol=1e-2, atol=1e-2)
    @test isapprox(g[2].beta, cpu(d_g[2].beta); rtol=1e-2, atol=1e-4)
    @test isapprox(g[3], cpu(d_g[3]); rtol=1e-2, atol=1e-4)

end

@testset "cuda: RNN" begin
    Random.seed!(128);

    # vanilla RNN
    m = RNN(10 => 5); x_seq = ones(10, 4, 10); h = init_hidden(m, 4)
    _, g = grad((m, x_seq, h) -> begin h_all, h = m(x_seq, h); sum(h_all) end, m, x_seq, h)
    d_m, d_x_seq, d_h = map(device, (m, x_seq, h))
    _, d_g = grad((m, x_seq, h) -> begin h_all, h = m(x_seq, h); sum(h_all) end, d_m, d_x_seq, d_h)
    @test g[2].cell.W_ih ≈ cpu(d_g[2].cell.W_ih)

    # LSTM
    m = LSTM(3 => 2); x_seq = ones(3, 2, 3); h, c = init_hidden(m, 2)
    _, g = grad((m, x_seq, h, c) -> begin h_all, h, c = m(x_seq, h, c); sum(h) end, m, x_seq, h, c)
    d_m, d_x_seq, d_h, d_c = map(device, (m, x_seq, h, c))
    _, d_g = grad((m, x_seq, h, c) -> begin h_all, h, c = m(x_seq, h, c); sum(h) end, d_m, d_x_seq, d_h, d_c)
    @test g[2].cell.W_ih ≈ cpu(d_g[2].cell.W_ih)

    # GRU
    m = GRU(3 => 2); x_seq = ones(3, 2, 3); h = init_hidden(m, 2)
    _, g = grad((m, x_seq, h) -> begin h_all, h = m(x_seq, h); sum(h_all) end, m, x_seq, h)
    d_m, d_x_seq, d_h = map(device, (m, x_seq, h))
    _, d_g = grad((m, x_seq, h) -> begin h_all, h = m(x_seq, h); sum(h_all) end, d_m, d_x_seq, d_h)
    @test g[2].cell.W_ih ≈ cpu(d_g[2].cell.W_ih)

end

@testset "cuda: activations" begin
    x = rand(Float32, 5, 5);
    d_x = device(x);

    @test grad(x -> sum(logistic.(x)), x)[2][2] ≈ grad(x -> sum(logistic.(x)), d_x)[2][2] |> cpu
    @test grad(x -> sum(softplus.(x)), x)[2][2] ≈ grad(x -> sum(softplus.(x)), d_x)[2][2] |> cpu
    @test grad(x -> sum(softsign.(x)), x)[2][2] ≈ grad(x -> sum(softsign.(x)), d_x)[2][2] |> cpu
    @test grad(x -> sum(relu.(x)), x)[2][2] ≈ grad(x -> sum(relu.(x)), d_x)[2][2] |> cpu
    @test grad(x -> sum(leakyrelu.(x, 0.01)), x)[2][2] ≈
        grad(x -> sum(leakyrelu.(x, 0.01)), d_x)[2][2] |> cpu
    # ELU on CUDA results in scalar operations warning followed by segfault, disabling it for now
    # @test grad(x -> sum(elu.(x, 1.0)), x)[2][1] ≈ grad(x -> sum(elu.(x, 1.0)), d_x)[2][1]

    g = grad(x -> sum(softmax(x)), x)[2][2]
    d_g = grad(x -> sum(softmax(x)), d_x)[2][2]
    @test isapprox(g, cpu(d_g), rtol = 1e-5, atol = 1e-5)

    @test grad(x -> sum(logsoftmax(x)), x)[2][2] ≈ grad(x -> sum(logsoftmax(x)), d_x)[2][2] |> cpu
end

@testset "cuda: losses" begin
    x = rand(5, 4); x = log.(x ./ sum(x; dims=1)); c = [3, 2, 1, 4]
    d_x = device(x); d_c = device(c)
    @test grad((x, c) -> nllloss(x, c), x, c)[2][2] ≈
        grad((x, c) -> nllloss(x, c), d_x, d_c)[2][2] |> cpu

    x = rand(5, 4); c = [3, 2, 1, 4]
    d_x = device(x); d_c = device(c)
    @test (grad((x, c) -> crossentropyloss(x, c), x, c)[2][2] ≈
           grad((x, c) -> crossentropyloss(x, c), d_x, d_c)[2][2] |> cpu)

    x = rand(5, 4); x_target = rand(5, 4)
    d_x = device(x); d_x_target = device(x_target)
    @test (grad((x, x_target) -> mseloss(x, x_target), x, x_target)[2][2] ≈
           grad((x, x_target) -> mseloss(x, x_target), d_x, d_x_target)[2][2] |> cpu)
end

@testset "cuda: optim" begin

    # not every parameter update will lead to descreased loss
    # but at least we can check that parameters are actually changed
    m = MyModel(Linear(5, 4)) |> device; x = rand(5, 10) |> device;
    old_m = deepcopy(m); old_x = deepcopy(x)
    _, g = grad(my_model_loss, m, x)

    # SGD
    update!(SGD(0.1; momentum=0.5), m, g[2])
    @test old_m.linear.W != m.linear.W
    update!(SGD(0.1; momentum=0.5), x, g[3])
    @test old_x != x

    # RMSprop
    update!(RMSprop(), m, g[2])
    @test old_m.linear.W != m.linear.W
    update!(RMSprop(), x, g[3])
    @test old_x != x

    # Adam
    update!(Adam(), m, g[2])
    @test old_m.linear.W != m.linear.W
    update!(Adam(), x, g[3])
    @test old_x != x

end


text:
# run with `julia scripts/precook.jl example.toml` to download precomputed spectra

using TOML
configfile = ARGS[1]
config = TOML.parsefile(configfile)

##
using Downloads
using CSV, DataFrames

run_name = config["general"]["name"]
mapids = [k for k in keys(config["map"])]
spectrapath = joinpath(config["scratch"], "rawspectra")
mkpath(spectrapath)

function download_if_necessary(url, dest; verbose=true)
    if isfile(dest) == false
        verbose && println("Downloading ", dest)
        @time Downloads.download(url, dest)
    else
        verbose && println("Extant, skip ", dest)
    end
end

spec_url_root = "https://github.com/xzackli/PSPipePlanckRender.jl/releases/download/v0.1/"
spec_file = "rawspectra.tar.gz"
spec_url = spec_url_root * spec_file
spec_dest = joinpath(config["scratch"], spec_file)
download_if_necessary(spec_url, spec_dest)

run(`tar -xzvf $(spec_dest) --overwrite -C $(config["scratch"])`)

max_forks_repo_name: sgiardie/PSpipe


text:
using OrthogonalPolynomials
using Test, Base.Math

# Laguerre
# Abramowitz and Stegun pg. 799
const TriangleCoeffsLaguerre = [
    1 1 1 2 6 24 120 720 5040 40320 362880 3628800 39916800 479001600;
    1 1 1 -4 -18 -96 -600 -4320 -35280 -322560 -3265920 -36288000 -439084800 -6748019200;
    2 2 -4 2 18 144 1200 10800 105840 1128960 13063680 163296000 2195424000 31614106600;
    6 6 -18 9 6 -96 -1200 -14400 -176400 -2257920 -30481920 -435456000 -6586272000 -105380352000;
    24 24 -96 72 -16 24 600 10800 176400 2822400 45722880 762048000 13172544000 237105792000;
    120 120 -600 600 -200 25 120 -4320 -106840 -2257920 -45722880 -914457600 -18441561600 -379369267200;
    720 720 -4320 5400 -2400 450 -36 720 35280 1128960 30481920 762048000 18441561600 442597478400;
    5040 5040 -35280 52920 -29400 7350 -882 49 5040 -322560 -13063680 -435456000 -13172544000 -379369267200;
    40320 40320 -322560 564480 -376320 117600 -18816 1568 -64 40320 3269520 163296000 6586272000 237105792000;
    362880 362880 -3265920 6531840 -5080320 1905120 -381024 42336 -2592 81 362880 -36288000 -2195424000 -105380352000;
    3623800 3628800 -36288000 81648000 -72576000 31752000 -7620480 1058400 -86400 4050 -100 3628800 439084800 31614106600;
    39916800 39916800 -439084800 1097712000 -1097712000 648806000 -103679680 25613280 -2613600 163350 -6050 121 39916800 -5748019200;
    479001600 479001600 -5748019200 15807052800 -17563392000 9879408000 -3161410560 614718720 -75271680 5880600 -290400 8712 -144 479001600]
L0(x) = 1.0
    L1(x) = -x + 1.0
    L2(x) = .5 * (x^2 - 4x + 2)
    L3(x) = (1/6)*(-x^3 +9x^2 -18x + 6)
    L4(x) = (1/24)*(x^4 - 16x^3 + 72x^2 - 96x + 24)
    L5(x) = (1/120)*(-x^5 + 25x^4 -200x^3 + 600x^2 - 600x + 120)
    L6(x) = (1/720)*(x^6 - 36x^5 + 450x^4 - 2400x^3 + 5400x^2 - 4320x + 720)
    L7(x) = (1/5040)*(-x^7 + 49x^6 - 882x^5 + 7350x^4 - 29_400x^3 + 52_920x^2 - 35_280x + 5040)
    L8(x) = (1/40320)*(x^8 - 64x^7 + 1568x^6 - 18_816x^5 + 117_600x^4 - 376_320x^3 + 564_480x^2 - 322_560x + 40320)
    L9(x) = (1/362880)*(-x^9 + 81x^8 - 2592x^7 + 42_336x^6 - 381_024x^5 + 1_905_120x^4 - 5_080_320x^3 + 6_531_840x^2 - 3_265_920x + 362_880)
    L10(x) = (1/3628800)*(x^10 - 100x^9 + 4_050x^8 - 86_400x^7 + 1_058_400x^6 - 7_620_480x^5 + 31_752_000x^4 - 72_576_000x^3  + 81_648_000x^2 - 36_288_000x + 3_628_800)
    L11(x) = (1/39916800)*(-x^11 + 121x^10 - 6_050x^9 + 163_350x^8 - 2_613_600x^7 + 25_613_280x^6 - 153_679_680x^5 + 548_856_000x^4 - 1_097_712_000x^3 + 1_097_712_000x^2 - 439_084_800x + 39_916_800)
    L12(x) = (1/479001600)*(x^12 - 144x^11 + 8_712x^10 - 290_400x^9 + 5_880_600x^8 - 75_271_680x^7 + 614_718_720x^6 - 3_161_410_560x^5 + 9_879_408_000x^4 - 17_563_392_000x^3 + 15_807_052_800x^2 - 5_748_019_200x + 479_001_600)

    Ls = [L0, L1, L2, L3, L4, L5, L6, L7, L8, L9, L10, L11, L12]

    # Credit to @yingboma
    function laguerre(n, x)

        n == 0 ? (return 1) : nothing

        l0, l1 = 1, 1-x
        for k in 1:n-1
            l0, l1 = l1, ( (2k+1-x)*l1 - k*l0 )/(k+1)
        end
        l1
    end
#n\z is the first column - Abramowitz and Stegun p. 800
const testvalslaguerre = [0. .5 1 3 5 10];
const TestTableLaguerre = [
    0 1.0000000000 1.0000000000 1.0000000000 1.0000000000 1.0000000000;
    1 0.5000000000 0.0000000000 -2.0000000000 -4.0000000000 -9.0000000000;
    2 0.1250000000 -0.5000000000 -0.5000000000 3.5000000000 31.0000000000;
    3 -0.1458333333 -0.6666666667 1.0000000000 2.6666666667 -45.6666666667;
    4 -0.3307291667 -0.6250000000 1.3750000000 -1.2916666667 11.0000000000;
    5 -0.4455729167 -0.4666666667 0.8500000000 -3.1666666667 34.3333333333;
    6 -0.5041449653 -0.2569444444 -0.0125000000 -2.0902777778 -3.4444444444;
    7 -0.5183392237 -0.0404761905 -0.7464285714 0.3253968254 -30.9047619048;
    8 -0.4983629984 0.1539930556 -1.1087053571 2.2357390873 -16.3015873016;
    9 -0.4529195204 0.3097442681 -1.0611607143 2.6917438272 14.7918871252;
    10 -0.3893744141 0.4189459325 -0.7000223214 1.7562761795 27.9841269841;
    11 -0.3139072988 0.4801341791 -0.1807995130 0.1075436909 14.5369568703;
    12 -0.2316496389 0.4962122235 0.3403546063 -1.4486042948 -9.9037464593];

# Hermite - Obtained via the HTable function in utils.jl
const TriangleCoeffsHermite =[
    1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
    0 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
    -2 0 4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
    0 -12 0 8 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
    12 0 -48 0 16 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
    0 120 0 -160 0 32 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; # H5
    -120 0 720 0 -480 0 64 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
    0 -1680 0 3360 0 -1344 0 128 0 0 0 0 0 0 0 0 0 0 0 0 0;
    1680 0 -13440 0 13440 0 -3584 0 256 0 0 0 0 0 0 0 0 0 0 0 0;
    0 30240 0 -80640 0 48384 0 -9216 0 512 0 0 0 0 0 0 0 0 0 0 0;
    -30240 0 302400 0 -403200 0 161280 0 -23040 0 1024 0 0 0 0 0 0 0 0 0 0; #H10
    0 -665280 0 2217600 0 -1774080 0 506880 0 -56320 0 2048 0 0 0 0 0 0 0 0 0;
    665280 0 -7983360 0 13305600 0 -7096320 0 1520640 0 -135168 0 4096 0 0 0 0 0 0 0 0; # H12
    0 17297280 0 -69189120 0 69189120 0 -263577600 0 4392960 0 -319488 0 8192 0 0 0 0 0 0 0;
    -17297280 0 242161920 0 -484323840 0 322882560 0 -92252160 0 12300288 0 -745472 0 16384 0 0 0 0 0 0;
    0 -5189184000 0 2421619200 0 -2905943040 0 1383782400 0 -307507200 0 33546240 0 -1720320 0 32768 0 0 0 0 0;
    518918400 0 -8302694400 0 19372953600 0 -15498362880 0 5535129600 0 -984023040 0 89456640 0 -3932160 0 65536 0 0 0 0;
    0 17643225600 0 -94097203200 0 131736084480 0 -75277762560 0 20910489600 0 -3041525760 0 233963520 0 -8912896 0 131072 0 0 0;
    -17643225600 0 317578060800 0 -846874828800 0 790416506880 0 -338749931520 0 75277762560 0 -9124577280 0 601620480 0 -20054016 0 262144 0 0;
    0 -670442572800 0 4022655436800 0 -6436248698880 0 4290832465920 0 -1430277488640 0 260050452480 0 -26671841280 0 1524105216 0 -44826624 0 524288 0;
    670442572800 0 -13408851456000 0 40226554368000 0 -42908324659200 0 21454162329600 0 -5721109954560 0 866834841600 0 -76205260800 0 3810263040 0 -99614720 0 1048576]
H0(x) = 1.0
    H1(x) = 2x
    H2(x) = -2 + 4x^2
    H3(x) = -12x + 8x^3
    H4(x) = 12 - 48x^2 + 16x^4
    H5(x) = 120x -160x^3 + 32x^5
    H6(x) = -120 + 720x^2 - 480x^4 + 64x^6
    H7(x) = -1680x + 3360x^3 - 1344x^5 + 128x^7
    H8(x) =1680 - 13440x^2 + 13440x^4 - 3584x^6 + 256x^8
    H9(x) =30240x - 80640x^3 + 48384x^5 - 9216x^7 + 512x^9
    H10(x) =-30240 + 302400x^2 - 403200x^4 + 161280x^6 - 23040x^8 + 1024x^10
    H11(x) =-665280x + 2217600x^3 - 1774080x^5 + 506880x^7 - 56320x^9 + 2048x^11
    H12(x) =665280 - 7983360x^2 + 13305600x^4 - 7096320x^6 + 1520640x^8 - 135168x^10 + 4096x^12

    Hs = [H0, H1, H2, H3, H4, H5, H6, H7, H8, H9, H10, H11, H12]

    function hermite(x,n)
        if n == 0
            return 1
        elseif n == 1
            return 2x
        else
            return 2x*hermite(x,n-1) - 2*(n-1)*hermite(x,n-2)
        end
    end

# A and S p. 802
const testvalshermite = [0 0.5 1 3 5 10]
const TestTableHermite = [
    0 1.00000 1.00000 1.0000000 1.0000000000 1.0000000000;
    1 1.00000 2.00000 6.0000000 1.0000000000e1 2.0000000000e1
    2 -1.00000 2.00000 3.4000000e1 9.8000000000e1 3.9800000000e2;
    3 -5.00000 -4.00000 1.8000000e2 9.4000000000e2 7.8800000000e3;
    4 1.00000 -2.00000e1 8.7600000e2 8.8120000000e3 1.5521200000e5;
    5 4.10000e1 -8.00000 3.8160000e3 8.0600000000e4 3.0412000000e6;
    6 3.10000e1 1.84000e2 1.4136000e4 7.1788000000e5 5.9271880000e7;
    7 -4.61000e2 4.64000e2 3.9024000e4 6.2116000000e6 1.1489432000e9;
    8 -8.95000e2 -1.64800e3 3.6240000e4 5.2065680000e7 2.2149057680e10;
    9 6.48100e3 -1.07200e4 -4.0694400e5 4.2127120000e8 4.2459806240e11;
    10 2.25910e4 8.22400e3 -3.0939840e6 3.2755297600e9 8.0932782098e12;
    11 -1.07029e5 2.30848e5 -1.0425024e7 2.4329873600e10 1.5337360295e14;
    12 -6.04031e5 2.80768e5 5.5175040e6 1.7123708128e11 2.8894199383e15];

# Legendre - Table made with script from utils.jl
const TriangleCoeffsLegendre = [
    1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
    2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
    6 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
    20 12 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
    70 60 6 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
    252 280 60 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; # P5
    924 1260 420 20 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
    3432 5544 2520 280 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
    12870 24024 13860 2520 70 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
    48620 102960 72072 18480 1260 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
    184756 437580 360360 120120 13860 252 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; #P10
    705432 1847560 1750320 720720 120120 5544 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
    2704156 7759752 8314020 4084080 900900 72072 924 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
    10400600 32449872 38798760 22170720 6126120 720720 24024 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
    40116600 135207800 178474296 11639628038798760 6126120 360360 3432 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
    155117520 561632400 811246800 594914320 232792560 46558512 4084080 102960 0 0 0 0 0 0 0 0 0 0 0 0 0; # P15
    601080390 2326762800 3650610600 2974571600 1338557220 325909584 38798760 1750320 12870 0 0 0 0 0 0 0 0 0 0 0 0;
    23336062209617286240 16287339600 14602442400 7436429000 2141691552 325909584 22170720 437580 0 0 0 0 0 0 0 0 0 0 0 0 0 ;
    9075135300 39671305740 72129646800 70578471600 40156716600 13385572200 2498640144 232792560 8314020 48620 0 0 0 0 0 0 0 0 0 0 0;
    35345263800 163352435400 317370445920 336605018400 211735414800 80313433200 17847429600 2141691552116396280 1847560 0 0 0 0 0 0 0 0 0 0 0 0 ;
    137846528820 671560012200 1388495700900 1586852229600 1093966309800 465817912560 120470149800 17847429600 1338557220 38798760 184756 0 0 0 0 0 0 0 0 0 0]
P0(x) = 1
    P1(x) =  x
    P2(x) = .5*(3x^2 -1)
    P3(x) = .5*(5x^3 -3x)
    P4(x) = (1/8)*(35x^4 - 30x^2 + 3)
    P5(x) = (1/8)*(63x^5 - 70x^3 + 15x)
    P6(x) = (1/16)*(231x^6 - 315x^4 + 105x^2 - 5)
    P7(x) = (1/16)*(429x^7 - 693x^5 + 315x^3 - 35x)
    P8(x) = (1/128)*(6435x^8 - 12012x^6 + 6930x^4 - 1260x^2 + 35)
    P9(x) = (1/128)*(12155x^9 - 25740x^7 + 18018x^5 - 4620x^3 + 315x)
    P10(x) = (1/256)*(46189x^10 - 109395x^8 + 90090x^6- 30030x^4 + 3465x^2 - 63)
    P11(x) = (1/256)*(2^11)*(705432x^11 - 1847560x^9 + 1750320x^7 - 720720x^5 + 120120x^3 - 5544x)
    P12(x) = (1/512)*(2^12)*(2704156x^12 - 7759752x^10 + 8314020x^8 - 4084080x^6 + 900900x^4 - 72072x^2 + 924)

    Ps = [P1, P2, P3, P4, P5, P6, P7, P8, P9, P10, P11, P12]
    testPs = [P2, P3, P9, P10]

    function legendre(x, n)
        if n == 0
            return 1
        elseif n == 1
            return x
        else
            return ((2n-1)*x*legendre(x,n-1) - (n - 1)*legendre(x, n-2))/n
        end
    end

# Credit: Abramowitz and Stegun, Table 8.1, p. 342
# X arccos x P2(x) P3(x) P9(x) P10(x)
const TestTableLegendre =[
    0.00 90.0000000 -0.50000  0.0000000 0.00000000 -0.2460937;
    0.01 89.4270326 -0.49985 -0.0149975 0.02457330 -0.2447414;
    0.02 88.8540080 -0.49940 -0.0299800 0.04893045 -0.2406984;
    0.03 88.2808687 -0.49865 -0.0449325 0.07285701 -0.2340069;
    0.04 87.7075572 -0.49760 -0.0598400 0.09614188 -0.2247364;
    0.05 87.1340160 -0.49625 -0.0746875 0.11857899 -0.2129835;
    0.06 86.5601872 -0.49460 -0.0894600 0.13996890 -0.1988711;
    0.07 85.9860128 -0.49265 -0.1041425 0.16012040 -0.1825468;
    0.08 85.4114343 -0.49040 -0.1187200 0.17885206 -0.1641820;
    0.09 84.8363929 -0.48785 -0.1331775 0.19599366 -0.1439702;
    0.10 84.2608295 -0.48500 -0.1475000 0.21138764 -0.1221250;
    0.11 83.6846844 -0.48185 -0.1616725 0.22489042 -0.0988786;
    0.12 83.1078974 -0.47840 -0.1756800 0.23637363 -0.0744793;
    0.13 82.5304077 -0.47465 -0.1895075 0.24572526 -0.0491890;
    0.14 81.9521537 -0.47060 -0.2031400 0.25285070 -0.0232812;
    0.15 81.3730734 -0.46625 -0.2165625 0.25767367 0.0029618;
    0.16 80.7931038 -0.46160 -0.2297600 0.26013706 0.0292520;
    0.17 80.2121810 -0.45665 -0.2427175 0.26020358 0.0552981;
    0.18 79.6302402 -0.45140 -0.2554200 0.25785632 0.0808085;
    0.19 79.0472158 -0.44585 -0.2678525 0.25309918 0.1054942;
    0.20 78.4630410 -0.44000 -0.2800000 0.24595712 0.1290720;
    0.21 77.8776477 -0.43385 -0.2918475 0.23647631 0.1512674;
    0.22 77.2909670 -0.42740 -0.3033800 0.22472407 0.1718175;
    0.23 76.7029282 -0.42065 -0.3145825 0.21078870 0.1904736;
    0.24 76.1134596 -0.41360 -0.3254400 0.19477914 0.2070049;
    0.25 75.5224878 -0.40625 -0.3359375 0.17682442 0.2212002;
    0.26 74.9299379 -0.39860 -0.3460600 0.15707305 0.2328714;
    0.27 74.3357331 -0.39065 -0.3557925 0.13569215 0.2418552;
    0.28 73.7397953 -0.38240 -0.3651200 0.11286642 0.2480162;
    0.29 73.1420440 -0.37385 -0.3740275 0.08879707 0.2512481;
    0.30 72.5423969 -0.36500 -0.3825000 0.06370038 0.2514763;
    0.31 71.9407695 -0.35585 -0.3905225 0.03780634 0.2486591;
    0.32 71.3370751 -0.34640 -0.3980800 0.01135691 0.2427889;
    0.33 70.7312245 -0.33665 -0.4051575 -0.01539566 0.2338937;
    0.34 70.1231259 -0.32660 -0.4117400 -0.04219085 0.2220373;
    0.35 69.5126849 -0.31625 -0.4178125 -0.06876185 0.2073200;
    0.36 68.8998039 -0.30560 -0.4233600 -0.09483780 0.1898783;
    0.37 68.2843827 -0.29465 -0.4283675 -0.12014608 0.1698848;
    0.38 67.6663173 -0.28340 -0.4328200 -0.14441472 0.1475472;
    0.39 67.0455006 -0.27185 -0.4367025 -0.16737489 0.1231073;
    0.40 66.4218215 -0.26000 -0.4400000 -0.18876356 0.0968391;
    0.41 65.7951652 -0.24785 -0.4426975 -0.20832609 0.0690471;
    0.42 65.1654125 -0.23540 -0.4447800 -0.22581900 0.0400639;
    0.43 64.5324399 -0.22265 -0.4462325 -0.24101269 0.0102469;
    0.44 63.8961188 -0.20960 -0.4470400 -0.25369426 -0.0200245;
    0.45 63.2563161 -0.19625 -0.4471875 -0.26367022 -0.0503530;
    0.46 62.6128925 -0.18260 -0.4466600 -0.27076932 -0.0803272;
    0.47 61.9657035 -0.16865 -0.4454425 -0.27484521 -0.1095264;
    0.48 61.3145980 -0.15440 -0.4435200 -0.27577908 -0.1375251;
    0.49 60.6594184 -0.13985 -0.4408775 -0.27348225 -0.1638987;
    0.50 60.0000000 -0.12500 -0.4375000 -0.26789856 -0.1882286;
    0.51 59.3361703 -0.10985 -0.4333725 -0.25900667 -0.2101083;
    0.52 58.6677485 -0.09440 -0.4284800 -0.24682215 -0.2291492;
    0.53 57.9945451 -0.07865 -0.4228075 -0.23139939 -0.2449873;
    0.54 57.3163611 -0.06260 -0.4163400 -0.21283321 -0.2572892;
    0.55 56.6329870 -0.04625 -0.4090625 -0.19126025 -0.2657585;
    0.56 55.9442022 -0.02960 -0.4009600 -0.16686000 -0.2701428;
    0.57 55.2497742 -0.01265 -0.3920175 -0.13985552 -0.2702397;
    0.58 54.5494574  0.00460 -0.3822200 -0.11051366 -0.2659030;
    0.59 53.8429918 0.02215 -0.3715525 -0.07914497 -0.2570492;
    0.60 53.1301024 0.04000 -0.3600000 -0.04610304 -0.2436627;
    0.61 52.4104970 0.05815 -0.3475475 -0.01178332 -0.2258016;
    0.62 51.6838655 0.07660 -0.3341800 0.02337862 -0.2036019;
    0.63 50.9498775 0.09535 -0.3198825 0.05890951 -0.1772816;
    0.64 50.2081805 0.11440 -0.3046400 0.09430141 -0.1471441;
    0.65 49.4583981 0.13375 -0.2884375 0.12901554 -0.1135805;
    0.66 48.7001272 0.15340 -0.2712600 0.16248693 -0.0770701;
    0.67 47.9329352 0.17335 -0.2530925 0.19412981 -0.0381808;
    0.68 47.1563569 0.19360 -0.2339200 0.22334410 0.0024330;
    0.69 46.3698911 0.21415 -0.2137275 0.24952270 0.0440337;
    0.70 45.5729960 0.23500 -0.1925000 0.27205993 0.0858058;
    0.71 44.7650847 0.25615 -0.1702225 0.29036111 0.1268631;
    0.72 43.9455196 0.27760 -0.1468800 0.30385323 0.1662589;
    0.73 43.1136059 0.29935 -0.1224575 0.31199698 0.2029976;
    0.74 42.2685844 0.32140 -0.0969400 0.31430004 0.2360508;
    0.75 41.4096221 0.34375 -0.0703125 0.31033185 0.2643745;
    0.76 40.5358021 0.36640 -0.0425600 0.29973981 0.2869319;
    0.77 39.6461111 0.38935 -0.0136675 0.28226712 0.3027179;
    0.78 38.7394246 0.41260 0.0163800 0.25777224 0.3107893;
    0.79 37.8144885 0.43615 0.0475975 0.22625012 0.3102979;
    0.80 36.8698976 0.46000 0.0800000 0.18785528 0.3005298;
    0.81 35.9040686 0.48415 0.1136025 0.14292678 0.2809487;
    0.82 34.9152062 0.50860 0.1484200 0.09201529 0.2512452;
    0.83 33.9012620 0.53335 0.1844675 0.03591226 0.2113919;
    0.84 32.8598804 0.55840 0.2217600 -0.02431874 0.1617050;
    0.85 31.7883306 0.58375 0.2603125 -0.08730820 0.1029123;
    0.86 30.6834171 0.60940 0.3001400 -0.15134456 0.0362291;
    0.87 29.5413605 0.63535 0.3412575 -0.21433544 -0.0365586;
    0.88 28.3576366 0.66160 0.3836800 -0.27376627 -0.1130029;
    0.89 27.1267531 0.68815 0.4274225 -0.32665610 -0.1898929;
    0.90 25.8419328 0.71500 0.4725000 -0.36951049 -0.2631456;
    0.91 24.4946485 0.74215 0.5189275 -0.39827146 -0.3276858;
    0.92 23.0739181 0.76960 0.5667200 -0.40826421 -0.3773158;
    0.93 21.5651850 0.79735 0.6158925 -0.39414060 -0.4045743;
    0.94 19.9484436 0.82540 0.6664600 -0.34981919 -0.4005829;
    0.95 18.1948723 0.85375 0.7184375 -0.26842182 -0.3548803;
    0.96 16.2602047 0.88240 0.7718400 -0.14220642 -0.2552433; #this last 3 is a 4 in the book :(
    0.97 14.0698677 0.91135 0.8266825 0.03750397 -0.0874940;
    0.98 11.4783409 0.94060 0.8829800 0.28039609 0.1647081;
    0.99  8.1096144 0.97015 0.9407475 0.59724553 0.5200890;
    1.00  0.0000000 1.00000 1.0000000 1.00000000 1.0000000;]

# Chebyshev T First kind - Credit to Michael Somos, https://oeis.org/A049310
const TriangleCoeffsChebyshevT = float([
    1    0    0     0     0      0      0      0      0      0      0       0     0; # T0
    0    1    0     0     0      0      0      0      0      0      0       0     0;
   -1    0    2     0     0      0      0      0      0      0      0       0     0;
    0   -3    0     4     0      0      0      0      0      0      0       0     0;
    1    0   -8     0     8      0      0      0      0      0      0       0     0;
    0    5    0   -20     0     16      0      0      0      0      0       0     0; # T5
   -1    0   18     0   -48      0     32      0      0      0      0       0     0;
    0   -7    0    56     0   -112      0     64      0      0      0       0     0;
    1    0  -32     0   160      0   -256      0    128      0      0       0     0;
    0    9    0  -120     0    432      0   -576      0    256      0       0     0;
   -1    0   50     0  -400      0   1120      0  -1280      0    512       0     0; # T10
    0  -11    0   220     0  -1232      0   2816      0  -2816      0    1024     0;
    1    0  -72     0   840      0  -3584      0   6912      0  -6144       0  2048;])
T0(x) = 1
    T1(x) = x
    T2(x) = -1 + 2x^2
    T3(x) = -3x + 4x^3
    T4(x) = 1 - 8x^2 + 8x^4
    T5(x) = 5x - 20x^3 + 16x^5
    T6(x) = -1 + 18x^2 - 48x^4 + 32x^6
    T7(x) = -7x + 56x^3 - 112x^5 + 64x^7
    T8(x) = 1 - 32x^2 + 160x^4 - 256x^6 + 128x^8
    T9(x) = 9x - 120x^3 + 432x^5 - 576x^7 + 256x^9
    T10(x) = -1 + 50x^2 - 400x^4 + 1120x^6 - 1280x^8 + 512x^10
    T11(x) = -11x + 220x^3 - 1232x^5 + 2816x^7 - 2816x^9 + 1024x^11
    T12(x) = 1 - 72x^2 + 840x^4 - 3584x^6 + 6912x^8 - 6144x^10 + 2048x^12

    Ts = [T0 T1 T2 T3 T4 T5 T6 T7 T8 T9 T10 T11 T12]

    function chebyshevfirst(x, n)
        if n == 0
            1
        elseif n == 1
            x
        else
            2x*chebyshevfirst(x,n-1) - chebyshevfirst(x, n-2)
        end
    end

# Credit: Abramowitz and Stegun
const testvalschebyshev = [0.2 0.4 0.6 0.8 1.0]
const TestTableChebyshevT = float([
    0 1.0000000000 1.0000000000 1.0000000000 1.0000000000 1;
    1 0.2000000000 0.4000000000 0.6000000000 0.8000000000 1;
    2 -0.920000000 -0.680000000 -0.280000000 0.2800000000 1;
    3 -0.568000000 -0.944000000 -0.936000000 -0.352000000 1;
    4 0.6928000000 -0.075200000 -0.843200000 -0.843200000 1;
    5 0.8451200000 0.8838400000 -0.075840000 -0.997120000 1;
    6 -0.354752000 0.7822720000 0.7521920000 -0.752192000 1;
    7 -0.987020800 -0.258022400 0.9784704000 -0.206387200 1;
    8 -0.040056320 -0.988689920 0.4219724800 0.4219724800 1;
    9 0.9709982720 -0.532929536 -0.4721034240 0.881543168 1;
    10 0.4284556288 0.5623462912 -0.9884965888 0.9884965888 1;
    11 -0.7996160205 0.9828065690 -0.7140924826 0.7000513741 1;
    12 -0.7483020370 0.2238989640 0.1315856097 0.1315856097 1])

# Chebyshev U Second kind
const TriangleCoeffsChebyshevU = [
    1    0    0     0     0      0      0      0      0      0       0     0     0; # T0
    0    2    0     0     0      0      0      0      0      0       0     0     0;
    -1   0    4     0     0      0      0      0      0      0       0     0     0;
    0   -4    0     8     0      0      0      0      0      0       0     0     0;
    1    0  -12     0    16      0      0      0      0      0       0     0     0;
    0    6    0   -32     0     32      0      0      0      0       0     0     0; # T5
    -1   0   24     0   -80      0     64      0      0      0       0     0     0;
    0   -8    0    80     0   -192      0    128      0      0       0     0     0;
    1    0  -40     0   240      0   -448      0    256      0       0     0     0;
    0   10    0  -160     0    672      0  -1024      0    512       0     0     0;
    -1   0   60     0  -560      0   1792      0  -2304      0    1024     0     0; # T10
    0  -12    0   280     0  -1792      0   4608      0  -5120       0  2048     0;
    1    0  -84     0  1120      0  -5376      0  11520      0  -11264     0  4096;]
U0(x) = 1
    U1(x) = 2x
    U2(x) = 4x^2 - 1
    U3(x) = 8x^3 - 4x
    U4(x) = 16x^4 - 12x^2 + 1
    U5(x) = 32x^5 - 32x^3 + 6x
    U6(x) = 64x^6  - 80x^4 + 24x^2 - 1
    U7(x) = 128x^7 - 192x^5 + 80x^3 - 8x
    U8(x) = 256x^8 - 448x^6 + 240x^4 - 40x^2 + 1
    U9(x) = 512x^9 - 1024x^7 + 672x^5 - 160x^3 + 10x
    U10(x) = 1024x^10 - 2304x^8 + 1792x^6 - 560x^4 + 60x^2 - 1
    U11(x) = 2048x^11 - 5120x^9 + 4608x^7 - 1792x^5 + 280x^3 - 12x
    U12(x) = 4096x^12 - 11264x^10 + 11520x^8 - 5376x^6 + 1120x^4 - 84x^2 + 1

    Us = [ U0 U1 U2 U3 U4 U5 U6 U7 U8 U9 U10 U11 U12]

    function chebyshevsecond(x, n)
        if n == 0
            1
        elseif n == 1
            2x
        else
            2x*chebyshevsecond(x,n-1) - chebyshevsecond(x, n-2)
        end
    end

# Credit: Abramowitz and Stegun
const TestTableChebyshevU = [
    0 1.0000000000 1.0000000000 1.0000000000 1.0000000000 1;
    1 0.4000000000 0.8000000000 1.2000000000 1.6000000000 2;
    2 -0.840000000 -0.360000000 0.4400000000 1.5600000000 3;
    3 -0.736000000 -1.088000000 -0.672000000 0.8960000000 4;
    4 0.5456000000 -0.510400000 -1.246400000 -0.126400000 5;
    5 0.9542400000 0.6796800000 -0.823680000 -1.098240000 6;
    6 -0.163904000 1.0541440000 0.2579840000 -1.630784000 7;
    7 -1.019801600 0.1636352000 1.1332608000 -1.511014400 8;
    8 -0.244016640 -0.923235840 1.1019289600 -0.786839040 9;
    9 0.9221949440 -0.902223872 0.1890539520 0.2520719360 10;
    10 0.6128946176 0.2014567424 -0.8750642176 1.1901541376 11;
    11 -0.6770370970 1.0633892659 -1.2391310131 1.6521746842 12;
    12 -0.8837094564 0.6492546703 -0.6118929981 1.4533253571 13;]

# Gegenbauer
function poch(x,n)
    if n == 0
        return 1
    end
    res = x
    for i in 1:n-1
        res *= x+i
    end
    return res
end
TriangleCoeffsGegenbauer = []
Cα0(x, α) = 1
    Cα1(x, α) = 2α*x
    Cα2(x, α) = -α + 2*poch(α, 2)*x^2
    Cα3(x, α) = (1/3) * (-6poch(α, 2)x + 4poch(α, 3)x^3)
    Cα4(x, α) = (1/6) * (3poch(α, 2) - 12poch(α, 3)x^2 + 4poch(α, 4)x^4)
    Cα5(x, α) = (1/15) * (15poch(α, 3)x - 20poch(α, 4)x^3 + 4poch(α, 5)x^5)
    Cα6(x, α) = (1/90) * (-15poch(α, 3) + 90poch(α, 4)x^2 - 60poch(α, 5)x^4 + 8poch(α, 6)x^6)

    Cαs = [Cα0 Cα1 Cα2 Cα3 Cα4 Cα5 Cα6]

    # Abramowitz and Stegun, example 3, Table 22.2, p. 789
    # For C(.5,n+1)(2.5) == [5(n+.25)*C(.5,n)(2.5) - (n-.5)*C(.5, n-1)(2.5)]*(1/(n+1))

# Obtained from ApproxFun with Ultraspherical(λ == 2)
# using ApproxFun
# cs = [Fun(Ultraspherical(2), [zeros(i);1]) for i in 0:12]
# [cs[i](j) for i in 1:7, j in -1:.125:1]
const TestTableGegenbauer = [
      1.0    1.0       1.0       1.0       1.0   1.0       1.0       1.0       1.0   1.0       1.0       1.0       1.0   1.0       1.0       1.0       1.0;
     -4.0   -3.5      -3.0      -2.5      -2.0  -1.5      -1.0      -0.5       0.0   0.5       1.0       1.5       2.0   2.5       3.0       3.5       4.0;
     10.0    7.1875    4.75      2.6875    1.0  -0.3125   -1.25     -1.8125   -2.0  -1.8125   -1.25     -0.3125    1.0   2.6875    4.75      7.1875   10.0;
    -20.0  -10.9375   -4.5      -0.3125    2.0   2.8125    2.5       1.4375    0.0  -1.4375   -2.5      -2.8125   -2.0   0.3125    4.5      10.9375   20.0;
     35.0   13.1445    1.3125   -3.54297  -4.0  -2.16797   0.3125    2.26953   3.0   2.26953   0.3125   -2.16797  -4.0  -3.54297   1.3125   13.1445   35.0;
    -56.0  -12.291     3.9375    5.75195   2.0  -1.98633  -3.6875   -2.69336   0.0   2.69336   3.6875    1.98633  -2.0  -5.75195  -3.9375   12.291    56.0;
     84.0    7.56812  -8.64063  -3.66431   3.0   4.62866   1.73438  -2.24048  -4.0  -2.24048   1.73438   4.62866   3.0  -3.66431  -8.64063   7.56812  84.0;]

# Jacobi
# Table 22.1 of Abramowitz and Stegun, p. 793
Pαβ0(α , β, x) = 1
    Pαβ1(α , β, x) = .5 * (2*(α + 1) + (α + β + 2)*(x - 1))
    Pαβ2(α , β, x) = .125 * (4*poch(α + 1, 2) + 4*(α + β + 3)*poch(α + 2, 1)*(x - 1) + poch(α + β + 3, 2)*(x - 1)^2)
    Pαβ3(α , β, x) = (1/48) * (8*poch(α + 1, 3) + 12*(α + β + 4)*poch(α + 2, 2)*(x - 1) + 6*poch(α + β + 4, 2)*poch(α + 3, 1)*(x - 1)^2 + poch(α + β + 4, 3)*(x - 1)^3)
    Pαβ4(α , β, x) = (1/384) * (16*poch(α + 1, 4) + 32*(α + β + 5)*poch(α + 2, 3)*(x - 1) + 24*poch(α + β + 5, 2)*poch(α + 3, 2)*(x - 1)^2 + 8*poch(α + β + 5, 3)*poch(α + 4, 1)*(x - 1)^3 + poch(α + β + 5,4)*(x - 1)^4)
    Pαβ5(α , β, x) = (1/3840) * (32*poch(α + 1, 5) + 80*(α + β + 6)*poch(α + 2, 4)*(x - 1) + 80*poch(α + β + 6, 2)*poch(α + 3, 3)*(x - 1)^2 + 40*poch(α + β + 6, 3)*poch(α + 4, 2)*(x - 1)^3 + 10*poch(α + β + 6, 4)*poch(α + 5, 1)*(x - 1)^4 + poch(α + β + 6, 5)*(x-1)^5)
    Pαβ6(α , β, x) = (1/46080) * (64*poch(α + 1, 6) + 192*(α + β + 7)*poch(α + 2, 5)*(x - 1) + 240*poch(α + β + 7, 2)*poch(α + 3, 4)*(x - 1)^2 + 160*poch(α + β + 7, 3)*poch(α + 4, 3)*(x - 1)^3 + 60*poch(α + β + 7, 4)*poch(α + 5, 2)*(x - 1)^4 + 12*poch(α + β + 7, 5)*(α + 6)*(x - 1)^5 + poch(α + β + 7, 6)*(x-1)^6)

    Pαβs = [Pαβ0 Pαβ1 Pαβ2 Pαβ3 Pαβ4 Pαβ5 Pαβ6]

# Obtained from ApproxFun.jl for Jacobi(2 == β, 3 == α)
# Note that we do α THEN β whereas they do it the other way around.
# using ApproxFun
# ps = [Fun(Jacobi(2,3), [zeros(i); 1]) for i in 0:12]
# [ps[i](j) for i in 1:7, j in -1:.125:1]
const TestTableJacobi = [
      1.0   1.0       1.0        1.0         1.0        1.0        1.0         1.0        1.0       1.0        1.0        1.0        1.0        1.0        1.0        1.0       1.0;
     -3.0  -2.5625   -2.125     -1.6875     -1.25      -0.8125    -0.375       0.0625     0.5       0.9375     1.375      1.8125     2.25       2.6875     3.125      3.5625    4.0;
      6.0   4.14063   2.5625     1.26563     0.25      -0.484375  -0.9375     -1.10938   -1.0      -0.609375   0.0625     1.01563    2.25       3.76563    5.5625     7.64063  10.0;
    -10.0  -5.21362  -1.94336    0.0524902   1.01563    1.18774    0.810547    0.125732  -0.625    -1.19995   -1.35742   -0.855713   0.546875   3.09204    7.02148   12.5769   20.0;
     15.0   5.4113    0.487061  -1.31355    -1.26953   -0.398026   0.545654    1.06804    0.9375    0.18425   -0.899658  -1.76033   -1.58203    0.712814   6.46362   17.2717   35.0;
    -21.0  -4.59868   1.14267    1.60243     0.280273  -0.859185  -1.04269    -0.344714   0.65625   1.22132    0.794891  -0.651664  -2.27637   -2.02295    3.72281   20.4768   56.0;
     28.0   2.90217  -2.18814   -0.735861    0.943359   1.01195   -0.0296936  -0.952181  -0.875     0.18675    1.30673    1.16287   -0.916016  -3.3725    -0.311676  21.1388   84.0]

@testset "OP Table Cases" begin
    @testset "Laguerre" begin
        @testset "Hand Implemenation" begin
            @test all(f(0) ≈ 1.0 for f in Ls)
            @test all(Ls[i](testvalslaguerre[j]) ≈ TestTableLaguerre[i,j] for i in 1:13, j in 2:6)
            @test all([isapprox.(Ls[i](testvalslaguerre[j]),TestTableLaguerre[i,j]) for i in 1:13, j in 2:6])
        end
        @testset "Naive Implementation" begin
            @test all(laguerre(i,0.0) ≈ 1.0 for i in 0:12)
            @test all([laguerre(i-1,testvalslaguerre[j]) ≈ TestTableLaguerre[i,j] for i in 1:13, j in 2:6])
        end
        @testset "Stubborn Way" begin
            α = 0
            @test_broken all(eval.([a(0,i) for i in 1:12])  .≈ 1)
            @test_broken all(eval.([a(testvalslaguerre[j],i) for i in 0:12, j in 2:6]) .≈ TestTableLaguerre[:,2:6])
        end
    end

    @testset "Hermite" begin
        @testset "Hand Implemented" begin
            @test all(Hs[i](testvalshermite[j]) ≈ TestTableHermite[i,j] for i in 1:13, j in 2:6)
        end
        @testset "Naive Implementation" begin
            @test all(hermite(i,0.0) ≈ 1.0 for i in 1:12)
            @test all([hermite(testvalshermite[j], i-1) ≈ TestTableHermite[i,j] for i in 1:13, j in 2:6])
        end
        @testset "Stubborn Way" begin
            @test_broken 1 == 2
        end
    end

    @testset "Legendre" begin
        @testset "Hand Implemented" begin
            @test all([isapprox.(testPs[j-2](TestTableLegendre[i,1]),TestTableLegendre[i,j]; atol= 6.9e-8) for i in 1:101, j in 3:6]) # Do not the godawful atol keyword argument to shimmy all tests into passing.
        end
        @testset "Naive Recursion" begin
            @test all([legendre(i,0) for i in 1:10] .≈ 1)
            @test all([isapprox.(legendre(TestTableLegendre[i,1], 2), TestTableLegendre[i,3]; atol=5.e-6) for i in 1:101])
            @test all([isapprox.(legendre(TestTableLegendre[i,1], 3), TestTableLegendre[i,4]; atol=5.e-6) for i in 1:101])
            @test all([isapprox.(legendre(TestTableLegendre[i,1], 9), TestTableLegendre[i,5]; atol=5.e-6) for i in 1:101])
            @test all([isapprox.(legendre(TestTableLegendre[i,1], 10), TestTableLegendre[i,6]; atol=5.e-6) for i in 1:101])
        end
        @testset "Stubborn Way" begin
            @test_broken 1 ==2
        end
    end

    @testset "Chebyshev T First Kind" begin
        @testset "Hand Implemented" begin
            @test all(Ts[i](testvalschebyshev[j]) ≈ TestTableChebyshevT[i,j+1] for i in 1:13, j in 1:5)
        end
        @testset "Naive Recursion" begin
            @test all([isapprox.(chebyshevfirst(testvalschebyshev[j-1], i-1), TestTableChebyshevT[i,j]) for i in 1:13, j in 2:6])
        end
        @testset "Stubborn Way" begin
            @test_broken 1 == 2
        end
    end

    @testset "Chebyshev U Second Kind" begin
        @testset "Hand Implemented" begin
            @test all(Us[i](testvalschebyshev[j-1]) ≈ TestTableChebyshevU[i,j] for i in 1:13, j in 2:6)
        end
        @testset "Naive Recursion" begin
            @test_broken 1 == 2
        end
        @testset "Stubborn Way" begin
            @test_broken 1 == 2
        end
    end

    @testset "Gegenbauer/Ultraspherical" begin
        @testset "Hand Implementation" begin
            testcs() = [Cαs[i](j, 2) for i in 1:7, j in -1:.125:1]
            @test isapprox.(testcs(), TestTableGegenbauer; atol = 4e-5) |> all
        end
        # Abramowitz and Stegun, example 3, Table 22.2, p. 789
        @testset "Naive Recursion" begin
        end
        @testset "Stubborn Way" begin
            @test_broken 1 == 2
        end
    end

    @testset "Jacobi" begin
        @testset "Hand Implementation" begin
            test() = [Pαβs[i](3,2,j) for i in 1:7, j in -1:.125:1]
            @test isapprox.(test(), TestTableJacobi; atol = 5e-5) |> all
        end
        @testset "Naive Recursion" begin
            @test_broken 1==2
        end
        @testset "Stubborn Way" begin
            @test_broken 1 == 2
        end
    end

        j = Jacobi(.5,1.5,8)
            ge = Gegenbauer(3,4)
            go = Gegenbauer(4,5)
            ue = ChebyshevSecondKind(4)
            uo = ChebyshevSecondKind(5)
            te = ChebyshevFirstKind(4)
            to = ChebyshevFirstKind(5)
            le = Legendre(4)
            lo = Legendre(5)
            he = Hermite(4)
            ho = Hermite(5)
            lae = Laguerre(2,4)
            lao = Laguerre(3,5)

    @testset "Constructors and Parity" begin
        @test j == Jacobi{.5,1.5,8}()
        @test ge == GegenbauerEven{3,2}()
        @test go == GegenbauerOdd{4,2}()
        @test ue == ChebyshevSecondKindEven{2}()
        @test uo == ChebyshevSecondKindOdd{2}()
        @test te == ChebyshevFirstKindEven{2}()
        @test to == ChebyshevFirstKindOdd{2}()
        @test le == Legendre{4}()
        @test lo == Legendre{5}()
        @test he == HermiteEven{2}()
        @test ho == HermiteOdd{2}()
        @test lae == Laguerre{2,4}()
        @test lao == Laguerre{3,5}()
    end

    # Abramowitz pg. 790, Example 2.
    @testset "Jacobi{.5, 1.5, 8}(2) for a, b, c, d, f" begin
        @test_broken all([a(j,2,i) for i in 8:-1:0] .≈ [1, 1.132353, 1.366667, 1.841026, 3.008392, 6.849651, 26.44156, 223.1091, 6545.533])
        @test all([b(j,i) for i in 8:-1:0] .== [18.0, 34.0, 48.0, 60.0, 70.0, 78.0, 84.0, 88.0, 90.0])
        @test all([c(j,i) for i in 8:-1:0] .== [136.0, 105.0, 78.0, 55.0, 36.0, 21.0, 10.0, 3.0, 0.0])
        @test d(j,2) ≈ 3.338470458984375
        @test f(j,2) == -1
        @test a(j,2) ≈ 21852.07 atol=.01
    end

    # Abramowitz pg. 790, Example 3
    @testset "Gegenbauer{.5, 0:6}(2.5)" begin
        g0 = Gegenbauer(.5, 0)
        g1 = Gegenbauer(.5, 1)
        g2 = Gegenbauer(.5, 2)
        g3 = Gegenbauer(.5, 3)
        g4 = Gegenbauer(.5, 4)
        g5 = Gegenbauer(.5, 5)
        g6 = Gegenbauer(.5, 6)
        @test isapprox(a(g0,2.5), 1; atol=5e-6)
        @test_broken isapprox(a(g1,2.5), 1.25; atol=5e-6)
        @test_broken isapprox(a(g2,2.5) , 3.65625; atol=5e-6)
        @test_broken isapprox(a(g3,2.5) , 13.08594; atol=5e-6)
        @test_broken isapprox(a(g4,2.5) , 50.87648; atol=5e-6)
        @test_broken isapprox(a(g5,2.5) , 207.0649; atol=5e-6)
        @test_broken isapprox(a(g6,2.5), 867.7516; atol=5e-6)

    end

end


text:
l = Array{Float64}(undef, 20, 3)
for t in 1:20
    l[t,1] = (1. + 0.015)^t *2000
    l[t,2] = (1. + 0.02)^t * 1250
    l[t,3] = (1. + 0.03)^t * 1700
end

tfp = Array{Float64}(undef, 20, 3)
for t in 1:20
    tfp[t,1] = (1 + 0.06)^t * 3.2
    tfp[t,2] = (1 + 0.03)^t * 1.8
    tfp[t,3] = (1 + 0.05)^t * 2.5
end

s = Array{Float64}(undef, 20, 3)
for t in 1:20
    s[t,1] = 0.21
    s[t,2] = 0.15
    s[t,3] = 0.28
end

depk = [0.11, 0.135 ,0.15]
k0   = [50.5, 22., 33.5]

sigma = Array{Float64}(undef, 20, 3)
for t in 1:20
    sigma[t,1] = (1. - 0.05)^t * 0.58
    sigma[t,2] = (1. - 0.04)^t * 0.5
    sigma[t,3] = (1. - 0.045)^t * 0.6
end


text:
"""
    module Domains

for constructing domains to be fed to Simulation
"""
module Domains

# const SURFACE_BUFFER_FACTOR = 0/4 # abandoned for now

export Domain,
Cavity,
Resonator,
Void,
Dielectric,
Waveguide,
Lead

using ...Defaults
using ..Boundaries
using ..DielectricFunctions
using ..Lattices
using ..Shapes
using LinearAlgebra
using RecipesBase

"""
    Cavity(boundary,lattice,[dielectric,pump];align=false,[name]) -> domain
    Cavity(lattice,boundary,[dielectric,pump];align=false,[name]) -> domain
"""
Cavity(args...;kwargs...) = Domain(:Cavity,args...;kwargs...)
"""
    Resonator(boundary,lattice,[dielectric,pump];align=false,[name]) -> domain
    Resonator(lattice,boundary,[dielectric,pump];align=false,[name]) -> domain
"""
Resonator(args...;kwargs...) = Domain(:Resonator,args...;kwargs...)
"""
    Void(boundary,lattice;align=false,[name]) -> domain
    Void(lattice,boundary;align=false,[name]) -> domain
"""
Void(args...; align=false, kwargs...) = Domain(:Void,args...; align=align, kwargs...)
"""
    Dielectric(lattice,boundary,[ε=1,F=0];align=false,[name]) -> domain
    Dielectric(boundary,lattice,[ε=1,F=0];align=false,[name]) -> domain
"""
Dielectric(args...;kwargs...) = Domain(:Dielectric,args...;kwargs...)
Dielectric(bnd::Union{Boundary,Lattice},lat::Union{Boundary,Lattice},ε::Number,args...;kwargs...) = Domain(:Dielectric,bnd,lat,DielectricFunction(ε),args...;kwargs...)
Dielectric(bnd::Union{Boundary,Lattice},lat::Union{Boundary,Lattice},ε,F::Number;kwargs...) = Domain(:Dielectric,bnd,lat,ε,PumpFunction(F);kwargs...)
"""
    Resonator(boundary,lattice,[dielectric,pump];align=true,[name]) -> domain
    Resonator(lattice,boundary,[dielectric,pump];align=true,[name]) -> domain
"""
Waveguide(args...; kwargs...) = Domain(:Waveguide,args...; align=true, kwargs...)
"""
    Lead(boundary,lattice,[dielectric,pump];align=true,[name]) -> domain
    Lead(lattice,boundary,[dielectric,pump];align=true,[name]) -> domain
"""
Lead(args...; kwargs...) = Domain(:Lead,args...; align=true, kwargs...)

"""
    struct Domain

    Domain(type,boundary,lattice,[dielectric,pump];align=false,[name]) -> domain
    Domain(type,lattice,boundary,[dielectric,pump];align=false,[name]) -> domain

combines `boundary` with `lattice` (in either order) to generate a list of sites.
`T` just labels the kinds of domain (e.g. Cavity, or Waveguide)
each site is labeled by being in the interior, or containing the interior (`domain.surface`),
or being a corner.

----------------------
    (::Domain)(args...;align=false) -> dom

construct a new domain with modified parameters in `args`
This is for "updating" non-geometric fields of the immutable `domain`.
For example, to change the polarity of the boundary layers, do
    `new_dom = old_dom(conj(old_dom.boundary))`

For geometric parameters it recomputes the whole domain, so this command is no
more efficient than explicitly constructing a new domain. It might save some typing, however.
"""
struct Domain
    type::Symbol
    name::Symbol
    boundary::Boundary
    lattice::Lattice
    dielectric::DielectricFunction
    pump::PumpFunction

    x::Array{Float64,1}
    y::Array{Float64,1}
    r::Array{Float64,1}
    θ::Array{Float64,1}
    ij::Array{CartesianIndex{2},1}
    ε::Array{ComplexF64,1}
    F::Array{Float64,1}

    interior::BitArray{1}
    bulk::BitArray{1}
    surface::BitArray{1}
    corner::BitArray{1}

    nnxm::Array{Int,1}
    nnxp::Array{Int,1}
    nnym::Array{Int,1}
    nnyp::Array{Int,1}

    function Domain(type::Symbol,lat::Lattice,bnd::Boundary,dielectric::TDF=DielectricFunction(),pump::TPF=PumpFunction();kwargs...) where {TDF,TPF}
        return Domain(type,bnd,lat,dom.dielectric,dom.pump,dom.name;kwargs...)
    end
    function Domain(type::Symbol,
        boundary::Boundary,
        lattice::Lattice,
        dielectric::TDF=DielectricFunction(),
        pump::TPF=PumpFunction();
        name::Symbol=:anonymous,
        fit::Bool=false,
        align::Bool=fit
        ) where
        {TDF<:DielectricFunction,
        TPF<:PumpFunction}

        bnd = boundary
        lat = align ? lattice(θ=boundary.shape.θ,x0=boundary.shape.x0,y0=boundary.shape.y0) : lattice
        if fit
            if lat.type==:Cartesian
                if !(typeof(bnd.shape)<:AbstractQuadrilateral)
                    @warn "`fit` with Cartesian lattice only applies to Squares or Rectangles. ignoring."
                else
                    r = Rectangle(bnd.shape)
                    nx = ceil(Int,r.a/lat.dx)
                    ny = ceil(Int,r.b/lat.dy)
                    isodd(nx) ? nothing : nx += 1
                    isodd(ny) ? nothing : ny += 1
                    lat=lat(dx=r.a/nx,dy=r.b/ny)
                end
            elseif lat.type==:Polar
                if !(typeof(bnd.shape)<:Union{Circle,Annulus})
                    @warn "`fit` with Polar lattice only applies to Circle or Annulus. ignoring. "
                else
                    nr = ceil(Int,bnd.shape.R/lat.dr)
                    nθ = ceil(Int,2π/lat.dθ)
                    lat=lat(dr=bnd.shape.R/nr,dθ=2π/nθ)
                end
            else
                throw(LatticeError(lat))
            end
        end

        ijminmax = generate_spanning_indices(bnd,lat)
        x,y = generate_xy(bnd,lat,ijminmax...)
        interior = generate_interior(bnd,lat,x,y)
        t1,t2,t3,t4 = generate_translations(lat,ijminmax...)
        surface = generate_surface!(interior,lat,bnd,t1,t2,t3,t4,x,y)
        nnxm,nnxp,nnym,nnyp = generate_neighbor_indices(interior,lat)
        corner = generate_corner(x,y,bnd,lat,surface,interior)

        inds_keep = findall(interior)
        x = x[inds_keep]
        y = y[inds_keep]
        interior = interior[inds_keep]
        surface = surface[inds_keep]
        bulk = interior .& .!surface
        corner = corner[inds_keep]
        nnxm = nnxm[inds_keep]
        nnxp = nnxp[inds_keep]
        nnym = nnym[inds_keep]
        nnyp = nnyp[inds_keep]
        inds_keep = map(z->z+CartesianIndex(ijminmax[1]-1,ijminmax[3]-1),inds_keep)
        r = map(z->lat.r0 + z[1]*lat.dr, inds_keep)
        θ = map(z->lat.θ  + z[2]*lat.dθ, inds_keep)

        return new(
            type,name,bnd,lat,
            dielectric,pump,
            x,y,r,θ,inds_keep,dielectric.(x,y),pump.(x,y),
            interior,bulk,surface,corner,
            nnxm,nnxp,nnym,nnyp
            )
    end

    # for modifying non-geometric parameters of an already existing domain
    function (dom::Domain)(args...)
        for a ∈ args
            dom = dom(a)
        end
        return dom
    end
    function (dom::Domain)(dielectric::DielectricFunction)
        return new(dom.type,
        dom.name,dom.boundary,dom.lattice,dielectric,dom.pump,
        dom.x,dom.y,dom.r,dom.θ,dom.ij,dom.ε,dom.F,dom.interior,dom.bulk,
        dom.surface,dom.corner,dom.nnxm,dom.nnxp,dom.nnym,dom.nnyp)
    end
    function (dom::Domain)(pump::PumpFunction)
        return new(dom.type,
        dom.name,dom.boundary,dom.lattice,dom.dielectric,pump,
        dom.x,dom.y,dom.r,dom.θ,dom.ij,dom.ε,dom.F,dom.interior,dom.bulk,
        dom.surface,dom.corner,dom.nnxm,dom.nnxp,dom.nnym,dom.nnyp)
    end
    function (dom::Domain)(bnd::Boundary)
        return new(dom.type,
        dom.name,bnd,dom.lattice,dom.dielectric,dom.pump,
        dom.x,dom.y,dom.r,dom.θ,dom.ij,dom.ε,dom.F,dom.interior,dom.bulk,
        dom.surface,dom.corner,dom.nnxm,dom.nnxp,dom.nnym,dom.nnyp)
    end

    function Base.show(io::IO,dom::Domain)
        println(io,"Domain ",dom.type,": ",dom.name)
        println(io,"Number of sites: ",length(dom.x))
        println(io,"Number of bulk points: ",sum(dom.bulk))
        println(io,"Number of surface points: ",sum(dom.surface))
        println(io,"Number of corner points: ",sum(dom.corner))
        println(io)
        println(io,dom.dielectric)
        println(io)
        println(io,dom.pump)
        println(io)
        println(io,dom.boundary)
        println(io)
        print(io,dom.lattice)
    end
end

function (dom::Domain)(lat::Lattice;kwargs...)
    return Domain(dom.type,dom.boundary,lat,dom.dielectric,dom.pump,dom.name;kwargs...)
end

# cartesian spanning indices
function generate_spanning_indices(bnd::Boundary,lat::Lattice)
    if lat.type==:Cartesian
        return generate_spanning_indices_cartesian(bnd,lat)
    elseif lat.type==:Polar
        return generate_spanning_indices_polar(bnd,lat)
    else
        throw("unrecognized lattice $(lat.type)")
    end
end
function generate_spanning_indices_cartesian(bnd::Boundary, lat::Lattice)
    i,j = get_lattice_index(lat,bnd.shape.x0,bnd.shape.y0)
    imin, imax = floor(Int,i), ceil(Int,i)
    jmin, jmax = floor(Int,j), ceil(Int,j)

    flag_xmin = flag_xmax = flag_ymin = flag_ymax = true
    init_flag = false
    while flag_ymin || flag_ymax || flag_xmin|| flag_xmax
        I = imin:imax
        J = jmin:jmax
        XYBool = BitArray(undef,2length(I)+2length(J))
        count1 = 0
        count2 = 0
        temp = Array{Float64}(undef,2)
        for i ∈ I
            count2 += 1
            X, Y = lat(i,J[1])
            XYBool[1+count1+i-I[1]] = bnd.shape(X,Y)
        end
        !init_flag || any(XYBool[(count1+1):count2]) ? (jmin -= 1; flag_ymin=true) : flag_ymin=false
        !init_flag ? init_flag = any(XYBool[(count1+1):count2]) : nothing
        count1 = count2
        for i ∈ I
            count2 += 1
            X,Y = lat(i,J[end])
            XYBool[1+count1+i-I[1]] = bnd.shape(X,Y)
        end
        !init_flag || any(XYBool[(count1+1):count2]) ? (jmax += 1; flag_ymax=true) : flag_ymax=false
        !init_flag ? init_flag = any(XYBool[(count1+1):count2]) : nothing
        count1 = count2
        for j ∈ J
            count2 += 1
            X,Y = lat(I[1],j)
            XYBool[1+count1+j-J[1]] = bnd.shape(X,Y)
        end
        !init_flag || any(XYBool[(count1+1):count2]) ? (imin -= 1; flag_xmin=true) : flag_xmin=false
        !init_flag ? init_flag = any(XYBool[(count1+1):count2]) : nothing
        count1 = count2
        for j ∈ J
            count2 += 1
            X,Y = lat(I[end],j)
            XYBool[2length(I) + length(J) + j+1-J[1]] = bnd.shape(X,Y)
        end
        !init_flag || any(XYBool[(count1+1):count2]) ? (imax += 1; flag_xmax=true) : flag_xmax=false
        !init_flag ? init_flag = any(XYBool[(count1+1):count2]) : nothing
    end
    return imin-2, imax+2, jmin-2, jmax+2
end

# polar spanning indices
function generate_spanning_indices_polar(bnd::Boundary, lat::Lattice)
    @assert typeof(bnd.shape)<:Union{Circle,Annulus} "polar coordinates only for circular outermost region"
    @assert lat.x0==bnd.shape.x0 && lat.y0==bnd.shape.y0 "polar lattice origin $((lat.x0,lat.y0)) must be same as circle origin $((bnd.shape.x0,bnd.shape.y0))"
    # i,j = get_lattice_index(lat,bnd.shape.x0,bnd.shape.y0)
    # imin, imax = floor(Int,i), ceil(Int,i)
    # jmin, jmax = floor(Int,j), ceil(Int,j)
    # imin ≤ 0 ? imin = 0 : nothing
    # jmin ≤ 0 ? jmin = 0 : nothing
    imin = 0
    imax = 1
    jmin = 0
    jmax = floor(Int,2π/lat.dθ)

    flag_rmin = flag_rmax = flag_θmin = flag_θmax = true
    init_flag = false
    while flag_rmin || flag_rmax || flag_θmin|| flag_θmax
        I = imin:imax
        J = jmin:jmax
        XYBool = BitArray(undef,2length(I)+2length(J))
        count1 = 0
        count2 = 0
        temp = Array{Float64}(undef,2)
        for i ∈ I
            count2 += 1
            X, Y = lat(i,J[1])
            XYBool[1+count1+i-I[1]] = bnd.shape(X,Y)
        end
        (!init_flag || any(XYBool[(count1+1):count2])) && jmin≥1 ? (jmin -= 1; flag_θmin=true) : flag_θmin=false
        !init_flag ? init_flag = any(XYBool[(count1+1):count2]) : nothing
        count1 = count2
        for i ∈ I
            count2 += 1
            X,Y = lat(i,J[end])
            XYBool[1+count1+i-I[1]] = bnd.shape(X,Y)
        end
        (!init_flag || any(XYBool[(count1+1):count2])) && (jmax+1)*lat.dθ≤2π ? (jmax += 1; flag_θmax=true) : flag_θmax=false
        !init_flag ? init_flag = any(XYBool[(count1+1):count2]) : nothing
        count1 = count2
        for j ∈ J
            count2 += 1
            X,Y = lat(I[1],j)
            XYBool[1+count1+j-J[1]] = bnd.shape(X,Y)
        end
        (!init_flag || any(XYBool[(count1+1):count2])) && imin≥1 ? (imin -= 1; flag_rmin=true) : flag_rmin=false
        !init_flag ? init_flag = any(XYBool[(count1+1):count2]) : nothing
        count1 = count2
        for j ∈ J
            count2 += 1
            X,Y = lat(I[end],j)
            XYBool[2length(I) + length(J) + j+1-J[1]] = bnd.shape(X,Y)
        end
        !init_flag || any(XYBool[(count1+1):count2]) ? (imax += 1; flag_rmax=true) : flag_rmax=false
        !init_flag ? init_flag = any(XYBool[(count1+1):count2]) : nothing
    end
    return imin, imax, jmin, jmax
end

function generate_xy(bnd::Boundary,lat::Lattice,imin,imax,jmin,jmax)
    I, J = imin:imax, jmin:jmax
    X = Array{Float64}(undef,length(I),length(J))
    Y = Array{Float64}(undef,length(I),length(J))
    for k ∈ CartesianIndices(X)
        i,j = k[1],k[2]
        X[k],Y[k] = lat(I[i],J[j])
    end
    return X, Y
end

function generate_interior(bnd::Boundary,lat::Lattice,x::Array,y::Array)
    interior = BitArray(undef,size(x)...)
    for k ∈ CartesianIndices(x)
        interior[k] = bnd.shape(x[k],y[k])
    end
    return interior
end

function generate_translations(lat::Lattice,imin,imax,jmin,jmax)
    I, J = imin:imax, jmin:jmax
    x1 = Array{Float64}(undef,length(I),length(J))
    x2 = Array{Float64}(undef,length(I),length(J))
    x3 = Array{Float64}(undef,length(I),length(J))
    x4 = Array{Float64}(undef,length(I),length(J))
    y1 = Array{Float64}(undef,length(I),length(J))
    y2 = Array{Float64}(undef,length(I),length(J))
    y3 = Array{Float64}(undef,length(I),length(J))
    y4 = Array{Float64}(undef,length(I),length(J))
    for k ∈ CartesianIndices(x1)
        i,j = k[1],k[2]
        x1[k],y1[k] = lat(I[i]-1,J[j])
        x2[k],y2[k] = lat(I[i]+1,J[j])
        x3[k],y3[k] = lat(I[i],J[j]-1)
        x4[k],y4[k] = lat(I[i],J[j]+1)
    end
    return (x1,y1),(x2,y2),(x3,y3),(x4,y4)
end

function generate_surface!(interior::BitArray,lat::Lattice,bnd::Boundary,t1,t2,t3,t4,x,y)
    surface = falses(size(interior)...)
    surface_temp = BitArray(undef,size(interior)...)
    for k ∈ CartesianIndices(interior)
        surface_temp[k] = interior[k] &&
            (!bnd.shape(t1[1][k],t1[2][k]) ||
             !bnd.shape(t2[1][k],t2[2][k]) ||
             !bnd.shape(t3[1][k],t3[2][k]) ||
             !bnd.shape(t4[1][k],t4[2][k])
            )
        if surface_temp[k]
            # n,t,d = normal_distance(bnd.shape,x[k],y[k])
            # if d[1]<min(lat.dx*SURFACE_BUFFER_FACTOR,lat.dy*SURFACE_BUFFER_FACTOR)
                # interior[k] = false
                # surface[k] = false
                # surface[k+CartesianIndex( 1, 0)]=true
                # surface[k+CartesianIndex(-1, 0)]=true
                # surface[k+CartesianIndex( 0, 1)]=true
                # surface[k+CartesianIndex( 0,-1)]=true
            # else
                surface[k] = true
            # end
        end
    end
    return surface
end

function generate_neighbor_indices(interior,lat::Lattice)
    if lat.type==:Cartesian
        return generate_neighbor_indices_cartesian(interior)
    elseif lat.type==:Polar
        return generate_neighbor_indices_polar(interior)
    else
        throw(LatticeError(lat))
    end
end

function generate_neighbor_indices_cartesian(interior)
    NN = LinearIndices(interior) - reshape(cumsum(.!interior[:]),size(interior)...)
    nnxm = vcat(zeros(Int,1,size(interior,2)),NN[1:end-1,:])
    nnxp = vcat(NN[2:end,:],zeros(Int,1,size(interior,2)))
    nnym = hcat(zeros(Int,size(interior,1),1),NN[:,1:end-1])
    nnyp = hcat(NN[:,2:end],zeros(Int,size(interior,1),1))
    return nnxm,nnxp,nnym,nnyp
end

function generate_neighbor_indices_polar(interior)
    NN = LinearIndices(interior) - reshape(cumsum(.!interior[:]),size(interior)...)
    nnrm = vcat(zeros(Int,1,size(interior,2)),NN[1:end-1,:])
    nnrp = vcat(NN[2:end,:],zeros(Int,1,size(interior,2)))
    nnθm = hcat(NN[:,end],NN[:,1:end-1])
    nnθp = hcat(NN[:,2:end],NN[:,1])
    return nnrm,nnrp,nnθm,nnθp
end

function generate_corner(x,y,bnd::Boundary,lat::Lattice,surface,interior)
    if lat.type==:Cartesian
        return generate_corner_cartesian(x,y,bnd,lat,surface,interior)
    elseif lat.type==:Polar
        return generate_corner_polar(x,y,bnd,lat,surface,interior)
    else
        throw(LatticeError(lat))
    end
end

function generate_corner_cartesian(x,y,bnd::Boundary,lat::Lattice,surface,interior)
    corner = falses(size(surface)...)
    c = bnd.shape.corners
    dx = hypot(lat.dx,lat.dy)#*(1+SURFACE_BUFFER_FACTOR)
    if length(c)>0
        for k ∈ CartesianIndices(surface)
            if surface[k]
                d = map(z->hypot(x[k]-z[1],y[k]-z[2]),c)
                corner[k] = minimum(d) ≤ dx+eps(dx)
                X1, Y1 = x[k+CartesianIndex( 2, 0)], y[k+CartesianIndex( 2, 0)]
                X2, Y2 = x[k+CartesianIndex(-2, 0)], y[k+CartesianIndex(-2, 0)]
                X3, Y3 = x[k+CartesianIndex( 0, 2)], y[k+CartesianIndex( 0, 2)]
                X4, Y4 = x[k+CartesianIndex( 0,-2)], y[k+CartesianIndex( 0,-2)]
                count = sum(
                            (!bnd.shape(X1,Y1) && !interior[k+CartesianIndex( 1, 0)],
                             !bnd.shape(X2,Y2) && !interior[k+CartesianIndex(-1, 0)],
                             !bnd.shape(X3,Y3) && !interior[k+CartesianIndex( 0, 1)],
                             !bnd.shape(X4,Y4) && !interior[k+CartesianIndex( 0,-1)],)
                            )
                corner[k] = corner[k] && (count>1)
            end
        end
    end
    return corner
end

function generate_corner_polar(x,y,bnd,lat,surface,interior)
    return falses(size(surface)...)
end

@recipe function f(d::Domain)
    aspect_ratio --> 1
    legend --> false
    @series begin
        seriestype --> :scatter
        markersize --> MARKERSIZE_SCALE/sqrt(length(d.x))
        markerstrokealpha --> 0
        shape --> MARKERSHAPE
        (d.x[d.bulk],d.y[d.bulk])
    end
    @series begin
        seriestype --> :scatter
        markersize --> MARKERSIZE_SCALE/sqrt(length(d.x))
        markerstrokealpha --> 0
        shape --> MARKERSHAPE
        (d.x[d.surface],d.y[d.surface])
    end
    @series begin
        seriestype --> :scatter
        markersize --> MARKERSIZE_SCALE/sqrt(length(d.x))
        markerstrokealpha --> 0
        shape --> MARKERSHAPE
        (d.x[d.corner],d.y[d.corner])
    end
    @series d.boundary
end

end # module