text:
# ---
# title: 662. Maximum Width of Binary Tree
# id: problem662
# author: Indigo
# date: 2021-06-28
# difficulty: Medium
# categories: Tree
# link: <https://leetcode.com/problems/maximum-width-of-binary-tree/description/>
# hidden: true
# ---
#
# Given a binary tree, write a function to get the maximum width of the given
# tree. The maximum width of a tree is the maximum width among all levels.
#
# The width of one level is defined as the length between the end-nodes (the
# leftmost and right most non-null nodes in the level, where the `null` nodes
# between the end-nodes are also counted into the length calculation.
#
# It is **guaranteed** that the answer will in the range of 32-bit signed
# integer.
#
# **Example 1:**
#
#
#
#     Input:
#
#                1
#              /   \
#             3     2
#            / \     \
#           5   3     9
#
#     Output: 4
#     Explanation: The maximum width existing in the third level with the length 4 (5,3,null,9).
#
#
# **Example 2:**
#
#
#
#     Input:
#
#               1
#              /
#             3
#            / \
#           5   3
#
#     Output: 2
#     Explanation: The maximum width existing in the third level with the length 2 (5,3).
#
#
# **Example 3:**
#
#
#
#     Input:
#
#               1
#              / \
#             3   2
#            /
#           5
#
#     Output: 2
#     Explanation: The maximum width existing in the second level with the length 2 (3,2).
#
#
# **Example 4:**
#
#
#
#     Input:
#
#               1
#              / \
#             3   2
#            /     \
#           5       9
#          /         \
#         6           7
#     Output: 8
#     Explanation: The maximum width existing in the fourth level with the length 8 (6,null,null,null,null,null,null,7).
#
#
#
#
# **Constraints:**
#
#   * The given binary tree will have between `1` and `3000` nodes.
#
#
## @lc code=start
using LeetCode

width_of_binary_tree(root::Nothing) = 0

function width_of_binary_tree(root::TreeNode{Int})::Int
    q1 = Tuple{TreeNode{Int}, Int}[]
    q2 = copy(q1)
    push!(q1, (root, 1))
    res = 1
    while !isempty(q1)
        res = max(res, q1[end][2] - q1[1][2] + 1)
        for (nd, idx) in q1
            !isnothing(nd.left) && push!(q2, (nd.left, idx << 1))
            !isnothing(nd.right) && push!(q2, (nd.right, (idx << 1) + 1))
        end
        q1, q2 = q2, q1
        empty!(q2)
    end
    return res
end
## @lc code=end


text:
using PyCall
@pyimport matplotlib.colors as matcolors
@pyimport matplotlib.colorbar as matcolorbar

module PlotFunctions
	import  ..MasslistFunctions
	import  ..DatasetFunctions
	using PyPlot
	using CSV
	using DataFrames
	using Dates

	export massDefectPlot, bananaPlot, addClickToggle, plotStages

	function massDefectPlot(masses, compositions, concentrations, colors, plotTitle, colorCodeTitle, ax; dotSize = 10, dotbase = 2, maxMass = 450, maxDefect = 0.25, maxDIVmin = 100, sumformulas = false, negativeMD = true, minConc = 0, maxConc = 0, colormap = "haline", norm = 0)
	  if (minConc == 0) & (maxConc == 0)
	  	maxConc = DatasetFunctions.nanmax(concentrations)
	  	minConc = maxConc/maxDIVmin
	  else
		maxDIVmin = maxConc/minConc
          end

	  #figure()
	  h2o = MasslistFunctions.createCompound(H=2,O=1, Hplus=0)
	  o = MasslistFunctions.createCompound(O=1, Hplus=0)

	  for i=1:length(masses)
	    m = masses[i]
	    adduct = h2o
	    if (MasslistFunctions.inCompositions(compositions[:,i] + adduct[4], compositions))
	      m1 = m + adduct[1]
	      #ax.plot([m, m1], [m-round(m), m1 - round(m1)], color="lightblue", zorder=1)
	    end
	    adduct = o
	    if (MasslistFunctions.inCompositions(compositions[:,i] + adduct[4], compositions))
	      m1 = m + adduct[1]
	      #ax.plot([m, m1], [m-round(m), m1 - round(m1)], color="red", zorder=1)
	    end
	    if sumformulas text(masses[i],masses[i]-round(masses[i]),sumFormulaStringFromCompositionArray(compositions[:,i]), color="grey", clip_on=true, verticalalignment="center", size=10, zorder=100) end
	  end
	  if norm != 0
		  if negativeMD == true
		    h = ax.scatter(masses, masses-round.(masses),dotSize*log.(concentrations./minConc), colors, zorder=10, linewidths=0.5, cmap=PyPlot.cm[colormap], norm=norm, alpha=0.9)
		    ax.set_ylim(-maxDefect,maxDefect)
		  else
		    h = ax.scatter(masses, masses-floor.(masses),dotSize*log.(concentrations./minConc), colors, zorder=10, linewidths=0.5, cmap=PyPlot.cm[colormap], norm=norm, alpha=0.9)
		    ax.set_ylim(0,maxDefect)
		  end
	  else
		  if negativeMD == true
		    h = ax.scatter(masses, masses-round.(masses),dotSize*log.(concentrations./minConc), colors, zorder=10, linewidths=0.5, cmap=PyPlot.cm[colormap], alpha=0.9)
		    ax.set_ylim(-maxDefect,maxDefect)
		  else
		    h = ax.scatter(masses, masses-floor.(masses),dotSize*log.(concentrations./minConc), colors, zorder=10, linewidths=0.5, cmap=PyPlot.cm[colormap], alpha=0.9)
		    ax.set_ylim(0,maxDefect)
		  end
	  end
	  #cb=colorbar(h, ax = ax)
	  #cb["ax"]["set_ylabel"](colorCodeTitle)
	  #ax.set_xlabel("Mass [amu]")
	  #ax.set_ylabel("Kendrick Mass Defect")
	  #ax.set_title(plotTitle)
	  #ax.grid("on")

	  maxsize = dotSize*log(maxDIVmin)
  	  msizes = dotSize.*log.(dotbase .^ (range(1e-5*maxDIVmin,stop=log(dotbase, maxDIVmin)+maxDIVmin*1e-5)))
  	  mvalues = round.(minConc .* (dotbase .^ (range(1e-5*maxDIVmin,stop=log(dotbase,maxDIVmin)+maxDIVmin*1e-5))), sigdigits = 2)
  	  for i=1:size(msizes)[1]
  	    ax.scatter([],[],s=msizes[i], color = "k", label=mvalues[i])
  	  end
  	  ax.legend()
	  return h
	end

function massDefectPlotVolatility(masses, compositions, concentrations, plotTitle, ax; dotSize = 10, dotbase = 2, maxMass = 450, maxDefect = 0.25, maxDIVmin = 100, sumformulas = false, negativeMD = true, minConc = 0, maxConc = 0, cmap = cmap, norm = norm)

	# calculate log10Cstar300K
	Nc = 25.0 	#basically const. which c-nr to start from
	bc = 0.475	#
	bo = 2.3
	bco = -0.3
	nc = compositions[1,:]
	no = compositions[6,:]

	badd = Float64.(copy(nc))
	for i = 1:length(nc)
	  if 5<=nc[i]<=15
	      badd[i] = 0.9
	  elseif nc[i] > 15
	      badd[i] = 1.13
	  end
	end

	log10Cstar300K = ((Nc.-nc).*bc .- no.*(bo.-badd) .- 2 .* bco .*(nc.*no)./(nc.+no)) #log10(ug/m³)  nonlinear approx of sat gas concentration at 300K

	if (minConc == 0) & (maxConc == 0)
		maxConc = DatasetFunctions.nanmax(concentrations)
		minConc = maxConc/maxDIVmin
	else
		maxDIVmin = maxConc/minConc
	end

	  #figure()
	h2o = MasslistFunctions.createCompound(H=2,O=1, Hplus=0)
	o = MasslistFunctions.createCompound(O=1, Hplus=0)

	  for i=1:length(masses)
	    m = masses[i]
	    adduct = h2o
	    if (MasslistFunctions.inCompositions(compositions[:,i] + adduct[4], compositions))
	      m1 = m + adduct[1]
	      #ax.plot([m, m1], [m-round(m), m1 - round(m1)], color="lightblue", zorder=1)
	    end
	    adduct = o
	    if (MasslistFunctions.inCompositions(compositions[:,i] + adduct[4], compositions))
	      m1 = m + adduct[1]
	      #ax.plot([m, m1], [m-round(m), m1 - round(m1)], color="red", zorder=1)
	    end
	    if sumformulas text(masses[i],masses[i]-round(masses[i]),sumFormulaStringFromCompositionArray(compositions[:,i]), color="grey", clip_on=true, verticalalignment="center", size=10, zorder=100) end
	  end
	  if negativeMD == true
	    h = ax.scatter(masses, masses-round.(masses),dotSize*log.(concentrations./minConc), log10Cstar300K, zorder=10, linewidths=0.5, cmap=cmap, norm=norm, alpha=0.9)
	    ax.set_ylim(-maxDefect,maxDefect)
	  else
	    h = ax.scatter(masses, masses-floor.(masses),dotSize*log.(concentrations./minConc), log10Cstar300K, zorder=10, linewidths=0.5, cmap=cmap, norm=norm, alpha=0.9)
	    ax.set_ylim(0,maxDefect)
	  end
	  #cb=colorbar(h, ax = ax)
	  #cb["ax"]["set_ylabel"](colorCodeTitle)
	  #ax.set_xlabel("Mass [amu]")
	  #ax.set_ylabel("Kendrick Mass Defect")
	  #ax.set_title(plotTitle)
	  #ax.grid("on")

	  maxsize = dotSize*log(maxDIVmin)
  	  msizes = dotSize.*log.(dotbase .^ (range(1e-5*maxDIVmin,stop=log(dotbase, maxDIVmin)+maxDIVmin*1e-5)))
  	  mvalues = round.(minConc .* (dotbase .^ (range(1e-5*maxDIVmin,stop=log(dotbase,maxDIVmin)+maxDIVmin*1e-5))), sigdigits = 2)
  	  for i=1:size(msizes)[1]
  	    ax.scatter([],[],s=msizes[i], color = "k", label=mvalues[i])
  	  end
  	  ax.legend()
	  return h
	end

	function bananaPlot(xbins,ybins,meshdataXY;subplotAx=0)
	    println("Not implemented yet!!!")
	    #if subplotAx == 0
	    #    figure()
	    #    subplotAx = subplot(1,1,1)
	    #end
	    pcolormesh([DateTime(2018,05,05),DateTime(2018,05,06),DateTime(2018,05,07)],[1,7,9],[1 2 3; 4 5 6; 7 8 9][1:end-1,1:end-1])
	    #yscale("log")
	end

	function addClickToggle(plotAxis)
	    f = plotAxis.get_figure()
	    plotLines = plotAxis.get_lines()
	    legLines = plotAxis.get_legend().get_lines()
	    for ll in legLines
	        ll.set_picker(5)
	    end
	    linesDict = Dict(zip(legLines,plotLines))
	    function onLegendPick(event)
	        legLine = event.artist
	        plotLine = linesDict[legLine]
	        isVisible = plotLine.get_visible()
	        plotLine.set_visible(! isVisible)
	        if ! isVisible
	            legLine.set_alpha(1.0)
	        else
	            legLine.set_alpha(0.1)
	        end
	        f.canvas.draw()
	    end
	    f.canvas.mpl_connect("pick_event",onLegendPick)
	end
	function getVisibleTraces(axis)
		return [ axis.get_lines()[i].get_visible() for i in 1:length(axis.get_lines())]
	end

	function plotStages(file, axes, starttime, endtime; headerrow = 1)
		data = DataFrame(CSV.File(file, header = headerrow))
		if ("year" in names(data)) && ("month" in names(data)) && ("day" in names(data)) && ("hour" in names(data)) && ("minute" in names(data))
			datetime = Dates.DateTime.(data.year, data.month, data.day, data.hour, data.minute)
		elseif ("datetime" in names(data))
			println(data.datetime[1])
			println("What is the timestringformat of this datetime stamp?")
			formatString = readline()
			datetime = Dates.DateTime.(data.datetime, formatString)
		elseif ("unixtime" in names(data))
			datetime = Dates.unix2datetime.(data.unixtime)
		end
		show = starttime .< datetime .< endtime
		vlines.(datetime[show], 0, 1)
		axvline.(datetime[show])
		text.(datetime[show]+Dates.Minute(5), axes.get_ylim()[2]*0.75 , string.(data.comments[show], "   "), rotation=90)
	end
end

module CalibrationFunctions
	import ..InterpolationFunctions
	import PyPlot

	export generateCalibFactorTrace, generateBgTraces, interpolateBgTraces
	function getMeanOfQuantile(samples,quant)
	    q = quantile(samples,quant)
	    return mean(samples[samples.<q])
	end
	function generateBgTraces(times, traces; slices=10, quant=0.05)
	    dt = false
	    if typeof(times[1]) == DateTime
		times = Dates.datetime2unix.(times)
		dt = true
	    end
	    l=size(traces,1)
	    w=size(traces,2)
	    m=slices
	    bgtraces = Array{Float64}(size(traces,1), size(traces,2))
	    bgtimes = [mean(times[Int(floor((n-1)*l/m+1)):Int(floor(n*l/m)),1]) for n = 1:m]
	    for i=1:w
		bgtraces[:,i] = InterpolationFunctions.interpolate(times, bgtimes, [getMeanOfQuantile(traces[Int(floor((n-1)*l/m+1)):Int(floor(n*l/m)),i],0.08) for n = 1:m])
	    end
	    if dt
		bgtimes = Dates.unix2datetime.(bgtimes)
	    end
	    return bgtraces
	end

	function interpolateBgTraces(times, bgtimes, bgvalues)
	    if typeof(times[1]) == DateTime
		times = Dates.datetime2unix.(times)
	    end
	    if typeof(bgtimes[1]) == DateTime
		bgtimes = Dates.datetime2unix.(bgtimes)
	    end

	    nMasses = size(bgvalues,2)
	    bgtraces = Array{Float64}(size(times,1), nMasses)

	    for i=1:nMasses
		bgtraces[:,i] = InterpolationFunctions.interpolate(times, bgtimes, bgvalues[:,i])
	    end

	    return bgtraces
	end

	function generateCalibFactorTrace(traceTimes, calibTimes, calibValues, transitionTimes)
	    if typeof(traceTimes[1]) == DateTime
		println("Convertung traceTimes to unixtime")
		traceTimes = Dates.datetime2unix(traceTimes)
	    end
	    if typeof(calibTimes[1]) == DateTime
		println("Convertung calibTimes to unixtime")
		calibTimes = Dates.datetime2unix(calibTimes)
	    end
	    if typeof(transitionTimes[1]) == DateTime
		println("Convertung transitionTimes to unixtime")
		transitionTimes = Dates.datetime2unix(transitionTimes)
	    end
	  # Stuff might be unsorted
	  traceTimesSorted = sort(traceTimes)
	  sortedCalibIndices = sortperm(calibTimes)
	  calibTimesSorted = calibTimes[sortedCalibIndices]
	  calibValuesSorted = calibValues[sortedCalibIndices]
	  transitionTimesSorted = sort(transitionTimes)

	  piecewiseTimes = Array{Float64}(0)
	  piecewiseValues = Array{Float64}(0)
	  currentValue = 0
	  currentCalibIndex = 1
	  currentTransitionIndex = 1
	  currentPiecewiseIndex = 1

	  # Drop transitions before start of trace
	  while (transitionTimesSorted[currentTransitionIndex] < traceTimes[1]) && (currentTransitionIndex < length(transitionTimesSorted))
	    currentTransitionIndex += 1
	  end
	  # Drop transitions before first calibration
	  while (transitionTimesSorted[currentTransitionIndex] < calibTimes[1]) && (currentTransitionIndex < length(transitionTimesSorted))
	    currentTransitionIndex += 1
	  end
	  println("First transition within data range: $currentTransitionIndex - $(transitionTimesSorted[currentTransitionIndex])")

	  # If no calib was done before start of trace, insert copy of first calib
	  println("No Calib before start of data range found, copying first successive calib")
	  push!(piecewiseTimes, minimum([traceTimes[1] calibTimes[1]]))
	  push!(piecewiseValues, calibValuesSorted[1])
	  currentPiecewiseIndex += 1

	  # Iterate over calibs and transitions
	  while (currentCalibIndex <= length(calibTimesSorted))
	    if currentTransitionIndex < length(transitionTimesSorted)
	      println("CurrCalib: $currentCalibIndex ($(calibTimesSorted[currentCalibIndex])), CurrTransition: $currentTransitionIndex ($(transitionTimesSorted[currentTransitionIndex]))")
	    end
	    # define piecewise values depending on what comes next, calib or transition
	    if ((currentTransitionIndex > length(transitionTimesSorted)) || (calibTimesSorted[currentCalibIndex] < transitionTimesSorted[currentTransitionIndex])) && (currentCalibIndex <= length(calibTimesSorted))
	      # add one point for a calib
	      println("Adding calib point  $(calibValuesSorted[currentCalibIndex]) at $(calibTimesSorted[currentCalibIndex])")
	      push!(piecewiseTimes, calibTimesSorted[currentCalibIndex])
	      push!(piecewiseValues, calibValuesSorted[currentCalibIndex])
	      currentPiecewiseIndex += 1
	      currentCalibIndex += 1
	    else
	      # skip all following transitions that come before a calibration, since there are no calibs in between
	      while (currentTransitionIndex < length(transitionTimesSorted)) && (transitionTimesSorted[currentTransitionIndex+1] < calibTimesSorted[currentCalibIndex])
		println("Skipping transition $currentTransitionIndex, there are no calibs in between!")
		currentTransitionIndex +=1
	      end
	      # add one point with last value before transition
	      println("Adding pre-transition point $(piecewiseValues[end]) at $(transitionTimesSorted[currentTransitionIndex])")
	      push!(piecewiseTimes, transitionTimesSorted[currentTransitionIndex] - 1e-99)
	      push!(piecewiseValues, piecewiseValues[end])
	      currentPiecewiseIndex += 1

	      # add one point with next value after transition
	      println("Adding post-transition point $(calibValuesSorted[currentCalibIndex]) at $(transitionTimesSorted[currentTransitionIndex])")
	      push!(piecewiseTimes, transitionTimesSorted[currentTransitionIndex] + 1e-99)
	      push!(piecewiseValues, calibValuesSorted[currentCalibIndex])
	      currentPiecewiseIndex += 1

	      currentTransitionIndex += 1
	    end
	  end
	  PyPlot.figure()
	  PyPlot.plot(piecewiseTimes, piecewiseValues, "x-")

	  return InterpolationFunctions.interpolate(traceTimes, piecewiseTimes, piecewiseValues)
	end
end

module ExportFunctions
	using ..InterpolationFunctions
	using Dates
	import ..MasslistFunctions

	export exportTracesCSV, exportTracesCSVLossCorr, toMatlabTime, fromMatlabTime

	function exportTracesCSV(saveFolderPath, elementNames, compositions, times, traces; average=0)
	  sumformulas = MasslistFunctions.sumFormulaStringListFromCompositionArrayList(compositions)
	  f = open("$saveFolderPath/ptr3compositions.txt", "w")
	  writedlm(f, hcat(["Mass" "SumFormula"],reshape(elementNames,(1,length(elementNames)))))
	  writedlm(f, hcat(MasslistFunctions.massFromCompositionArrayList(compositions),sumformulas , compositions'))
	  close(f)
	  f = open("$saveFolderPath/ptr3traces.csv", "w")
	  writedlm(f, hcat(["Time"], reshape(sumformulas,(1,length(sumformulas)))))
	  if (average==0)
	    writedlm(f, hcat(times ,traces))
	  else
	    writedlm(f, hcat(averageSamples(times,average) ,averageSamples(traces,average)))
	  end
	  close(f)
	end

	function exportTracesCSVLossCorr(saveFolderPath, elementNames, compositions, times, traces, lossfactor, lossfactorerr, corrfactor, corrfactorerr, corrnotes; average=0)
	  sumformulas = MasslistFunctions.sumFormulaStringListFromCompositionArrayList(compositions)
	  f = open("$saveFolderPath/ptr3compositions.txt", "w")
	  writedlm(f, hcat(["Mass"	"SumFormula"],reshape(elementNames,(1,length(elementNames))),["LossFactor"	"LossFactorError"	"CorrFactor"	"CorrFactorErr"	"CorrNotes"]))
	  writedlm(f, hcat(MasslistFunctions.massFromCompositionArrayList(compositions),sumformulas , compositions', lossfactor, lossfactorerr, corrfactor, corrfactorerr, corrnotes))
	  close(f)
	  f = open("$saveFolderPath/ptr3tracesInletLossCorr.csv", "w")
	  writedlm(f, hcat(["Time"], reshape(sumformulas,(1,length(sumformulas)))))
	  if (average==0)
	    writedlm(f, hcat(times , (corrfactor.*traces' )' ))
	  else
	    writedlm(f, hcat(averageSamples(times,average) ,(corrfactor.*(averageSamples(traces,average))' )' ))
	  end
	  close(f)
	end

	################ EXAMPLE from Matlab: 737551.673515479 should be 06-May-2019 16:09:51 #############
	function toMatlabTime(t::Dates.DateTime)
	    timespan = ((t+Dates.Day(1)) - Dates.DateTime(0,1,1,0,0,0))
	    millis::Float64 = Dates.value(timespan)
	    return millis/24/3600/1000
	end

	function fromMatlabTime(timestamp::Number)
	    days=Int(floor(timestamp))
	    millisecondsRemainder = Int(round((timestamp-days)*24*3600*1000))
	    return Dates.DateTime(0,1,1,0,0,0)+Dates.Day(days)+Dates.Millisecond(millisecondsRemainder)-Dates.Day(1)
	end
end

module ImportFunctions
	using Dates
	using CSV
	using DataFrames

	export importFromCSV

	function importFromCSV(file; headerrow = 1)
		data = DataFrame(CSV.File(file, header = headerrow))
		if ("year" in names(data)) && ("month" in names(data)) && ("day" in names(data)) && ("hour" in names(data)) && ("minute" in names(data))
			datetime = Dates.DateTime.(data.year, data.month, data.day, data.hour, data.minute)
		elseif ("datetime" in names(data))
			println(data.datetime[1])
			println("What is the timestringformat of this datetime stamp?")
			formatString = readline()
			datetime = Dates.DateTime.(data.datetime, formatString)
		elseif ("unixtime" in names(data))
			datetime = Dates.unix2datetime.(data.unixtime)
		end
		return (datetime, data)
	end
end

ext: jl


text:
module PowSysMod

include("PowSysMod_body.jl")

export load_OPFproblems
export build_Problem!
export build_globalpb!
export MatpowerInput, MatpowerSimpleInput, IIDMInput, GOCInput

end

max_stars_repo_name: Gimaju/GOC

max_issues_repo_name: Gimaju/GOC

max_forks_repo_name: Gimaju/GOC


text:
module Harminv

using Compat

const SMALL = 1e-8
const NPOW = 8
const UNITY_THRESH = 1e-4
const SINGULAR_THRESHOLD = 1e-5

struct HarminvData{C,F}
    c::Vector{C}
    n::Int
    K::Int
    J::Int
    nfreqs::Ref{Int}
    fmin::F
    fmax::F
    z::Vector{C}
    U0::Matrix{C}
    U1::Matrix{C}
    G0::Vector{C}
    G0_M::Vector{C}
    D0::Vector{C}
    B::Ref{Matrix{C}}
    u::Ref{Vector{C}}
    amps::Vector{C}
    errs::Vector{F}
end

harmonics(fmin, fmax, nf) = @. exp(-im * 2π * fmin + (fmin-fmax)/(nf-1) * (0:nf-1))

function HarminvData(signal::Vector{C}, fmin, fmax, nf, nfreqs = Ref(-1), z = harmonics(fmin,fmax,nf)) where C
    J = nf
    n = length(signal)

    @assert nf > 1 || nfreqs[] != -1
    @assert n > 0
    @assert !isempty(signal)
    @assert fmin < fmax

    K = div(n,2) - 1
    U0 = Matrix{C}(undef, J, J)
    U1 = Matrix{C}(undef, J, J)
    amps = Vector{C}(undef, J)
    errs = Vector{Float64}(undef, J)
    G0, G0_M, D0 = generate_U(U0, U0, 0, signal, n, K, J, J, z, z, nothing, nothing, nothing)
    HarminvData(signal, n, K, nf, nfreqs, fmin, fmax, z, U0, U1,
                G0, G0_M, D0, Ref{Matrix{C}}(), Ref{Vector{C}}(), C[], Float64[])
end

# for iterating on a previous solution, e.g. in solve_again
function HarminvData(d::HarminvData)
    newd = HarminvData(d.c, d.fmin, d.fmax, d.nfreqs[], Ref(0), d.u[][:])
end

function generate_U(U::AbstractMatrix{C}, U1, p, c::AbstractArray{C},
                    n, K, J, J2, z::AbstractArray{C}, z2::AbstractArray{C},
                    G0 , G0_M, D0) where C
    M = K - 1
    @assert n >= 2K + p

    if G0 != nothing
        G = G0#[:]
        G_M = G0_M#[:]
        D = D0#[:]
    else
        G = zeros(C, J)
        G_M = zeros(C, J)
        D = zeros(C, J)
    end

    z_m = ones(C, J)
    z_inv = 1.0 ./ z
    z_M = z.^-M
    if z !== z2
        z2_inv = 1 ./ z2
        z2_m = ones(C, J2)
        z2_M = z2.^-M
        G2 = zeros(C, J2)
        G2_M = zeros(C, J2)
    end

    for m = 1:M
        c1 = c[m + p]
        c2 = c[m + p + M + 1]
        d = float(m)
        d2 = float(M - m - 1)
        if G0 == nothing
            for i = 1:J
                x1 = z_m[i] * c1
                x2 = z_m[i] * c2
                G[i] += x1
                G_M[i] += x2
                D[i] += x1*d + x2*d2*z_M[i] * z_inv[i]
                if m % NPOW == 0
                    z_m[i] = z_inv[i]^m
                else
                    z_m[i] *= z_inv[i]
                end
            end
        end
        if z !== z2
            for i = 1:J2
                G2[i] = z2_m[i] * c1
                G2_M[i] = z2_m[i] * c2
                if m % NPOW == 0
                    z2_m[i] = z2_inv[i]^m
                else
                    z2_m[i] *= z2_inv[i]
                end
            end
        end
    end

    if z !== z2
        for i = 1:J
            for j = 1:J2
                if abs(z[i] - z2[j]) < SMALL
                    U[i,j] = D[i]
                else
                    U[i,j] = ( z[i] * G2[j] - z2[j] * G[i] +
                                        z2_M[j] * G_M[i] - z_M[i] * G2_M[j] ) /
                                      (z[i] - z2[j])
                end
            end
        end
        if U1 !== nothing
            for i = 1:J
                for j = 1:J2
                    if abs(z[i] - z2[j]) < SMALL
                        U1[i,j] = z[i] * (D[i] - G[i]) +
                                           z_M[i] * G_M[i]
                    else
                        U1[i,j] = ( z[i] * z2[j] * (G2[j] - G[i]) +
                                             z2_M[j] * z[i] * G_M[i] -
                                             z_M[i] * z2[j] * G2_M[j] ) /
                                           (z[i] - z2[j])
                    end
                end
            end
        end
    else
        for i = 1:J2
            U[i,i] = D[i]
            for j = (i+2):J
                U[i,j] = ( z[i] * G[j] - z[j] * G[i] +
                           z_M[j] * G_M[i] - z_M[i] * G_M[j] ) /
                         (z[i] - z[j])
            end
        end
        if U1 !== nothing
            for i = 1:J2
                U1[i,i] = z[i] * (D[i] - G[i]) + z_M[i] * G_M[i]
                for j = (i+2):J
                    U1[i,j] = ( z[i] * z[j] * (G[j] - G[i]) +
                                z_M[j] * z[i] * G_M[i] -
                                z_M[i] * z[j] * G_M[j] ) /
                              (z[i] - z[j])
                end
            end
        end
        if z === z2
            U .= Symmetric(U, :U)
            if U1 !== nothing
                U1 .= Symmetric(U1, :U)
            end
        end
    end

    return G, G_M, D
end

function symmetric_normalize!(vect)
    symm_norm = sqrt(sum(vect.*vect))
    vect ./= symm_norm
end

function solve_eigenvects(A0)
    v0, mm = eig(A0)
    nvals = length(v0)
    # symmetric-normalized eigenvectors
    @inbounds for col = 1:nvals
        symmetric_normalize!(@view mm[:,col])
    end
    v0, mm
end

function solve_once!(d::HarminvData{C,F}) where {C,F}
    zone = one(C)
    zzero = zero(C)

    J = d.J
    v0, V0 = solve_eigenvects(d.U0)
    abs_v0 = abs.(v0)
    max_v0 = maximum(abs_v0)
    c0 = zero(eltype(v0))
    threshold = SINGULAR_THRESHOLD * abs(max_v0)

    d.nfreqs[] = J
    for i = 1:J
        if abs_v0[i] < threshold
            v0[i] = c0
            d.nfreqs[] -= 1
        else
            j = findfirst(v0, c0)
            if j > 0 && j < i
                V0[:, j] = V0[:, i]
                v0[j] = v0[i]
                v0[i] = c0
            else
                j = i
            end
            scale!(V0[:,j], 1/sqrt(v0[j]))
        end
    end
    # the good eigenvectors of V0
    V = @view V0[:,1:d.nfreqs[]]
    # doesn't seem to agree with the comments
    # in @stevengj's harminv.c but does with the
    # actual code + ZGEMM BLAS reference.
    d.B[] = d.U1 * V
    H1 = V' * d.B[]
    d.u[], V1 = solve_eigenvects(H1)
    d.B[] = V * V1
end

function solve_again!(d::HarminvData)
    if d.nfreqs == 0
        solve_once!(d)
    end
    println(d.nfreqs[])

    mode_ok = BitArray(d.nfreqs[])
    new_d = HarminvData(d)
    solve_once!(d)
end

# TODO: actually use ok
function solve_ok_modes!(d::HarminvData, ok::Function, ok_d)
    solve_once!(d)
    cur_nf = d.nfreqs[]
    while true
        prev_nf = cur_nf
        solve_again!(d, ok, ok_d)
        cur_nf = d.nfreqs[]
        if cur_nf >= prev_nf
            break
        end
    end

end

solve!(d::HarminvData) = solve_ok_modes!(d, nothing, nothing)

# frequency_errors!(d::HarminvData) =

function u_near_unity(u, n)
    nlgabsu = n * log(abs(u))
    nlgabsu > log(UNITY_THRESH) && nlgabsu < -log(UNITY_THRESH)
end

function compute_amplitudes(d::HarminvData{C,F}) where {C,F}
    a = Array{C}(undef, d.nfreqs[])
    u = Array{C}(undef, d.nfreqs[])
    ku = 1

    for k = 1:d.nfreqs[]
        if u_near_unity(d.u[][k], d.n)
            u[ku] = d.u[][k]
            ku += 1
        end
    end

    nu = ku - 1
    Uu = Matrix{C}(undef, d.J, nu)
    generate_U(Uu, nothing, 0, d.c, d.n, d.K, d.J, nu, d.z, u, d.G0, d.G0_M, d.D0)

    for k = 1:d.nfreqs[]
        asum = zero(C)
        if u_near_unity(d.u[][k], d.n)
            for j = 1:d.J
                asum += d.B[][j,k] * Uu[j, ku]
            end
            asum /= d.K
            ku += 1
        else
            for j = 1:d.J
                asum += d.B[][j,k] * d.G0[j]
            end
        end
        a[k] = asum^2
    end
    return a
end

end


text:
export WallBoundaries, FreeSlip, NoSlip

abstract type SlipComponent{N} end

struct FreeSlip{N} <: SlipComponent{N}
    FreeSlip(::Dim{N}, ::Integer) where {N} = new{N}()
end

Base.show(io::IO, ::Type{FreeSlip{N}}) where {N} = print(io, "FreeSlip")

mutable struct NoSlip{N} <: SlipComponent{N}
    virtual_velocities::Vector{SVectorF{N}}

    NoSlip(::Dim{N}, number_of_particles::Integer) where {N} =
        new{N}([zeros(SVectorF{N}) for _ = 1:number_of_particles])
end

Base.show(io::IO, ::Type{NoSlip{N}}) where {N} = print(io, "NoSlip")

mutable struct WallBoundaries{N,S} <: BoundaryComponent{N}
    positions::Vector{SVectorF{N}}
    velocities::Vector{SVectorF{N}}
    acceleration::SVectorF{N}
    virtual_pressures::Vector{Float}
    slip_component::S

    function WallBoundaries(
        positions::Vector{SVectorF{N}},
        velocities::Vector{SVectorF{N}},
        acceleration::SVectorF{N},
        slip_component::Type{S} = FreeSlip,
    ) where {N,S<:SlipComponent{N}}
        number_of_particles = length(positions)
        @checkeq(length(velocities), number_of_particles)
        new{N,S}(
            positions,
            velocities,
            acceleration,
            fill(0.0, number_of_particles),
            slip_component(Dim{N}(), number_of_particles),
        )
    end
end

Base.show(io::IO, ::Type{WallBoundaries{N,S}}) where {N,S} =
    print(io, "WallBoundaries{$(S)}")

get_positions(wall::WallBoundaries) = wall.positions
get_velocities(wall::WallBoundaries) = wall.velocities
get_pressures(wall::WallBoundaries) = wall.virtual_pressures

function compute_mass_density(
    wall::WallBoundaries,
    energy_component::EnergyComponent,
    eos::EquationOfState,
    idx_fluid::Unsigned,
    idx_boundary::Unsigned,
)
    mass_density_from(
        eos,
        wall.virtual_pressures[idx_boundary],
        get_specific_energy(energy_component, idx_fluid),
    )
end

function updateboundaries!(
    wall::WallBoundaries{N},
    positions::AbstractVector{SVectorF{N}},
    velocities::Velocities{N},
    standard_mass::StandardMass{ConstantWidthKernel{N}},
    pressures::Pressures,
) where {N}
    update_wall_properties!(
        wall,
        positions,
        get_velocities(velocities),
        get_mass_densities(standard_mass),
        get_pressures(pressures),
        get_kernel(standard_mass),
    )
end

function evolveboundaries!(wall::WallBoundaries, Δt::Number)
    wall.positions .+= Δt * wall.velocities
end

function update_wall_properties!(
    wall::WallBoundaries{N},
    fluid_positions::AbstractVector{SVectorF{N}},
    fluid_velocities::AbstractVector{SVectorF{N}},
    fluid_mass_densities::AbstractVector{Float},
    fluid_pressures::AbstractVector{Float},
    W::ConstantWidthKernel{N},
) where {N}
    h = get_kernel_width(W)
    h⁻¹ = get_inverse_kernel_width(W)
    for (i, (rᵢ, vᵢ)) in enumerate(zip(wall.positions, wall.velocities))
        ΣPⱼWᵢⱼ::Float = 0.0
        ΣρⱼrᵢⱼWᵢⱼ = zeros(SVectorF{N})
        ΣvⱼWᵢⱼ = zeros(SVectorF{N})
        ΣWᵢⱼ::Float = 0.0
        for (rⱼ, vⱼ, ρⱼ, Pⱼ) in zip(
            fluid_positions,
            fluid_velocities,
            fluid_mass_densities,
            fluid_pressures,
        )
            rᵢⱼ = rⱼ - rᵢ
            qᵢⱼ = norm(rᵢⱼ) * h⁻¹
            if nonzeroat(W, qᵢⱼ)
                Wᵢⱼ = W(qᵢⱼ, h)
                ΣPⱼWᵢⱼ += Pⱼ * Wᵢⱼ
                ΣρⱼrᵢⱼWᵢⱼ .+= (ρⱼ * Wᵢⱼ) * rᵢⱼ
                accumulate_virtual_velocity(
                    ΣvⱼWᵢⱼ,
                    vⱼ,
                    Wᵢⱼ,
                    wall.slip_component,
                )
                ΣWᵢⱼ += W_rᵢⱼ_h
            end
        end
        ΣWᵢⱼ⁻¹ = 1 / ΣWᵢⱼ
        wall.virtual_pressures[i] =
            (ΣPⱼWᵢⱼ - wall.acceleration ⋅ ΣρⱼrᵢⱼWᵢⱼ) * ΣWᵢⱼ⁻¹
        set_virtual_velocity(wall.slip_component, vᵢ, ΣvⱼWᵢⱼ, ΣWᵢⱼ⁻¹, i)
    end
end

function accumulate_virtual_velocity(
    ::SVectorF{N},
    ::SVectorF{N},
    ::Number,
    ::FreeSlip{N},
) where {N} end

function set_virtual_velocity(
    ::FreeSlip{N},
    ::SVectorF{N},
    ::SVectorF{N},
    ::Number,
    ::Integer,
) where {N} end

function accumulate_virtual_velocity(
    ΣvⱼWᵢⱼ::SVectorF{N},
    vⱼ::SVectorF{N},
    Wᵢⱼ::Number,
    no_slip::NoSlip{N},
) where {N}
    ΣvⱼWᵢⱼ .+= vⱼ * Wᵢⱼ
end

function set_virtual_velocity(
    no_slip::NoSlip{N},
    vᵢ::SVectorF{N},
    ΣvⱼWᵢⱼ::SVectorF{N},
    ΣWᵢⱼ⁻¹::Number,
    i::Integer,
) where {N}
    ṽᵢ = ΣvⱼWᵢⱼ * ΣWᵢⱼ⁻¹
    no_slip.virtual_velocities[i] = 2 * vᵢ - ṽᵢ
end

max_stars_repo_name: lars-frogner/Whirl

max_issues_repo_name: lars-frogner/Whirl

max_forks_repo_name: lars-frogner/Whirl


text:
using Test, MPI
import GaussQuadrature
using ClimateMachine
using ClimateMachine.ConfigTypes
using ClimateMachine.Mesh.Topologies
using ClimateMachine.Mesh.Grids
using ClimateMachine.Mesh.Geometry
using ClimateMachine.Mesh.Interpolation
using ClimateMachine.Writers
using StaticArrays
using GPUifyLoops

using ClimateMachine.VariableTemplates
#------------------------------------------------
using ClimateMachine.DGmethods
using ClimateMachine.DGmethods.NumericalFluxes
using ClimateMachine.MPIStateArrays
using ClimateMachine.ODESolvers
using ClimateMachine.GenericCallbacks
using ClimateMachine.Atmos
using ClimateMachine.VariableTemplates
using ClimateMachine.MoistThermodynamics
using ClimateMachine.TicToc
using LinearAlgebra
using StaticArrays
using Logging, Printf, Dates
using ClimateMachine.VTK

using ClimateMachine.Atmos: vars_state_conservative, vars_state_auxiliary

using CLIMAParameters
using CLIMAParameters.Planet: R_d, planet_radius, grav, MSLP
struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

using Random
using Statistics
const seed = MersenneTwister(0)

const ArrayType = ClimateMachine.array_type()
#-------------------------------------
function Initialize_Brick_Interpolation_Test!(
    bl,
    state::Vars,
    aux::Vars,
    (x, y, z),
    t,
)
    FT = eltype(state)
    # Dummy variables for initial condition function
    state.ρ = FT(0)
    state.ρu = SVector{3, FT}(0, 0, 0)
    state.ρe = FT(0)
    state.moisture.ρq_tot = FT(0)
end
#------------------------------------------------
struct TestSphereSetup{DT}
    p_ground::DT
    T_initial::DT
    domain_height::DT

    function TestSphereSetup(
        p_ground::DT,
        T_initial::DT,
        domain_height::DT,
    ) where {DT <: AbstractFloat}
        return new{DT}(p_ground, T_initial, domain_height)
    end
end
#----------------------------------------------------------------------------
function (setup::TestSphereSetup)(bl, state, aux, coords, t)
    # callable to set initial conditions
    FT = eltype(state)
    _grav::FT = grav(param_set)
    _R_d::FT = R_d(param_set)

    z = altitude(bl, aux)

    scale_height::FT = _R_d * setup.T_initial / _grav
    p::FT = setup.p_ground * exp(-z / scale_height)
    e_int = internal_energy(bl.param_set, setup.T_initial)
    e_pot = gravitational_potential(bl.orientation, aux)

    # TODO: Fix type instability: typeof(setup.T_initial) == typeof(p) fails
    state.ρ = air_density(bl.param_set, FT(setup.T_initial), p)
    state.ρu = SVector{3, FT}(0, 0, 0)
    state.ρe = state.ρ * (e_int + e_pot)
    return nothing
end
#----------------------------------------------------------------------------
function run_brick_interpolation_test()
    ClimateMachine.init()
    for FT in (Float32, Float64)
        ArrayType = ClimateMachine.array_type()
        DA = ClimateMachine.array_type()
        mpicomm = MPI.COMM_WORLD
        root = 0
        pid = MPI.Comm_rank(mpicomm)
        npr = MPI.Comm_size(mpicomm)

        xmin, ymin, zmin = 0, 0, 0                   # defining domain extent
        xmax, ymax, zmax = 2000, 400, 2000
        xres = [FT(10), FT(10), FT(10)] # resolution of interpolation grid

        Ne = (20, 4, 20)

        polynomialorder = 5 #8# 5 #8 #4
        #-------------------------
        _x, _y, _z = ClimateMachine.Mesh.Grids.vgeoid.x1id,
        ClimateMachine.Mesh.Grids.vgeoid.x2id,
        ClimateMachine.Mesh.Grids.vgeoid.x3id
        #-------------------------

        brickrange = (
            range(FT(xmin); length = Ne[1] + 1, stop = xmax),
            range(FT(ymin); length = Ne[2] + 1, stop = ymax),
            range(FT(zmin); length = Ne[3] + 1, stop = zmax),
        )
        topl = StackedBrickTopology(
            mpicomm,
            brickrange,
            periodicity = (true, true, false),
        )
        grid = DiscontinuousSpectralElementGrid(
            topl,
            FloatType = FT,
            DeviceArray = DA,
            polynomialorder = polynomialorder,
        )
        model = AtmosModel{FT}(
            AtmosLESConfigType,
            param_set;
            ref_state = NoReferenceState(),
            turbulence = ConstantViscosityWithDivergence(FT(0)),
            source = (Gravity(),),
            init_state_conservative = Initialize_Brick_Interpolation_Test!,
        )

        dg = DGModel(
            model,
            grid,
            RusanovNumericalFlux(),
            CentralNumericalFluxSecondOrder(),
            CentralNumericalFluxGradient(),
        )

        Q = init_ode_state(dg, FT(0))
        #------------------------------
        x1 = @view grid.vgeo[:, _x, :]
        x2 = @view grid.vgeo[:, _y, :]
        x3 = @view grid.vgeo[:, _z, :]

        fcn(x, y, z) = sin.(x) .* cos.(y) .* cos.(z) # sample function
        #----calling interpolation function on state variable # st_idx--------------------------
        nvars = size(Q.data, 2)

        for vari in 1:nvars
            Q.data[:, vari, :] = fcn(x1 ./ xmax, x2 ./ ymax, x3 ./ zmax)
        end

        xbnd = Array{FT}(undef, 2, 3)

        xbnd[1, 1] = FT(xmin)
        xbnd[2, 1] = FT(xmax)
        xbnd[1, 2] = FT(ymin)
        xbnd[2, 2] = FT(ymax)
        xbnd[1, 3] = FT(zmin)
        xbnd[2, 3] = FT(zmax)
        #----------------------------------------------------------
        x1g = collect(range(xbnd[1, 1], xbnd[2, 1], step = xres[1]))
        nx1 = length(x1g)
        x2g = collect(range(xbnd[1, 2], xbnd[2, 2], step = xres[2]))
        nx2 = length(x2g)
        x3g = collect(range(xbnd[1, 3], xbnd[2, 3], step = xres[3]))
        nx3 = length(x3g)

        filename = "test.nc"
        varnames = ("ρ", "ρu", "ρv", "ρw", "e", "other")

        intrp_brck = InterpolationBrick(grid, xbnd, x1g, x2g, x3g)        # sets up the interpolation structure
        iv = DA(Array{FT}(undef, intrp_brck.Npl, nvars))                  # allocating space for the interpolation variable
        if pid == 0
            fiv = DA(Array{FT}(undef, nx1, nx2, nx3, nvars))    # allocating space for the full interpolation variables accumulated on proc# 0
        else
            fiv = DA(Array{FT}(undef, 0, 0, 0, 0))
        end
        interpolate_local!(intrp_brck, Q.data, iv)                    # interpolation
        accumulate_interpolated_data!(intrp_brck, iv, fiv)      # write interpolation data to file
        #------------------------------
        err_inf_dom = zeros(FT, nvars)

        x1g = intrp_brck.x1g
        x2g = intrp_brck.x2g
        x3g = intrp_brck.x3g

        if pid == 0
            nx1 = length(x1g)
            nx2 = length(x2g)
            nx3 = length(x3g)
            x1 = Array{FT}(undef, nx1, nx2, nx3)
            x2 = similar(x1)
            x3 = similar(x1)

            fiv_cpu = Array(fiv)

            for k in 1:nx3, j in 1:nx2, i in 1:nx1
                x1[i, j, k] = x1g[i]
                x2[i, j, k] = x2g[j]
                x3[i, j, k] = x3g[k]
            end
            fex = fcn(x1 ./ xmax, x2 ./ ymax, x3 ./ zmax)

            for vari in 1:nvars
                err_inf_dom[vari] =
                    maximum(abs.(fiv_cpu[:, :, :, vari] .- fex[:, :, :]))
            end
        end

        MPI.Bcast!(err_inf_dom, root, mpicomm)

        if FT == Float64
            toler = 1.0E-9
        elseif FT == Float32
            toler = 1.0E-6
        end

        if maximum(err_inf_dom) > toler
            if pid == 0
                println("err_inf_domain = $(maximum(err_inf_dom)) is larger than prescribed tolerance of $toler")
            end
            MPI.Barrier(mpicomm)
        end
        @test maximum(err_inf_dom) < toler
    end
    return nothing
    #----------------
end #function run_brick_interpolation_test
#----------------------------------------------------------------------------

#----------------------------------------------------------------------------
# Cubed sphere, lat/long interpolation test
#----------------------------------------------------------------------------
function run_cubed_sphere_interpolation_test()
    ClimateMachine.init()
    for FT in (Float32, Float64) #Float32 #Float64
        DA = ClimateMachine.array_type()
        device = ClimateMachine.array_type() <: Array ? CPU() : CUDA()
        mpicomm = MPI.COMM_WORLD
        root = 0
        pid = MPI.Comm_rank(mpicomm)
        npr = MPI.Comm_size(mpicomm)

        domain_height = FT(30e3)

        polynomialorder = 5
        numelem_horz = 6
        numelem_vert = 4

        #-------------------------
        _x, _y, _z = ClimateMachine.Mesh.Grids.vgeoid.x1id,
        ClimateMachine.Mesh.Grids.vgeoid.x2id,
        ClimateMachine.Mesh.Grids.vgeoid.x3id
        _ρ, _ρu, _ρv, _ρw = 1, 2, 3, 4
        #-------------------------
        _planet_radius::FT = planet_radius(param_set)

        vert_range = grid1d(
            _planet_radius,
            FT(_planet_radius + domain_height),
            nelem = numelem_vert,
        )

        lat_res = FT(1) # 1 degree resolution
        long_res = FT(1) # 1 degree resolution
        nel_vert_grd = 20 #100 #50 #10#50
        rad_res = FT((vert_range[end] - vert_range[1]) / FT(nel_vert_grd)) #1000.00    # 1000 m vertical resolution
        #----------------------------------------------------------
        _MSLP::FT = MSLP(param_set)
        setup = TestSphereSetup(_MSLP, FT(255), FT(30e3))

        topology = StackedCubedSphereTopology(mpicomm, numelem_horz, vert_range)

        grid = DiscontinuousSpectralElementGrid(
            topology,
            FloatType = FT,
            DeviceArray = DA,
            polynomialorder = polynomialorder,
            meshwarp = ClimateMachine.Mesh.Topologies.cubedshellwarp,
        )

        model = AtmosModel{FT}(
            AtmosLESConfigType,
            param_set;
            orientation = SphericalOrientation(),
            ref_state = NoReferenceState(),
            turbulence = ConstantViscosityWithDivergence(FT(0)),
            moisture = DryModel(),
            source = nothing,
            init_state_conservative = setup,
        )

        dg = DGModel(
            model,
            grid,
            RusanovNumericalFlux(),
            CentralNumericalFluxSecondOrder(),
            CentralNumericalFluxGradient(),
        )

        Q = init_ode_state(dg, FT(0))

        device = typeof(Q.data) <: Array ? CPU() : CUDA()
        #------------------------------
        x1 = @view grid.vgeo[:, _x, :]
        x2 = @view grid.vgeo[:, _y, :]
        x3 = @view grid.vgeo[:, _z, :]

        xmax = _planet_radius
        ymax = _planet_radius
        zmax = _planet_radius

        fcn(x, y, z) = sin.(x) .* cos.(y) .* cos.(z) # sample function

        nvars = size(Q.data, 2)

        for i in 1:nvars
            Q.data[:, i, :] .= fcn(x1 ./ xmax, x2 ./ ymax, x3 ./ zmax)
        end
        #------------------------------
        lat_min, lat_max = FT(-90.0), FT(90.0)            # inclination/zeinth angle range
        long_min, long_max = FT(-180.0), FT(180.0)     # azimuthal angle range
        rad_min, rad_max = vert_range[1], vert_range[end] # radius range

        lat_grd = collect(range(lat_min, lat_max, step = lat_res))
        n_lat = length(lat_grd)
        long_grd = collect(range(long_min, long_max, step = long_res))
        n_long = length(long_grd)
        rad_grd = collect(range(rad_min, rad_max, step = rad_res))
        n_rad = length(rad_grd)

        _ρu, _ρv, _ρw = 2, 3, 4

        filename = "test.nc"
        varnames = ("ρ", "ρu", "ρv", "ρw", "e")
        projectv = true

        intrp_cs = InterpolationCubedSphere(
            grid,
            collect(vert_range),
            numelem_horz,
            lat_grd,
            long_grd,
            rad_grd,
        ) # sets up the interpolation structure
        iv = DA(Array{FT}(undef, intrp_cs.Npl, nvars))             # allocating space for the interpolation variable
        if pid == 0
            fiv = DA(Array{FT}(undef, n_rad, n_lat, n_long, nvars))    # allocating space for the full interpolation variables accumulated on proc# 0
        else
            fiv = DA(Array{FT}(undef, 0, 0, 0, 0))
        end

        interpolate_local!(intrp_cs, Q.data, iv)                   # interpolation
        project_cubed_sphere!(intrp_cs, iv, (_ρu, _ρv, _ρw))         # project velocity onto unit vectors along rad, lat & long
        accumulate_interpolated_data!(intrp_cs, iv, fiv)           # accumulate interpolated data on to proc# 0
        #----------------------------------------------------------
        # Testing
        err_inf_dom = zeros(FT, nvars)
        rad = Array(intrp_cs.rad_grd)
        lat = Array(intrp_cs.lat_grd)
        long = Array(intrp_cs.long_grd)
        fiv_cpu = Array(fiv)
        if pid == 0
            nrad = length(rad)
            nlat = length(lat)
            nlong = length(long)
            x1g = Array{FT}(undef, nrad, nlat, nlong)
            x2g = similar(x1g)
            x3g = similar(x1g)

            fex = zeros(FT, nrad, nlat, nlong, nvars)

            for vari in 1:nvars
                for k in 1:nlong, j in 1:nlat, i in 1:nrad
                    x1g_ijk = rad[i] * cosd(lat[j]) * cosd(long[k]) # inclination -> latitude; azimuthal -> longitude.
                    x2g_ijk = rad[i] * cosd(lat[j]) * sind(long[k]) # inclination -> latitude; azimuthal -> longitude.
                    x3g_ijk = rad[i] * sind(lat[j])

                    fex[i, j, k, vari] =
                        fcn(x1g_ijk / xmax, x2g_ijk / ymax, x3g_ijk / zmax)
                end
            end

            if projectv
                for k in 1:nlong, j in 1:nlat, i in 1:nrad
                    fex[i, j, k, _ρu] =
                        fex[i, j, k, _ρ] * cosd(lat[j]) * cosd(long[k]) +
                        fex[i, j, k, _ρ] * cosd(lat[j]) * sind(long[k]) +
                        fex[i, j, k, _ρ] * sind(lat[j])

                    fex[i, j, k, _ρv] =
                        -fex[i, j, k, _ρ] * sind(lat[j]) * cosd(long[k])
                    -fex[i, j, k, _ρ] * sind(lat[j]) * sind(long[k]) +
                    fex[i, j, k, _ρ] * cosd(lat[j])

                    fex[i, j, k, _ρw] =
                        -fex[i, j, k, _ρ] * sind(long[k]) +
                        fex[i, j, k, _ρ] * cosd(long[k])
                end
            end

            for vari in 1:nvars
                err_inf_dom[vari] =
                    maximum(abs.(fiv_cpu[:, :, :, vari] .- fex[:, :, :, vari]))
            end
        end

        MPI.Bcast!(err_inf_dom, root, mpicomm)

        if FT == Float64
            toler = 1.0E-7
        elseif FT == Float32
            toler = 1.0E-6
        end

        if maximum(err_inf_dom) > toler
            if pid == 0
                println("err_inf_domain = $(maximum(err_inf_dom)) is larger than prescribed tolerance of $toler")
            end
            MPI.Barrier(mpicomm)
        end
        @test maximum(err_inf_dom) < toler
    end
    return nothing
end
#----------------------------------------------------------------------------
@testset "Interpolation tests" begin
    run_brick_interpolation_test()
    run_cubed_sphere_interpolation_test()
end
#------------------------------------------------

max_stars_repo_name: sandreza/CLIMA

max_issues_repo_name: sandreza/CLIMA

max_forks_repo_name: sandreza/CLIMA


text:
@testset "$TEST $G" begin

# TODO some tests here intermittently fail on travis-ci
# using @test_skip for the time being

if !@isdefined(nonbacktrack_embedding_dense)
    #= Spectral embedding of the non-backtracking matrix of `g`
    (see [Krzakala et al.](http://www.pnas.org/content/110/52/20935.short)).

    `g`: imput Graph
    `k`: number of dimensions in which to embed

    return : a matrix ϕ where ϕ[:,i] are the coordinates for vertex i.
    =#
    function nonbacktrack_embedding_dense(g::AGraph, k::Int)
        B, edgeid = nonbacktracking_matrix(g)
        λ, eigv = Arpack.eigs(B, nev=k+1, v0=ones(Float64, size(B,1)))
        ϕ = zeros(ComplexF64, k-1, nv(g))
        # TODO decide what to do with the stationary distribution ϕ[:,1]
        # this code just throws it away in favor of eigv[:,2:k+1].
        # we might also use the degree distribution to scale these vectors as is
        # common with the laplacian/adjacency methods.
        E = Edge{vertextype(g)}
        for n=1:k-1
            v= eigv[:,n+1]
            for i=1:nv(g)
                for e in in_edges(g, i)
                    j = src(e)
                    u = edgeid[e]
                    ϕ[n,i] += v[u]
                end
            end
        end
        return ϕ
    end
end

n = 10; k = 5
pg = PathGraph(n, G)
ϕ1 = nonbacktrack_embedding(pg, k)'
#check that we can use the implicit matvec in nonbacktrack_embedding
ϕ2 = nonbacktrack_embedding_dense(pg, k)'
@test size(ϕ2) == size(ϕ1)

# TODO remove if Nonbacktracking type won't be reintroduced again
# nbt = Nonbacktracking(pg)
# B, emap = nonbacktracking_matrix(pg)
# Bs = sparse(nbt)
# @test sparse(B) == Bs

# # check that matvec works
# x = ones(Float64, nbt.m)
# y = nbt * x
# z = B * x
# @test norm(y-z) < 1e-8

# #check that matmat works and Matrix(nbt) == B
# @test norm(nbt*eye(nbt.m) - B) < 1e-8

#check that this recovers communities in the path of cliques
n=10
g10 = CompleteGraph(n, G)
z = copy(g10)
for k=2:5
    z = blockdiag(z, g10)
    add_edge!(z, (k-1)*n, k*n)

    c = community_detection_nback(z, k)
    @test sort(union(c)) == [1:k;]
    a = collect(n:n:k*n)
    #@test_skip length(c[a]) == length(unique(c[a]))
    for i=1:k
        for j=(i-1)*n+1:i*n
            #@test_skip c[j] == c[i*n]
        end
    end

    c = community_detection_bethe(z, k)
    @test sort(union(c)) == [1:k;]
    a = collect(n:n:k*n)
    #@test_skip length(c[a]) == length(unique(c[a]))

    for i=1:k
        for j=(i-1)*n+1:i*n
            @test c[j] == c[i*n]
        end
    end

    c = community_detection_bethe(z)
    @test sort(union(c)) == [1:k;]
    a = collect(n:n:k*n)
    #@test_skip length(c[a]) == length(unique(c[a]))
    for i=1:k
        for j=(i-1)*n+1:i*n
            @test c[j] == c[i*n]
        end
    end
end

n=10
g10 = CompleteGraph(n, G)
z = copy(g10)
for k=2:5
    z = blockdiag(z, g10)
    add_edge!(z, (k-1)*n, k*n)
    c, ch = label_propagation(z)
    a = collect(n:n:k*n)
    a = Int[div(i-1,n)+1 for i=1:k*n]
    # check the number of community
    # TODO fix tests
    # @test length(unique(a)) == length(unique(c))
    #@test_skip length(unique(a))== length(unique(c))
    # check the partition
    # @test a == c
end

end # testset


text:
export logdensity

function logdensity(m::JointDistribution{A0,A,B,M},x) where {A0,A,B,M}
    _logdensity(M, m.model, m.args, x)
end

function logdensity(m::JointDistribution{A0,A,B,M},x, ::typeof(logdensity)) where {A0,A,B,M}
    _logdensity(M, m.model, m.args, x)
end

@gg M function _logdensity(_::Type{M}, _m::Model, _args, _data) where M <: TypeLevel{Module}
    Expr(:let,
        Expr(:(=), :M, from_type(M)),
        type2model(_m) |> sourceLogdensity() |> loadvals(_args, _data))
end

function sourceLogdensity()
    function(_m::Model)
        proc(_m, st :: Assign)     = :($(st.x) = $(st.rhs))
        proc(_m, st :: Sample)     = :(_ℓ += logdensity($(st.rhs), $(st.x)))
        proc(_m, st :: Return)     = nothing
        proc(_m, st :: LineNumber) = nothing

        wrap(kernel) = @q begin
            _ℓ = 0.0
            $kernel
            return _ℓ
        end

        buildSource(_m, proc, wrap) |> flatten
    end
end


text:
# =========================================================================== #
## Block - Primitives
# =========================================================================== #

abstract type AbstractBlock end
abstract type AbstractSource <: AbstractBlock end
abstract type AbstractSink <: AbstractBlock end

abstract type AbstractUnaryTransformer <: AbstractBlock end
abstract type AbstractBinaryTransformer <: AbstractBlock end

abstract type AbstractIntegrator <: AbstractBlock end

include("primitives_source.jl")
include("primitives_unary.jl")
include("primitives_binary.jl")
include("primitives_integrator.jl")

function Base.show(io::IO, block::T) where {T <: AbstractBlock}
    println(T)
    fnames = fieldnames(T)
    for fn in fnames
        print("\t")
        print(fn)
        print(": ")
        print("$(getfield(block, fn))")
        print("\n")
    end
end

function Base.show(io::IO, block::T) where {T <: AbstractInputPin}
    println(T)
    println("\tvalue: $(block.value)")
    # containing_object
end
function Base.show(io::IO, block::T) where {T <: AbstractOutputPin}
    println(T)
    println("\tvalue: $(block.value)")
    println("\tprotocol: $(block.protocol)")

    # containing_object
end

function is_integrator(block::T) where {T <: AbstractBlock}
    return false
end

function is_integrator(block::Integrator)
    return true
end

function has_input_pin(block::T) where {T <: AbstractBlock}
    return true
end

function has_input_pin(block::T) where {T <: AbstractSource}
    return false
end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
module Configurations

export @option, from_dict, from_kwargs, from_toml, to_dict, to_toml

using TOML
using ExproniconLite
using OrderedCollections

export no_default,
    Maybe,
    Reflect,
    # reflection
    field_default,
    # field_alias,
    type_alias,
    # traits
    is_option,
    # parse
    from_dict,
    from_kwargs,
    from_toml,
    from_toml_if_exists,
    # serialize
    to_dict,
    DuplicatedFieldError,
    DuplicatedAliasError,
    InvalidKeyError,
    TOMLStyle,
    YAMLStyle,
    JSONStyle

@static if VERSION < v"1.1"
    function fieldtypes(T::Type)
        ntuple(fieldcount(T)) do idx
            fieldtype(T, idx)
        end
    end
end

include("errors.jl")
include("types.jl")
include("convert.jl")
include("reflection.jl")
include("utils.jl")
include("codegen.jl")

include("from_dict.jl")
include("from_kwargs.jl")
include("from_toml.jl")
# include("parse.jl")
include("serialize.jl")

end


text:
using DSGE, Test, DelimitedFiles, JLD2, ModelConstructors, Random, FileIO
import DSGE: klein_transition_matrices, n_model_states, n_backward_looking_states

# What do you want to do?
check_steady_state = true
check_jacobian = true
check_solution = true
check_irfs = true
check_steady_state_calibrate = true
write_steady_state_calibrate = false
check_likelihood = true
write_likelihood = false

path = dirname(@__FILE__)

addtl_save_str = VERSION >= v"1.5" ? "_1.5" : ""

m = HetDSGEGovDebt(testing_gamma = true, ref_dir = HETDSGEGOVDEBT)
m <= Setting(:steady_state_only, true)
m <= Setting(:reduce_ell, true)
m <= Setting(:binsize, 10)

# Steady-state computation
if check_steady_state
    steadystate!(m)

    file = jldopen("$path/reference/steady_state_reduce_ell.jld2", "r")
    saved_Rk   = read(file, "Rk")
    saved_ω    = read(file, "omega")
    saved_kl   = read(file, "kl")
    saved_k    = read(file, "k")
    saved_y    = read(file, "y")
    saved_T    = read(file, "T")
    saved_ell  = read(file, "ell")
    saved_c    = read(file, "c")
    saved_μ    = read(file, "mu")
    saved_β    = read(file, "beta")
    close(file)

    @testset "Check steady state outputs" begin
        @test saved_Rk ≈ m[:Rkstar].value
        @test saved_ω  ≈ m[:ωstar].value
        @test saved_kl ≈ m[:klstar].value
        @test saved_k  ≈ m[:kstar].value
        @test saved_y  ≈ m[:ystar].value
        @test saved_T  ≈ m[:Tstar].value

        @test isapprox(saved_ell, m[:lstar].value, atol = .01)
        @test isapprox(saved_c, m[:cstar].value, atol = .01)
        @test isapprox(saved_μ, m[:μstar].value, atol = .01)
        # Tolerance of convergence of β is 1e-5
        @test isapprox(saved_β, m[:βstar].value, atol = .01)
    end
end

if check_jacobian
    m.testing = true # Ensures it tests against the unnormalized Jacobian
    m <= Setting(:steady_state_only, true)
    steadystate!(m)
    JJ = DSGE.jacobian(m)

    file = jldopen("$path/reference/jacobian_reduce_ell.jld2", "r")
    saved_JJ  = read(file, "JJ")
    KFP = read(file, "KFP")
    KKP = read(file, "KKP")
    LIIP = read(file, "LIIP")
    LYP = read(file, "LYP")
    LWP = read(file, "LWP")
    LXP = read(file, "LXP")
    BP = read(file, "BP")
    GP = read(file, "GP")
    ZP = read(file, "ZP")
    MUP = read(file, "MUP")
    LAMWP = read(file, "LAMWP")
    LAMFP = read(file, "LAMFP")
    MONP = read(file, "MONP")
    ELLP = read(file, "ELLP")
    CP = read(file, "CP")
    RRP = read(file, "RRP")
    IIP = read(file, "IIP")
    TTP = read(file, "TTP")
    WP = read(file, "WP")
    HHP = read(file, "HHP")
    PIP = read(file, "PIP")
    PIWP = read(file, "PIWP")
    LAMP = read(file, "LAMP")
    YP = read(file, "YP")
    XP = read(file, "XP")
    MCP = read(file, "MCP")
    QP = read(file, "QP")
    RKP = read(file, "RKP")
    KF = read(file, "KF")
    KK = read(file, "KK")
    LII = read(file, "LII")
    LY = read(file, "LY")
    LW = read(file, "LW")
    LX = read(file, "LX")
    B = read(file, "B")
    G = read(file, "G")
    Z = read(file, "Z")
    MU = read(file, "MU")
    LAMW = read(file, "LAMW")
    LAMF = read(file, "LAMF")
    MON = read(file, "MON")
    ELL = read(file, "ELL")
    C = read(file, "C")
    RR = read(file, "RR")
    II = read(file, "II")
    TT = read(file, "TT")
    W = read(file, "W")
    HH = read(file, "HH")
    PI = read(file, "PI")
    PIW = read(file, "PIW")
    LAM = read(file, "LAM")
    Y = read(file, "Y")
    X = read(file, "X")
    MC = read(file, "MC")
    Q = read(file, "Q")
    RK = read(file, "RK")
    TG = read(file, "TG")
    TGP = read(file, "TGP")
    BG = read(file, "BG")
    BGP = read(file, "BGP")
    F1 = read(file, "F1")
    F2 = read(file, "F2")
    F5 = read(file, "F5")
    F6 = read(file, "F6")
    F7 = read(file, "F7")
    F8 = read(file, "F8")
    F9 = read(file, "F9")
    F10 = read(file, "F10")
    F11 = read(file, "F11")
    F12 = read(file, "F12")
    F13 = read(file, "F13")
    F14 = read(file, "F14")
    F15 = read(file, "F15")
    F16 = read(file, "F16")
    F17 = read(file, "F17")
    F18 = read(file, "F18")
    F20 = read(file, "F20")
    F24 = read(file, "F24")
    F31 = read(file, "F31")
    F32 = read(file, "F32")
    F33 = read(file, "F33")
    F34 = read(file, "F34")
    F35 = read(file, "F35")
    F36 = read(file, "F36")
    F37 = read(file, "F37")
    F38 = read(file, "F38")
    F39 = read(file, "F39")
    F40 = read(file, "F40")
    F41 = read(file, "F41")
    F42 = read(file, "F42")
    close(file)
    endo = DSGE.augment_model_states(m.endogenous_states_original,
                                     DSGE.n_model_states_original(m))

    eq = m.equilibrium_conditions
    @testset "Check indices" begin
        @testset "Function blocks which output a function" begin
            @test KFP == endo[:kf′_t]  ## lagged ell function
        end
        @testset "Endogenous scalar-valued states:" begin
            @test KKP   == first(endo[:k′_t])
            @test LIIP  == first(endo[:i′_t1])
            @test LYP   == first(endo[:y′_t1])
            @test LWP   == first(endo[:w′_t1])
            @test LXP   == first(endo[:I′_t1])
            @test BG == first(endo[:bg_t])
            @test BGP == first(endo[:bg′_t])
        end
        @testset "Exogenous scalar-valued states" begin
            @test BP    == first(endo[:b′_t])
            @test GP    == first(endo[:g′_t])
            @test ZP    == first(endo[:z′_t])
            @test MUP   == first(endo[:μ′_t])
            @test LAMWP == first(endo[:λ_w′_t])
            @test LAMFP == first(endo[:λ_f′_t])
            @test MONP  == first(endo[:rm′_t])
        end
        @testset "Function-valued jumps" begin
            @test ELLP  == endo[:l′_t]
        end
        @testset "Scalar-valued jumps" begin
            @test CP   == first(endo[:C′_t])
            @test RRP  == first(endo[:R′_t])
            @test IIP  == first(endo[:i′_t])
            @test TTP  == first(endo[:t′_t])
            @test WP   == first(endo[:w′_t])
            @test HHP  == first(endo[:L′_t])
            @test PIP  == first(endo[:π′_t])
            @test PIWP == first(endo[:π_w′_t])
            @test LAMP == first(endo[:margutil′_t])
            @test YP   == first(endo[:y′_t])
            @test XP   == first(endo[:I′_t])
            @test MCP  == first(endo[:mc′_t])
            @test QP   == first(endo[:Q′_t])
            @test RKP  == first(endo[:capreturn′_t])
            @test KF   == endo[:kf_t]
            @test KK   == first(endo[:k_t])
            @test LII  == first(endo[:i_t1])
            @test LY   == first(endo[:y_t1])
            @test LW   == first(endo[:w_t1])
            @test LX   == first(endo[:I_t1])
            @test B    == first(endo[:b_t])
            @test G    == first(endo[:g_t])
            @test Z    == first(endo[:z_t])
            @test MU   == first(endo[:μ_t])
            @test LAMW == first(endo[:λ_w_t])
            @test LAMF == first(endo[:λ_f_t])
            @test MON  == first(endo[:rm_t])
            @test ELL  == endo[:l_t]
            @test C    == first(endo[:C_t])
            @test RR   == first(endo[:R_t])
            @test II   == first(endo[:i_t])
            @test TT   == first(endo[:t_t])
            @test W    == first(endo[:w_t])
            @test HH   == first(endo[:L_t])
            @test PI   == first(endo[:π_t])
            @test PIW  == first(endo[:π_w_t])
            @test LAM  == first(endo[:margutil_t])
            @test Y    == first(endo[:y_t])
            @test X    == first(endo[:I_t])
            @test MC   == first(endo[:mc_t])
            @test Q    == first(endo[:Q_t])
            @test RK   == first(endo[:capreturn_t])
            @test TG   == first(endo[:tg_t])
            @test TGP  == first(endo[:tg′_t])
        end
        @testset "Function blocks" begin
            @test F1  == eq[:eq_euler]
            @test F2  == eq[:eq_kolmogorov_fwd]
        end
        @testset "Function blocks which map functions to scalars" begin
            @test F5  == eq[:eq_agg_consumption]
            @test F6  == eq[:eq_lambda]
        end
        @testset "Scalar blocks involving endogenous variables" begin
            @test F7  == eq[:eq_transfers]
            @test F8  == eq[:eq_investment]
            @test F9  == eq[:eq_tobin_q]
            @test F10 == eq[:eq_capital_accumulation]
            @test F11 == eq[:eq_wage_phillips]
            @test F12 == eq[:eq_price_phillips]
            @test F13 == eq[:eq_marginal_cost]
            @test F14 == eq[:eq_gdp]
            @test F15 == eq[:eq_optimal_kl]
            @test F16 == eq[:eq_taylor]
            @test F17 == eq[:eq_fisher]
            @test F18 == eq[:eq_nominal_wage_inflation]
        end
        @testset "Lagged variables" begin
            @test F20 == eq[:LI]
            @test F24 == eq[:LY]
            @test F31 == eq[:LW]
            @test F32 == eq[:LX]
        end
        @testset "Shocks" begin
            @test F33 == eq[:eq_b]
            @test F34 == eq[:eq_g]
            @test F35 == eq[:eq_z]
            @test F36 == eq[:eq_μ]
            @test F37 == eq[:eq_λ_w]
            @test F38 == eq[:eq_λ_f]
            @test F39 == eq[:eq_rm]
        end
        @testset "New equations" begin
            @test F40 == eq[:eq_fiscal_rule]
            @test F41 == eq[:eq_g_budget_constraint]
            @test F42 == eq[:eq_resource_constraint]
        end
    end

    @testset "Check Jacobian outputs" begin
        @testset "Euler Equation" begin
            @test saved_JJ[F1,ELLP] ≈ JJ[F1, ELLP]
            @test saved_JJ[F1,ZP] ≈ JJ[F1,ZP]
            @test saved_JJ[F1,WP] ≈ JJ[F1,WP]
            @test saved_JJ[F1,HHP] ≈ JJ[F1,HHP]
            @test saved_JJ[F1,TTP] ≈ JJ[F1,TTP]
            @test saved_JJ[F1,B] ≈ JJ[F1,B]
            @test saved_JJ[F1,ELL] ≈ JJ[F1,ELL]
            @test saved_JJ[F1,RR] ≈ JJ[F1,RR]
        end
        @testset "Check Kolmogorov Fwd Equation" begin
            @test saved_JJ[F2,KFP] ≈ JJ[F2,KFP]
            @test saved_JJ[F2,KF] ≈ JJ[F2,KF]
            @test saved_JJ[F2,ELL] ≈ JJ[F2,ELL]
            @test saved_JJ[F2,RR] ≈ JJ[F2,RR]
            @test saved_JJ[F2,Z] ≈ JJ[F2,Z]
            @test saved_JJ[F2,W] ≈ JJ[F2,W]
            @test saved_JJ[F2,HH] ≈ JJ[F2,HH]
            @test saved_JJ[F2,TT] ≈ JJ[F2,TT]
        end
        @testset "Check aggregate consumption" begin
            @test saved_JJ[F5,C]   ≈ JJ[F5,C]
            @test saved_JJ[F5,ELL] ≈ JJ[F5,ELL]
            @test saved_JJ[F5,KF]  ≈ JJ[F5,KF]
            @test saved_JJ[F5,Z]   ≈ JJ[F5,Z]
            @test saved_JJ[F5,W]   ≈ JJ[F5,W]
            @test saved_JJ[F5,HH]  ≈ JJ[F5,HH]
            @test saved_JJ[F5,TT]  ≈ JJ[F5,TT]
        end
        @testset "Check lambda = average marginal utility" begin
            @test saved_JJ[F6,LAM] ≈ JJ[F6,LAM]
            @test saved_JJ[F6,KF]  ≈ JJ[F6,KF]
            @test saved_JJ[F6,RR]  ≈ JJ[F6,RR]
            @test saved_JJ[F6,Z]   ≈ JJ[F6,Z]
            @test saved_JJ[F6,W]   ≈ JJ[F6,W]
            @test saved_JJ[F6,HH]  ≈ JJ[F6,HH]
            @test saved_JJ[F6,TT]  ≈ JJ[F6,TT]
            @test saved_JJ[F6,ELL] ≈ JJ[F6,ELL]
        end
        @testset "transfer equation" begin
            @test saved_JJ[F7,TT] ≈ JJ[F7,TT]
            @test saved_JJ[F7,RK] ≈ JJ[F7,RK]
            @test saved_JJ[F7,KK] ≈ JJ[F7,KK]
            @test saved_JJ[F7,Z]  ≈ JJ[F7,Z]
            @test saved_JJ[F7,X]  ≈ JJ[F7,X]
            @test saved_JJ[F7,MC] ≈ JJ[F7,MC]
            @test saved_JJ[F7,Y]  ≈ JJ[F7,Y]
            @test saved_JJ[F7,G]  ≈ JJ[F7,G]
        end
        @testset "investment" begin
            @test saved_JJ[F8,Q]  ≈ JJ[F8,Q]
            @test saved_JJ[F8,MU] ≈ JJ[F8,MU]
            @test saved_JJ[F8,XP] ≈ JJ[F8,XP]
            @test saved_JJ[F8,ZP] ≈ JJ[F8,ZP]
            @test saved_JJ[F8,X]  ≈ JJ[F8,X]
            @test saved_JJ[F8,LX] ≈ JJ[F8,LX]
            @test saved_JJ[F8,Z]  ≈ JJ[F8,Z]
        end
        @testset "tobins q" begin
            @test saved_JJ[F9,RR]  ≈ JJ[F9,RR]
            @test saved_JJ[F9,Q]   ≈ JJ[F9,Q]
            @test saved_JJ[F9,RKP] ≈ JJ[F9,RKP]
            @test saved_JJ[F9,QP]  ≈ JJ[F9,QP]
        end
        @testset "capital accumulation" begin
            @test saved_JJ[F10,KKP] ≈ JJ[F10,KKP]
            @test saved_JJ[F10,KK] ≈ JJ[F10,KK]
            @test saved_JJ[F10,Z] ≈ JJ[F10,Z]
            @test saved_JJ[F10,MU] ≈ JJ[F10,MU]
            @test saved_JJ[F10,X] ≈ JJ[F10,X]
        end
        @testset "wage phillips curve" begin
            @test saved_JJ[F11,PIW] ≈ JJ[F11,PIW]
            @test saved_JJ[F11,LAMW] ≈ JJ[F11,LAMW]
            @test saved_JJ[F11,HH] ≈ JJ[F11,HH]
            @test saved_JJ[F11,LAM] ≈ JJ[F11,LAM]
            @test saved_JJ[F11,W] ≈ JJ[F11,W]
            @test saved_JJ[F11,PIP] ≈ JJ[F11,PIP]
        end
        @testset "price phillips curve" begin
            @test saved_JJ[F12,PI] ≈ JJ[F12,PI]
            @test saved_JJ[F12,MC] ≈ JJ[F12,MC]
            @test saved_JJ[F12,LAMF] ≈ JJ[F12,LAMF]
            @test saved_JJ[F12,PIP] ≈ JJ[F12,PIP]
        end
        @testset "marginal cost" begin
            @test saved_JJ[F13,MC] ≈ JJ[F13,MC]
            @test saved_JJ[F13,W] ≈ JJ[F13,W]
            @test saved_JJ[F13,RK] ≈ JJ[F13,RK]
        end
        @testset "gdp" begin
            @test saved_JJ[F14,Y] ≈ JJ[F14,Y]
            @test saved_JJ[F14,Z] ≈ JJ[F14,Z]
            @test saved_JJ[F14,KK] ≈ JJ[F14,KK]
            @test saved_JJ[F14,HH] ≈ JJ[F14,HH]
        end
        @testset "optimal k/l ratio" begin
            @test saved_JJ[F15,RK] ≈ JJ[F15,RK]
            @test saved_JJ[F15,W] ≈ JJ[F15,W]
            @test saved_JJ[F15,HH] ≈ JJ[F15,HH]
            @test saved_JJ[F15,KK] ≈ JJ[F15,KK]
            @test saved_JJ[F15,Z] ≈ JJ[F15,Z]
        end
        @testset "taylor rule" begin
            @test saved_JJ[F16,II] ≈ JJ[F16,II]
            @test saved_JJ[F16,LII] ≈ JJ[F16,LII]
            @test saved_JJ[F16,PI] ≈ JJ[F16,PI]
            @test saved_JJ[F16,Y] ≈ JJ[F16,Y]
            @test saved_JJ[F16,LY] ≈ JJ[F16,LY]
            @test saved_JJ[F16,Z] ≈ JJ[F16,Z]
            @test saved_JJ[F16,MON] ≈ JJ[F16,MON]
        end
        @testset "fisher eqn" begin
            @test saved_JJ[F17,RR] ≈ JJ[F17,RR]
            @test saved_JJ[F17,PIP] ≈ JJ[F17,PIP]
            @test saved_JJ[F17,II] ≈ JJ[F17,II]
        end
       @testset "wage inflation" begin
           @test saved_JJ[F18,PIW] ≈ JJ[F18,PIW]
           @test saved_JJ[F18,PI] ≈ JJ[F18,PI]
           @test saved_JJ[F18,Z] ≈ JJ[F18,Z]
           @test saved_JJ[F18,W] ≈ JJ[F18,W]
           @test saved_JJ[F18,LW] ≈ JJ[F18,LW]
       end
       @testset "fiscal rule" begin
           @test saved_JJ[F40,TG] ≈ JJ[F40,TG]
           @test saved_JJ[F40,RR] ≈ JJ[F40,RR]
           @test saved_JJ[F40,BG] ≈ JJ[F40,BG]
           @test saved_JJ[F40,Z] ≈ JJ[F40,Z]
           @test saved_JJ[F40,Y] ≈ JJ[F40,Y]
           @test saved_JJ[F40,G] ≈ JJ[F40,G]
       end
       @testset "govt budget constraint" begin
           @test saved_JJ[F41,BGP] ≈ JJ[F41,BGP]
           @test saved_JJ[F41,RR] ≈ JJ[F41,RR]
           @test saved_JJ[F41,BG] ≈ JJ[F41,BG]
           @test saved_JJ[F41,Z] ≈ JJ[F41,Z]
           @test saved_JJ[F41,Y] ≈ JJ[F41,Y]
           @test saved_JJ[F41,G] ≈ JJ[F41,G]
           @test saved_JJ[F41,TG] ≈ JJ[F41,TG]
       end
       @testset "update lagged variables" begin
            @test saved_JJ[F20,LIIP] ≈ JJ[F20,LIIP]
            @test saved_JJ[F20,II] ≈ JJ[F20,II]
            @test saved_JJ[F24,LYP] ≈ JJ[F24,LYP]
            @test saved_JJ[F24,Y] ≈ JJ[F24,Y]
            @test saved_JJ[F31,LWP] ≈ JJ[F31,LWP]
            @test saved_JJ[F31,W] ≈ JJ[F31,W]
            @test saved_JJ[F32,LXP] ≈ JJ[F32,LXP]
            @test saved_JJ[F32,X] ≈ JJ[F32,X]
        end
        @testset "discount factor shock" begin
            @test saved_JJ[F33,BP] ≈ JJ[F33,BP]
            @test saved_JJ[F33,B] ≈ JJ[F33,B]
        end
        @testset "g/y shock" begin
            @test saved_JJ[F34,GP] ≈ JJ[F34,GP]
            @test saved_JJ[F34,G] ≈ JJ[F34,G]
        end
        @testset " tfp growth shock" begin
            @test saved_JJ[F35,ZP] ≈ JJ[F35,ZP]
            @test saved_JJ[F35,Z] ≈ JJ[F35,Z]
        end
        @testset "investment shock" begin
            @test saved_JJ[F36,MUP] ≈ JJ[F36,MUP]
            @test saved_JJ[F36,MU] ≈ JJ[F36,MU]
        end
        @testset "wage mkup shock" begin
            @test saved_JJ[F37,LAMWP] ≈ JJ[F37,LAMWP]
            @test saved_JJ[F37,LAMW] ≈ JJ[F37,LAMW]
        end
        @testset "price mkup shock" begin
            @test saved_JJ[F38,LAMFP] ≈ JJ[F38,LAMFP]
            @test saved_JJ[F38,LAMF] ≈ JJ[F38,LAMF]
        end
        @testset "monetary policy shock" begin
            @test saved_JJ[F39,MONP] ≈ JJ[F39,MONP]
            @test saved_JJ[F39,MON] ≈ JJ[F39,MON]
        end
    end
end

if check_solution
    if !check_steady_state
        steadystate!(m)
    end
    nx_save = get_setting(m, :nx)
    m.testing = false
    m <= Setting(:nx, nx_save)
    gx, hx, _ = klein(m)

    file = JLD2.jldopen("$path/reference/solve_reduce_ell.jld2", "r")
    saved_gx   = read(file, "gx")
    saved_hx   = read(file, "hx")
    close(file)

    @testset "Check solve outputs" begin
        @test saved_gx  ≈ gx
        @test saved_hx  ≈ hx
    end
end

if check_irfs
    steadystate!(m)
    file = jldopen("$path/reference/irfs_reduce_ell.jld2", "r")

    IRFkf = read(file, "IRFkf")
    IRFZ = read(file, "IRFZ")
    IRFsh = read(file, "IRFsh")
    IRFell = read(file, "IRFell")
    IRFR = read(file, "IRFR")
    IRFi = read(file, "IRFi")
    IRFW = read(file, "IRFW")
    IRFPI = read(file, "IRFPI")
    IRFT = read(file, "IRFT")
    IRFY = read(file, "IRFY")
    IRFPIW = read(file, "IRFPIW")
    IRFX = read(file, "IRFX")
    IRFK = read(file, "IRFK")
    IRFQ = read(file, "IRFQ")
    IRFRK = read(file, "IRFRK")
    IRFH = read(file, "IRFH")
    IRFc = read(file, "IRFc")
    IRFC = read(file, "IRFC")
    IRFMC = read(file, "IRFMC")
    IRFLAM = read(file, "IRFLAM")
    IRFm = read(file, "IRFm")
    close(file)

    sys = compute_system(m)
    states, obs, pseudo = impulse_responses(m, sys, use_augmented_states = true)
    endo = m.endogenous_states_original

    @testset "Check IRFs" begin
        #Last entry is 3 because it's the third shock, the Z shock
        @test IRFkf ≈ states[endo[:kf′_t], 1:20, 3]
        @test IRFsh ≈ vec(states[endo[:z′_t], 1:20, 3])
        @test IRFell ≈ states[endo[:l′_t], 1:20, 3]
        @test IRFR ≈ vec(states[endo[:R′_t], 1:20, 3])
        @test IRFi ≈ vec(states[endo[:i′_t], 1:20, 3])
        @test IRFW ≈ vec(states[endo[:w′_t], 1:20, 3])
        @test IRFPI ≈ vec(states[endo[:π′_t], 1:20, 3])
        @test IRFT ≈ vec(states[endo[:t′_t], 1:20, 3])
        @test IRFY ≈ vec(states[endo[:y′_t], 1:20, 3])
        @test IRFPIW ≈ vec(states[endo[:π_w′_t], 1:20, 3])
        @test IRFX ≈ vec(states[endo[:I′_t], 1:20, 3])
        @test IRFK ≈ vec(states[endo[:k′_t], 1:20, 3])
        @test IRFQ ≈ vec(states[endo[:Q′_t], 1:20, 3])
        @test IRFRK ≈ vec(states[endo[:capreturn′_t], 1:20, 3])
        @test IRFH ≈ vec(states[endo[:L′_t], 1:20, 3])
        #@test IRFC ≈ vec(obs[m.observables[:obs_consumption], 1:20, 3])
        @test IRFMC ≈ vec(states[endo[:mc′_t], 1:20, 3])
        @test IRFLAM ≈ vec(states[endo[:margutil′_t], 1:20, 3])
        #@test IRFm ≈ vec(states[endo[:margutil′_t], 1:20, 3])
    end
end

m <= Setting(:steady_state_only, false)

# Steady-state computation
if check_steady_state_calibrate || write_steady_state_calibrate
    Random.seed!(0)
    steadystate!(m)

    if write_steady_state_calibrate
        JLD2.jldopen("$path/reference/steady_state_calibration$(addtl_save_str).jld2", "w") do file
            file["c"]    = m[:cstar].value
            file["m"]    = m[:μstar].value
            file["mpc"]  = m[:mpc].value
            file["pc0"]  = m[:pc0].value
            file["beta"] = m[:βstar].value
        end
    end

    if check_steady_state_calibrate
        file = jldopen("$path/reference/steady_state_calibration$(addtl_save_str).jld2", "r")
        saved_c    = read(file, "c")
        saved_μ    = read(file, "m")
        saved_mpc  = read(file, "mpc")
        saved_β    = read(file, "beta")
        saved_pc0  = read(file, "pc0")
        close(file)

        @testset "Check steady state outputs" begin
            @test saved_c   ≈ m[:cstar].value
            @test saved_μ   ≈ m[:μstar].value
            @test saved_pc0 ≈ m[:pc0].value
            @test saved_mpc ≈ m[:mpc].value
            # Tolerance of convergence of β is 1e-5
            @test saved_β ≈ m[:βstar].value
        end
    end
end

if check_likelihood || write_likelihood
    data = Matrix{Float64}(readdlm("$path/../../../reference/YY.txt")')

    if write_likelihood
        println("Overwriting likelihood...")
        JLD2.jldopen("$path/reference/likelihood_reduce_ell$(addtl_save_str).jld2", true, true,
                     true, IOStream) do file
            file["likelihood"] = likelihood(m, data)
        end
    end
    if check_likelihood
        lik_save = load("$path/reference/likelihood_reduce_ell$(addtl_save_str).jld2", "likelihood")

        @testset "Checking Likelihood output" begin
            @test likelihood(m, data) ≈ lik_save
        end
    end
end


text:
using Test, NaturallyUnitful
using NaturallyUnitful: c, ħ, ħc, kB, ϵ0

@testset "tests" begin
    @test natural(1u"m") ≈ uconvert(u"eV^-1", 1u"m"/(ħc))

    @test unnatural(u"m", uconvert(u"eV^-1", 1u"m"/(ħc))) ≈ 1.0u"m"

    @test natural(1e8u"m/s") ≈ convert(Float64, 1e8u"m/s"/c)

    @test unnatural(u"m/s", convert(Float64, 1e8u"m/s"/c)) ≈ 1e8u"m/s"

    @test natural(1u"kg") ≈ uconvert(u"GeV", 1u"kg"*c^2)

    @test natural(1u"m", base=u"GeV") ≈ uconvert(u"GeV^-1", 1u"m"/(ħc))

    @test unnatural(u"K", 1u"eV") ≈ uconvert(u"K", 1u"eV"/kB)
end


text:
immutable Environment
    clock_resolution::Float64
    clock_cost::Float64
end

Base.print(env::Environment) = let res  = env.clock_resolution
                                   cost = env.clock_cost
	print("Environment(\n  clock_resolution:$res ns\n  clock_cost:$cost ns\n)")
end

const SEED = 100 :: Int64

function resolution(k::Integer)
    times = zeros(Float64, k)
    for i in 1:k
        times[i] = time_sec()
    end
    return nonzeros(times[2:end] - times[1:end-1])
end

function cost(timelimit::Float64)
    time_clock(1)
    (_, iters, elapsed) = run_for_atleast(0.25, SEED, time_clock)
    ntimes = ceil(timelimit / elapsed)
    [time_clock(iters) / iters  for _ in 1:ntimes]
end

function analyze_mean(sample, niter)
    u = mean(sample)
    @printf("mean is %s (%d iterations)\n", time_str(u), niter)
    return u
end

function analyze_mean(sample)
    analyze_mean(sample, length(sample))
end

function measure_environment()
    Base.gc()
    @printf("Warming up...\n")
    (_, seed, _) = run_for_atleast(0.25, SEED, resolution)
    Base.gc()
    @printf("Estimating clock resolution... ")
    (_, _, tdiffs) = run_for_atleast(0.25, seed, resolution)
    clock_res = analyze_mean(tdiffs)
    @printf("Estimating cost of clock call... ")
    time_limit = min(1.0e5 * clock_res, 1.0)
    clock_cost = analyze_mean(cost(time_limit), SEED)
    return Environment(clock_res, clock_cost)
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct ActiveDirectoryProperties <: SwaggerModel
    domainName::Any # spec type: Union{ Nothing, String } # spec name: domainName
    netBiosDomainName::Any # spec type: Union{ Nothing, String } # spec name: netBiosDomainName
    forestName::Any # spec type: Union{ Nothing, String } # spec name: forestName
    domainGuid::Any # spec type: Union{ Nothing, String } # spec name: domainGuid
    domainSid::Any # spec type: Union{ Nothing, String } # spec name: domainSid
    azureStorageSid::Any # spec type: Union{ Nothing, String } # spec name: azureStorageSid

    function ActiveDirectoryProperties(;domainName=nothing, netBiosDomainName=nothing, forestName=nothing, domainGuid=nothing, domainSid=nothing, azureStorageSid=nothing)
        o = new()
        validate_property(ActiveDirectoryProperties, Symbol("domainName"), domainName)
        setfield!(o, Symbol("domainName"), domainName)
        validate_property(ActiveDirectoryProperties, Symbol("netBiosDomainName"), netBiosDomainName)
        setfield!(o, Symbol("netBiosDomainName"), netBiosDomainName)
        validate_property(ActiveDirectoryProperties, Symbol("forestName"), forestName)
        setfield!(o, Symbol("forestName"), forestName)
        validate_property(ActiveDirectoryProperties, Symbol("domainGuid"), domainGuid)
        setfield!(o, Symbol("domainGuid"), domainGuid)
        validate_property(ActiveDirectoryProperties, Symbol("domainSid"), domainSid)
        setfield!(o, Symbol("domainSid"), domainSid)
        validate_property(ActiveDirectoryProperties, Symbol("azureStorageSid"), azureStorageSid)
        setfield!(o, Symbol("azureStorageSid"), azureStorageSid)
        o
    end
end # type ActiveDirectoryProperties

const _property_map_ActiveDirectoryProperties = Dict{Symbol,Symbol}(Symbol("domainName")=>Symbol("domainName"), Symbol("netBiosDomainName")=>Symbol("netBiosDomainName"), Symbol("forestName")=>Symbol("forestName"), Symbol("domainGuid")=>Symbol("domainGuid"), Symbol("domainSid")=>Symbol("domainSid"), Symbol("azureStorageSid")=>Symbol("azureStorageSid"))
const _property_types_ActiveDirectoryProperties = Dict{Symbol,String}(Symbol("domainName")=>"String", Symbol("netBiosDomainName")=>"String", Symbol("forestName")=>"String", Symbol("domainGuid")=>"String", Symbol("domainSid")=>"String", Symbol("azureStorageSid")=>"String")
Base.propertynames(::Type{ ActiveDirectoryProperties }) = collect(keys(_property_map_ActiveDirectoryProperties))
Swagger.property_type(::Type{ ActiveDirectoryProperties }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ActiveDirectoryProperties[name]))}
Swagger.field_name(::Type{ ActiveDirectoryProperties }, property_name::Symbol) =  _property_map_ActiveDirectoryProperties[property_name]

function check_required(o::ActiveDirectoryProperties)
    (getproperty(o, Symbol("domainName")) === nothing) && (return false)
    (getproperty(o, Symbol("netBiosDomainName")) === nothing) && (return false)
    (getproperty(o, Symbol("forestName")) === nothing) && (return false)
    (getproperty(o, Symbol("domainGuid")) === nothing) && (return false)
    (getproperty(o, Symbol("domainSid")) === nothing) && (return false)
    (getproperty(o, Symbol("azureStorageSid")) === nothing) && (return false)
    true
end

function validate_property(::Type{ ActiveDirectoryProperties }, name::Symbol, val)
end


text:
using MultivariateOrthogonalPolynomials, DifferentialEquations, Plots
pyplot() # pyplot supports disks

Z = Zernike(1)
W = Weighted(Z)
xy = axes(W,1)
x,y = first.(xy),last.(xy)
Δ = Z \ Laplacian(xy) * W
S = Z \ W

# initial condition is (1-r^2) * exp(-(x-0.1)^2 - (y-0.2)^2)

K = Block.(Base.OneTo(11))

Δₙ = Δ[K,K]
Sₙ = S[K,K]
Zₙ = Z[:,K]
Wₙ = W[:,K]
c₀ = Zₙ \ @.(exp(-(x-0.1)^2 - (y-0.2)^2))

diskheat(c, (Δₙ, Sₙ), t) = Sₙ \ (Δₙ * c)
u = solve(ODEProblem(diskheat, c₀, (0.,1.), (Δₙ, Sₙ)), Tsit5(), reltol=1e-8, abstol=1e-8)

surface(Wₙ * u(1.0))

diskheat(c, (Δₙ, Sₙ), t) = Δₙ * c
u = solve(ODEProblem(ODEFunction(diskheat; jac=(u, (Δₙ, Sₙ), t) -> Δₙ, mass_matrix=Sₙ), c₀, (0.,1.), (Δₙ, Sₙ)), reltol=1e-8, abstol=1e-8)