text:
# using Omega
# using OmegaTestModels
# import ForwardDiff
# using Flux
# import Tracker

# const ΩTs = [LinearΩ{Vector{Int}, UnitRange{Int64}, Vector{Float64}},
#              Omega.SimpleΩ{Vector{Int}, Float64},
#              LinearΩ{Vector{Int}, UnitRange{Int64}, Tracker.TrackedArray{Float64, 1, Array{Float64,1}}}]

# const gradalgs = [Omega.TrackerGrad,
#                   Omega.ForwardDiffGrad,
#                   # Omega.ZygoteGrad
#                   ]

# function testgrad(OT = Omega.SimpleΩ{Vector{Int}, Float64})
#   μ = uniform(0.0, 1.0)
#   x = normal(μ, 1.0) + normal(μ, 1.0)
#   y = logerr(x ==ₛ 1.0)
#   ω = OT()
#   y(ω)
#   Omega.lineargradient(y, ω, Omega.ForwardDiffGrad)
# end

# testgrad(LinearΩ{Vector{Int}, UnitRange{Int64}, Vector{Float64}})

# # Test them all
# istracked(::Type{LinearΩ{I, SEG, T}}) where {I, SEG, T <: Tracker.TrackedArray} = true
# istracked(ΩT) = false

# function testgrad2(; modelcond, ΩT, gradalg)
#   @show ΩT
#   println()
#   Omega.lineargradient(logerr(modelcond), ΩT(), gradalg)
# end

# # Filter to differentiable models
# models = filter(hascond ∧ isdiff, allmodels)

# # Dont use TrackerGrad on wrong Omegas, etc
# f1(m, ΩT, gradalg) = gradalg == Omega.TrackerGrad ?  istracked(ΩT) : true
# f2(m, ΩT, gradalg) = gradalg == Omega.ForwardDiffGrad ? !istracked(ΩT) : true

# function runtests()
#   for model in models, gradalg in gradalgs, ΩT in ΩTs
#     if f1(model, ΩT, gradalg) && f2(model, ΩT, gradalg)
#       @show model.name
#       @show ΩT
#       @show gradalg
#       testgrad2(; modelcond = model.cond, ΩT = ΩT, gradalg = gradalg)
#     end
#   end
# end

# runtests()

# function zygotetest()
#   x = normal(0, 1, (10,))
#   y = sum(x)
#   lineargradient(y, rand(defΩ()), ZygoteGrad)
# end

# zygotetest()

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
using CollisionDetection

"""
Returns a mesh on the same vertexbuffer as the input mesh. The submesh
will be a mesh of dimension k containing all the k-cells that are in mesh
and that fulfill the predicate pred.

`pred` is a function with signature `pred(cell) -> Bool` returning true if
the simplex is to be added to the submesh under construction.
"""
function submesh(pred, mesh::Mesh)

    kcells = similar(mesh.faces)

    j = 1
    for i in 1:length(kcells)

        kcell = mesh.faces[i]
        if pred(kcell)
            kcells[j] = kcell
            j += 1
        end

    end

    kcells = kcells[1:j-1]

    Mesh(mesh.vertices, kcells)
end

"""
Store the k-cells of a mesh in an octree.

    function octree{U,D,T}(mesh::Mesh{U,D,T}, kcells::Array{Int,2})
"""
function octree(mesh::AbstractMesh)

    nverts = dimension(mesh) + 1 # number of vertices per cell
    ncells = numcells(mesh)

    P = vertextype(mesh)
    T = coordtype(mesh)

    points = zeros(P, ncells)
    radii = zeros(T, ncells)

    for (i,cl) in enumerate(cells(mesh))
        ch = chart(mesh, cl)
        points[i] = cartesian(center(ch))
        radii[i] = maximum(norm(v-points[i]) for v in ch.vertices)
    end

    return Octree(points, radii)
end

"""
Returns the boundingbox of a patch in terms of its center and halfsize.

    function boundingbox{U,D,C,N,T}(p::Simplex{U,D,C,N,T}) -> center, halfsize
"""
function boundingbox(p::Simplex{U,D,C,N,T}) where {U,D,C,N,T}

    # ll = minimum(p.vertices)
    # ur = maximum(p.vertices)

    ll = first(p.vertices); for v in p.vertices; ll = min.(v, ll); end
    ur = first(p.vertices); for v in p.vertices; ll = max.(v, ll); end

    center = (ll + ur) / 2
    halfsize = maximum(ur - center)

    return center, halfsize
end

function boundingbox(V)

    # ll = minimum(V)
    # ur = maximum(V)

    ll = first(V); for v ∈ V; ll = min.(v, ll); end
    ur = first(V); for v ∈ V; ur = max.(v, ur); end

    c = (ll + ur) / 2
    s = maximum(ur - c)

    return c, s
end

boundingbox(p::SVector{N,T}) where {N,T<:Number} = p, zero(eltype(p))

"""
    overlap_gpredicate(γ::Mesh) -> (patch -> Bool)

Create a predicate that for a given patch determinees if it overlaps with the provided target mesh `γ`.
"""
function overlap_gpredicate(γ::AbstractMesh)

    if numcells(γ) == 0
        return simplex -> false
    end

    # build an octree with the k-cells of smallmesh
    tree = octree(γ)

    function pred(p1)

        # find all simplices of the small mesh that potentially
        # collide with this patch
        c1, s1 = boundingbox(p1)
        for box in boxes(tree, (c,s)->boxesoverlap(c,s,c1,s1))
            for i in box
                # p2 = chart(γ, γ.faces[i])
                p2 = chart(γ, cells(γ)[i])
                overlap(p1,p2) && return true
            end
        end

        return false
    end

    return pred
end

"""
Geometric predicate for determining in log(N) complexity if a the image of a chart is in
the closure of mesh `γ`.
"""
function inclosure_gpredicate(γ::AbstractMesh)

    if numcells(γ) == 0
        return simplex -> false
    end

    # build an octree with the k-cells of smallmesh
    tree = octree(γ)

    function pred(p1)

        # find all simplices of the small mesh that potentially
        # collide with this patch
        c1, s1 = boundingbox(p1)
        for box in boxes(tree, (c,s)->boxesoverlap(c,s,c1,s1))
            for i in box
                # p2 = simplex(vertices(γ, γ.faces[i]))
                p2 = chart(γ, cells(γ)[i])
                isinclosure(p2, p1) && return true
            end
        end

        return false
    end

    return pred
end

"""
Creates a predicate that can be used to check wheter an edge is interior to
a surface (true) or on its boundary (false). This predicate is based on combinatorics.
In particular it expects as argument a tuple of indices pointing into the vertex
buffer of `mesh`.
"""
function interior_tpredicate(mesh::Mesh)

    @assert dimension(mesh) == 2

    vtoc, vton = vertextocellmap(mesh)

    function pred(simplex::SVector{2,Int})
        v = simplex[1]
        n = vton[v]
        a = vtoc[v,1:n]

        v = simplex[2]
        n = vton[v]
        b = vtoc[v,1:n]

        ab = a ∩ b
        return length(ab) == 2
    end

    return pred
end

"""
Creates a predicate that can be used to check wheter a vertex is interior to
a surface (true) or on its boundary (false).
In particular it expects as argument an index pointing into the vertex buffer of `mesh`.
"""
function interior_vpredicate(mesh::AbstractMesh)
    #TODO: update to accept simplex as argument
    @assert dimension(mesh) == 2
    skel = skeleton(mesh,1)
    edges = cells(skel)
    vtoc,vton = vertextocellmap(skel)
    pred = interior_tpredicate(mesh)
    function vpred(vidx)
        for j = 1:vton[vidx]
            pred(edges[vtoc[vidx,j]]) ? nothing : return false;
        end
        return true
    end
    return vpred
end

"""
    submesh(selection, mesh)

Create a submesh from `mesh` comprising those elements that overlap with elements from `selection`. It is assumed that `selection` and `mesh` have the same dimension.
"""
function submesh(sm::Mesh, bm::Mesh)
    overlaps = overlap_gpredicate(sm)
    in_smallmesh = c -> overlaps(simplex(vertices(bm,c)))
    submesh(in_smallmesh, bm)
end

ext: jl


text:
export tree_greedy, MinSpaceOut, MinSpaceDiff
export timespace_complexity

struct MinSpaceOut end
struct MinSpaceDiff end

struct LegInfo{ET}
    l1::Vector{ET}
    l2::Vector{ET}
    l12::Vector{ET}
    l01::Vector{ET}
    l02::Vector{ET}
    l012::Vector{ET}
end

"""
    tree_greedy(incidence_list, log2_sizes)

Compute greedy order, and the time and space complexities, the rows of the `incidence_list` are vertices and columns are edges.
`log2_sizes` are defined on edges.
"""
function tree_greedy(incidence_list::IncidenceList{VT,ET}, log2_edge_sizes; method=MinSpaceOut()) where {VT,ET}
    incidence_list = copy(incidence_list)
    n = nv(incidence_list)
    #tree_dict = collect(Any, 1:n)
    log2_tcs = Float64[] # time complexity
    log2_scs = Float64[]

    tree = Dict{VT,Any}([v=>v for v in vertices(incidence_list)])
    while true
        cost_values = evaluate_costs(method, incidence_list, log2_edge_sizes)
        local pair
        vmin = Inf
        for (k,v) in cost_values
            if v < vmin
                pair = k
                vmin = v
            end
        end
        log2_tc_step, sc, code = contract_pair!(incidence_list, pair..., log2_edge_sizes)
        push!(log2_tcs, log2_tc_step)
        push!(log2_scs, space_complexity(incidence_list, log2_edge_sizes))
        if nv(incidence_list) > 1
            tree[pair[1]] = ContractionTree(tree[pair[1]], tree[pair[2]])
        else
            return ContractionTree(tree[pair[1]], tree[pair[2]]), log2_tcs, log2_scs
        end
    end
end

function contract_pair!(incidence_list, vi, vj, log2_edge_sizes)
    log2dim(legs) = sum(l->log2_edge_sizes[l], legs, init=0)
    @assert vj > vi
    # compute time complexity and output tensor
    legsets = analyze_contraction(incidence_list, vi, vj)
    D12,D01,D02,D012 = log2dim.(getfield.(Ref(legsets),3:6))
    tc = D12+D01+D02+D012  # dangling legs D1 and D2 do not contribute

    # einsum code
    eout = legsets.l01 ∪ legsets.l02 ∪ legsets.l012
    code = (edges(incidence_list, vi), edges(incidence_list, vj)) => eout
    sc = log2dim(eout)

    # change incidence_list
    delete_vertex!(incidence_list, vj)
    change_edges!(incidence_list, vi, eout)
    for e in eout
        replace_vertex!(incidence_list, e, vj=>vi)
    end
    remove_edges!(incidence_list, legsets.l1 ∪ legsets.l2 ∪ legsets.l12)
    return tc, sc, code
end

function evaluate_costs(method, incidence_list::IncidenceList{VT,ET}, log2_edge_sizes) where {VT,ET}
    # initialize cost values
    cost_values = Dict{Tuple{VT,VT},Float64}()
    for vi = vertices(incidence_list)
        for vj in neighbors(incidence_list, vi)
            if vj > vi
                cost_values[(vi,vj)] = greedy_loss(method, incidence_list, log2_edge_sizes, vi, vj)
            end
        end
    end
    return cost_values
end

function analyze_contraction(incidence_list::IncidenceList{VT,ET}, vi::VT, vj::VT) where {VT,ET}
    ei = edges(incidence_list, vi)
    ej = edges(incidence_list, vj)
    leg012,leg12,leg1,leg2,leg01,leg02 = ET[], ET[], ET[], ET[], ET[], ET[]
    # external legs
    for leg in ei ∪ ej
        isext = leg ∈ incidence_list.openedges || !all(x->x==vi || x==vj, vertices(incidence_list, leg))
        if isext
            if leg ∈ ei
                if leg ∈ ej
                    push!(leg012, leg)
                else
                    push!(leg01, leg)
                end
            else
                push!(leg02, leg)
            end
        else
            if leg ∈ ei
                if leg ∈ ej
                    push!(leg12, leg)
                else
                    push!(leg1, leg)
                end
            else
                push!(leg2, leg)
            end
        end
    end
    return LegInfo(leg1, leg2, leg12, leg01, leg02, leg012)
end

function greedy_loss(::MinSpaceOut, incidence_list, log2_edge_sizes, vi, vj)
    log2dim(legs) = sum(l->log2_edge_sizes[l], legs, init=0)
    legs = analyze_contraction(incidence_list, vi, vj)
    log2dim(legs.l01)+log2dim(legs.l02)+log2dim(legs.l012) - 0.01*(log2dim(legs.l12)+log2dim(legs.l1)+log2dim(legs.l2))  # only counts external legs
end

function greedy_loss(::MinSpaceDiff, incidence_list, log2_edge_sizes, vi, vj)
    log2dim(legs) = sum(l->log2_edge_sizes[l], legs, init=0)
    legs = analyze_contraction(incidence_list, vi, vj)
    D1,D2,D12,D01,D02,D012 = log2dim.(getfield.(Ref(legs), 1:6))
    exp2(D01+D02+D012) - exp2(D1+D01+D12) - exp2(D2+D02+D12)  # out - in
end

function space_complexity(incidence_list, log2_sizes)
    sc = 0.0
    for v in vertices(incidence_list)
        for e in edges(incidence_list, v)
            sc += log2_sizes[e]
        end
    end
    return sc
end

function contract_tree!(incidence_list::IncidenceList, tree::ContractionTree, log2_edge_sizes, tcs, scs)
    vi = tree.left isa ContractionTree ? contract_tree!(incidence_list, tree.left, log2_edge_sizes, tcs, scs) : tree.left
    vj = tree.right isa ContractionTree ? contract_tree!(incidence_list, tree.right, log2_edge_sizes, tcs, scs) : tree.right
    tc, sc, code = contract_pair!(incidence_list, vi, vj, log2_edge_sizes)
    push!(tcs, tc)
    push!(scs, sc)
    return vi
end

function timespace_complexity(incidence_list::IncidenceList, tree::ContractionTree, log2_edge_sizes)
    tcs, scs = Float64[], Float64[]
    contract_tree!(copy(incidence_list), tree, log2_edge_sizes, tcs, scs)
    return log2sumexp2(tcs), maximum(scs)
end


text:
#=
startup.jl is a startup script which does several things:

1) Parses command line arguments
2) Loads parameters from given parameter file and constructs an array of parameters
3) Displays some messages

After running the script, the following variables are available:

const N_th::Int64 - the number of threads
const N_S::Int64 - number of aggregate states (parametrisations)
parfile::String - name of the parameter file
par::Array{pars{Float64,Int64},1} - parameters for each steady state
=#

#Get command line arguments, save them in parsed_args named tuple (global variable). Also save parfile - the name of the parameter file
using ArgParse #for parsing command line arguments.
#Load command line arguments. Move this to an included file later:
#Set up table of arguments passed by command line - see doc of ArgParse for details
s = ArgParseSettings()
@add_arg_table! s begin
    "--parfile"
        help = "parameter file name"
        arg_type = String
        default = ""
    "--numworkers"
        help = "number of workers"
        arg_type = Int
        default = 1
    "--slurm"
        help = "Use ClusterManagers for SLURM. 1 yes, 0 no (default)."
        arg_type = Int
        default = 0
end
#now get the parsed arguments
parsed_args = parse_args(ARGS, s)

#Copy the often used variables from the dictionary
parfile = parsed_args["parfile"]

#Intro message
const N_th = Threads.nthreads()
println("**************************************")
println("****** Brexit investment model *******")
println("**************************************\n")

println("Number of available threads: $N_th\n")
if N_th == 1
    println("Consider increasing the number of threads by 'export JULIA_NUM_THREADS=#'")
end

#Load parameters from file
if length(parfile) == 0
    println("No parameter file given. Using default parameter file baseline.jl.")
    parfile="baseline.jl"
end
include("../parameters/$parfile") #directory above
println("Parameters loaded from file parameters/$parfile\n")

#Number of aggregate states of the economy/parametrisations
const N_S = length(par_diff)

#Generate array par of parameters for different long-run equilibria
#using values of parameters from the parameter file.
expr = "par=["
for i=1:N_S
    global expr
   expr*="pars("
   if length(par_comm)>0 expr *= par_comm*"," end
   if length(par_diff)>0 expr *= par_diff[i] end
    expr*="),"
end
expr*="]"

#parse the string as a command and evaluate
expr2 = Meta.parse(expr)
eval(expr2)
expr = nothing;expr2 = nothing

#Check parameters (element-wise)
check_par.(par,N_S)

#Initialise results folder:
foldername = "results/"*genFolderName(parfile)
#Create results folder if it does not already exist
if !isdir("results")
    mkdir("results")
end
#Create the results/foldername folder. If it already exists, an error will occur
#(This could happen if the program is executed multiple times within a single second so the folder names are the same).
mkdir(foldername)
#Copy patemeter file in the results folder (for reproducibility)
cp("parameters/$parfile","$foldername/$parfile")

#Print a message about parameters:
println("*********Summary of parameters:*************")
println("All baseline parameters (pre-Brexit):")
@show par[1] #par[1] is the baseline (pre-Brexit)
println("_______________________")
println("Parameters that differ between stationary equilibria:")
@show par_diff
println(" ")
println("********************************************")


text:
function ComputeShiftold(s,l,m,n,a; ϵ = 10^(-10))
    #println("computing qnm")
    ψ = qnmfunctionnew(s,l,m,n,a)
    #println("Onto Defining Functions")
    F₊₁ = DoubleIntegrand1(ψ,δT,ψ)
    F₊₂ = DoubleIntegrand2(ψ,δT,ψ)
    Fₒ = CircleIntegrand(ψ,δT,ψ)
    G₊₁ = DoubleIntegrand1(ψ,∂ωT,ψ)
    G₊₂ = DoubleIntegrand2(ψ,∂ωT,ψ)
    Gₒ = CircleIntegrand(ψ,∂ωT,ψ)
    #println("THe value of F₊ is")
    #println(F₊(2.1,0.2))
    numerator1 = hcubature(t-> F₊₁(t[1],t[2]), [(0.0+ϵ),-(1.0-ϵ)], [(1.0-ϵ),(1.0-ϵ)])
    numerator2 = hcubature(t-> F₊₂(t[1],t[2]), [(0.0+ϵ),-(1.0-ϵ)], [(1.0-ϵ),(1.0-ϵ)])
    numerator3 = hcubature(t-> Fₒ(t[1],t[2]), [(0.0+ϵ),-(1.0-ϵ)], [(1.0-ϵ),(1.0-ϵ)])
    denomenator1 = hcubature(t-> G₊₁(t[1],t[2]), [(0.0+ϵ),-(1.0-ϵ)], [(1.0-ϵ),(1.0-ϵ)])
    denomenator2 = hcubature(t-> G₊₂(t[1],t[2]), [(0.0+ϵ),-(1.0-ϵ)], [(1.0-ϵ),(1.0-ϵ)])
    denomenator3 = hcubature(t-> Gₒ(t[1],t[2]), [(0.0+ϵ),-(1.0-ϵ)], [(1.0-ϵ),(1.0-ϵ)])
    δω =  -(numerator1[1]+numerator2[1]+numerator3[1])/(denomenator1[1]+denomenator2[1]+denomenator3[1])
    ψ.ω,δω
end

function Converter(f1s::Function)
function f2(x,v)
    kk = size(x)
    thesize = length(kk) == 1 ? 1 : kk[2]
    for i ∈ 1:(thesize)
        a = f1s(x[1,i],x[2,i])
        v[1,i] = real(a)
        v[2,i] = imag(a)
    end
end
f2
end

function ComputeShift(s,l,m,n,a,δT::Function,∂ωT::Function; ϵ = 10^(-10), abstol = 10^(-10))
    #println("computing qnm")
    println("newfunc")
    ψ = qnmfunctionnew(s,l,m,n,a)
    #println("Onto Defining Functions")
    #IntegralCHanger
    Δr = (ψ.R.r₊-ψ.R.r₋)
    wheretosendit = 0.8
    β = Δr*(1/wheretosendit - 1)
    F₊₁ = DoubleIntegrand1(ψ,δT,ψ; β=β)
    F₊₂ = DoubleIntegrand2(ψ,δT,ψ; β=β)
    Fₒ = CircleIntegrand(ψ,δT,ψ)
    G₊₁ = DoubleIntegrand1(ψ,∂ωT,ψ; β=β)
    G₊₂ = DoubleIntegrand2(ψ,∂ωT,ψ; β=β)
    Gₒ = CircleIntegrand(ψ,∂ωT,ψ)
    F₊₁(0.4,0.3) |> println
    F₊₂(0.4,0.3) |> println
    Fₒ(0.4,0.3) |> println
    G₊₁(0.4,0.3) |> println
    G₊₂(0.4,0.3) |> println
    Gₒ(0.4,0.3) |> println
    FF₊₁ = F₊₁ |> Converter;
    FF₊₂ = F₊₂ |> Converter;
    FFₒ = Fₒ  |> Converter;
    GG₊₁ = G₊₁ |> Converter;
    GG₊₂ = G₊₂ |> Converter;
    GGₒ = Gₒ |> Converter;
    #println("THe value of F₊ is
    #println(F₊(2.1,0.2))
    println("starting")
    numerator1 = pcubature_v(2, FF₊₁, [(0.0+ϵ),-(1.0-ϵ)], [(1.0-ϵ),(1.0-ϵ)]; error_norm = Cubature.PAIRED, abstol=abstol)
    print("First Done")
    numerator2 = pcubature_v(2, FF₊₂, [(0.0+ϵ),-(1.0-ϵ)], [(1.0-ϵ),(1.0-ϵ)]; error_norm = Cubature.PAIRED, abstol=abstol)
    print("Second Done")
    tovec = (x -> ([real(x[1]),imag(x[1])],x[2]))
    numerator3 = pcubature_v(2, FFₒ, [(0.0+ϵ),-(1.0-ϵ)], [(1.0-ϵ),(1.0-ϵ)]; error_norm = Cubature.PAIRED, abstol=abstol)
    print("3 Done")
    denomenator1 = pcubature_v(2, GG₊₁, [(0.0+ϵ),-(1.0-ϵ)], [(1.0-ϵ),(1.0-ϵ)]; error_norm = Cubature.PAIRED, abstol=abstol)
    print("4 Done")
    denomenator2 = pcubature_v(2, GG₊₂, [(0.0+ϵ),-(1.0-ϵ)], [(1.0-ϵ),(1.0-ϵ)]; error_norm = Cubature.PAIRED, abstol=abstol)
    print("5 Done")
    denomenator3 = pcubature_v(2, GGₒ, [(0.0+ϵ),-(1.0-ϵ)], [(1.0-ϵ),(1.0-ϵ)]; error_norm = Cubature.PAIRED, abstol=abstol)
    print("6 Done")
    tonum = (x -> x[1]+im*x[2])
    numer = tonum(numerator1[1]+numerator2[1]+numerator3[1]);
    δnumer = sum((numerator1[2]+numerator2[2]+numerator3[2]).^2) |> sqrt;
    denom = tonum(denomenator1[1]+denomenator2[1]+denomenator3[1]);
    δdenom = sum((denomenator1[2]+denomenator2[2]+denomenator3[2]).^2) |> sqrt;
    δω =  -tonum(numerator1[1]+numerator2[1]+numerator3[1])/tonum(denomenator1[1]+denomenator2[1]+denomenator3[1])
    δω_error = sqrt((abs(1/denom)*δnumer)^2 + (abs(numer/(denom^2))*δdenom)^2)
    ψ.ω,δω,δω_error
end

#NEW SPEEDUP TRY
function DoubleIntegrand1(ψL,δTs::Function,ψ; β=0.5)
    ψ1 = ∂r(ψ)
    ψ2 = ∂θ(ψ)
    ψ11 = ∂r(ψ1);
    ψ12 = ∂r(ψ2);
    ψ22 = ∂θ(ψ2);
    FF = let r₊ = ψ.R.r₊, r₋ = ψ.R.r₋, a = ψ.a, m = ψ.m ,ω = ψ.ω ,s = ψ.s, ψ=ψ, ψL=ψL, δTs=δTs, ψ11 = ψ11,ψ12 = ψ12,ψ22 = ψ22,ψ1 = ψ1,ψ2 = ψ2
    function F(r,z)
        δTs(r,z,a,m,ω,s,ψ11,ψ12,ψ22,ψ1,ψ2,ψ,ψL)*((r-r₊)*(r-r₋))^s
    end
    end
    Δr = (ψ.R.r₊ - ψ.R.r₋)*0.4
    let r₊ = ψ.R.r₊, Δr = Δr
        (t,z) -> β*im*((FF(r₊ - Δr*(im+1.0) + im*β*t/(1-t),z))/((1-t)^2))
    end
end

function DoubleIntegrand2(ψL,δTs::Function,ψ; β=0.5)
    ψ1 = ∂r(ψ)
    ψ2 = ∂θ(ψ)
    ψ11 = ∂r(ψ1);
    ψ12 = ∂r(ψ2);
    ψ22 = ∂θ(ψ2);
    FF = let r₊ = ψ.R.r₊, r₋ = ψ.R.r₋, a = ψ.a, m = ψ.m ,ω = ψ.ω ,s = ψ.s, ψ=ψ, ψL=ψL, δTs=δTs, ψ11 = ψ11,ψ12 = ψ12,ψ22 = ψ22,ψ1 = ψ1,ψ2 = ψ2
    function F(r,z)
        δTs(r,z,a,m,ω,s,ψ11,ψ12,ψ22,ψ1,ψ2,ψ,ψL)*((r-r₊)*(r-r₋))^s
    end
    end
    Δr = (ψ.R.r₊ - ψ.R.r₋)*0.4
    let r₊ = ψ.R.r₊, Δr = Δr
        (t,z) -> β*im*((-FF(r₊ - Δr*(im-1.0) + im*β*t/(1-t),z))/((1-t)^2))
    end
end

function CircleIntegrand(ψL,δTs::Function,ψ)
    ψ1 = ∂r(ψ)
    ψ2 = ∂θ(ψ)
    ψ11 = ∂r(ψ1);
    ψ12 = ∂r(ψ2);
    ψ22 = ∂θ(ψ2);
    FF = let r₊ = ψ.R.r₊, r₋ = ψ.R.r₋, a = ψ.a, m = ψ.m ,ω = ψ.ω ,s = ψ.s, ψ=ψ, ψL=ψL, δTs=δTs, ψ11 = ψ11,ψ12 = ψ12,ψ22 = ψ22,ψ1 = ψ1,ψ2 = ψ2
    function F(r,z)
        δTs(r,z,a,m,ω,s,ψ11,ψ12,ψ22,ψ1,ψ2,ψ,ψL)*((r-r₊)*(r-r₋))^s
    end
    end
    #(t,z) -> (-2*Δr*FF(r₊ + Δr*(2*t-(im+1.0)),z))
    Δr = (ψ.R.r₊ - ψ.R.r₋)*0.4
    let r₊ = ψ.R.r₊, Δr = Δr
        (t,z) -> (-2*Δr*FF(r₊ + Δr*(1.0-2.0*t-im),z))
    end
end

max_stars_repo_name: Potatoasad/KerrQNMShifts

max_issues_repo_name: Potatoasad/KerrQNMShifts

max_forks_repo_name: Potatoasad/KerrQNMShifts


text:
using NiLang
using Test

@testset "instructs.jl" begin
    include("instructs.jl")
end

@testset "functions.jl" begin
    include("functions.jl")
end

@testset "autodiff" begin
    include("autodiff/autodiff.jl")
end


text:
using RecipesBase

# Custom show methods
title(x::BalancedPanel{SingleUnitTreatment{Continuous}}) = "single treated unit, continuous treatment"
title(x::BalancedPanel{SingleUnitTreatment{Discontinuous}}) = "single treated unit, discontinuous treatment"
title(x::BalancedPanel{MultiUnitTreatment{Simultaneous{Continuous}}}) = "multiple treated units, simultaneous continuous treatment"
title(x::BalancedPanel{MultiUnitTreatment{Simultaneous{Discontinuous}}}) = "multiple treated units, simultaneous discontinuous treatment"
title(x::BalancedPanel{MultiUnitTreatment{Staggered{Continuous}}}) = "multiple treated units, staggered continuous treatment"
title(x::BalancedPanel{MultiUnitTreatment{Staggered{Discontinuous}}}) = "multiple treated units, staggered discontinuous treatment"

function Base.show(io::IO, mime::MIME"text/plain", x::BalancedPanel{SingleUnitTreatment{Continuous}})
    println("Balanced Panel - $(title(x))")
    println("    Treated unit: $(treated_labels(x))")
    println("    Number of untreated units: $(size(x.Y, 1) - 1)")
    println("    First treatment period: $(first_treated_period_labels(x))")
    println("    Number of pretreatment periods: $(length_T₀(x))")
    println("    Number of treatment periods: $(length_T₁(x))")
end

#!# TO DO - add show methods for other panel types

# Plotting recipe
@recipe function f(bp::BalancedPanel; kind = "treatment")

    xguide := "Period"
    #xticks --> bp.ts

    if kind == "treatment"

        legend --> :none
        yguide := "Unit"
        yticks --> (1:bp.N, bp.is)

        #@series begin
        #    seriestype := :heatmap
        #    x := bp.W
        #end

        for (i, r) ∈ enumerate(eachrow(bp.W))
            if sum(r) > 0
                @series begin
                    seriestype := :scatter
                    alpha --> [x ? 1.0 : 0.3 for x ∈ r]
                    markerstrokewidth --> 0.01
                    bp.ts, fill(i, length(r))
                end
            else
                @series begin
                    seriestype := :scatter
                    color --> "grey"
                    markerstrokewidth --> 0.01
                    bp.ts, fill(i, length(r))
                end
           end
       end

    elseif kind == "outcome"

        yguide := "Outcome value"
        legend --> :outertopright

        for i ∈ 1:bp.N

            if in(treated_ids(bp))(i)
                @series begin
                    label --> bp.is[i]
                    color --> i
                    bp.ts[1:findfirst(bp.W[i, :])], bp.Y[i, 1:findfirst(bp.W[i, :])]
                end

                @series begin
                    seriestype := :scatter
                    color --> i
                    markersize --> 5
                    markerstrokewidth --> 0.01
                    markershape --> :square
                    label --> ""
                    [bp.ts[findfirst(bp.W[i, :])]], [bp.Y[i, findfirst(bp.W[i, :])]]
                end

                @series begin
                    label --> ""
                    color --> i
                    alpha --> 0.7
                    #linestyle --> :dash
                    bp.ts[findfirst(bp.W[i, :]):end], bp.Y[i, findfirst(bp.W[i, :]):end]
                end
            else
                @series begin
                    label --> bp.is[i]
                    color --> i
                    alpha --> 0.3
                    bp.ts, bp.Y[i, :]
                end
            end
        end
    end

    nothing
end


text:
"""
    read_circ(file_input)

Takes circ file with name `file_input` and returns arrays
of radii, x, and y coordinates.
"""
function read_circ(file_input::String)
    # read in circ file
    circfile = string(file_input)
    f         = open(circfile,"r")
    lines     = readlines(f)
    num_lines = length(lines)
    close(f)

    # obtain pre-processing parameters (# bodies, radii, x/y centers)
    nbodies = parse(Int,lines[2])
    radiusArr = zeros(nbodies)
    xArr = zeros(nbodies)
    yArr = zeros(nbodies)

    tind = 1
    for ti = 4:3:num_lines-1
        radiusArr[tind] = parse(Float64,lines[ti])
        xArr[tind] = parse(Float64,lines[ti+1])
        yArr[tind] = parse(Float64,lines[ti+2])
        tind += 1
    end
    return radiusArr, xArr, yArr
end

"""
    write_circ(file_name,radiusArr,xArr,yArr)

Writes `circ` file (text file) with name `file_name`, using arrays
for radius, x and y positions
"""
function write_circ(file_name::String,radiusArr,xArr,yArr)
    filepath = string(file_name)
    open(filepath, "w") do io
        println(io, "# nbods")
        println(io, ps.param.n_objects)
        println(io, "# data below: radius, xc, yc for each body.")
        for ti=1:ps.param.n_objects
            println(io, radiusArr[ti])
            println(io, xArr[ti])
            println(io, yArr[ti])
        end
    end
    nothing
end
"""
    write_circ(file_name,ps)

Takes vector of particles `pList` and writes `circ` file
(text file) with name `file_name`
"""
function write_circ(file_name::String,pList)
    filepath = string(file_name)
    open(filepath, "w") do io
        println(io, "# nbods")
        println(io, ps.param.n_objects)
        println(io, "# data below: radius, xc, yc for each body.")
        for ti=1:ps.param.n_objects
            println(io, pList[ti].radius)
            println(io, pList[ti].pos.x)
            println(io, pList[ti].pos.y)
        end
    end
    nothing
end

"""
    read_in_circ2(file_name)
"""
function read_circ2(file_name::String)
    f         = open(file_name,"r")
    lines     = readlines(f)
    num_lines = length(lines)
    close(f)

    # obtain pre-processing parameters (# bodies, radii, x/y centers, materialID)
    nbodies   = parse(Int,lines[2])
    radiusArr = zeros(nbodies)
    xArr      = zeros(nbodies)
    yArr      = zeros(nbodies)
    matID     = zeros(Int,nbodies)

    tind = 1
    for ti = 4:4:num_lines-1
        radiusArr[tind] = parse(Float64,lines[ti])
        xArr[tind] = parse(Float64,lines[ti+1])
        yArr[tind] = parse(Float64,lines[ti+2])
        matID[tind] = parse(Int,lines[ti+3])
        tind += 1
    end

    return radiusArr, xArr, yArr, matID
end
"""
    write_circ2(file_name,pList)
"""
function write_circ2(file_name,radiusArr,xArr,yArr,matIDArr)
    filepath = string(file_name)
    open(filepath, "w") do io
        println(io, "# nbods")
        println(io, length(radiusArr))
        println(io, "# data below: radius, xc, yc, materialID for each body.")
        for ti=1:length(radiusArr)
            println(io, radiusArr[ti])
            println(io, xArr[ti])
            println(io, yArr[ti])
            println(io, matIDArr[ti])
        end
    end
end
function write_circ2(file_name,pList)
    filepath = string(file_name)
    open(filepath, "w") do io
        println(io, "# nbods")
        println(io, length(pList))
        println(io, "# data below: radius, xc, yc, materialID for each body.")
        for ti=1:length(pList)
            println(io, pList[ti].radius)
            println(io, pList[ti].pos.x)
            println(io, pList[ti].pos.y)
            println(io, pList[ti].materialID)
        end
    end
end

"""
    read_wall(file_name)
"""
function read_wall(file_name)
    f         = open(file_name,"r")
    lines     = readlines(f)
    num_lines = length(lines)
    close(f)

    nwalls  = parse(Int,lines[2])
    wList   = Array{AbstractWall}(undef,nwalls)

    tind = 1
    for ti = 4:8:num_lines-1
        linetype = lines[ti]
        x0 = parse(Float64,lines[ti+1])
        y0 = parse(Float64,lines[ti+2])
        if linetype == ":Line"
            x1 = parse(Float64,lines[ti+3])
            y1 = parse(Float64,lines[ti+4])
            thickness = parse(Float64,lines[ti+5])
            wList[tind] = LineWall(Point2D(x0,y0),Point2D(x1,y1),thickness)
        elseif linetype == ":Circle"
            radius = parse(Float64,lines[ti+3])
            thickness = parse(Float64,lines[ti+4])
            wList[tind] = CircleWall(Point2D(x0,y0),radius,thickness)
        elseif linetype == ":Arc"
            x1 = parse(Float64,lines[ti+3])
            y1 = parse(Float64,lines[ti+4])
            x2 = parse(Float64,lines[ti+5])
            y2 = parse(Float64,lines[ti+6])
            thickness = parse(Float64,lines[ti+7])
            wList[tind] = ArcWall(Point2D(x0,y0),Point2D(x1,y1),Point2D(x2,y2),thickness)
        else
            error("Linetype $(linetype) not recognized.")
        end
        tind +=1
    end
    return wList
end

"""
    write_wall(file_name,wList)
"""
function write_wall(file_name,wList)
    n_walls = length(wList)

    filepath = string(file_name)
    open(filepath, "w") do io
        println(io, "# walls")
        println(io, n_walls)
        println(io, "# data below: depends on linetype.")
        for ti=1:n_walls
            # print wall type
            if typeof(wList[ti]) <: LineWall
                println(io,":Line")
                println(io,wList[ti].nodes[1].x)
                println(io,wList[ti].nodes[1].y)
                println(io,wList[ti].nodes[2].x)
                println(io,wList[ti].nodes[2].y)
                println(io,wList[ti].thickness)
                println(io,"n/a")
                println(io,"n/a")
            elseif typeof(wList[ti]) <: CircleWall
                println(io,":Circle")
                println(io,wList[ti].center.x)
                println(io,wList[ti].center.y)
                println(io,wList[ti].radius)
                println(io,wList[ti].thickness)
                println(io,"n/a")
                println(io,"n/a")
                println(io,"n/a")
            elseif typeof(wList[ti]) <: ArcWall
                println(io,":Arc")
                println(io,wList[ti].nodes[1].x)
                println(io,wList[ti].nodes[1].y)
                println(io,wList[ti].nodes[2].x)
                println(io,wList[ti].nodes[2].y)
                println(io,wList[ti].nodes[3].x)
                println(io,wList[ti].nodes[3].y)
                println(io,wList[ti].thickness)
            end
        end
    end
end

"""
    combine_circ(output_name,input_files...)

Combines circ and circ2 files into one large circ (or circ2) file named
"output_name". If any of the input files are circ2, then all circ files
are converted to circ2 with materialID equal to zero.
"""
function combine_circ(output_name,input_files...)
    n_files = length(input_files)
    n_particles = 0

    # check if any extensions are circ2
    circ2_bool_arr = Array{Bool}(undef,n_files)
    for ti=1:n_files
        circ2_bool_arr[ti] = is_circ2_extension(input_files[ti])
    end
    any_circ2_bool = any(circ2_bool_arr)

    # get number of total points with 1st pass through
    for ti=1:n_files
        if !circ2_bool_arr[ti]
            tempradArr, _, _ = read_circ(input_files[ti])
        elseif circ2_bool_arr[ti]
            tempradArr, _, _, _ = read_circ2(input_files[ti])
        end
        n_particles += length(tempradArr)
    end

    # initialize full arrays
    radiusArr = zeros(n_particles)
    xArr      = zeros(n_particles)
    yArr      = zeros(n_particles)
    matIDArr  = zeros(Int,n_particles)

    # fill in arrays with 2nd pass through
    tk = 1
    for ti=1:n_files
        if !(circ2_bool_arr[ti])
            tempradArr, tempxArr, tempyArr = read_circ(input_files[ti])
        elseif circ2_bool_arr[ti]
            tempradArr, tempxArr, tempyArr, tempMatIDArr = read_circ2(input_files[ti])
        end
        n_particles = length(tempradArr)
        for tj=1:n_particles
            radiusArr[tk] = tempradArr[tj]
            xArr[tk] = tempxArr[tj]
            yArr[tk] = tempyArr[tj]
            if circ2_bool_arr[ti]
                matIDArr[tk] = tempMatIDArr[ti]
            end
            tk += 1
        end
    end

    if any_circ2_bool
        write_circ2(output_name,radiusArr,xArr,yArr,matIDArr)
    else
        write_circ(output_name,radiusArr,xArr,yArr)
    end

    nothing
end

function is_circ_extension(filename::String)
    _,ext = splitext(filename)
    return (ext == ".circ" ? true : false)
end

function is_circ2_extension(filename::String)
    _,ext = splitext(filename)
    return (ext == ".circ2" ? true : false)
end


text:
module BLAST

export
    blastn,
    blastp,
    readblastXML,
    BLASTResult

import Bio.Align: AlignedSequence
import Bio.Seq: BioSequence, DNASequence, AminoAcidSequence
import EzXML

include("blastcommandline.jl")

end # module BLAST


text:
# # Bound on Global Extremum

#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/generated/Polynomial Optimization/bound_on_blobal_extremum.ipynb)
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/generated/Polynomial Optimization/bound_on_global_extremum.ipynb)
# **Adapted from**: SOSTOOLS' SOSDEMO3 (See Section 4.3 of [SOSTOOLS User's Manual](http://sysos.eng.ox.ac.uk/sostools/sostools.pdf))

using Test #src
using DynamicPolynomials
@polyvar x1 x2

# The Goldstein-Price function $f(x)$ is defined as follows:

f1 = x1 + x2 + 1
f2 = 19 - 14x1 + 3x1^2 - 14x2 + 6x1*x2 + 3x2^2
f3 = 2x1 - 3x2
f4 = 18 - 32x1 + 12x1^2 + 48x2 - 36x1*x2 + 27x2^2
f = (1 + f1^2 * f2) * (30 + f3^2 * f4)

# We need to pick an SDP solver, see [here](https://jump.dev/JuMP.jl/v0.21.6/installation/#Supported-solvers) for a list of the available choices.
# We use `SOSModel` instead of `Model` to be able to use the `>=` syntax for Sum-of-Squares constraints.

using SumOfSquares
using CSDP
solver = optimizer_with_attributes(CSDP.Optimizer, MOI.Silent() => true)
model = SOSModel(solver);

# We create the decision variable $\gamma$ that will be the lower bound to the Goldstein-Price function.
# We maximize it to have the highest possible lower bound.

@variable(model, γ)
@objective(model, Max, γ)

# We constrain $\gamma$ to be a lower bound with the following constraint
# that ensures that $f(x_1, x_2) \ge \gamma$ for all $x_1, x_2$.

@constraint(model, f >= γ)

JuMP.optimize!(model)

# We verify that the solver has found a feasible solution:

JuMP.primal_status(model)
@test JuMP.primal_status(model) == MOI.FEASIBLE_POINT || JuMP.primal_status(model) == MOI.NEARLY_FEASIBLE_POINT #src

# We can now obtain the lower bound either with `value(γ)` or `objective_value(model)`:

objective_value(model)
@test objective_value(model) ≈ 3 rtol=1e-2 #src

max_stars_repo_name: blegat/SumOfSquares

max_issues_repo_name: blegat/SumOfSquares

max_forks_repo_name: blegat/SumOfSquares


text:
using Test
using ClimateMachine
using ClimateMachine.Mesh.Grids
using ClimateMachine.Mesh.Topologies
using ClimateMachine.MPIStateArrays
import ClimateMachine.Abstractions: DiscontinuousSpectralElementGrid
using Impero, Printf, MPI, LinearAlgebra

using ImperoHooks
#import ClimateMachine.Mesh.Grids: DiscontinuousSpectralElementGrid
ClimateMachine.init()
const ArrayType = ClimateMachine.array_type()
const mpicomm = MPI.COMM_WORLD
Ω = Circle(0,2π) × Circle(0,2π) # Ω = S¹(0,2π) × Interval(-1,1) × Interval(-2,2), Ω = Earth()

# error messages
DiscontinuousSpectralElementGrid(Ω)
DiscontinuousSpectralElementGrid(Ω, elements = (10,10,10))
DiscontinuousSpectralElementGrid(Ω, elements = (10,10,10), polynomialorder = (3,3,3))
DiscontinuousSpectralElementGrid(Ω, elements = (10,10), polynomialorder = (3,4))
DiscontinuousSpectralElementGrid(Ω×Ω, elements = (10,10,10,10), polynomialorder = (3,3,3,3))
# functional 2D
grid = DiscontinuousSpectralElementGrid(Ω, elements = (10,10), polynomialorder = (4,4), array = ArrayType)
# functional 3D
Ω = Circle(0,2π) × Circle(0,2π) × Interval(-1,1)
grid = DiscontinuousSpectralElementGrid(Ω, elements = (2,2,2), polynomialorder = (4,4,4), array = ArrayType)
##
@testset "domain test" begin
    # Previous
    dim = 3
    FT = Float64
    Nh = 3
    Nv = 3
    # Defining grid structure
    periodicity = ntuple(j -> false, dim)
    brickrange = (
        ntuple(
            j -> range(FT(-1); length = Nh + 1,
                stop = 1),
            dim - 1,
        )...,
        range(FT(-5); length = Nv + 1, stop = 5),
    )
    topl = StackedBrickTopology(
                                mpicomm,
                                brickrange;
                                periodicity = periodicity,
                                boundary = (
                                    ntuple(j -> (1, 2), dim - 1)...,
                                    (3, 4),
                                )
    )
    N = 1
    grid = DiscontinuousSpectralElementGrid(
        topl,
        FloatType = FT,
        DeviceArray = ArrayType,
        polynomialorder = N,
    )
    # Impero
    Ω = Interval(-1,1) × Interval(-1,1) × Interval(-5,5)
    igrid = DiscontinuousSpectralElementGrid(
        Ω,
        elements = (Nh,Nh,Nv),
        polynomialorder = (N,N,N),
        array = ArrayType,
        )
    @test norm(grid.D .- igrid.D) ≈ 0.0
    @test norm(grid.vgeo - igrid.vgeo) ≈ 0.0
    @test norm(grid.sgeo - igrid.sgeo) ≈ 0.0

    # Multiplying by 1.0 because norm(CuArray{Int}) gives error
    # TODO: update packages to fix this error
    @test norm((grid.elemtobndy - igrid.elemtobndy)*1.0) ≈ 0.0
    # Previous
    dim = 3
    FT = Float64
    Nh = 3
    Nv = 3
    # Defining grid structure
    periodicity = (false, true, true)
    brickrange = (
        ntuple(
            j -> range(FT(-1); length = Nh + 1,
                stop = 1),
            dim - 1,
        )...,
        range(FT(-5); length = Nv + 1, stop = 5),
    )
    topl = StackedBrickTopology(
                                mpicomm,
                                brickrange;
                                periodicity = periodicity,
                                boundary = (
                                    ntuple(j -> (1, 2), dim - 1)...,
                                    (3, 4),
                                )
    )
    N = 1
    grid = DiscontinuousSpectralElementGrid(
        topl,
        FloatType = FT,
        DeviceArray = ArrayType,
        polynomialorder = N,
    )
    # Impero
    Ω = Interval(-1,1) × Circle(-1,1) × Circle(-5,5)
    igrid = DiscontinuousSpectralElementGrid(
        Ω,
        elements = (Nh,Nh,Nv),
        polynomialorder = (N,N,N),
        array = ArrayType,
        )
    @test norm(grid.D .- igrid.D) ≈ 0.0
    @test norm(grid.vgeo - igrid.vgeo) ≈ 0.0
    @test norm(grid.sgeo - igrid.sgeo) ≈ 0.0

    # Multiplying by 1.0 because norm(CuArray{Int}) gives error
    # TODO: update packages to fix this error
    @test norm((grid.elemtobndy - igrid.elemtobndy)*1.0) ≈ 0.0
end

max_stars_repo_name: CliMA/ImperoHooks

max_issues_repo_name: CliMA/ImperoHooks

max_forks_repo_name: CliMA/ImperoHooks


text:
using TemporalGPs:
    dtc,
    dtcify,
    DTCSeparable,
    RectilinearGrid,
    RegularInTime,
    get_time,
    get_space,
    Separable,
    approx_posterior_marginals

@testset "pseudo_point" begin

    rng = MersenneTwister(123456)

    @testset "dtcify" begin
        z = randn(rng, 3)
        k_sep = Separable(SEKernel(), Matern32Kernel())
        @test dtcify(z, k_sep) isa DTCSeparable
        @test dtcify(z, 0.5 * k_sep) isa ScaledKernel{<:DTCSeparable}
        @test dtcify(z, transform(k_sep, 0.5)) isa TransformedKernel{<:DTCSeparable}
        @test dtcify(z, k_sep + k_sep) isa KernelSum{<:Tuple{DTCSeparable, DTCSeparable}}
    end

    # A couple of "base" kernels used as components in more complicated kernels below.
    separable_1 = Separable(SEKernel(), Matern12Kernel())
    separable_2 = Separable(SEKernel(), Matern52Kernel())

    # The various spatio-temporal kernels to try out.
    kernels = [

        (name="separable-1", val=separable_1),
        (name="separable-2", val=separable_2),

        (name="scaled-separable", val=0.5 * Separable(Matern52Kernel(), Matern32Kernel())),
        (
            name="stretched-separable",
            val=Separable(SEKernel(), transform(Matern12Kernel(), 1.3)),
        ),

        (name="sum-separable-1", val=separable_1 + separable_2),
        (name="sum-separable-2", val=1.3 * separable_1 + 0.95 * separable_2),
    ]

    # Input locations.
    xs = [
        (
            name="rectilinear",
            val=RectilinearGrid(randn(rng, 2), RegularSpacing(0.0, 0.3, 3)),
        ),
        (
            name="regular-in-time",
            val=RegularInTime(
                RegularSpacing(0.0, 0.1, 11),
                [randn(rng, 3) for _ in 1:11],
            ),
        ),
    ]

    # Spatial-locations of pseudo-inputs and predictions.

    z_r = randn(rng, 2)
    x_pr_r = randn(rng, 10)

    @testset "kernel=$(k.name), x=$(x.name)" for k in kernels, x in xs

        # Compute pseudo-input locations. These have to share time points with `x`.
        t = get_time(x.val)
        z = RectilinearGrid(z_r, t)
        z_naive = collect(z)

        # Construct naive GP.
        f_naive = GP(k.val)
        fx_naive = f_naive(collect(x.val), 0.1)
        y = rand(rng, fx_naive)

        # Construct state-space GP.
        f = to_sde(f_naive)
        fx = f(x.val, 0.1)

        # Verify dimensions of the LGSSM constructed to compute the DTC. This catches a
        # surprisingly large number of bugs during development.
        fx_dtc = TemporalGPs.dtcify(z_r, fx)
        lgssm = TemporalGPs.build_lgssm(fx_dtc)
        @test sum(map(dim_out, lgssm.emissions)) == length(y)
        validate_dims(lgssm)

        # The two approaches to DTC computation should be equivalent up to roundoff error.
        dtc_naive = dtc(fx_naive, y, f_naive(z_naive))
        dtc_sde = dtc(fx, y, z_r)
        @test dtc_naive ≈ dtc_sde rtol=1e-6

        elbo_naive = elbo(fx_naive, y, f_naive(z_naive))
        elbo_sde = elbo(fx, y, z_r)
        @test elbo_naive ≈ elbo_sde rtol=1e-6

        adjoint_test(
            (y, z_r) -> elbo(fx, y, z_r), (y, z_r);
            rtol=1e-7,
            context=Zygote.Context(),
            check_infers=false,
        )

        # Compute approximate posterior marginals naively.
        f_approx_post_naive = approx_posterior(VFE(), fx_naive, y, f_naive(z_naive))
        x_pr = RectilinearGrid(x_pr_r, get_time(x.val))
        naive_approx_post_marginals = marginals(f_approx_post_naive(collect(x_pr)))

        # This is a horrible interface, but it's the best that I can do on short notice.
        approx_post_marginals = approx_posterior_marginals(dtc, fx, y, z_r, x_pr_r)

        @test mean.(naive_approx_post_marginals) ≈ mean.(approx_post_marginals) rtol=1e-7
        @test std.(naive_approx_post_marginals) ≈ std.(approx_post_marginals) rtol=1e-7
        @show norm(std.(naive_approx_post_marginals) - std.(approx_post_marginals))

        # Similarly awful interface, make predictions for each point separately.
        approx_post_marginals_individual = map(eachindex(t)) do t
            approx_posterior_marginals(dtc, fx, y, z_r, x_pr_r, t)
        end

        approx_post_marginals_vec = reduce(vcat, approx_post_marginals_individual)

        @test mean.(approx_post_marginals) ≈ mean.(approx_post_marginals_vec) rtol=1e-7
        @test std.(approx_post_marginals) ≈ std.(approx_post_marginals_vec) rtol=1e-7

        # Do the RegularInTime one and compare it against RectilinearGrid.
        x_pr_rit = RegularInTime(get_time(x_pr), [get_space(x_pr) for _ in get_time(x.val)])
        approx_post_marginals_rit = approx_posterior_marginals(dtc, fx, y, z_r, x_pr_rit)

        @test mean.(approx_post_marginals) ≈ mean.(approx_post_marginals_rit) rtol=1e-7
        @test std.(approx_post_marginals) ≈ std.(approx_post_marginals_rit) rtol=1e-7

        @testset "missings" begin

            # Construct missing data.
            y_missing = Vector{Union{eltype(y), Missing}}(undef, size(y))
            y_missing .= y
            missing_idx = sort(randperm(rng, length(y))[1:Int(floor(length(y) / 3))])
            missing_idx = eachindex(y)
            y_missing[missing_idx] .= missing

            # Construct naive missing inputs and outputs.
            present_idx = setdiff(eachindex(y), missing_idx)
            naive_inputs_missings = collect(x.val)[present_idx]
            naive_y_missings = y[present_idx]
            fx_naive = f_naive(naive_inputs_missings, 0.1)

            # Compute DTC using both approaches.
            dtc_naive = dtc(fx_naive, naive_y_missings, f_naive(z_naive))
            dtc_sde = dtc(fx, y_missing, z_r)
            @test dtc_naive ≈ dtc_sde rtol=1e-7 atol=1e-7

            elbo_naive = elbo(fx_naive, naive_y_missings, f_naive(z_naive))
            elbo_sde = elbo(fx, y_missing, z_r)
            @test elbo_naive ≈ elbo_sde rtol=1e-7 atol=1e-7

            # Compute approximate posterior marginals naively with missings.
            f_approx_post_naive = approx_posterior(
                VFE(), fx_naive, naive_y_missings, f_naive(z_naive),
            )
            naive_approx_post_marginals = marginals(f_approx_post_naive(collect(x_pr)))

            # Compute approximate posterior marginals using the state-space approximation.
            approx_post_marginals = approx_posterior_marginals(
                dtc, fx, y_missing, z_r, x_pr_r,
            )

            @test mean.(naive_approx_post_marginals) ≈ mean.(approx_post_marginals) rtol=1e-7
            @test std.(naive_approx_post_marginals) ≈ std.(approx_post_marginals) rtol=1e-7
        end
    end
end


text:
module UnregisteredTests
import Unregistered
import Test

Unregistered.greet()

end


text:
"""
    CKANConnection(url::String)
    CKANConnection(url::String, apiversion::Int)

Creates a CKANConnection object with an API url,  by first checking that the
url is valid (i.e. return a status value of 200). An API version number (Int)
can be specified. If not specified, the "connection" will use the latest
version of the target API.
"""
struct CKANConnection
    url::String
    apiversion::Union{Int,Nothing}

    function CKANConnection(url::String; useversion::Bool=true)
        if useversion
            url = make_api_url(url, "")
        end
        test_url(url) ? new(url, nothing) : error("Invalid url")
    end

    function CKANConnection(url::String, apiversion::Int)
        url = make_api_url(url, apiversion)
        test_url(url) ? new(url, apiversion) : error("Invalid url")
    end

end

"""
    CKANResource(connection::CKANConnection, id::String, datasetid::String, metadata::Dict

Creates a CKANResource objetcs. Ressources are usually containes in CKAN datasets
(i.e. packages).
"""
struct CKANResource
    connection::CKANConnection
    id::String
    datasetid::String
    metadata::Dict
end

"""

    CKANDataset(connection::CKANConnection, id::String, metadata::Dict, resources::Union{Array{CKANResource}, Array{String}})

Creates a CKANDataset object. Datasets (i.e. packages) are the main objects in CKAN.
CKANDataset objects usually contains one or more CKANResource objects.
"""
struct CKANDataset
    connection::CKANConnection
    id::String
    metadata::Dict
    resources::Union{Array{CKANResource}, Array{String}}
end


text:
"""
    _get_warmup_points(cbm; random_objective=false, numstop=Inf)

Generate warmup points for all reactions in the model that are not fixed.
Assumes that the input JuMP model in `cbm` is already constrained.

The warmup points are sampled randomly from all possibilities until `numstop`
is reached; by default all points are generated.

By default, the warmup points are generated as in FVA by minimizing and
maximizing all reactions; `random_objective` switches this to completely random
objectives.
"""
function _get_warmup_points(cbm; random_objective = false, numstop = Inf)
    v = cbm[:x]
    ubs = cbm[:ubs]
    lbs = cbm[:lbs]
    # determine which rxns should be max/min-ized (non fixed fluxes)
    fixed_rxns = Int[]
    for i in eachindex(v)
        ub_val = normalized_rhs(ubs[i])
        lb_val = normalized_rhs(lbs[i])
        if ub_val >= 0 && ub_val == -lb_val
            push!(fixed_rxns, i)
        elseif ub_val < 0 && -ub_val == lb_val
            push!(fixed_rxns, i)
        end
    end

    # determine number of warmup points (if numstops is less than the number of non-fixed reactions.)
    var_rxn_inds = shuffle!(filter(x -> !(x in fixed_rxns), 1:length(v))) # shuffle the optimization points
    NN = numstop > length(var_rxn_inds) ? length(var_rxn_inds) : numstop
    wpoints = zeros(length(v), 2 * NN)

    # generate warmup points.
    for (i, ii) in zip(1:length(var_rxn_inds), 1:2:(2*length(var_rxn_inds)))
        i > NN && break

        if random_objective
            @objective(cbm, Max, sum(rand() * v[iii] for iii in var_rxn_inds))
        else
            @objective(cbm, Max, v[var_rxn_inds[i]])
        end

        optimize!(cbm)
        for j = 1:size(wpoints, 1)
            wpoints[j, ii] = value(v[j])
        end

        if random_objective
            @objective(cbm, Min, sum(rand() * v[iii] for iii in var_rxn_inds))
        else
            @objective(cbm, Min, v[var_rxn_inds[i]])
        end

        optimize!(cbm)
        for j = 1:size(wpoints, 1)
            wpoints[j, ii+1] = value(v[j])
        end
    end

    return wpoints
end


text:
# ------------------------------------------------------------------------------
# Preconditioner abstract base class
# ------------------------------------------------------------------------------

"""
Finite element operator preconditioner
"""
abstract type AbstractPreconditioner end

# ------------------------------------------------------------------------------


text:
#== # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Tests related to the position of the Sun.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
#
#   [1] Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
#       Microcosm Press, Hawthorn, CA, USA.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ==#

# File: ./src/sun/sun_position.jl
# ===============================

# Function sun_position_i
# -----------------------

################################################################################
#                                 Test Results
################################################################################
#
# Example 5-1: Finding the Sun position vector [1, p. 280]
#
#   Using JD_UTC ≈ JD_UT1 = 2453827.5, one gets
#
#           | 0.9771945 |
#   r_MOD = | 0.1924424 | AU
#           | 0.0834308 |
#
################################################################################

@testset "Function sun_position_i" begin
    JD_UTC = 2453827.5
    JD_UT1 = JD_UTC

    s_MOD = sun_position_i(JD_UT1)/au2m

    @test s_MOD[1] ≈ 0.9771945 atol = 2e-6
    @test s_MOD[2] ≈ 0.1924424 atol = 2e-6
    @test s_MOD[3] ≈ 0.0834308 atol = 2e-6
end


text:
# Interface for any diffusion model

"""
  Abstract type for diffusion tensors.  Used mostly for shock capturing.
  This type specifies the diffusion tensor.  It must implement the functions
  below (see the documentation for these specific function signatures
  for more detail:

  ```
    applyDiffusionTensor(obj::AbstractDiffusion, w::AbstractMatrix,
                    i::Integer, nodes::AbstractVector,
                    dx::Abstract3DArray, flux::Abstract3DArray)
  ```
"""
abstract type AbstractDiffusion end

#------------------------------------------------------------------------------
# Functions each diffusion model must implement
# Users should not generally call these functions, use the API functions below
"""
  Applies the diffusion tensor to a given array, for an element

  More specifically,

  ```
  for i in nodes
    for d=1:mesh.dim
      for d2=1:mesh.dim
        flux[:, i, d] = C[:, :, d1, d2]*dx[:, i, d2]
      end
    end
  end
  ```
  where C[:, :, :, :] is the diffusion tensor computed at state w[:, i].
  Note that only the nodes of the element specified by the `nodes` vector
  should be computed.  The values of `dx` not in `nodes` are undefined.

  **Inputs**

   * obj: the [`AbstractDiffusion`](@ref) object
   * sbp: `AbstractOperator`
   * params: `AbstractParamType`
   * q: the `numDofPerNode` x `numNodesPerElement` array of conservative
        variables for the element
   * w: the `numDofPerNode` x `numNodesPerElement` array of entropy variables
        for the element
   * coords: the `dim` x `numNodesPerElement` array of the xyz coordinates of
             the element nodes
   * dxidx: the `dim` x `dim` x `numNodesPerElement` array of the (scaled)
            mapping jacobian determinant
   * jac: the `numNodesPerElement` array of the mapping jacobian determinant.
   * i: the element number.  This is useful for diffusions where the coeffients
        are precomputed and stored in an array
   * nodes: a vector of integers describing which nodes of `flux` need to be
            populated.
   * dx: the values to multiply against, `numDofPerNode` x `numNodesPerElement`
         x `dim`

  **Inputs/Outputs**

   * flux: array to overwrite with the result, same size as `dx`

"""
function applyDiffusionTensor(obj::AbstractDiffusion, sbp::AbstractOperator,
                    params::AbstractParamType, q::AbstractMatrix,
                    w::AbstractMatrix,
                    coords::AbstractMatrix, dxidx::Abstract3DArray,
                    jac::AbstractVector,
                    i::Integer, nodes::AbstractVector,
                    dx::Abstract3DArray, flux::Abstract3DArray{Tres}
                   ) where {Tres}

  error("abstract fallback for applyDiffusionTensor() reached.  Did you forget to extend applyDiffusionTensor with a new method?")
end

"""
  More specifically,

  ```
  for i in nodes
    for d=1:mesh.dim
      for d2=1:mesh.dim
        # primal calculation
        # flux[:, i, d] = C[:, :, d1, d2]*dx[:, i, d2]
        dx_bar[:, i, d1] += C[:, :, d1, d2].'*flux_bar[:, i, d]
        C_bar[:, :, d1, d2] += flux_bar[:, i, d1]*dx[:, i, d1].'  # outer product
      end
    end
  end

  **Inputs**

   * obj: the [`AbstractDiffusion`](@ref) object
   * sbp: `AbstractOperator`
   * params: `AbstractParamType`
   * q: the `numDofPerNode` x `numNodesPerElement` array of conservative
        variables for the element
   * w: the `numDofPerNode` x `numNodesPerElement` array of entropy variables
        for the element
   * coords: the `dim` x `numNodesPerElement` array of the xyz coordinates of
             the element nodes
   * dxidx: the `dim` x `dim` x `numNodesPerElement` array of the (scaled)
            mapping jacobian determinant
   * jac: the `numNodesPerElement` array of the mapping jacobian determinant.
   * i: the element number.  This is useful for diffusions where the coeffients
        are precomputed and stored in an array
   * nodes: a vector of integers describing which nodes of `flux` need to be
            populated.
   * dx: the values to multiply against, `numDofPerNode` x `numNodesPerElement`
         x `dim`
   * flux_bar: adjoint part of `flux`, same size

  **Inputs/Outputs**

   * flux: array to overwrite with the result, same size as `dx`
   * q_bar: array to be updated with the result of back propgiation, same size
            as `q`
   * w_bar: array to be updated with the result of back propgiation, same size
            as `w`
   * dx_bar: array to be updated with result of back propigation, same size
             as `dx`

"""
function applyDiffusionTensor_revq(obj::AbstractDiffusion,
                    sbp::AbstractOperator,
                    params::AbstractParamType,
                    q::AbstractMatrix, q_bar::AbstractMatrix,
                    w::AbstractMatrix, w_bar::AbstractMatrix,
                    coords::AbstractMatrix, dxidx::Abstract3DArray,
                    jac::AbstractVector,
                    i::Integer, nodes::AbstractVector,
                    dx::Abstract3DArray, dx_bar::Abstract3DArray,
                    flux::Abstract3DArray, flux_bar::Abstract3DArray
                   )

  error("abstract fallback for applyDiffusionTensor_revq() reached.  Did you forget to extend applyDiffusionTensor_revq with a new method?")
end

"""
  Similar to [`applyDiffusionTensor_revq`](@ref), but computes the reverse mode
  with respect to the metrics rather than `q`.

  **Inputs**

   * obj: the [`AbstractDiffusion`](@ref) object
   * sbp: `AbstractOperator`
   * params: `AbstractParamType`
   * q: the `numDofPerNode` x `numNodesPerElement` array of conservative
        variables for the element
   * w: the `numDofPerNode` x `numNodesPerElement` array of entropy variables
        for the element
   * coords: the `dim` x `numNodesPerElement` array of the xyz coordinates of
             the element nodes
   * dxidx: the `dim` x `dim` x `numNodesPerElement` array of the (scaled)
            mapping jacobian determinant
   * jac: the `numNodesPerElement` array of the mapping jacobian determinant.
   * i: the element number.  This is useful for diffusions where the coeffients
        are precomputed and stored in an array
   * nodes: a vector of integers describing which nodes of `flux` need to be
            populated.
   * dx: the values to multiply against, `numDofPerNode` x `numNodesPerElement`
         x `dim`
   * flux_bar: adjoint part of `flux`, same size

  **Inputs/Outputs**

   * flux: array to overwrite with the result, same size as `dx`
   * dxidx_bar: array to update with result, same size as `dxidx`
   * jac_bar: array to update with result, same size as `jac`.
   * dx_bar: array to be updated with result of back propigation, same size
             as `dx`
"""
function applyDiffusionTensor_revm(obj::AbstractDiffusion,
                    sbp::AbstractOperator,
                    params::AbstractParamType,
                    q::AbstractMatrix, w::AbstractMatrix,
                    coords::AbstractMatrix,
                    dxidx::Abstract3DArray, dxidx_bar::Abstract3DArray,
                    jac::AbstractVector, jac_bar::AbstractVector,
                    i::Integer, nodes::AbstractVector,
                    dx::Abstract3DArray, dx_bar::Abstract3DArray,
                    flux::Abstract3DArray, flux_bar::Abstract3DArray
                   )

  error("abstract fallback for applyDiffusionTensor_revm() reached.  Did you forget to extend applyDiffusionTensor_revm with a new method?")
end

"""
  Differentiated version of [`ApplyDiffusionTensor`](@ref)

  **Inputs**

   * obj: the [`AbstractDiffusion`](@ref) object
   * sbp: `AbstractOperator`
   * params: `AbstractParamType`
   * q: the `numDofPerNode` x `numNodesPerElement` array of conservative
        variables for the element
   * w: the `numDofPerNode` x `numNodesPerElement` array of entropy variables
        for the element
   * coords: the `dim` x `numNodesPerElement` array of the xyz coordinates of
             the element nodes
   * dxidx: the `dim` x `dim` x `numNodesPerElement` array of the (scaled)
            mapping jacobian determinant
   * jac: the `numNodesPerElement` array of the mapping jacobian determinant.
   * i: the element number.  This is useful for diffusions where the coeffients
        are precomputed and stored in an array
   * nodes: vector of integers identifying  nodes at which to compute the
            output (4th dimension of `t1_dot` and `t2_dot`)
   * t1: the values to multiply against, `numDofPerNode` x `numNodesPerElement`
         x `dim`
   * t1_dot: the derivative of t1, `numDofPerNode` x `numDofPerNode` x `dim`
            x `numNodesPerElement` x `numNodesPerElement`

  **Inputs/Outputs**

   * t2_dot: the derivative of t1, `numDofPerNode` x `numDofPerNode` x `dim`
            x `numNodesPerElement` x `numNodesPerElement`, overwritten

"""
function applyDiffusionTensor_diff(obj::AbstractDiffusion,
                    sbp::AbstractOperator,
                    params::AbstractParamType,
                    q_el::AbstractMatrix, w_el::AbstractMatrix,
                    coords::AbstractMatrix, dxidx::Abstract3DArray,
                    jac::AbstractVector,
                    i::Integer, nodes::AbstractVector, t1::Abstract3DArray,
                    t1_dot::Abstract5DArray, t2_dot::Abstract5DArray)

  error("abstract fallback for applyDiffusionTensor_diff() reached.  Did you forget to extend applyDiffusionTensor_diff() with a new method?")

end

"""
  Method for the case where t1 is a function of 2 states, qL and qR, ie.

  t2_dotL = Lambda * t1_dotL + dLambda/dq * t1
  t2_dotR = Lambda * t1_dotR

  Note that the chain rule term dLambda/dq is added to t2_dotL, so if using
  this function for elementR of a given interface, the argument `t1_dotL`
  should really be the derivative of `t1` with respect to qR (and similarly
  `t2_dotL` and `t2_dotR` should be reversed.

  **Inputs**

   * obj: the [`AbstractDiffusion`](@ref) object
   * sbp: `AbstractOperator`
   * params: `AbstractParamType`
   * q: the `numDofPerNode` x `numNodesPerElement` array of conservative
        variables for the element
   * w: the `numDofPerNode` x `numNodesPerElement` array of entropy variables
        for the element
   * coords: the `dim` x `numNodesPerElement` array of the xyz coordinates of
             the element nodes
   * dxidx: the `dim` x `dim` x `numNodesPerElement` array of the (scaled)
            mapping jacobian determinant
   * jac: the `numNodesPerElement` array of the mapping jacobian determinant.
   * i: element number
   * nodes: vector of integers identifying  nodes at which to compute the
            output (4th dimension of `t1_dot` and `t2_dot`)
   * t1: the values Lambda is multiplied against for the primal operation,
         `numDofPerNode` x `numNodesPerElement` x `dim`
   * t1_dotL: 5D Jacobian containing the derivative of t1 wrt the left element
   * t1_dotR: 5D Jacobian containing the derivative of t1 wrt the right element

  **Inputs/Outputs**

   * t2_dotL: 5D Jacobian overwritten with result for elementL
   * t2_dotR: 5D Jacobian overwritten with result for elementR
"""
function applyDiffusionTensor_diff(obj::AbstractDiffusion, sbp::AbstractOperator,
                    params::AbstractParamType,
                    q_el::AbstractMatrix, w_el::AbstractMatrix,
                    coords::AbstractMatrix, dxidx::Abstract3DArray,
                    jac::AbstractVector,
                    i::Integer, nodes::AbstractVector, t1::Abstract3DArray,
                    t1_dotL::Abstract5DArray, t1_dotR::Abstract5DArray,
                    t2_dotL::Abstract5DArray, t2_dotR::Abstract5DArray)

  error("abstract fallback for applyDiffusionTensor_diff() reached.  Did you forget to extend applyDiffusionTensor_diff() with a new method?")
end

#------------------------------------------------------------------------------
# API functions, common to all diffusion models (implemented in terms of the
# above)

"""
  This function applies the diffusion tensor to only the nodes in the stencil of
  R for a given face.  This is useful when `flux` will be interpolated to the
  face in the next step, and can be significantly faster than applying the
  diffusion tensor to all nodes of the element.

  **Inputs**

   * obj: the [`AbstractDiffusion`](@ref) object
   * sbp: `AbstractOperator`
   * params: `AbstractParamType`
   * q: the `numDofPerNode` x `numNodesPerElement` array of conservative
        variables for the element
   * w: the `numDofPerNode` x `numNodesPerElement` array of entropy variables
        for the element
   * coords: the `dim` x `numNodesPerElement` array of the xyz coordinates of
             the element nodes
   * dxidx: the `dim` x `dim` x `numNodesPerElement` array of the (scaled)
            mapping jacobian determinant
   * jac: the `numNodesPerElement` array of the mapping jacobian determinant.
   * i: the element number.  This is useful for diffusions where the coeffients
        are precomputed and stored in an array
   * sbpface: an `AbstractFace` object that describes the stencil of R
   * face: integer describing which face of the element will be used.
   * dx: the values to multiply against, `numDofPerNode` x `numNodesPerElement`
         x `dim`

  **Inputs/Outputs**

   * flux: array to overwrite with the result, same size as `dx`

"""
function applyDiffusionTensor(obj::AbstractDiffusion, sbp::AbstractOperator,
                    params::AbstractParamType,
                    q::AbstractMatrix, w::AbstractMatrix,
                    coords::AbstractMatrix, dxidx::Abstract3DArray,
                    jac::AbstractVector,
                    i::Integer, sbpface::AbstractFace, face::Integer,
                    dx::Abstract3DArray, flux::Abstract3DArray)

  nodes = sview(sbpface.perm, :, face)
  applyDiffusionTensor(obj, sbp, params, q, w, coords, dxidx, jac, i, nodes,
                       dx, flux)

  return nothing
end

"""
  This method applies the diffusion tensor to all nodes of the element.

  **Inputs**

   * obj: the [`AbstractDiffusion`](@ref) object
   * sbp: `AbstractOperator`
   * params: `AbstractParamType`
   * q: the `numDofPerNode` x `numNodesPerElement` array of conservative
        variables for the element
   * w: the `numDofPerNode` x `numNodesPerElement` array of entropy variables
        for the element
   * coords: the `dim` x `numNodesPerElement` array of the xyz coordinates of
             the element nodes
   * dxidx: the `dim` x `dim` x `numNodesPerElement` array of the (scaled)
            mapping jacobian determinant
   * jac: the `numNodesPerElement` array of the mapping jacobian determinant.
   * i: the element number.  This is useful for diffusions where the coeffients
        are precomputed and stored in an array
   * dx: the values to multiply against, `numDofPerNode` x `numNodesPerElement`
         x `dim`

  **Inputs/Outputs**

   * flux: array to overwrite with the result, same size as `dx`
"""
function applyDiffusionTensor(obj::AbstractDiffusion, sbp::AbstractOperator,
                    params::AbstractParamType,
                    q::AbstractMatrix, w::AbstractMatrix,
                    coords::AbstractMatrix, dxidx::Abstract3DArray,
                    jac::AbstractVector,
                    i::Integer,
                    dx::Abstract3DArray, flux::Abstract3DArray)

  numNodesPerElement = size(dx, 2)
  nodes = 1:numNodesPerElement
  applyDiffusionTensor(obj, sbp, params, q, w, coords, dxidx, jac, i, nodes,
                       dx, flux)

  return nothing
end

#------------------------------------------------------------------------------
# diff

"""
  This method computes the Jacobian of applying the diffusion tensor to
  only the nodes on the stencil of R.

  **Inputs**

   * obj: the [`AbstractDiffusion`](@ref) object
   * sbp: `AbstractOperator`
   * params: `AbstractParamType`
   * q: the `numDofPerNode` x `numNodesPerElement` array of conservative
        variables for the element
   * w: the `numDofPerNode` x `numNodesPerElement` array of entropy variables
        for the element
   * coords: the `dim` x `numNodesPerElement` array of the xyz coordinates of
             the element nodes
   * dxidx: the `dim` x `dim` x `numNodesPerElement` array of the (scaled)
            mapping jacobian determinant
   * jac: the `numNodesPerElement` array of the mapping jacobian determinant.
   * i: the element number.  This is useful for diffusions where the coeffients
        are precomputed and stored in an array
   * sbpface: the `AbstractFace` that identifies the stencil of R
   * face: the local face number
   * t1: the values to multiply against, `numDofPerNode` x `numNodesPerElement`
         x `dim`
   * t1_dot: the derivative of t1, `numDofPerNode` x `numDofPerNode` x `dim`
            x `numNodesPerElement` x `numNodesPerElement`

  **Inputs/Outputs**

   * t2_dot: the derivative of t1, `numDofPerNode` x `numDofPerNode` x `dim`
            x `numNodesPerElement` x `numNodesPerElement`, overwritten
"""
function applyDiffusionTensor_diff(obj::AbstractDiffusion,
                    sbp::AbstractOperator, params::AbstractParamType,
                    q::AbstractMatrix, w::AbstractMatrix,
                    coords::AbstractMatrix, dxidx::Abstract3DArray,
                    jac::AbstractVector,
                    i::Integer, sbpface::AbstractFace, face::Integer,
                    t1::Abstract3DArray,
                    t1_dot::Abstract5DArray, t2_dot::Abstract5DArray)

  nodes = sview(sbpface.perm, :, face)
  applyDiffusionTensor_diff(obj, sbp, params, q, w, coords, dxidx, jac, i,
                            nodes, t1, t1_dot, t2_dot)

  return nothing
end

"""
  This method computes the Jacobian of applying the diffusion tensor to
  all nodes of the element

  **Inputs**

   * obj: the [`AbstractDiffusion`](@ref) object
   * sbp: `AbstractOperator`
   * params: `AbstractParamType`
   * q: the `numDofPerNode` x `numNodesPerElement` array of conservative
        variables for the element
   * w: the `numDofPerNode` x `numNodesPerElement` array of entropy variables
        for the element
   * coords: the `dim` x `numNodesPerElement` array of the xyz coordinates of
             the element nodes
   * dxidx: the `dim` x `dim` x `numNodesPerElement` array of the (scaled)
            mapping jacobian determinant
   * jac: the `numNodesPerElement` array of the mapping jacobian determinant.
   * i: the element number.  This is useful for diffusions where the coeffients
        are precomputed and stored in an array
   * t1: the values to multiply against, `numDofPerNode` x `numNodesPerElement`
         x `dim`
   * t1_dot: the derivative of t1, `numDofPerNode` x `numDofPerNode` x `dim`
            x `numNodesPerElement` x `numNodesPerElement`

  **Inputs/Outputs**

   * t2_dot: the derivative of t1, `numDofPerNode` x `numDofPerNode` x `dim`
            x `numNodesPerElement` x `numNodesPerElement`, overwritten

"""
function applyDiffusionTensor_diff(obj::AbstractDiffusion,
                    sbp::AbstractOperator, params::AbstractParamType,
                    q::AbstractMatrix, w::AbstractMatrix,
                    coords::AbstractMatrix, dxidx::Abstract3DArray,
                    jac::AbstractVector,
                    i::Integer,
                    t1::Abstract3DArray,
                    t1_dot::Abstract5DArray, t2_dot::Abstract5DArray)

  nodes = 1:size(t1, 2)
  applyDiffusionTensor_diff(obj, sbp, params, q,  w, coords, dxidx, jac, i,
                            nodes, t1, t1_dot, t2_dot)

  return nothing
end

"""
  Method for computing the Jacobian of applying the diffusion tensor to
  only the nodes in the stencil of R, for the case where `t1` depends
  on both `qL` and `qR`.

  **Inputs**

   * obj: the [`AbstractDiffusion`](@ref) object
   * sbp: `AbstractOperator`
   * params: `AbstractParamType`
   * q: the `numDofPerNode` x `numNodesPerElement` array of conservative
        variables for the element
   * w: the `numDofPerNode` x `numNodesPerElement` array of entropy variables
        for the element
   * coords: the `dim` x `numNodesPerElement` array of the xyz coordinates of
             the element nodes
   * dxidx: the `dim` x `dim` x `numNodesPerElement` array of the (scaled)
            mapping jacobian determinant
   * jac: the `numNodesPerElement` array of the mapping jacobian determinant.
   * i: element number
   * sbpface: the `AbstractFace` that describes the stencil of R
   * face: the local face number
   * t1: the values Lambda is multiplied against for the primal operation,
         `numDofPerNode` x `numNodesPerElement` x `dim`
   * t1_dotL: 5D Jacobian containing the derivative of t1 wrt the left element
   * t1_dotR: 5D Jacobian containing the derivative of t1 wrt the right element

  **Inputs/Outputs**

   * t2_dotL: 5D Jacobian overwritten with result for elementL
   * t2_dotR: 5D Jacobian overwritten with result for elementR
"""
function applyDiffusionTensor_diff(obj::AbstractDiffusion,
                    sbp::AbstractOperator, params::AbstractParamType,
                    q::AbstractMatrix, w::AbstractMatrix,
                    coords::AbstractMatrix, dxidx::Abstract3DArray,
                    jac::AbstractVector,
                    i::Integer, sbpface::AbstractFace, face::Integer,
                    t1::Abstract3DArray,
                    t1_dotL::Abstract5DArray, t1_dotR::Abstract5DArray,
                    t2_dotL::Abstract5DArray, t2_dotR::Abstract5DArray)

  nodes = 1:size(t1_dotL, 4)
  #nodes = sview(sbpface.perm, :, face)
  applyDiffusionTensor_diff(obj, sbp, params, q, w, coords, dxidx, jac, i,
                            nodes, t1, t1_dotL, t1_dotR, t2_dotL, t2_dotR)

  return nothing
end

"""
  Method for computing the Jacobian of applying the diffusion tensor to
  all the nodes of the element, for the case where `t1` depends
  on both `qL` and `qR`.

  **Inputs**

   * obj: the [`AbstractDiffusion`](@ref) object
   * sbp: `AbstractOperator`
   * params: `AbstractParamType`
   * q: the `numDofPerNode` x `numNodesPerElement` array of conservative
        variables for the element
   * w: the `numDofPerNode` x `numNodesPerElement` array of entropy variables
        for the element
   * coords: the `dim` x `numNodesPerElement` array of the xyz coordinates of
             the element nodes
   * dxidx: the `dim` x `dim` x `numNodesPerElement` array of the (scaled)
            mapping jacobian determinant
   * jac: the `numNodesPerElement` array of the mapping jacobian determinant.
   * i: element number
   * t1: the values Lambda is multiplied against for the primal operation,
         `numDofPerNode` x `numNodesPerElement` x `dim`
   * t1_dotL: 5D Jacobian containing the derivative of t1 wrt the left element
   * t1_dotR: 5D Jacobian containing the derivative of t1 wrt the right element

  **Inputs/Outputs**

   * t2_dotL: 5D Jacobian overwritten with result for elementL
   * t2_dotR: 5D Jacobian overwritten with result for elementR

"""
function applyDiffusionTensor_diff(obj::AbstractDiffusion,
                    sbp::AbstractOperator, params::AbstractParamType,
                    q::AbstractMatrix, w::AbstractMatrix,
                    coords::AbstractMatrix, dxidx::Abstract3DArray,
                    jac::AbstractVector,
                    i::Integer,
                    t1::Abstract3DArray,
                    t1_dotL::Abstract5DArray, t1_dotR::Abstract5DArray,
                    t2_dotL::Abstract5DArray, t2_dotR::Abstract5DArray)

  nodes = 1:size(t1_dotL, 4)
  applyDiffusionTensor_diff(obj, sbp, params, q, w, coords, dxidx, jac, i, nodes,
                            t1, t1_dotL, t1_dotR, t2_dotL, t2_dotR)

  return nothing
end

#------------------------------------------------------------------------------
# revq

function applyDiffusionTensor_revq(obj::AbstractDiffusion,
                    sbp::AbstractOperator,
                    params::AbstractParamType,
                    q::AbstractMatrix, q_bar::AbstractMatrix,
                    w::AbstractMatrix, w_bar::AbstractMatrix,
                    coords::AbstractMatrix, dxidx::Abstract3DArray,
                    jac::AbstractVector,
                    i::Integer,
                    dx::Abstract3DArray, dx_bar::Abstract3DArray,
                    flux::Abstract3DArray{Tres}, flux_bar::Abstract3DArray
                   ) where {Tres}

  nodes = 1:size(q, 2)
  applyDiffusionTensor_revq(obj, sbp, params, q, q_bar, w, w_bar, coords,
                            dxidx, jac, i, nodes, dx, dx_bar, flux, flux_bar)
end

function applyDiffusionTensor_revq(obj::AbstractDiffusion,
                    sbp::AbstractOperator,
                    params::AbstractParamType,
                    q::AbstractMatrix, q_bar::AbstractMatrix,
                    w::AbstractMatrix, w_bar::AbstractMatrix,
                    coords::AbstractMatrix, dxidx::Abstract3DArray,
                    jac::AbstractVector,
                    i::Integer, sbpface::AbstractFace, face::Integer,
                    dx::Abstract3DArray, dx_bar::Abstract3DArray,
                    flux::Abstract3DArray{Tres}, flux_bar::Abstract3DArray
                   ) where {Tres}

  nodes = sview(sbpface.perm, :, face)
  applyDiffusionTensor_revq(obj, sbp, params, q, q_bar, w, w_bar, coords,
                            dxidx, jac, i, nodes, dx, dx_bar, flux, flux_bar)

  return nothing
end

#------------------------------------------------------------------------------
# revm

function applyDiffusionTensor_revm(obj::AbstractDiffusion,
                    sbp::AbstractOperator,
                    params::AbstractParamType,
                    q::AbstractMatrix, w::AbstractMatrix,
                    coords::AbstractMatrix,
                    dxidx::Abstract3DArray, dxidx_bar::Abstract3DArray,
                    jac::AbstractVector, jac_bar::AbstractVector,
                    i::Integer,
                    dx::Abstract3DArray, dx_bar::Abstract3DArray,
                    flux::Abstract3DArray, flux_bar::Abstract3DArray
                   )

  nodes = 1:size(q, 2)
  applyDiffusionTensor_revm(obj, sbp, params, q, w, coords, dxidx, dxidx_bar,
                            jac, jac_bar, i, nodes, dx, dx_bar, flux, flux_bar)

  return nothing
end

function applyDiffusionTensor_revm(obj::AbstractDiffusion,
                    sbp::AbstractOperator,
                    params::AbstractParamType,
                    q::AbstractMatrix, w::AbstractMatrix,
                    coords::AbstractMatrix,
                    dxidx::Abstract3DArray, dxidx_bar::Abstract3DArray,
                    jac::AbstractVector, jac_bar::AbstractVector,
                    i::Integer, sbpface::AbstractFace, face::Integer,
                    dx::Abstract3DArray, dx_bar::Abstract3DArray,
                    flux::Abstract3DArray, flux_bar::Abstract3DArray
                   )

  nodes = sview(sbpface.perm, :, face)
  applyDiffusionTensor_revm(obj, sbp, params, q, w, coords, dxidx, dxidx_bar,
                            jac, jac_bar, i, nodes, dx, dx_bar, flux, flux_bar)

  return nothing
end


text:
# Separable sum, using tuples of arrays as variables

export SeparableSum

"""
**Separable sum of functions**

    SeparableSum(f_1, ..., f_k)

Given functions `f_1` to `f_k`, returns their separable sum, that is
```math
g(x_1, ..., x_k) = \\sum_{i=1}^k f_i(x_i).
```
The object `g` constructed in this way can be evaluated at `Tuple`s of length
`k`. Likewise, the `prox` and `prox!` methods for `g` operate with
(input and output) `Tuple`s of length `k`.

Example:

    f = SeparableSum(NormL1(), NuclearNorm()); # separable sum of two functions
    x = randn(10); # some random vector
    Y = randn(20, 30); # some random matrix
    f_xY = f((x, Y)); # evaluates f at (x, Y)
    (u, V), f_uV = prox(f, (x, Y), 1.3); # computes prox at (x, Y)
"""
struct SeparableSum{T <: Tuple} <: ProximableFunction
    fs::T
end

SeparableSum(fs::Vararg{ProximableFunction}) = SeparableSum((fs...,))

is_prox_accurate(f::SeparableSum) = all(is_prox_accurate.(f.fs))
is_convex(f::SeparableSum) = all(is_convex.(f.fs))
is_set(f::SeparableSum) = all(is_set.(f.fs))
is_singleton(f::SeparableSum) = all(is_singleton.(f.fs))
is_cone(f::SeparableSum) = all(is_cone.(f.fs))
is_affine(f::SeparableSum) = all(is_affine.(f.fs))
is_smooth(f::SeparableSum) = all(is_smooth.(f.fs))
is_quadratic(f::SeparableSum) = all(is_quadratic.(f.fs))
is_generalized_quadratic(f::SeparableSum) = all(is_generalized_quadratic.(f.fs))
is_strongly_convex(f::SeparableSum) = all(is_strongly_convex.(f.fs))

function (f::SeparableSum)(x::TupleOfArrays{R}) where R <: Real
    sum = R(0)
    for k in eachindex(x)
        sum += f.fs[k](x[k])
    end
    return sum
end

function prox!(ys::TupleOfArrays{R}, fs::Tuple, xs::TupleOfArrays{R}, gamma::R=R(1)) where R <: Real
    sum = R(0)
    for k in eachindex(xs)
        sum += prox!(ys[k], fs[k], xs[k], gamma)
    end
    return sum
end

function prox!(ys::TupleOfArrays{R}, fs::Tuple, xs::TupleOfArrays{R}, gamma::Tuple) where R <: Real
    sum = R(0)
    for k in eachindex(xs)
        sum += prox!(ys[k], fs[k], xs[k], gamma[k])
    end
    return sum
end

prox!(ys::TupleOfArrays{R}, f::SeparableSum, xs::TupleOfArrays{R}, gamma=R(1)) where R <: Real = prox!(ys, f.fs, xs, gamma)

function gradient!(grad::TupleOfArrays{R}, fs::Tuple, x::TupleOfArrays{R}) where R <: Real
    val = R(0)
    for k in eachindex(fs)
        val += gradient!(grad[k], fs[k], x[k])
    end
    return val
end

gradient!(grad::TupleOfArrays, f::SeparableSum, x::TupleOfArrays) = gradient!(grad, f.fs, x)

fun_name(f::SeparableSum) = "separable sum"
fun_dom(f::SeparableSum) = "n/a"
fun_expr(f::SeparableSum) = "(x₁, …, xₖ) ↦ f₁(x₁) + … + fₖ(xₖ)"
fun_params(f::SeparableSum) = "n/a"

function prox_naive(f::SeparableSum, xs::TupleOfArrays{R}, gamma::Union{R, Tuple}=R(1)) where R <: Real
    fys = R(0)
    ys = []
    for k in eachindex(xs)
        y, fy = prox_naive(f.fs[k], xs[k], typeof(gamma) <: Real ? gamma : gamma[k])
        fys += fy
        append!(ys, [y])
    end
    return Tuple(ys), fys
end


text:
module Trebekian

using Flux
using Random
Random.seed!(1234);

simple_rnn = Chain(Flux.RNN(1, 10), Dense(10, 1))

function loss(x, y)
  l = Flux.mse(simple_rnn.(x)[end], y)
  Flux.reset!(simple_rnn)
  return l
end

# generate some data for training
# from the shared dict on your computer
# and return a training and testing batch
# TODO: specify how much data to use for both training and testing
function generate_data()
	words = readlines("/usr/share/dict/words")
	lens = length.(words)
	words = shuffle!((v -> Float32.(v)).(collect.(words)))

	# split into train and test batches - just use the last 10 as
	# my test
	num_in_test_set = Int(round(length(words) / 10))
	num_in_training_set = length(words) - num_in_test_set

	start_index = 1
	break_index = start_index + num_in_training_set
	end_index = break_index + num_in_test_set

	words[start_index:break_index - 1], lens[start_index:break_index - 1], words[break_index:end_index - 1], lens[break_index:end_index - 1]
end

function train_simplernn(train_data, train_labels, test_data, test_labels)
	simple_rnn = Chain(Flux.RNN(1, 10), Dense(10, 1))

	ps = Flux.params(simple_rnn)

	# just use gradient descent
	opt = Flux.Descent(0.01)

	# eval callback function
	evalcb() = @show(sum(loss.(test_data, test_labels)))
	evalcb()

	Flux.train!(loss, ps, zip(train_data, train_labels), opt, cb = Flux.throttle(evalcb, 1))

	# print the training loss at the end
	evalcb()
end

end # module


text:
#
# Regression test
# (each julia script of test can also be run independently)
#

tests = [
  # Unit tests
  "unit_ahs.jl",
  "unit_wca.jl",
  "unit_lwca.jl",
  "unit_nfe.jl",
  "unit_nfetb.jl",
  "unit_hhtb.jl",

  # Investigational calculations
  "tm_optim.jl",
  "tm_binary.jl",
  # "tm_energy.jl",
  "tm_hhtb.jl"
]

for test in tests
  @time include(test)
end


text:
module TestStacks

using Test
import ThreadPools: QueuePool, LoggedQueuePool
using ThreadPools

include("util.jl")

function fn(x)
    x.data
end

@testset "stack tests" begin

    @testset "QueuePool foreground" begin

        objs = [TestObj(x) for x in 1:64]
        output = []
        pool = QueuePool(1)

        stack = Channel{TestObj}(1024) do stack
            for item in stack
                put!(fn, pool, item)
            end
            close(pool)
        end

        for item in objs
            put!(stack, item)
        end

        for result in poolresults(pool)
            push!(output, result)
            if result % 3 == 0
                put!(stack, TestObj(11))
            end
            if !isready(stack)
                close(stack)
            end
        end

        @test length(output) == 85

    end

    @testset "QueuePool background" begin

        objs = [TestObj(x) for x in 1:64]
        output = []
        pool = QueuePool(2)

        stack = Channel{TestObj}(1024) do stack
            for item in stack
                put!(fn, pool, item)
            end
            close(pool)
        end

        for item in objs
            put!(stack, item)
        end

        for result in poolresults(pool)
            push!(output, result)
            if result % 3 == 0
                put!(stack, TestObj(11))
            end
            if !isready(stack)
                close(stack)
            end
        end

        @test length(output) == 85

    end

    @testset "LoggedQueuePool foreground" begin

        objs = [TestObj(x) for x in 1:64]
        output = []
        pool = LoggedQueuePool(1)

        stack = Channel{TestObj}(1024) do stack
            for item in stack
                put!(fn, pool, item)
            end
            close(pool)
        end

        for item in objs
            put!(stack, item)
        end

        for result in poolresults(pool)
            push!(output, result)
            if result % 3 == 0
                put!(stack, TestObj(11))
            end
            if !isready(stack)
                close(stack)
            end
        end

        @test length(output) == 85
        @test length(collect(keys(pool.log))) == Threads.nthreads()
        @test sum(length, values(pool.log)) == 85

    end

    @testset "LoggedQueuePool background" begin

        objs = [TestObj(x) for x in 1:64]
        output = []
        pool = LoggedQueuePool(2)

        stack = Channel{TestObj}(1024) do stack
            for item in stack
                put!(fn, pool, item)
            end
            close(pool)
        end

        for item in objs
            put!(stack, item)
        end

        for result in poolresults(pool)
            push!(output, result)
            if result % 3 == 0
                put!(stack, TestObj(11))
            end
            if !isready(stack)
                close(stack)
            end
        end

        @test length(output) == 85
        @test length(collect(keys(pool.log))) == max(1, Threads.nthreads()-1)
        @test sum(length, values(pool.log)) == 85

    end

end

end # module


text:
#nb # %% A slide [markdown] {"slideshow": {"slide_type": "subslide"}}
# ## Symmetric Strip Grating
# This example consists of a symmetric strip grating, i.e. a grating where the strip width
# is half the unit cell period ``P``:
#
# ![diagram](./assets/symmetric_strip_diagram.png)
#
# Only three of the infinite number of strips in the grating are shown, and they extend infinitely to the left and right.
# The grating lies in the ``z=0`` plane with free space on both sides. The shaded areas represent metalization.
# The dashed lines show two possible choices for the unit cell location: "J" for a formulation in terms of electric
# surface currents, and "M" for magnetic surface currents.
#
# For normal incidence there is a closed-form solution due to Weinstein,
# but for a more recent reference one can find the solution in Problem 10.6 of R. E. Collin,
# *Field Theory of Guided Waves, Second Ed.*,
# IEEE Press, 1991.  Here is the code for computing the exact solution:
"""
    grating(kP, nterms=30) -> (Γ, T)

Compute the normal incidence reflecton and transmission coefficients of a symmetric grid of
zero-thickness conducting strips.  The product of the period of the strips and the incident
electric field wavenumber is `kP` (dimensionless).  The incident electric field is
perpendicular to the direction along the axis of the strips.  The series have been
accelerated by applying a Kummer's transformation, using the first two terms in the Maclaurin
series for the inverse sin function.  `kP` must be in the half-open interval [0,1). The
default number of summed terms `nterms` yields better than 10 digits of accuracy over the
interval [0.01,0.99].
"""
function grating(kP; nterms=30)
    sum1 = 1.3862943611198906 # \sum_{n=1}^{\infty} 1/(n-1/2) - 1/n = log(4)
    sum3 = 7.2123414189575710 # \sum_{n=1}^{\infty} (n-1/2)^{-3} - n^{-3} = 6 * \zeta(3)
    x = kP/(4π)
    θ = x*sum1 + x^3/6 * sum3
    for n = 1:nterms
        xonmhalf = x/(n - 0.5)
        xon = x/n
        term = asin(xonmhalf) - (xonmhalf + (xonmhalf)^3/6) -
              (asin(xon) - (xon + xon^3/6))
        θ += term
    end
    Γ = sin(θ) * cis(-π/2 - θ)
    T = 1 + Γ
    return (Γ, T)
end

#  Note that using the extension of
# [Babinet's Principle for electromagnetic fields](http://kirkmcd.princeton.edu/examples/babinet.pdf)
# this also provides the solution (upon appropriate interchange and sign change of the coefficients) for
# the case where the incident wave polarization is parallel to the direction of the strips.

# Here is the PSSFSS code to analyze this structure using electric currents as the unknowns.  We will
# scale the geometry so that the frequency in GHz is numerically equal to the period of the strips
# measured in wavelengths.

using Plots, PSSFSS
c = 11.802852677165355 # light speed [inch*GHz]
period = c  # so the period/wavelength = freq in GHz
Py = period
Ly = period/2
Px = Lx = Ly/10 # Infinite in x direction so this can be anything
Ny = 60
Nx = round(Int, Ny*Lx/Ly)
sheet = rectstrip(;Px, Py, Lx, Ly, Nx, Ny, units=inch)
flist = 0.02:0.02:0.98
steering = (θ=0, ϕ=0)
strata = [Layer()
          sheet
          Layer()]
results_j = analyze(strata, flist, steering, showprogress=false,
                    resultfile=devnull, logfile=devnull);
p1 = plot(sheet)
p2 = plot(sheet, unitcell=true)
ptitle = plot(title = "Symmetric Strip Triangulation",
             grid = false, showaxis = false, xtick=[], ytick=[],
             bottom_margin = -50Plots.px)
plot(ptitle, p1, p2, layout = @layout([A{0.09h}; [B C]]))
#md savefig("symstrip1.png"); nothing  # hide
#-
#md # ![](symstrip1.png)

# Note that setting `Lx = Px` causes the strip to fully occupy the x-extent
# of the unit cell.  PSSFSS automatically ensures that the triangle edges at these unit
# cell boundaries define basis functions that satisfy the Floquet (phase shift) boundary
# conditions, so that currents are free to flow across these unit cell boundaries.

# We can also analyze the same structure using magnetic currents in the areas free of
# metalization as the unknowns:
sheet = rectstrip(;class='M', Px, Py, Lx, Ly, Nx, Ny, units=inch)
strata = [Layer()
          sheet
          Layer()]
results_m = analyze(strata, flist, steering, showprogress=false,
                    resultfile=devnull, logfile=devnull);

#
# Each 50-frequency run of `analyze` takes about 14 seconds
# for this geometry of 720 triangles on my machine.
# More detailed timing information is available in the log file
# (which is omitted for generating this documentation).

# We will compare the PSSFSS results to the analytic solution:
## Generate exact results:
rt = grating.(2π*flist)
rperp_exact = first.(rt)
tperp_exact = last.(rt)
rpar_exact = -tperp_exact
tpar_exact = -rperp_exact;

# Obtain PSSFSS results for electric and magnetic currents:
outrequest = @outputs s11(v,v) s21(v,v) s11(h,h) s21(h,h)
rperp_j, tperp_j, rpar_j, tpar_j =
      collect.(eachcol(extract_result(results_j, outrequest)))
rperp_m, tperp_m, rpar_m, tpar_m =
      collect.(eachcol(extract_result(results_m, outrequest)));

# Generate the comparison plots:
angdeg(z) = rad2deg(angle(z)) # Convenience function

p1 = plot(title = "Perpendicular Reflection Magnitude",
          xlabel = "Period (wavelengths)",
          ylabel = "Coefficient Magnitude",
          legend=:topleft)
plot!(p1, flist, abs.(rperp_exact), ls=:dash, label="Exact")
plot!(p1, flist, abs.(rperp_j), label="PSSFSS J")
plot!(p1, flist, abs.(rperp_m), label="PSSFSS M")
#md savefig("symstrip2.png"); nothing  # hide
#-
#md # ![](symstrip2.png)
#-
p2 = plot(title = "Perpendicular Reflection Phase",
          xlabel = "Period (wavelengths)",
          ylabel = "Phase (deg)")
plot!(p2, flist, angdeg.(rperp_exact), ls=:dash, label="Exact")
plot!(p2, flist, angdeg.(rperp_j), label="PSSFSS J")
plot!(p2, flist, angdeg.(rperp_m), label="PSSFSS M")
#md savefig("symstrip3.png"); nothing  # hide
#-
#md # ![](symstrip3.png)
#-
p1 = plot(title = "Parallel Reflection Magnitude",
          xlabel = "Period (wavelengths)",
          ylabel = "Coefficient Magnitude")
plot!(p1, flist, abs.(rpar_exact), ls=:dash, label="Exact")
plot!(p1, flist, abs.(rpar_j), label="PSSFSS J")
plot!(p1, flist, abs.(rpar_m), label="PSSFSS M")
#md savefig("symstrip4.png"); nothing  # hide
#-
#md # ![](symstrip4.png)
#-
p2 = plot(title = "Parallel Reflection Phase",
          xlabel = "Period (wavelengths)",
          ylabel = "Phase (deg)")
plot!(p2, flist, angdeg.(rpar_exact), ls=:dash, label="Exact")
plot!(p2, flist, angdeg.(rpar_j), label="PSSFSS J")
plot!(p2, flist, angdeg.(rpar_m), label="PSSFSS M")
#md savefig("symstrip5.png"); nothing  # hide
#-
#md # ![](symstrip5.png)

# Now look at the transmission coefficients:
p1 = plot(title = "Perpendicular Transmission Magnitude",
          xlabel = "Period (wavelengths)",
          ylabel = "Coefficient Magnitude")
plot!(p1, flist, abs.(tperp_exact), ls=:dash, label="Exact")
plot!(p1, flist, abs.(tperp_j), label="PSSFSS J")
plot!(p1, flist, abs.(tperp_m), label="PSSFSS M")
#md savefig("symstrip6.png"); nothing  # hide
#-
#md # ![](symstrip6.png)
#-
p2 = plot(title = "Perpendicular Transmission Phase",
          xlabel = "Period (wavelengths)",
          ylabel = "Phase (deg)")
plot!(p2, flist, angdeg.(tperp_exact), ls=:dash, label="Exact")
plot!(p2, flist, angdeg.(tperp_j), label="PSSFSS J")
plot!(p2, flist, angdeg.(tperp_m), label="PSSFSS M")
#md savefig("symstrip7.png"); nothing  # hide
#-
#md # ![](symstrip7.png)
#-
p1 = plot(title = "Parallel Transmission Magnitude",
          xlabel = "Period (wavelengths)",
          ylabel = "Coefficient Magnitude", legend=:topleft)
plot!(p1, flist, abs.(tpar_exact), ls=:dash, label="Exact")
plot!(p1, flist, abs.(tpar_j), label="PSSFSS J")
plot!(p1, flist, abs.(tpar_m), label="PSSFSS M")
#md savefig("symstrip8.png"); nothing  # hide
#-
#md # ![](symstrip8.png)
#-
p2 = plot(title = "Parallel Transmission Phase",
          xlabel = "Period (wavelengths)",
          ylabel = "Phase (deg)")
plot!(p2, flist, angdeg.(tpar_exact), ls=:dash, label="Exact")
plot!(p2, flist, angdeg.(tpar_j), label="PSSFSS J")
plot!(p2, flist, angdeg.(tpar_m), label="PSSFSS M")
#md savefig("symstrip9.png"); nothing  # hide
#-
#md # ![](symstrip9.png)

# ### Conclusion
# Although good agreement is obtained, as expected the best agreement between
# all three results occurs for the lowest frequencies, where the triangles are
# smallest in terms of wavelength.  This emphasizes the fact that it is necessary for the
# user to check that enough triangles have been requested for good convergence
# over the frequency band of interest.  This example is an extremely demanding case
# in terms of bandwidth, as the ratio of maximum to minimum frequency here
# is ``0.98/0.02 = 49:1``