text:
"""
    statscholesky(xtx::Symmetric{T}, tol::Real=-1) where {T<:AbstractFloat}
Return a `CholeskyPivoted` object created from `xtx` where the pivoting scheme
retains the original order unless singularity is detected.  Columns that are
(computationally) linearly dependent on columns to their left are moved to the
right hand side in a left circular shift.
"""
function statscholesky(xtx::Symmetric{T}, tol::Real = -1) where {T<:AbstractFloat}
    n = size(xtx, 2)
    chpiv = cholesky(xtx, Val(true), tol = T(tol), check = false)
    chunp = cholesky(xtx, check = false);

    piv = LinearAlgebra.BlasInt.(collect(1:n))
    r = chpiv.rank

    if r < n

        k = chunp.info
        if k > r
            @warn """Fixed-effects matrix may be ill conditioned.
                     Left-circular shift may fail."""
        end

        nleft = n
        while nleft > r
            # the 0 lowerbound is for MKL compatibility
            # this arises when the unpivoted Cholesky succeeds but the pivoted
            # Chokesky estimates less than full rank (see #367)
            if 0 < k < nleft
                piv = piv[[1:k-1; k+1:n; k]]
                chunp = cholesky!(Symmetric(xtx[piv, piv]), check = false)
            end
            k = chunp.info
            nleft -= 1
        end

        for j = (r+1):n   # an MKL <-> OpenBLAS difference
            for i = (r+1):j
                chunp.factors[i, j] = zero(T)
            end
        end

    end

    CholeskyPivoted(chunp.factors, chunp.uplo, piv, r, tol, chpiv.info)
end

"""
    statsqr(x::Matrix{T}, ranktol::Real=1e-8) where {T<:AbstractFloat}
Return a `QRPivoted` object created from `x` where the pivoting scheme
retains the original order unless singularity is detected.
"""
function statsqr(x::Matrix{T}; ranktol=1e-8) where {T<:AbstractFloat}
    n = size(x, 2)

    qrpiv = qr(x, Val(true))
    rank =  searchsortedlast(abs.(diag(qrpiv.R)), ranktol, rev=true);
    if rank < n
        piv = qrpiv.p
    else
        piv = LinearAlgebra.BlasInt.(collect(1:n))
    end

    if rank < n && piv[1] ≠ 1
        # make sure the first column isn't moved
        # this is usually the intercept and it's desirable to avoid
        # pivoting the intercept out. note that this inflation may
        # not always be enough, but if it isn't, then you probably
        # have other problems

        x = deepcopy(x) # this is horrible
        x[:, 1] .*= 1e6
        qrpiv = qr(x, Val(true))
        piv = qrpiv.p
    end

    # preserve as much of the original column ordering as possible
    piv = [sort!(piv[1:rank]); sort!(piv[rank+1:n])];
    qrunp = qr(x[:, piv], Val(false));

    QRPivoted(qrunp.factors, qrpiv.τ, piv)
end


text:
using LinearAlgebra, StableRNGs, StatsModels, Test

import MixedModels: statsqr, statscholesky

xtx(X) = Symmetric(X'X, :U)  # creat the symmetric matrix X'X from X
LinearAlgebra.rank(F::LinearAlgebra.QRPivoted; tol=1e-8) = searchsortedlast(abs.(diag(F.R)), tol, rev=true);

const rng = StableRNG(4321234)

const simdat = (
    G = repeat('A':'T', inner=10),
    H = repeat('a':'e', inner=2, outer=20),
    U = repeat(0.:9, outer=20),
    V = repeat(-4.5:4.5, outer=20),
    Y = 0.1 * randn(rng, 200),
    Z = rand(rng, 200)
)

@testset "fullranknumeric" begin
    mm = modelmatrix(@formula(Y ~ 1 + U), simdat)
    st = statsqr(mm)
    @test st.p == 1:2
end

@testset "fullrankcategorical" begin
    mm = modelmatrix(@formula(Y ~ 1 + G*H), simdat)
    st = statsqr(mm)
    @test rank(st) == 100
    @test st.p == 1:100
end

@testset "dependentcolumn" begin
    mm = modelmatrix(@formula(Y ~ 1 + U + V + Z), simdat)
    st = statsqr(mm)
    perm = [1,2,4,3]
    @test rank(st) == 3
    @test st.p == perm
end

@testset "qr missing cells" begin
    mm = modelmatrix(@formula(Y ~ 1 + G*H), simdat)[5:end,:]
    st = statsqr(mm)
    @test rank(st) == 98
    # we no longer offer ordering guarantees besides preserving
    # relative order of linearly independent columns
    # and trying to keep the first column in the first position
    unpivoted = st.p[begin:rank(st)]
    @test unpivoted == sort(unpivoted)
end

@testset "cholesky missing cells" begin
    mm = modelmatrix(@formula(Y ~ 1 + G*H), simdat)
    # when a cell is missing, the indicator for it is always zero
    mm[:, 43] .= 0
    XtX = xtx(mm)
    ch = statscholesky(XtX)
    perm = [1:42; 44:100; 43]
    @test ch.rank == 99
    @test ch.piv == perm
    @test isapprox(xtx(ch.U), XtX[perm, perm], atol=0.00001)
end

ext: jl


text:
## ----------------------------------------------------------------------------
_load_ME_dat(iJR, method, exp) = ldat(iJR, "maxent_ep_dat", (;method, exp), ".jls")

## ----------------------------------------------------------------------------
# exchanges ave var correlations
function plot_ME_methods_comparizon()

    me_method1 = :ME_MAX_POL
    me_method2 = :ME_Z_EXPECTED_G_BOUNDED
    me_name1 = "ME^2"
    me_name2 = "ME^1"

    corr_av_p0 = plot(;
        # title = _textbf("Flux average"),
        xlabel = _textbf(me_name1, " flux ave."),
        ylabel = _textbf(me_name2, " flux ave."),
    )
    corr_va_p0 = plot(;
        # title = _textbf("Flux variance"),
        xlabel = _textbf(me_name1, " flux var. (log scale)"),
        ylabel = _textbf(me_name2, " flux var. (log scale)"),
    )
    corr_av_p, corr_va_p = deepcopy.([corr_av_p0, corr_va_p0])
    corr_ref_params = (;label = "", lw = 3, alpha = 0.8, ls = :dash)

    vas, avs = [], []
    for (iJR, Data) in [
        (ChK.iJR904, ChK.KayserData),
        (ChN.iJR904, ChN.NanchenData),
        (ChF.iJR904, ChF.FolsomData),
    ]
        src = string(nameof(Data))
        marker = (8, source_markers[Data])

        # load resume
        DAT = ChE.load_DAT(src)
        EXPS = DAT[:EXPS]

        lavs, lvas = [], []
        for exp in EXPS

            ## -------------------------------------------------------------
            # Common iders

            ME1_DAT = _load_ME_dat(iJR, me_method1, exp)
            model1 = ME1_DAT[:model]
            ME2_DAT = _load_ME_dat(iJR, me_method2, exp)
            model2 = ME2_DAT[:model]

            iders = intersect(model1.rxns, model2.rxns)

            ## -------------------------------------------------------------
            # method1
            epouts = ME1_DAT[:epouts]
            exp_beta = maximum(keys(epouts))
            epout = epouts[exp_beta]

            me1_av = ChU.av(model1, epout, iders)
            me1_va = ChU.va(model1, epout, iders)

            ## -------------------------------------------------------------
            # method2
            epouts = ME2_DAT[:epouts]
            exp_beta = maximum(keys(epouts))
            epout = epouts[exp_beta]

            me2_av = ChU.av(model2, epout, iders)
            me2_va = ChU.va(model2, epout, iders)

            ## -------------------------------------------------------------
            scatter!(corr_av_p, me1_av, me2_av; label = "", marker,
                color = :black, alpha = 0.5
            )

            scatter!(corr_va_p, log.(me1_va), log.(me2_va); label = "", marker,
                color = :black, alpha = 0.5
            )

            push!(lavs, me1_av..., me2_av...)
            push!(lvas, me1_va..., me2_va...)
            sort!(lavs); sort!(lvas)
        end

        ## -------------------------------------------------------------
        push!(avs, lavs...)
        push!(vas, lvas...)
        sort!(avs); sort!(vas)

    end

    plot!(corr_av_p, avs, avs; corr_ref_params...)
    plot!(corr_va_p, log.(vas), log.(vas); corr_ref_params...)

    # sfig(ChE, [corr_av_p, corr_va_p],
    #     @fileid, "corrs_me_max_vs_bounded", ".png",
    #     layout = (1, 2)
    # )

    return corr_av_p, corr_va_p
end

## ----------------------------------------------------------------------------
# glc, biom marginals
function plot_me_marginals()

    iJR = ChN.iJR904
    Data = ChN.NanchenData
    exp = 4
    src = nameof(Data)

    me_method1 = :ME_MAX_POL
    me_method2 = :ME_Z_EXPECTED_G_BOUNDED

    colors = Dict(
        me_method1 => :red,
        me_method2 => :blue,
    )

    labels = Dict(
        me_method1 => "ME^2",
        me_method2 => "ME^1",
    )

    rxns_map = iJR.load_rxns_map()

    @info("Doing", src, exp)

    ps = Plots.Plot[]
    for (sign, exp_ider, ylabel, xlim) in [
            (-1.0, "GLC", "u_g" , (0.8, 2.6)),
            (-1.0, "AC", "u_a", (-0.075, 0.015)),
            (1.0, "D", "z", (0.0, 0.2)),
        ]

        model_ider = rxns_map[exp_ider]

        p = plot(;
            xlabel = _textbf(ylabel),
            ylabel = _textbf("norm. pdf"),
        )

        for method in [me_method1, me_method2]

            ME_DAT = _load_ME_dat(iJR, method, exp)
            model = ME_DAT[:model]
            epouts = ME_DAT[:epouts]
            exp_beta = maximum(keys(epouts))
            epout = epouts[exp_beta]

            # uses absolute values

            lb_, ub_  = ChU.bounds(model, model_ider)
            @show model_ider, lb_, ub_
            @assert lb_ * ub_ >= 0 # both bounds with equal sign
            lb_, ub_ = sort!(sign .* [lb_, ub_])
            @show lb_, ub_

            μ_ = sign * ChU.μ(model, epout, model_ider)
            av_ = sign * ChU.av(model, epout, model_ider)
            σ_ = sqrt(ChU.σ(model, epout, model_ider))
            @show μ_, av_, σ_
            ChP.plot_marginal!(p, μ_, σ_, lb_, ub_;
                av = av_,
                color = colors[method],
                label = _textbf(labels[method]),
                xlim,
                lw = 5,
                legend=:topleft,
                normalize = true
            )

            println()
        end
        push!(ps, p)

    end # for exp_ider
    # sfig(ChE, ps,
    #     @fileid, "exchange_marginals", (;src, exp), ".png",
    #     layout = (3, 1)
    # )

    return ps

end


text:
import DifferentialEquations
using ReachabilityAnalysis: _solve_ensemble, ReachSolution

@testset "Simulation traces using DifferentialEquations.jl solvers" begin
    prob = @ivp(x' = 1.0x, x(0) ∈ [1/2])

    # `solve` extends CommonSolve.jl
    sol = DifferentialEquations.solve(prob, tspan=(0.0, 1.0))
    @test sol isa ReachSolution

    # ensemble problems API
    sol = _solve_ensemble(prob, tspan=(0.0, 1.0))
    @test sol isa DifferentialEquations.EnsembleSolution
    sol = _solve_ensemble(prob, T=1.0)
    @test sol isa DifferentialEquations.EnsembleSolution

    # inplace (default)
    sol = _solve_ensemble(prob, T=1.0, inplace=true)

    # out-of-place
    sol = _solve_ensemble(prob, T=1.0, inplace=false)
end


text:
#using .TensorTrains

"""
Example of the PPP Hamiltonian for N=6
E_mals ≃ -0.46752535023652453
"""

N = 6
H_tto = PPP_C_NH_N(N)
ψ_0 = tt_up_rks(half_filling(N),32;ϵ_wn=1e-2)
ψ_0 = orthogonalize(ψ_0)
E_mals, ψ_tt = mals_eigsolv(H_tto,ψ_0)

max_stars_repo_name: msdupuy/Tensor-Train-Julia

max_issues_repo_name: msdupuy/Tensor-Train-Julia

max_forks_repo_name: msdupuy/Tensor-Train-Julia


text:
using MolecularGaussians
using GaussianMixtureAlignment
using MolecularGraph
using MolecularGraph.Graph
using LinearAlgebra
using StaticArrays
using CoordinateTransformations
using Optim
using Test

using GaussianMixtureAlignment: distance

const MG = MolecularGaussians

@testset "Partial Charge" begin
    mol = sdftomol(joinpath(@__DIR__, "..", "data", "E1050_3d.sdf"))
    # adding partial charges to attributes
    @test_throws KeyError pc = mol.attributes[:partialcharges]
    partialcharges!(mol)
    @test typeof(mol.attributes[:partialcharges]) <: Dict
    # A sulphur atom should have the most positive partial charge
    @test nodeattr(mol,findmax(partialcharges!(mol))[2]).symbol == :S
    # An oxygen atom should have the most negative partial charge
    @test nodeattr(mol,findmin(partialcharges!(mol))[2]).symbol == :O
end

@testset "Van der Waals radii" begin
    mol = sdftomol(joinpath(@__DIR__, "..", "data", "E1050_3d.sdf"))
    # adding radii to attributes
    @test_throws KeyError pc = mol.attributes[:vdwradii]
    vdwradii!(mol)
    @test typeof(mol.attributes[:vdwradii]) <: Dict
    # All atoms should have the appropriate radius
    for idx in nodeset(mol)
        @test vdwradii!(mol)[idx] == MolecularGraph.ATOM_VANDERWAALS_RADII[atomnumber(nodeattr(mol,idx).symbol)]
    end
end

@testset "Gaussian Mixture Distance" begin
    mol = sdftomol(joinpath(@__DIR__, "..", "data", "E1050_3d.sdf"))
    gonane = sdftomol(joinpath(@__DIR__, "..", "data", "gonane5α.sdf"))
    # identical mixture models have no distance
    mol_gmm = MolGMM(mol)
    gonane_gmm = MolGMM(gonane)
    @test abs(distance(mol_gmm, mol_gmm)) < 1e-12
    @test abs(distance(gonane_gmm, gonane_gmm)) < 1e-12
    # different mixture models have some distance
    @test distance(mol_gmm, gonane_gmm) > 0.1
    # different transforms of a molecule have some distance
    gmm= MolGMM(mol)
    tform = AffineMap(10*rand(6)...)
    @test distance(tform(gmm), gmm) > 0.1
    # subgraphs of a molecule have some distance
    submol = nodesubgraph(mol, collect(nodeset(mol))[1:Int(floor(end/2))])
    sub_gmm = MolGMM(submol)
    @test distance(sub_gmm, gmm) > 0.1
end

@testset "MolGMM and PharmacophoreGMM alignment" begin
    ## MolGMM alignment
    mol1= sdftomol(joinpath(@__DIR__, "..", "data", "E1050_3d.sdf"))
    mol2 = sdftomol(joinpath(@__DIR__, "..", "data", "E1103_3d.sdf"))
    molgmm1 = MolGMM(mol1)
    molgmm2 = MolGMM(mol2)
    gonane = sdftomol(joinpath(@__DIR__, "..", "data", "gonane5α.sdf"))
    # No rotation/translation when aligning a mol to itself
    self_res = tiv_gogma_align(molgmm1, molgmm1, 0.5, 0.5; maxstagnant=1000)
    self_tform = AffineMap(self_res.tform_params...)
    @test norm(self_tform.translation) ≈ 0 atol=1e-6
    @test self_tform.linear ≈ I
    # Do you get similar distances when performing the alignment in both directions?
    f_ovrlp = tiv_gogma_align(molgmm1, molgmm2, 0.5, 0.5; maxstagnant=1000).upperbound
    b_ovrlp = tiv_gogma_align(molgmm2, molgmm1, 0.5, 0.5; maxstagnant=1000).upperbound
    @test abs(2*(f_ovrlp-b_ovrlp)/(f_ovrlp+b_ovrlp)) < 0.05     # 5% difference too forgiving?
    # Do you get the inverse transformation when swapping molfixed and molmoving (requires more evals)?
    # This may not be important for generating a distance matrix for many molecules
    # f_res= tiv_gogma_align(PharmacophoreGMM(mol1), PharmacophoreGMM(mol2); maxstagnant=1000)
    # f_tform = AffineMap(f_res.tform_params...)
    # b_res = tiv_gogma_align(PharmacophoreGMM(mol2), PharmacophoreGMM(mol1); maxstagnant=1000)
    # b_tform = AffineMap(b_res.tform_params...)
    # @test (b_tform ∘ f_tform).linear ≈ I atol=0.1
    # @test (b_tform ∘ f_tform).tranlsation ≈ zeros(3) atol=0.1
end


text:
using DataFrames
using Dates
using Printf
using SLiDE

"""
# 1A: Convert Type
See: src/utils/utils.jl

Convert between datatypes using one function. There are some julia functions
(see [signed](https://docs.julialang.org/en/v1/base/numbers/#Base.signed),
[parse](https://docs.julialang.org/en/v1/base/numbers/#Base.parse)) that can do this,
but different functions are required depending on the input type.

This function (and all of its methods) take a more stylistically-julia approach.
This is helpful when a user can specify a desired DataType.
"""

# Helper function for this example to print simple type conversions.
print_convert(x1, x2) = @printf("\t%10s %-9s --> %7s %-9s\n",
    x1, @sprintf("(%s)",typeof(x1)),
    x2, @sprintf("(%s)",typeof(x2)))

# ******************************************************************************************
# Now supports DataFrames (with one "value" column) to dictionary of keys =>value.
# This is more stylistically julia than specifying the conversion in the function name.
NYR = 3  # change if desired

yr = 2020-(NYR-1):2020
r = ["co","wi"]
NR = length(r)

df1 = DataFrame(yr = yr, value = Float64.(1:NYR))
df2 = sort(DataFrame(
    yr = repeat(yr, outer=[NR]),
    r = repeat(r, inner=[NYR]),
    value = Float64.(1:NYR*NR)))
df3 = edit_with(copy(df2), Add(:units, "millions of us dollars (usd)"))
df4 = edit_with(copy(df3), Rename(:value, :Val))

d1 = convert_type(Dict, df1)
d2 = convert_type(Dict, df2)
d3a = convert_type(Dict, df3; drop_cols = [:units])
d3b = convert_type(Dict, df3; drop_cols = :units)
d4 = convert_type(Dict, df4; drop_cols = :units, value_col = :Val)

# ******************************************************************************************
# Simple examples.
x = Date(2020)
print_convert(x, convert_type(Int64, x))

x = 1
print_convert(x, convert_type(Float64, x))
print_convert(x, convert_type(String, x))
print_convert(x, convert_type(Symbol, x))

x = 1.0
print_convert(x, convert_type(Int64, x))
print_convert(x, convert_type(String, x))
print_convert(x, convert_type(Symbol, x))

x = "1,234"
print_convert(x, convert_type(Int64, x))
print_convert(x, convert_type(Float64, x))
print_convert(x, convert_type(Symbol, x))

x = Symbol("1")
print_convert(x, convert_type(Int64, x))
print_convert(x, convert_type(Float64, x))
print_convert(x, convert_type(String, x))

x = missing
print_convert(x, convert_type(Int64, x))
print_convert(x, convert_type(Float64, x))
print_convert(x, convert_type(String, x))

max_stars_repo_name: matbesancon/SLiDE

max_issues_repo_name: matbesancon/SLiDE

max_forks_repo_name: matbesancon/SLiDE


text:
using DemoActions
using Documenter

makedocs(;
    modules=[DemoActions],
    authors="anand <anandj@uchicago.edu> and contributors",
    repo="https://github.com/anandijain/DemoActions.jl/blob/{commit}{path}#L{line}",
    sitename="DemoActions.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://anandijain.github.io/DemoActions.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/anandijain/DemoActions.jl",
    devbranch="main"
)


text:
module Restful
    using HTTP

    const methods = [:get, :post, :put, :patch, :delete, :head, :options]

    include("app.jl")
    include("route.jl")

    include("Default.jl")
    using .Default

    include("Json.jl")
    using .Json

    include("Template.jl")
    using .Template
end


text:
"""
    getindex(g)

Return graph metadata.
"""
Base.getindex(g::MetaGraph) = g.graph_data

"""
    getindex(g, label)

Return vertex metadata for `label`.
"""
Base.getindex(g::MetaGraph, label) = g.vertex_properties[label][2]

"""
    getindex(g, label_1, label_2)

Return edge metadata for the edge between `label_1` and `label_2`.
"""
Base.getindex(g::MetaGraph, label_1, label_2) = g.edge_data[arrange(g, label_1, label_2)]

"""
    haskey(g, label)

Determine whether a graph `g` contains the vertex `label`.
"""
Base.haskey(g::MetaGraph, label) = haskey(g.vertex_properties, label)

"""
    haskey(g, label_1, label_2)

Determine whether a graph `g` contains an edge from `label_1` to `label_2`.

The order of `label_1` and `label_2` only matters if `g` is a digraph.
"""
function Base.haskey(g::MetaGraph, label_1, label_2)
    return (
        haskey(g, label_1) &&
        haskey(g, label_2) &&
        haskey(g.edge_data, arrange(g, label_1, label_2))
    )
end

"""
    setindex!(g, data, label)

Set vertex metadata for `label` to `data`.
"""
function Base.setindex!(g::MetaGraph, data, label)
    if haskey(g, label)
        set_data!(g, label, data)
    else
        add_vertex!(g, label, data)
    end
    return nothing
end

"""
    setindex!(g, data, label_1, label_2)

Set edge metadata for `(label_1, label_2)` to `data`.
"""
function Base.setindex!(g::MetaGraph, data, label_1, label_2)
    if haskey(g, label_1, label_2)
        set_data!(g, label_1, label_2, data)
    else
        add_edge!(g, label_1, label_2, data)
    end
    return nothing
end

"""
    delete!(g, label)

Delete vertex `label`.
"""
function Base.delete!(g::MetaGraph, label)
    if haskey(g, label)
        v = code_for(g, label)
        _rem_vertex!(g, label, v)
    end
    return nothing
end

"""
    delete!(g, label_1, label_2)

Delete edge `(label_1, label_2)`.
"""
function Base.delete!(g::MetaGraph, label_1, label_2)
    v1, v2 = code_for(g, label_1), code_for(g, label_2)
    rem_edge!(g, v1, v2)
    return nothing
end

"""
    _copy_props!(oldg, newg, vmap)

Copy properties from `oldg` to `newg` following vertex map `vmap`.
"""
function _copy_props!(oldg::G, newg::G, vmap) where {G<:MetaGraph}
    for (newv, oldv) in enumerate(vmap)
        oldl = oldg.vertex_labels[oldv]
        _, data = oldg.vertex_properties[oldl]
        newg.vertex_labels[newv] = oldl
        newg.vertex_properties[oldl] = (newv, data)
    end
    for newe in edges(newg.graph)
        vertex_labels = newg.vertex_labels
        v1, v2 = Tuple(newe)
        label_1 = vertex_labels[v1]
        label_2 = vertex_labels[v2]
        newg.edge_data[arrange(newg, label_1, label_2, v1, v2)] = oldg.edge_data[arrange(
            oldg, label_1, label_2
        )]
    end
    return nothing
end

# TODO - It would be nice to be able to apply a function to properties.
# Not sure how this might work, but if the property is a vector,
# a generic way to append to it would be a good thing.


text:
@reset
@set printstyle none

Econv = [
-76.05293659641389
-56.20535441552476
-230.62430998333460
-154.09239259145326
-113.86482635166755
-279.11530153592219
-40.21342962572409
-378.34082589236215
-227.76309337209935
-108.95448240313179
]

Edf = [
-76.05293666271527
-56.20536076021406
-230.62430445834565
-154.09238732189544
-113.86482178236324
-279.11536397687240
-40.21342776121311
-378.34082725476162
-227.76307917196510
-108.95446350768205
]

@testset "RHF" begin
    @testset "Conventional" begin
        Fermi.Options.set("df", false)

        for i = eachindex(molecules)
            # Read molecule
            path = joinpath(@__DIR__, "xyz/"*molecules[i]*".xyz")
            mol = open(f->read(f,String), path)

            # Define options
            Fermi.Options.set("molstring", mol)
            Fermi.Options.set("basis", basis[i])

            Iu = Fermi.Integrals.IntegralHelper(eri_type=Fermi.Integrals.SparseERI())

            wf = @energy Iu => rhf
            @test isapprox(wf.energy, Econv[i], rtol=tol) # Energy from Psi4
        end
    end

    @testset "Density Fitted" begin
        Fermi.Options.set("df", true)
        Fermi.Options.set("jkfit", "cc-pvqz-jkfit")

        for i = eachindex(molecules)
            # Read molecule
            path = joinpath(@__DIR__, "xyz/"*molecules[i]*".xyz")
            mol = open(f->read(f,String), path)

            # Define options
            Fermi.Options.set("molstring", mol)
            Fermi.Options.set("basis", basis[i])

            wf = @energy rhf
            @test isapprox(wf.energy, Edf[i], rtol=tol) # Energy from Psi4
        end
    end
end


text:
# Advent of code for 2020

# Day 1

# Before you leave, the Elves in accounting just need you to fix your expense
# report (your puzzle input); apparently, something isn't quite adding up.

day_1_sample = [1721, 979, 366, 299, 675, 1456];
day_1_input = parse.(Int, readlines(open("input-01")))

function match2020(lines; total = 2020)
    expense_set = Set(lines);
    for line in lines
        if (total - line) in expense_set
            return line * (total - line)
        end
    end
    return -1
end

@assert 514579 == @show match2020(day_1_sample)
@assert 1010299 == @show match2020(day_1_input)

## Part 2
# In your expense report, what is the product of the three entries that sum to 2020?

function match2020three(lines)
    for line in lines
        product = match2020(lines; total = 2020 - line)
        if product >= 0
            return product * line
        end
    end
    return -1
end

@assert 241861950 == @show match2020three(day_1_sample)
@assert 42140160 == @show match2020three(day_1_input)

# Day 2

# Each line gives the password policy and then the password. The password policy
# indicates the lowest and highest number of times a given letter must appear
# for the password to be valid. For example, 1-3 a means that the password must
# contain a at least 1 time and at most 3 times.

function check_policy_v1(line)
    # TODO: Consider a regex for extraction instead.
    a = split(line, ":")
    b = split(a[1], " ")
    c = parse.(Int, split(b[1], "-"))
    password = strip(a[2])
    min_count = c[1]
    max_count = c[2]
    letter = b[2][1] # Turn the string into a char
    letter_count = count(i->(i==letter), password)
    min_count <= letter_count <= max_count
end

@assert true == @show check_policy_v1("1-3 a: abcde")
@assert false == @show check_policy_v1("1-3 b: cdefg")
@assert true == @show check_policy_v1("2-9 c: ccccccccc")

day_2_input = readlines(open("input-02"))

@assert 506 == @show count(check_policy_v1.(day_2_input))

# Each policy actually describes two positions in the password, where 1 means
# the first character, 2 means the second character, and so on. (Be careful;
# Toboggan Corporate Policies have no concept of "index zero"!) Exactly one of
# these positions must contain the given letter. Other occurrences of the letter
# are irrelevant for the purposes of policy enforcement.

function check_policy_v2(line)
    # TODO: Consider a regex for extraction instead.
    a = split(line, ":")
    b = split(a[1], " ")
    index = parse.(Int, split(b[1], "-"))
    password = strip(a[2])
    letter = b[2][1] # Turn the string into a char
    match_count = count(i -> (i == letter), getindex(password, index))
    1 == match_count
end

@assert true == @show check_policy_v2("1-3 a: abcde")
@assert false == @show check_policy_v2("1-3 b: cdefg")
@assert false == @show check_policy_v2("2-9 c: ccccccccc")

@assert 443 == @show count(check_policy_v2.(day_2_input))

# Day 3

# The toboggan can only follow a few specific slopes (you opted for a cheaper
# model that prefers rational numbers); start by counting all the trees you
# would encounter for the slope right 3, down 1:

day_3_sample = [
    "..##.......",
    "#...#...#..",
    ".#....#..#.",
    "..#.#...#.#",
    ".#...##..#.",
    "..#.##.....",
    ".#.#.#....#",
    ".#........#",
    "#.##...#...",
    "#...##....#",
    ".#..#...#.#"]
day_3_input = readlines(open("input-03"))

parse_trees(line) = map(c -> (c == '#'), collect(line))

@assert [false, false, true, true, false] == parse_trees("..##.")

function count_trees(input_map, delta)
    map = parse_trees.(input_map)
    dx, dy = delta
    x, y = 1, 1
    tree_count = 0
    while y <= length(map)
        if map[y][x]
            tree_count += 1
        end
        x += dx
        y += dy
        if x > length(map[1])
            # Assumption: dx < width(input_map)
            x -= length(map[1])
        end
    end
    tree_count
end

@assert 7 == @show count_trees(day_3_sample, (3, 1))
@assert 250 == @show count_trees(day_3_input, (3, 1))

## Part 2

# What do you get if you multiply together the number of trees encountered on each of the listed slopes?
# Right 1, down 1.
# Right 3, down 1. (This is the slope you already checked.)
# Right 5, down 1.
# Right 7, down 1.
# Right 1, down 2.

slopes = [(1, 1), (3, 1), (5, 1), (7, 1), (1, 2)]

@assert 336 == @show prod(count_trees.(Ref(day_3_sample),slopes))
@assert 1592662500 == @show prod(count_trees.(Ref(day_3_input),slopes))

# Day 4

# Count the number of valid passports - those that have all required fields.
# Treat cid as optional. In your batch file, how many passports are valid?

day_4_sample = [
    "ecl:gry pid:860033327 eyr:2020 hcl:#fffffd",
    "byr:1937 iyr:2017 cid:147 hgt:183cm",
    "",
    "iyr:2013 ecl:amb cid:350 eyr:2023 pid:028048884",
    "hcl:#cfa07d byr:1929",
    "",
    "hcl:#ae17e1 iyr:2013",
    "eyr:2024",
    "ecl:brn pid:760753108 byr:1931",
    "hgt:179cm",
    "",
    "hcl:#cfa07d eyr:2025 pid:166559648",
    "iyr:2011 ecl:brn hgt:59in" ]

day_4_input = readlines(open("input-04"))

function parse_passport_list(input)
    passports = []
    next_passport = Dict()
    for line in input
        if line == ""
            push!(passports, next_passport)
            next_passport = Dict()
        else
            for pair in split(line, " ")
                pair_split = split(pair, ":")
                push!(next_passport, pair_split[1] => pair_split[2])
            end
        end
    end
    push!(passports, next_passport)
    passports
end

@assert 4 == length( parse_passport_list(day_4_sample))

# byr (Birth Year)
# iyr (Issue Year)
# eyr (Expiration Year)
# hgt (Height)
# hcl (Hair Color)
# ecl (Eye Color)
# pid (Passport ID)
# cid (Country ID) # This is where we cheat

passport_cheat_fields = ["byr", "iyr", "eyr", "hgt", "hcl", "ecl", "pid"]

function is_passport_valid_cheat(passport)
    all(haskey.(Ref(passport), passport_cheat_fields))
end

@assert [true, false, true, false] == is_passport_valid_cheat.(parse_passport_list(day_4_sample))
@assert 2 == @show count(is_passport_valid_cheat.(parse_passport_list(day_4_sample)))
@assert 210 == @show count(is_passport_valid_cheat.(parse_passport_list(day_4_input)))

## Part 2
# The line is moving more quickly now, but you overhear airport security talking
# about how passports with invalid data are getting through. Better add some
# data validation, quick!

function is_field_valid(key::String, value)
    if key == "byr"
        occursin(r"^[0-9]{4}$", value) && 1920 <= parse(Int, value) <= 2002
    elseif key == "iyr"
        occursin(r"^[0-9]{4}$", value) && 2010 <= parse(Int, value) <= 2020
    elseif key == "eyr"
        occursin(r"^[0-9]{4}$", value) && 2020 <= parse(Int, value) <= 2030
    elseif key == "hgt"
        m = match(r"([0-9]+)(in|cm)", value)
        if m == Nothing()
            false
        elseif m[2] == "cm"
            150 <= parse(Int, m[1]) <= 193
        else
            59 <= parse(Int, m[1]) <= 76
        end
    elseif key == "hcl"
        occursin(r"^#[0-9a-f]{6}$", value)
    elseif key == "ecl"
        value in ["amb", "blu", "brn", "gry", "grn", "hzl", "oth"]
    elseif key == "pid"
        occursin(r"^[0-9]{9}$", value)
    end
end

# Test suite
@assert true == is_field_valid("byr", "2002" )
@assert false == is_field_valid("byr", "2003" )
@assert true == is_field_valid("hgt", "60in" )
@assert true == is_field_valid("hgt", "190cm" )
@assert false == is_field_valid("hgt", "190in" )
@assert false == is_field_valid("hgt", "190" )
@assert true == is_field_valid("hcl", "#123abc" )
@assert false == is_field_valid("hcl", "#123abz" )
@assert false == is_field_valid("hcl", "123abc" )
@assert true == is_field_valid("ecl", "brn" )
@assert false == is_field_valid("ecl", "wat" )
@assert true == is_field_valid("pid", "000000001" )
@assert false == is_field_valid("pid", "0123456789" )

function is_field_valid(passport::Dict, key::String)
    haskey(passport, key) && is_field_valid(key, passport[key])
end

function is_passport_valid_fields_cheat(passport)
    all(is_field_valid.(Ref(passport), passport_cheat_fields))
end

@assert 131 == @show count(is_passport_valid_fields_cheat.(parse_passport_list(day_4_input)))

# Day 5

# As a sanity check, look through your list of boarding passes. What is the
# highest seat ID on a boarding pass?

# Well that just sounds like binary with extra steps.

# "BFFFBBFRRR": row 70, column 7, seat ID 567.
# "FFFBBBFRRR": row 14, column 7, seat ID 119.
# "BBFFBBFRLL": row 102, column 4, seat ID 820.

seat_digits = Dict('B' => 1, 'F' => 0, 'R' => 1, 'L' => 0)

function parse_seat_id(seat_id::String)
    acc = 0
    for c in seat_id
        acc = 2 * acc + seat_digits[c]
    end
    acc
end

@assert 567 == parse_seat_id("BFFFBBFRRR")
@assert 119 == parse_seat_id("FFFBBBFRRR")
@assert 820 == parse_seat_id("BBFFBBFRLL")

day_5_input = parse_seat_id.(readlines(open("input-05")))

@assert 896 == @show maximum(day_5_input)

# Ding! The "fasten seat belt" signs have turned on. Time to find your seat.

function find_empty_seat(taken_seats)
    sorted_seats = sort(taken_seats)
    check_seat_id = sorted_seats[1]
    for seat_id in sorted_seats
        if seat_id != check_seat_id
            return check_seat_id
        else
            check_seat_id += 1
        end
    end
end

@assert 659 == @show find_empty_seat(day_5_input)

# Day 6

# As your flight approaches the regional airport where you'll switch to a much
# larger plane, customs declaration forms are distributed to the passengers.

day_6_input = readlines(open("input-06"))

function group_declarations_union( single_declarations )
    groups = []
    group = Set()
    for line in single_declarations
        if line == ""
            push!(groups, group)
            group = Set()
        else
            line_declarations = Set(collect(line))
            union!(group, line_declarations)
        end
    end
    push!(groups, group)
    groups
end

function group_declarations_intersect( single_declarations )
    groups = []
    group = Set(collect("abcdefghijklmnopqrstuvwxyz"))
    for line in single_declarations
        if line == ""
            push!(groups, group)
            group = Set(collect("abcdefghijklmnopqrstuvwxyz"))
        else
            line_declarations = Set(collect(line))
            intersect!(group, line_declarations)
        end
    end
    push!(groups, group)
    groups
end

@assert 6551 == @show sum(length.(group_declarations_union(day_6_input)))
@assert 3358 == @show sum(length.(group_declarations_intersect(day_6_input)))

# Day 7

# You land at the regional airport in time for your next flight. In fact, it
# looks like you'll even have time to grab some food: all flights are currently
# delayed due to issues in luggage processing.

day_7_sample = [
    "light red bags contain 1 bright white bag, 2 muted yellow bags.",
    "dark orange bags contain 3 bright white bags, 4 muted yellow bags.",
    "bright white bags contain 1 shiny gold bag.",
    "muted yellow bags contain 2 shiny gold bags, 9 faded blue bags.",
    "shiny gold bags contain 1 dark olive bag, 2 vibrant plum bags.",
    "dark olive bags contain 3 faded blue bags, 4 dotted black bags.",
    "vibrant plum bags contain 5 faded blue bags, 6 dotted black bags.",
    "faded blue bags contain no other bags.",
    "dotted black bags contain no other bags." ]
day_7_input = readlines(open("input-07"))

function parse_bag_rule(rules) :: Tuple{Dict{Any, Set}, Dict{Any, Vector}}
    reverse_map = Dict{Any, Set}()
    forward_map = Dict{Any, Vector}()
    for line in rules
        outside_color = match(r"(.+) bags contain", line)[1]
        has_no_content = occursin(r"no other bags", line)
        for m in eachmatch(r"([0-9]+) ([^,.]+) bags?", line)
            content_count, content_color = m.captures
            connect!(reverse_map, content_color, outside_color)
            connect!(forward_map, outside_color, (parse(Int, content_count), content_color))
            # push!(forward_map, outside_color => push!(get(forward_map, outside_color, Vector()), (content_count, content_color)))
        end
    end
    (reverse_map, forward_map)
end

function connect!(map :: Dict{Any, Set}, source, target)
    push!(map, source => push!(get(map, source, Set()), target))
end

function connect!(map :: Dict{Any, Vector}, source, target)
    push!(map, source => push!(get(map, source, Vector()), target))
end

# You have a shiny gold bag. If you wanted to carry it in at least one other
# bag, how many different bag colors would be valid for the outermost bag?
# (In other words: how many colors can, eventually, contain at least one shiny
# gold bag?)

# This is a graph question. Find all colors c s.th. there exists a
# "can contain" path from c to "shiny gold".

function reachable_from(directed_graph::Dict, starting_node) :: Set
    seen = Set()
    todo_list = Set(Ref(starting_node))
    while length(todo_list) > 0
        node = pop!(todo_list)
        targets = Set(get(directed_graph, node, Set())) # Copy, for mutation
        setdiff!(targets, seen) # we only care about never seen nodes
        union!(seen, targets)
        union!(todo_list, targets)
    end
    seen
end

@assert 4 == @show length(reachable_from(parse_bag_rule(day_7_sample)[1], "shiny gold"))
@assert 131 == @show length(reachable_from(parse_bag_rule(day_7_input)[1], "shiny gold"))

# Part 2

# How many individual bags are required inside your single shiny gold bag?

# This function assumes the graph described by forward_map is loop free.
# Otherwise it will run into an infinite loop.
function count_bags!(cache, forward_map, color)
    cached_value = get(cache, color, -1)
    if cached_value >= 0
        return cached_value
    end
    # Not in cache, so we need to calculate it.
    total = 0
    for content in get(forward_map, color, Vector())
        count, inner_color = content
        inner_count = count_bags!(cache, forward_map, inner_color)
        total = total + count * (inner_count + 1)
    end
    push!(cache, color => total)
    total
end

day_7_sample_2 = [
    "shiny gold bags contain 2 dark red bags.",
    "dark red bags contain 2 dark orange bags.",
    "dark orange bags contain 2 dark yellow bags.",
    "dark yellow bags contain 2 dark green bags.",
    "dark green bags contain 2 dark blue bags.",
    "dark blue bags contain 2 dark violet bags.",
    "dark violet bags contain no other bags."]

@assert 126 == @show count_bags!(Dict(), parse_bag_rule(day_7_sample_2)[2], "shiny gold")
@assert 11261 == @show count_bags!(Dict(), parse_bag_rule(day_7_input)[2], "shiny gold")

# Day 8 - Write a VM day :-)

@enum Opcode begin
    NoOp
    Accumulate
    Jump
end

# For now, an instruction is Tuple{Opcode, Int64}

day_8_sample = [
    "nop +0",
    "acc +1",
    "jmp +4",
    "acc +3",
    "jmp -3",
    "acc -99",
    "acc +1",
    "jmp -4",
    "acc +6"]
day_8_input = readlines(open("input-08"))

function parse_instruction(line :: String) :: Tuple{Opcode, Int64}
    m = match(r"([a-z]+) ((\+|-)[0-9]+)", line)
    (parse_opcode(m[1]), parse(Int64, m[2]))
end

function parse_opcode(code) :: Opcode
    if code == "nop"
        NoOp
    elseif code == "acc"
        Accumulate
    elseif code == "jmp"
        Jump
    else
        throw(ArgumentError("code $code not recognized"))
    end
end

mutable struct ElfVM
    code::Vector{Tuple{Opcode, Int64}}
    instructionPointer::Int
    accumulator::Int64
end

"""Initalizes an ElfVM with the instructionPointer set to the first line."""
function ElfVM(code)
    ElfVM(code, 1, 0)
end

function Base.copy(vm::ElfVM)::ElfVM
    ElfVM(copy(vm.code), vm.instructionPointer, vm.accumulator)
end

function step!(vm::ElfVM)
    instr, param = vm.code[vm.instructionPointer]
    if instr == NoOp
        vm.instructionPointer += 1
    elseif instr == Accumulate
        vm.accumulator += param
        vm.instructionPointer += 1
    elseif instr == Jump
        vm.instructionPointer += param
    end
end

function is_halted(vm::ElfVM) :: Bool
    vm.instructionPointer > length(vm.code)
end

function detect_cycle!(vm::ElfVM) :: Union{Int64, Nothing}
    seen = Set()
    while true
        if is_halted(vm)
            return Nothing()
        elseif vm.instructionPointer in seen
            return vm.accumulator
        end
        push!(seen, vm.instructionPointer)
        step!(vm)
    end
end

@assert 5 == @show detect_cycle!(ElfVM(parse_instruction.(day_8_sample)))
@assert 1675 == @show detect_cycle!(ElfVM(parse_instruction.(day_8_input)))

# Part 2

# Fix the program so that it terminates normally by changing exactly one jmp
# (to nop) or nop (to jmp). What is the value of the accumulator after the
# program terminates?

function switch_opcode(original::Opcode) :: Opcode
    if original == NoOp
        Jump
    elseif original == Jump
        NoOp
    else
        original
    end
end

function fix_program(original_vm::ElfVM) :: Int64
    for i in 1:length(original_vm.code)
        vm = copy(original_vm)

        # Try to fix the vm
        old_opcode, old_param = vm.code[i]
        vm.code[i] = (switch_opcode(old_opcode), old_param)

        cycle_result = detect_cycle!(vm)
        if cycle_result == Nothing()
            return vm.accumulator
        end
    end
    throw("No fix found")
end

@assert 8 == @show fix_program(ElfVM(parse_instruction.(day_8_sample)))
@assert 1532 == @show fix_program(ElfVM(parse_instruction.(day_8_input)))

# Day 9

# The data appears to be encrypted with the eXchange-Masking Addition System (XMAS)
# which, conveniently for you, is an old cypher with an important weakness.

day_9_sample = [35, 20, 15, 25, 47, 40, 62, 55, 65, 95, 102, 117, 150, 182, 127,
    219, 299, 277, 309, 576]
day_9_input = parse.(Int, readlines(open("input-09")))

function check_XMAS_code(input; preamble_length=15, offset=0)
    choices = input[1+offset:preamble_length+offset]
    target_sum = input[preamble_length+offset+1]
    for i in 1:preamble_length, j in i+1:preamble_length
        if choices[i] + choices[j] == target_sum
            return (i, j)
        end
    end
    Nothing()
end

function find_first_invalid_XMAS_code(input; preamble_length=25)
    for offset in 0:length(input)-preamble_length
        combination = check_XMAS_code(input; preamble_length= preamble_length, offset = offset)
        if combination == Nothing()
            return input[preamble_length+offset+1]
        end
    end
    Nothing()
end

@assert 127 == @show find_first_invalid_XMAS_code(day_9_sample; preamble_length=5)
@assert 18272118 == @show find_first_invalid_XMAS_code(day_9_input; preamble_length=25)

# Part 2

# The final step in breaking the XMAS encryption relies on the invalid number
# you just found: you must find a contiguous set of at least two numbers in your
# list which sum to the invalid number from step 1.

function find_continuous(input, invalid_number)
    # start_ptr >= end_ptr, imagine a snake game snake.
    start_ptr = 1
    end_ptr = 1
    sum = input[1]
    while sum != invalid_number
        if sum > invalid_number
            sum -= input[end_ptr]
            end_ptr += 1
        elseif sum < invalid_number
            start_ptr += 1
            sum += input[start_ptr]
        end
    end
    (start_ptr, end_ptr)
end

function find_min_max_sum(input; preamble_length=25)
    invalid_number = find_first_invalid_XMAS_code(input; preamble_length = preamble_length)
    start_ptr, end_ptr = find_continuous(input, invalid_number)
    # remember, start_ptr is the snake head and is greater than the tail end_ptr
    sum_range = input[end_ptr:start_ptr]
    minimum(sum_range) + maximum(sum_range)
end

@assert 62 == @show find_min_max_sum(day_9_sample; preamble_length=5)
@assert 2186361 == @show find_min_max_sum(day_9_input; preamble_length=25)

# Day 10, Joltage

day_10_sample_1 = [ 16, 10, 15, 5, 1, 11, 7, 19, 6, 12, 4 ]
day_10_sample_2 = [ 28, 33, 18, 42, 31, 14, 46, 20, 48, 47, 24, 23, 49, 45, 19,
    38, 39, 11, 1, 32, 25, 35, 8, 17, 7, 9, 4, 2, 34, 10, 3 ]
day_10_input = parse.(Int, readlines(open("input-10")))

function joltage_difference_list(chargers)
    chargers_sorted = copy(chargers)
    push!(chargers_sorted, 0)
    chargers_sorted = sort!(chargers_sorted)
    push!(chargers_sorted, chargers_sorted[end]+3)
    differences = [chargers_sorted[i] - chargers_sorted[i-1] for i in 2:length(chargers_sorted)]
end

function joltage_chain_product(chargers)
    differences = joltage_difference_list(chargers)
    @assert 0 == count(i -> i == 2, differences)
    count(i -> i == 1, differences) * count(i -> i == 3, differences)
end

@assert 35 == @show joltage_chain_product(day_10_sample_1)
@assert 220 == @show joltage_chain_product(day_10_sample_2)
@assert 1820 == @show joltage_chain_product(day_10_input)

# Part 2

# How many choices are there for [3, 1, 1, ..., 1, 3] with n-1 ones?
_bcc = [1, 1, 2]

for i in 4:100
    push!(_bcc, _bcc[i-1] + _bcc[i-2] + _bcc[i-3])
end

block_choices(n) = _bcc[n]

@assert block_choices(4) == 4
@assert block_choices(5) == 7

function total_pathways(chargers)
    differences = joltage_difference_list(chargers)
    product = 1
    len = 0
    for i in differences
        if i == 1
            len += 1
        elseif i == 3
            len += 1
            product *= block_choices(len)
            len = 0
        else
            throw("Unexpected difference between chargers")
        end
    end
    product
end

@assert 8 == @show total_pathways(day_10_sample_1)
@assert 19208 == @show total_pathways(day_10_sample_2)
@assert 3454189699072 == @show total_pathways(day_10_input)

# Day 11 -- Cellular Automata

day_11_sample = [
    "L.LL.LL.LL",
    "LLLLLLL.LL",
    "L.L.L..L..",
    "LLLL.LL.LL",
    "L.LL.LL.LL",
    "L.LLLLL.LL",
    "..L.L.....",
    "LLLLLLLLLL",
    "L.LLLLLL.L",
    "L.LLLLL.LL" ]
day_11_input = readlines(open("input-11"))

function parse_seat_layout(layout)
    is_seat.(hcat(collect.(layout)...))
end

is_seat(char::Char) = char == 'L'

function evolve_seat_occupation(layout, occupation = zeros(Int8, size(layout)))
    w, h = size(layout)
    count = count_nbhd.(Ref(occupation), 1:w, transpose(1:h))
    next_occupation.(count, occupation, layout)
end

function evolve_seat_occupation_part_2(layout, occupation = zeros(Int8, size(layout)))
    w, h = size(layout)
    count = count_view.(Ref(layout), Ref(occupation), 1:w, transpose(1:h))
    next_occupation.(count, occupation, layout; tollerance = 5)
end

function count_nbhd(occupation, x, y)
    w, h = size(occupation)
    count = Int8(0)
    for i in -1:1, j in -1:1
        if i == 0 && j == 0
            continue
        end
        _x = x + i
        _y = y + j
        if (1 <= _x <= w) && (1 <= _y <= h)
            count += occupation[_x, _y]
        end
    end
    count
end

# Now, instead of considering just the eight immediately adjacent seats,
# consider the first seat in each of those eight directions.
function count_view(layout, occupation, x, y)
    w, h = size(occupation)
    count = Int8(0)
    for i in -1:1, j in -1:1
        if i == 0 && j == 0
            continue
        end
        for t in 1:w
            _x = x + i * t
            _y = y + j * t
            if !((1 <= _x <= w) && (1 <= _y <= h))
                break
            elseif layout[_x, _y]
                if (1 <= _x <= w) && (1 <= _y <= h)
                    count += occupation[_x, _y]
                end
                break
            end
        end
    end
    count
end

function next_occupation(count::Int8, occupation::Int8, is_seat::Bool; tollerance = 4) :: Int8
    if is_seat && occupation == 0 && count == 0
        # If a seat is empty (L) and there are no occupied seats adjacent to it,
        # the seat becomes occupied.
        1
    elseif is_seat && occupation == 1 && count >= tollerance
        # If a seat is occupied (#) and four or more seats adjacent to it are
        # also occupied, the seat becomes empty.
        0
    else
        # Otherwise, the seat's state does not change.
        occupation
    end
end

function evolve_to_fixed_point(layout)
    occupation = evolve_seat_occupation(layout)
    new_occupation = evolve_seat_occupation(layout, occupation)
    while occupation != new_occupation
        occupation = new_occupation
        new_occupation = evolve_seat_occupation(layout, occupation)
    end
    occupation
end

@assert 37 == @show sum(evolve_to_fixed_point(parse_seat_layout(day_11_sample)))
@assert 2427 == @show sum(evolve_to_fixed_point(parse_seat_layout(day_11_input)))

function evolve_to_fixed_point_part_2(layout)
    occupation = evolve_seat_occupation_part_2(layout)
    new_occupation = evolve_seat_occupation_part_2(layout, occupation)
    while occupation != new_occupation
        occupation = new_occupation
        new_occupation = evolve_seat_occupation_part_2(layout, occupation)
    end
    occupation
end

@assert 26 == @show sum(evolve_to_fixed_point_part_2(parse_seat_layout(day_11_sample)))
@assert 2199 == @show sum(evolve_to_fixed_point_part_2(parse_seat_layout(day_11_input)))

# Day 12

day_12_sample = [
    "F10",
    "N3",
    "F7",
    "R90",
    "F11"]
day_12_input = readlines(open("input-12"))

function drive_ferry(instructions) :: Tuple
    # Coordinate system is a default 'math' coordinate system
    # x is right (east), y is up (north)
    position = (0, 0)
    # The ship starts by facing east.
    ship_orientation = 'E'
    for instruction in instructions
        (direction, distance) = decode_ferry_instruction(instruction)
        if direction == 'F'
            position = drive_ferry_compass(ship_orientation, distance, position...)
        elseif direction == 'R'
            for _ in 1:(distance/90)
                ship_orientation = turn_ferry_right[ship_orientation]
            end
        elseif direction == 'L'
            for _ in 1:(distance/90)
                ship_orientation = turn_ferry_left[ship_orientation]
            end
        else
            position = drive_ferry_compass(direction, distance, position...)
        end
    end
    position
end

function drive_ferry_compass(direction::Char, distance, x, y) :: Tuple
    if direction == 'N'
        (x, y + distance)
    elseif direction == 'E'
        (x + distance, y)
    elseif direction == 'S'
        (x, y - distance)
    elseif direction == 'W'
        (x - distance, y)
    end
end

turn_ferry_right = Dict('N' => 'E', 'E' => 'S', 'S' => 'W', 'W' => 'N')
turn_ferry_left = Dict('N' => 'W', 'W' => 'S', 'S' => 'E', 'E' => 'N')

decode_ferry_instruction(instruction::String) = (instruction[1], parse(Int, instruction[2:end]))

@assert 25 == @show sum(abs.(drive_ferry(day_12_sample)))
@assert 858 == @show sum(abs.(drive_ferry(day_12_input)))

# Part 2

function drive_ferry_waypoint(instructions) :: Tuple
    position = (0, 0)
    waypoint = (10, 1)
    ship_orientation = 'E'
    for instruction in instructions
        (direction, distance) = decode_ferry_instruction(instruction)

        if direction == 'F'
            # drive to the waypoint `distance` times.
            position = position .+ waypoint .* distance
        elseif direction == 'R'
            # rotate the waypoint around the ship
            for _ in 1:(distance/90)
                waypoint = (waypoint[2], -waypoint[1])
            end
        elseif direction == 'L'
            for _ in 1:(distance/90)
                waypoint = (-waypoint[2], waypoint[1])
            end
        else
            waypoint = drive_ferry_compass(direction, distance, waypoint...)
        end
    end
    position
end

@assert 286 == @show sum(abs.(drive_ferry_waypoint(day_12_sample)))
@assert 39140 == @show sum(abs.(drive_ferry_waypoint(day_12_input)))

# Day 13, Bus schedule

day_13_sample = (939, [7,13,Nothing(),Nothing(),59,Nothing(),31,19])
day_13_lines = readlines(open("input-13"))
parse_bus_number(word) = word == "x" ? Nothing() : parse(Int, word)
day_13_input = (parse(Int, day_13_lines[1]), parse_bus_number.(split(day_13_lines[2], ",")))

# What is the ID of the earliest bus you can take to the airport multiplied by
# the number of minutes you'll need to wait for that bus?

function minutes_of_wait(now, bus_number)
    mod_offset = now % bus_number
    (bus_number - mod_offset, bus_number)
end

function best_bus(input)
    wait_vector = minutes_of_wait.(input[1], filter(x -> x != Nothing(), input[2]))
    i = argmin(getindex.(wait_vector, 1))
    wait_time, bus_number = wait_vector[i]
    wait_time * bus_number
end

@assert 295 == @show best_bus(day_13_sample)
@assert 4135 == @show best_bus(day_13_input)

# Part 2
# The shuttle company is running a contest: one gold coin for anyone that can
# find the earliest timestamp such that the first bus ID departs at that time
# and each subsequent listed bus ID departs at that subsequent minute. (The
# first line in your input is no longer relevant.)

# Looks like this requires you to apply the chinese remainder theorem.

function chinese_bus(bus_numbers)
    modulus = []
    remainder = []
    for (index, value) in enumerate(bus_numbers)
        if value != Nothing()
            push!(modulus, value)
            # t + offset = 0 (mod modulus)
            # => t = modulus - offset (mod modulus)
            push!(remainder, value - (index - 1))
        end
    end
    # Now, we appy the chinese remainder theorem to these requirements.
    chineseremainder(Int128.(modulus), Int128.(remainder))
end

# See https://rosettacode.org/wiki/Chinese_remainder_theorem
# This implementation leaves the Int64 bounds for my puzzle input -.-
function chineseremainder(n::Array{Int128}, a::Array{Int128})
    full_product = prod(n)
    mod(sum(ai * invmod(full_product ÷ ni, ni) * full_product ÷ ni for (ni, ai) in zip(n, a)), full_product)
end

@assert 3417 == chinese_bus([17, nothing, 13, 19])
@assert 754018 == chinese_bus([67,7,59,61 ])
@assert 779210 == chinese_bus([67,nothing,7,59,61 ])
@assert 1261476 == chinese_bus([67,7,nothing,59,61 ])
@assert 1202161486 == chinese_bus([1789,37,47,1889 ])

@assert 1068781 == @show chinese_bus(day_13_sample[2])
@assert Int128(640856202464541) == @show chinese_bus(day_13_input[2])

# Day 14

# Doube bit masks and memory

day_14_sample =[
    "mask = XXXXXXXXXXXXXXXXXXXXXXXXXXXXX1XXXX0X",
    "mem[8] = 11",
    "mem[7] = 101",
    "mem[8] = 0" ]
day_14_input = readlines(open("input-14"))

zero_mask_dict = Dict('1' => 1, '0' => 0, 'X' => 1)
one_mask_dict = Dict('1' => 1, '0' => 0, 'X' => 0)

function binary_to_int(vec) :: Int64
    acc = 0
    for v in vec
        acc = 2 * acc + v
    end
    acc
end

@assert binary_to_int([0, 1, 1, 0]) == 6

# Mask where only the '1' char maps to a 1 bit. We | with this mask.
parse_one_mask(str) = binary_to_int(getindex.(Ref(one_mask_dict), collect(str)))

function parse_masks(str) :: Tuple{Int64, Int64}
    # Mask where only the '0' char maps to a 0 bit. We & with this mask.
    zero_mask = binary_to_int(getindex.(Ref(zero_mask_dict), collect(str)))
    one_mask = parse_one_mask(str)
    (zero_mask, one_mask)
end

use_masks(value::Int64, masks::Tuple{Int64, Int64}) :: Int64 = (value & masks[1]) | masks[2]

function run_bitmask_program(lines)
    mask = nothing
    mem = Dict()
    for line in lines
        mask_match = match(r"mask = ([01X]+)", line)
        mem_match = match(r"mem\[([0-9]+)\] = ([0-9]+)", line)
        if mask_match != nothing
            mask = parse_masks(mask_match[1])
        elseif mem_match != nothing
            index = parse(Int64, mem_match[1])
            value = parse(Int64, mem_match[2])
            mem[index] = use_masks(value, mask)
        else
            throw("Could not parse input: " + line)
        end
    end
    sum(values(mem))
end

@assert 165 == @show run_bitmask_program(day_14_sample)
@assert 15018100062885 == @show run_bitmask_program(day_14_input)

# Part 2

day_14_sample_2 = [
    "mask = 000000000000000000000000000000X1001X",
    "mem[42] = 100",
    "mask = 00000000000000000000000000000000X0XX",
    "mem[26] = 1" ]

# By running the regex (X[^X]*){9} and (X[^X]*){10} on my input I find that
# there are at most 9 floating bits ('X's) in any mask. This means just
# iterating the <= 512 combinations is ok and I don't need any clever math.

# If the bitmask bit is 0, the corresponding memory address bit is unchanged.
# If the bitmask bit is 1, the corresponding memory address bit is overwritten with 1.
# If the bitmask bit is X, the corresponding memory address bit is floating.

parse_floating_mask(str) = getindex.(Ref(Dict('1' => 0, '0' => 0, 'X' => 1)), collect(str))
function build_floating_mask_vector(str)
    parsed = parse_floating_mask(str)
    # Zeroes stuff, so bits are 1 by default.
    zero_masks = [Int64(0)]
    # Sets stuff to one, so bits are 0 by default.
    one_masks = [Int64(0)]
    for p in parsed
        if p == 0
            zero_masks .*= 2
            zero_masks .+= 1
            one_masks .*= 2
        else
            zero_masks .*= 2
            one_masks .*= 2
            append!(zero_masks, zero_masks .+ 1)
            append!(one_masks, one_masks .+ 1)
        end
    end
    (zero_masks, one_masks)
end

@assert ([68719476702, 68719476734, 68719476703, 68719476735], [0, 32, 1, 33]) ==
    build_floating_mask_vector("000000000000000000000000000000X1001X")

function write_dict!(d, k, v)
    d[k] = v
end

function run_bitmask_program_2(lines)
    one_mask = nothing
    z, o = nothing, nothing
    mem = Dict()
    for line in lines
        mask_match = match(r"mask = ([01X]+)", line)
        mem_match = match(r"mem\[([0-9]+)\] = ([0-9]+)", line)
        if mask_match != nothing
            one_mask = parse_one_mask(mask_match[1])
            z, o = build_floating_mask_vector(mask_match[1])
        elseif mem_match != nothing
            index = parse(Int64, mem_match[1])
            value = parse(Int64, mem_match[2])
            indices = ((index .& z) .| o) .| one_mask
            write_dict!.(Ref(mem), indices, value)
        else
            throw("Could not parse input: " + line)
        end
    end
    sum(values(mem))
end

@assert 208 == @show run_bitmask_program_2(day_14_sample_2)
@assert 5724245857696 == @show run_bitmask_program_2(day_14_input)

# Day 15 - The game ends when the Elves get sick of playing.

function number_diff_game(starting_sequence::Vector{Int64}, until::Int64)
    i = 1
    # Setup with the starting sequence
    last_spoken_dict = Dict{Int64, Int64}()
    for number in starting_sequence
        last_spoken_dict[number] = i
        i += 1
    end
    to_speak = 0
    while i < until
        if to_speak in keys(last_spoken_dict)
            j = last_spoken_dict[to_speak]
            last_spoken_dict[to_speak] = i
            to_speak = i - j
            i += 1
        else
            last_spoken_dict[to_speak] = i
            to_speak = 0
            i += 1
        end
    end
    to_speak
end

@assert 436 == number_diff_game([0,3,6], 2020)
@assert 1 == number_diff_game([1,3,2], 2020)
@assert 10 == number_diff_game([2,1,3], 2020)
@assert 27 == number_diff_game([1,2,3], 2020)
@assert 78 == number_diff_game([2,3,1], 2020)
@assert 438 == number_diff_game([3,2,1], 2020)
@assert 1836 == number_diff_game([3,1,2], 2020)

day_15_input = [2,0,1,7,4,14,18]
@assert 496 == @show number_diff_game(day_15_input, 2020)

# This takes a long time, so we are not recomputing it each time.
# @assert 175594 == number_diff_game([0,3,6], 30000000)
# @assert 2578 == number_diff_game([1,3,2], 30000000)
# @assert 3544142 == number_diff_game([2,1,3], 30000000)
# @assert 261214 == number_diff_game([1,2,3], 30000000)
# @assert 6895259 == number_diff_game([2,3,1], 30000000)
# @assert 18 == number_diff_game([3,2,1], 30000000)
# @assert 362 == number_diff_game([3,1,2], 30000000)

# @assert 883 == @show number_diff_game(day_15_input, 30000000)

# Day 16

day_16_sample = [
    "class: 1-3 or 5-7",
    "row: 6-11 or 33-44",
    "seat: 13-40 or 45-50",
    "",
    "your ticket:",
    "7,1,14",
    "",
    "nearby tickets:",
    "7,3,47",
    "40,4,50",
    "55,2,20",
    "38,6,12" ]
day_16_input = readlines(open("input-16"))

struct TicketConstraint
    name :: String
    min_1 :: Int
    max_1 :: Int
    min_2 :: Int
    max_2 :: Int
end

function parse_ticket_constraint(line)
    m = match(r"([^:]+): ([0-9]+)-([0-9]+) or ([0-9]+)-([0-9]+)", line)
    if m != nothing
        TicketConstraint(m[1], parse(Int, m[2]), parse(Int, m[3]), parse(Int, m[4]), parse(Int, m[5]))
    else
        nothing
    end
end

function fits_constraint(c::TicketConstraint, number::Int)::Bool
    (c.min_1 <= number <= c.max_1) || (c.min_2 <= number <= c.max_2)
end

function do_day_16_part_1(input)
    constraints = []
    while input[1] != ""
        push!(constraints, parse_ticket_constraint(popat!(input, 1)))
    end
    popat!(input, 1) # Remove the empty line that indicates now comes our ticket
    @assert "your ticket:" == popat!(input, 1)
    my_ticket = parse.(Int, split(popat!(input, 1), ","))
    @assert "" == popat!(input, 1)
    @assert "nearby tickets:" == popat!(input, 1)
    other_tickets = []
    for line in input
        push!(other_tickets, parse.(Int, split(line, ",")))
    end

    # Now do the actual computation:
    error_sum = 0
    for ticket in other_tickets, number in ticket
        if !any(fits_constraint.(constraints, number))
            error_sum += number
        end
    end
    error_sum
end

@assert 71 == @show do_day_16_part_1(copy(day_16_sample))
@assert 23925 == @show do_day_16_part_1(copy(day_16_input))

# Part 2:

# TODO: Missing for now.

# Day 17: Conway Cubes

# - If a cube is active and exactly 2 or 3 of its neighbors are also active, the
#   cube remains active. Otherwise, the cube becomes inactive.
# - If a cube is inactive but exactly 3 of its neighbors are active, the cube
#   becomes active. Otherwise, the cube remains inactive.

day_17_sample = [
    ".#.",
    "..#",
    "###"]

day_17_input = [
    "##.#....",
    "...#...#",
    ".#.#.##.",
    "..#.#...",
    ".###....",
    ".##.#...",
    "#.##..##",
    "#.####.."]

function parse_conway_cube(input)
    grid = Int8.('#' .== hcat(collect.(input)...))
    reshape(grid, size(grid)..., 1)
end

"""Grows the conway space and computes the next evolution step."""
function next_conway_cube(old_small)
    # Copy data into bigger space
    old = zeros(Int8, (size(old_small) .+ 2)...)
    w, h, l = size(old_small)
    for x in 1:w, y in 1:h, z in 1:l
        old[x+1, y+1, z+1] = old_small[x, y, z]
    end

    new = zeros(Int8, (size(old))...)
    w, h, l = size(new)
    for x in 1:w, y in 1:h, z in 1:l
        acc = Int8(0)
        for dx in -1:1, dy in -1:1, dz in -1:1
            if dx == dy == dz == 0
                continue
            elseif 1 <= x+dx <= w && 1 <= y+dy <= w && 1 <= z+dz <= l
                acc += old[x+dx, y+dy, z+dz]
                # if acc >= 4
                #     # no need to continue with the full cube
                #     break
                # end
            end
        end
        if acc == 3
            new[x, y, z] = 1
        elseif acc == 2 && old[x, y, z] == 1
            new[x, y, z] = 1
        end
    end
    new
end

iterate_conway_cube(cube, i=6) = i <= 0 ? cube : iterate_conway_cube(next_conway_cube(cube), i-1)

@assert 112 == @show sum(iterate_conway_cube(parse_conway_cube(day_17_sample)))
@assert 291 == @show sum(iterate_conway_cube(parse_conway_cube(day_17_input)))

# Part 2

# Well that just sounds like part 1 with extra dimensions.

function parse_conway_cube_4(input)
    grid = Int8.('#' .== hcat(collect.(input)...))
    reshape(grid, size(grid)..., 1, 1)
end

"""Grows the conway space and computes the next evolution step."""
function next_conway_cube_4(old_small)
    # Copy data into bigger space
    old = zeros(Int8, (size(old_small) .+ 2)...)
    wi, h, l, m = size(old_small)
    for x in 1:wi, y in 1:h, z in 1:l, w in 1:m
        old[x+1, y+1, z+1, w+1] = old_small[x, y, z, w]
    end

    new = zeros(Int8, (size(old))...)
    wi, h, l, m = size(new)
    for x in 1:wi, y in 1:h, z in 1:l, w in 1:m
        acc = Int8(0)
        for dx in -1:1, dy in -1:1, dz in -1:1, dw in -1:1
            if dx == dy == dz == dw == 0
                continue
            elseif 1 <= x+dx <= wi && 1 <= y+dy <= h && 1 <= z+dz <= l && 1 <= w+dw <= m
                acc += old[x+dx, y+dy, z+dz, w+dw]
            end
        end
        if acc == 3
            new[x, y, z, w] = 1
        elseif acc == 2 && old[x, y, z, w] == 1
            new[x, y, z, w] = 1
        end
    end
    new
end

iterate_conway_cube_4(cube, i=6) = i <= 0 ? cube : iterate_conway_cube_4(next_conway_cube_4(cube), i-1)

@assert 848 == @show sum(iterate_conway_cube_4(parse_conway_cube_4(day_17_sample)))
@assert 1524 == @show sum(iterate_conway_cube_4(parse_conway_cube_4(day_17_input)))


text:
function difference_of_psd(A)
if ishermitian(A)
    S = eigen(A);
    D = S.values;
    V = S.vectors;
    P = (D+abs.(D))/2;
    B = V*Diagonal(P)*V'; B = (B+B')/2;
    C = B-A;
    return(B,C)
end
end

max_stars_repo_name: cbingane/optigon

max_issues_repo_name: cbingane/optigon

max_forks_repo_name: cbingane/optigon


text:
##
function quitMenuItem_activate_cb(widgetptr::Ptr, user_data)
    #widget = convert(GtkMenuItem, widgetptr)
    main_window = user_data
    destroy(main_window)#TODO global
    return nothing
end
function newMenuItem_activate_cb(widgetptr::Ptr, user_data)
    main_window = user_data
    add_tab(main_window.editor)
    save(project)##FIXME this souldn't be here
    return nothing
end
function openMenuItem_activate_cb(widgetptr::Ptr, user_data)
    main_window = user_data
    openfile_dialog(main_window.editor)
    return nothing
end

function openREADME_MenuItem_activate_cb(widgetptr::Ptr, user_data)
    main_window = user_data
    open_in_new_tab(joinpath(HOMEDIR, "..", "README.md"), main_window.editor)
    nothing
end
function user_settings_MenuItem_activate_cb(widgetptr::Ptr, user_data)
    main_window = user_data
    open_in_new_tab(joinpath(HOMEDIR, "config", "user_settings.ini"), main_window.editor)
    nothing
end

function MainMenu(main_window::MainWindow)

    menubar = GtkMenuBar() |>
        (fileMenu = GtkMenuItem("_File")) |>
        (settingsMenu= GtkMenuItem("_Settings"))|>
        (helpMenu = GtkMenuItem("_Help"))

    buildmenu([
        MenuItem("New File", newMenuItem_activate_cb),
        MenuItem("Open File", openMenuItem_activate_cb),
        GtkSeparatorMenuItem,
        MenuItem("Quit", quitMenuItem_activate_cb),
        ],
        fileMenu,
        (main_window)
    )

    #Settings
    buildmenu([
        MenuItem("User Settings", user_settings_MenuItem_activate_cb),
        ],
        settingsMenu,
        (main_window)
    )

    # Help
    buildmenu([
        MenuItem("README", openREADME_MenuItem_activate_cb),
        ],
        helpMenu,
        (main_window)
    )

    return menubar
end


text:
using Test
using BinaryProvider, FileIO, Random, Pkg
using MakieGallery
using Makie, AbstractPlotting
using Statistics

database = (get(ENV, "MAKIEGALLERY_MINIMAL", "false") == "true") ? MakieGallery.load_tests() : MakieGallery.load_database()

# THese examples download additional data - don't want to deal with that!
# to_skip = ["WorldClim visualization", "Image on Geometry (Moon)", "Image on Geometry (Earth)"]
# # we directly modify the database, which seems easiest for now
# filter!(entry-> !(entry.title in to_skip), database)
tested_diff_path = joinpath(@__DIR__, "tested_different")
test_record_path = joinpath(@__DIR__, "test_recordings")
rm(tested_diff_path, force = true, recursive = true)
mkpath(tested_diff_path)
rm(test_record_path, force = true, recursive = true)
mkpath(test_record_path)
# rerecord = [
#     "DifferentialEquations path animation"
# ]
# filter!(x-> x.title in rerecord, database)
examples = MakieGallery.record_examples(test_record_path)
# MakieGallery.generate_preview(test_record_path, joinpath(homedir(), "Desktop", "index.html"))
# MakieGallery.generate_thumbnails(test_record_path)
# MakieGallery.gallery_from_recordings(test_record_path, joinpath(test_record_path, "index.html"))
MakieGallery.run_comparison(test_record_path, tested_diff_path)


text:
struct Regularize{N,F}

  "x values of points, normalized to grid index space"
  x :: Vector{Float64}

  "y values of points, normalized to grid index space"
  y :: Vector{Float64}

  "1/dV factor"
  overdv :: Float64

  "weights for each point (e.g. arclengths), divided by dV"
  wgt :: Vector{Float64}

  "Discrete Delta function"
  ddf :: AbstractDDF

  "Radius of the DDF"
  ddf_radius :: Float64

  "Symmetry flag"
  _issymmetric :: Bool

end

"""
    Regularize(x,y,dx,[ddftype=Yang3],[graddir=0],[I0=(1,1)], [weights=1.0], [filter=false],
                       [issymmetric=false])

Constructor to set up an operator for regularizing and interpolating data from/to
points immersed in the grid to/from fields on the grid itself. The supplied
`x` and `y` represent physical coordinates of the immersed points, and `dx`
denotes a uniform physical cell size of the grid. The separate arguments `x` and
`y` can be replaced by a single argument `X` of type `VectorData` holding the
coordinates.

The operations of regularization and interpolation are carried out with a discrete
delta function (ddf), which defaults to the type `Yang3`. Others are also possible,
such as `Roma`, `Goza` or `M3`. The optional argument `graddir`, if set to 1 or 2, will
generate an interpolation operator that evaluates the negative of the
respective component of the gradient of a grid field at the immersed points. The
default value of this argument is 0, which simply interpolates. Note that the
regularization form of this gradient type is also possible.

The optional tuple
`I0` represents the indices of the primary node that coincides with `(x,y) = (0,0)`.
This defaults to `(1,1)`, which leaves one layer of ghost (dual) cells and sets
the physical origin in the lower left corner of the grid of interior dual cells.

Another optional parameter, `weights`, sets the weight of each point in the
regularization. This would generally be set with, say, the differential arc
length for regularization of data on a curve. It can be a vector (of the same length
as x and y) or a scalar if uniform. It defaults to 1.0.

The optional Boolean parameter `filter` can be set to `true` if it is desired to
apply filtering (see Goza et al, J Comput Phys 2016) to the grid data before
interpolating. This is generally only used in the context of preconditioning
the solution for forces on the immersed points.

If the optional Boolean parameter `issymmetric` is set to `true`, then the
regularization and interpolation are constructed to be transposes of each other.
Note that this option overrides any supplied weights. The default of this
parameter is `false`.

The resulting operator can be used in either direction, regularization and
interpolation, with the first argument representing the *target* (the entity
to regularize/interpolate to), and the second argument
the *source* (the entity to regularize/interpolate from). The regularization
does not use the filtering option.

# Example

In the example below, we set up a 12 x 12 grid. Using the default value for `I0`
and setting `dx = 0.1`, the physical dimensions of the non-ghost part of the grid
are 1.0 x 1.0. Three points are set up in the interior, and a vector field is assigned
to them, with the x component of each of them set to 1.0. These data are regularized
to a field of primal edges on the grid, using the Roma DDF kernel.

```jldoctest
julia> x = [0.25,0.75,0.25]; y = [0.75,0.25,0.25];

julia> X = VectorData(x,y);

julia> q = Edges(Primal,(12,12));

julia> dx = 0.1;

julia> H = Regularize(x,y,dx;ddftype=Roma)
Regularization/interpolation operator with non-filtered interpolation
  DDF type CartesianGrids.Roma
  3 points in grid with cell area 0.01

julia> f = VectorData(X);

julia> fill!(f.u,1.0);

julia> H(q,f)
Edges{Primal,12,12,Float64} data
u (in grid orientation)
11×12 Array{Float64,2}:
 0.0  0.0  0.0       0.0     0.0      …  0.0       0.0     0.0      0.0  0.0
 0.0  0.0  0.0       0.0     0.0         0.0       0.0     0.0      0.0  0.0
 0.0  0.0  8.33333  33.3333  8.33333     0.0       0.0     0.0      0.0  0.0
 0.0  0.0  8.33333  33.3333  8.33333     0.0       0.0     0.0      0.0  0.0
 0.0  0.0  0.0       0.0     0.0         0.0       0.0     0.0      0.0  0.0
 0.0  0.0  0.0       0.0     0.0      …  0.0       0.0     0.0      0.0  0.0
 0.0  0.0  0.0       0.0     0.0         0.0       0.0     0.0      0.0  0.0
 0.0  0.0  8.33333  33.3333  8.33333     8.33333  33.3333  8.33333  0.0  0.0
 0.0  0.0  8.33333  33.3333  8.33333     8.33333  33.3333  8.33333  0.0  0.0
 0.0  0.0  0.0       0.0     0.0         0.0       0.0     0.0      0.0  0.0
 0.0  0.0  0.0       0.0     0.0      …  0.0       0.0     0.0      0.0  0.0
v (in grid orientation)
12×11 Array{Float64,2}:
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
```
"""
function Regularize(x::AbstractVector{T},y::AbstractVector{T},dx::T;
                    ddftype=Yang3,graddir::Int=0,
                    I0::Tuple{Int,Int}=(1,1),
                    weights::Union{T,Vector{T}}=1.0,
                    filter::Bool = false,
                    issymmetric::Bool = false) where {T<:Real}

  _issymmetric = (filter ? false : issymmetric)

  n = length(x)
  @assert length(y)==n
  if !_issymmetric
    if typeof(weights) == T
      wtvec = similar(x)
      fill!(wtvec,weights/(dx*dx))
    else
      @assert length(weights)==n
      wtvec = deepcopy(weights)./(dx*dx)
    end
  else
    # if the regularization and interpolation are symmetric, then the
    # weights are automatically set to be the cell area in order to cancel it
    # in the denominator of the regularization operator.
    wtvec = similar(x)
    fill!(wtvec,1.0)
  end

  baseddf = DDF(ddftype=ddftype,dx=1.0)
  if graddir == 0
    ddf = baseddf
  else
    ddf = GradDDF(graddir,ddftype=ddftype,dx=1.0)
  end

  Regularize{length(x),filter}(x./dx.+I0[1],y./dx.+I0[2],1.0/(dx*dx),
                      wtvec,ddf,_get_regularization_radius(baseddf),_issymmetric)
end

Regularize(x::T,y::T,a...;b...) where {T<:Real} = Regularize([x],[y],a...;b...)

Regularize(x::VectorData,a...;b...) = Regularize(x.u,x.v,a...;b...)

function Base.show(io::IO, H::Regularize{N,F}) where {N,F}
    filter = F ? "filtered" : "non-filtered"
    ddftype,_ = typeof(H.ddf).parameters
    op = H._issymmetric ? "Symmetric regularization/interpolation" : "Regularization/interpolation"
    println(io, "$op operator with $filter interpolation")
    println(io, "  DDF type $ddftype")
    println(io, "  $N points in grid with cell area $(sprint(show,1.0/H.overdv;context=:compact => true))")
end

function _get_regularization_radius(ddf::DDF)
  v = 1.0
  r = 0.0
  dr = 0.01
  while (v = abs(ddf(r))) > eps()
    r += dr
  end
  return r
end

_delta_block(radius,shift) = -radius+shift, radius+shift
_index_range(x,xmin,xmax,n,dn) = max(1,ceil(Int,x+xmin)), min(n-dn,floor(Int,x+xmax))
_distance_list(x,imin,imax,shift) = float(imin)-shift-x:float(imax)-shift-x

"""
    RegularizationMatrix(H::Regularize,f::PointData,u::CellData) -> Hmat

Construct and store a matrix representation of regularization associated with `H`
for data of type `f` to data of type `u`. The resulting matrix `Hmat` can then be
used to apply on point data of type `f` to regularize it to grid data of type `u`,
using `mul!(u,Hmat,f)`. It can also be used as just `Hmat*f`.

If `H` is a symmetric regularization and interpolation operator, then this
actually returns a tuple `Hmat, Emat`, where `Emat` is the interpolation matrix.
"""
struct RegularizationMatrix{TU,TF} <: AbstractMatrix{Float64}
  M :: SparseMatrixCSC{Float64,Int64}
end

"""
    InterpolationMatrix(H::Regularize,u::CellData,f::PointData) -> Emat

Construct and store a matrix representation of interpolation associated with `H`
for data of type `u` to data of type `f`. The resulting matrix `Emat` can then be
used to apply on grid data of type `u` to interpolate it to point data of type `f`,
using `mul!(f,Emat,u)`. It can also be used as just `Emat*u`.
"""
struct InterpolationMatrix{TU,TF} <: AbstractMatrix{Float64}
  M :: SparseMatrixCSC{Float64,Int64}
end

@wraparray RegularizationMatrix M 2
@wraparray InterpolationMatrix M 2

# ======  Regularization and interpolation operators of scalar types ======== #

pointtype = :ScalarData
for (gridtype,ctype,dnx,dny,shiftx,shifty) in @generate_scalarlist(SCALARLIST)

# Regularization
  @eval function (H::Regularize{N,F})(target::$gridtype{$ctype,NX,NY,T,DDT},source::$pointtype{N,S,DT}) where {N,F,NX,NY,S,T,DT,DDT}
        radius = H.ddf_radius
        fill!(target,0.0)
        xmin, xmax = _delta_block(radius,$shiftx)
        ymin, ymax = _delta_block(radius,$shifty)
        nzinds = findall(!iszero,source)
        @inbounds for pt in nzinds
            fact = source[pt]*H.wgt[pt]
            minx, maxx = _index_range(H.x[pt],xmin,xmax,NX,$dnx)
            miny, maxy = _index_range(H.y[pt],ymin,ymax,NY,$dny)
            prangex = _distance_list(H.x[pt],minx,maxx,$shiftx)
            prangey = _distance_list(H.y[pt],miny,maxy,$shifty)
            target[minx:maxx,miny:maxy] .+= fact.*H.ddf(prangex,prangey)
        end
        target
  end

# Interpolation
  @eval function (H::Regularize{N,false})(target::$pointtype{N,S,DT},source::$gridtype{$ctype,NX,NY,T,DDT}) where {N,NX,NY,S,T,DT,DDT}
        radius = H.ddf_radius
        fill!(target,0.0)
        xmin, xmax = _delta_block(radius,$shiftx)
        ymin, ymax = _delta_block(radius,$shifty)
        @inbounds for pt in eachindex(target)
            minx, maxx = _index_range(H.x[pt],xmin,xmax,NX,$dnx)
            miny, maxy = _index_range(H.y[pt],ymin,ymax,NY,$dny)
            prangex = _distance_list(H.x[pt],minx,maxx,$shiftx)
            prangey = _distance_list(H.y[pt],miny,maxy,$shifty)
            target[pt] += sum(source[minx:maxx,miny:maxy].*H.ddf(prangex,prangey))
        end
        target
  end

# Interpolation with filtering
  @eval function (H::Regularize{N,true})(target::$pointtype{N,S,DT},source::$gridtype{$ctype,NX,NY,T,DDT}) where {N,NX,NY,S,T,DT,DDT}
        tmp = typeof(source)()
        radius = H.ddf_radius
        fill!(target,0.0)
        xmin, xmax = _delta_block(radius,$shiftx)
        ymin, ymax = _delta_block(radius,$shifty)
        @inbounds for pt in eachindex(target)
            minx, maxx = _index_range(H.x[pt],xmin,xmax,NX,$dnx)
            miny, maxy = _index_range(H.y[pt],ymin,ymax,NY,$dny)
            prangex = _distance_list(H.x[pt],minx,maxx,$shiftx)
            prangey = _distance_list(H.y[pt],miny,maxy,$shifty)
            tmp[minx:maxx,miny:maxy] .+= H.wgt[pt].*H.ddf(prangex,prangey)
        end
        nzinds = findall(x -> abs(x) > eps(),tmp)
        tmp[nzinds] .= inv.(tmp[nzinds])

        @inbounds for pt in eachindex(target)
            minx, maxx = _index_range(H.x[pt],xmin,xmax,NX,$dnx)
            miny, maxy = _index_range(H.y[pt],ymin,ymax,NY,$dny)
            prangex = _distance_list(H.x[pt],minx,maxx,$shiftx)
            prangey = _distance_list(H.y[pt],miny,maxy,$shifty)
            target[pt] += sum(tmp[minx:maxx,miny:maxy].*
                              source[minx:maxx,miny:maxy].*
                              H.ddf(prangex,prangey))
            #end
        end
        target
  end

  # Construct regularization matrix
  @eval function RegularizationMatrix(H::Regularize{N,F},
    f::$pointtype{N,S,DT},
    u::$gridtype{$ctype,NX,NY,T,DDT}) where {N,F,NX,NY,S,T,DT,DDT}

    linI = LinearIndices(u)
    rad = H.ddf_radius
    xmin, xmax = _delta_block(rad,$shiftx)
    ymin, ymax = _delta_block(rad,$shifty)

    rows = Int64[]
    cols = Int64[]
    vals = eltype(u)[]

    for pt = 1:N
      minx, maxx = _index_range(H.x[pt],xmin,xmax,NX,$dnx)
      miny, maxy = _index_range(H.y[pt],ymin,ymax,NY,$dny)
      prangex = _distance_list(H.x[pt],minx,maxx,$shiftx)
      prangey = _distance_list(H.y[pt],miny,maxy,$shifty)

      I1 = vec(linI[minx:maxx,miny:maxy])
      append!(rows,I1)
      append!(cols,fill(pt,length(I1)))
      append!(vals,H.wgt[pt].*vec(H.ddf(prangex,prangey)))
    end
    Hmat = sparse(rows,cols,vals,length(u),length(f))
    if H._issymmetric
      # In symmetric case, these matrices are identical. (Interpolation is stored
      # as its transpose.)
      return RegularizationMatrix{typeof(u),typeof(f)}(Hmat),InterpolationMatrix{typeof(u),typeof(f)}(Hmat)
    else
      return RegularizationMatrix{typeof(u),typeof(f)}(Hmat)
    end
  end

  # Construct interpolation matrix
  @eval function InterpolationMatrix(H::Regularize{N,false},
    u::$gridtype{$ctype,NX,NY,T,DDT},
    f::$pointtype{N,S,DT}) where {N,NX,NY,S,T,DT,DDT}

    linI = LinearIndices(u)
    rad = H.ddf_radius
    xmin, xmax = _delta_block(rad,$shiftx)
    ymin, ymax = _delta_block(rad,$shifty)

    rows = Int64[]
    cols = Int64[]
    vals = eltype(u)[]

    for pt = 1:N
      minx, maxx = _index_range(H.x[pt],xmin,xmax,NX,$dnx)
      miny, maxy = _index_range(H.y[pt],ymin,ymax,NY,$dny)
      prangex = _distance_list(H.x[pt],minx,maxx,$shiftx)
      prangey = _distance_list(H.y[pt],miny,maxy,$shifty)

      I1 = vec(linI[minx:maxx,miny:maxy])
      append!(rows,I1)
      append!(cols,fill(pt,length(I1)))
      append!(vals,vec(H.ddf(prangex,prangey)))
    end
    Emat = sparse(rows,cols,vals,length(u),length(f))
    InterpolationMatrix{typeof(u),typeof(f)}(Emat)
  end

  # Construct interpolation matrix with filtering. Because we construct
  # the matrix by supplying unit vectors from the points, it is easier
  # to carry out the filtering here rather than call the filtered operator
  # we defined above
  @eval function InterpolationMatrix(H::Regularize{N,true},
    u::$gridtype{$ctype,NX,NY,T,DDT},
    f::$pointtype{N,S,DT}) where {N,NX,NY,S,T,DT,DDT}

    linI = LinearIndices(u)
    rad = H.ddf_radius
    xmin, xmax = _delta_block(rad,$shiftx)
    ymin, ymax = _delta_block(rad,$shifty)

    g = similar(f)
    v = similar(u)
    fill!(g,1.0)
    wt = sparsevec(H(v,g))
    wt.nzval .= inv.(wt.nzval)

    rows = Int64[]
    cols = Int64[]
    vals = eltype(u)[]
    for pt = 1:N
      minx, maxx = _index_range(H.x[pt],xmin,xmax,NX,$dnx)
      miny, maxy = _index_range(H.y[pt],ymin,ymax,NY,$dny)
      prangex = _distance_list(H.x[pt],minx,maxx,$shiftx)
      prangey = _distance_list(H.y[pt],miny,maxy,$shifty)

      I1 = vec(linI[minx:maxx,miny:maxy])
      append!(rows,I1)
      append!(cols,fill(pt,length(I1)))
      append!(vals,wt[I1].*vec(H.ddf(prangex,prangey)))
    end
    Emat = sparse(rows,cols,vals,length(u),length(f))

    InterpolationMatrix{typeof(u),typeof(f)}(Emat)
  end

end

# ===== Regularization and interpolation operators of vector data to edges ===== #
pointtype = :VectorData
for (gridtype,ctype,dunx,duny,dvnx,dvny,shiftux,shiftuy,shiftvx,shiftvy) in @generate_collectionlist(VECTORLIST)

# Regularization
  @eval function (H::Regularize{N,F})(target::$gridtype{$(ctype...),NX,NY,T,DDT},source::$pointtype{N,S,DT}) where {N,F,NX,NY,S,T,DT,DDT}
        H(target.u,source.u)
        H(target.v,source.v)
        target
  end

# Interpolation
  @eval function (H::Regularize{N,F})(target::$pointtype{N,S,DT},source::$gridtype{$(ctype...),NX,NY,T,DDT}) where {N,F,NX,NY,S,T,DT,DDT}
        H(target.u,source.u)
        H(target.v,source.v)
        target
  end

  # Construct regularization matrix
  @eval function RegularizationMatrix(H::Regularize{N,F},src::$pointtype{N,S,DT},target::$gridtype{$(ctype...),NX,NY,T,DDT}) where {N,F,NX,NY,S,T,DT,DDT}

    lenu = length(target.u)
    lenv = length(target.v)
    Hmat = spzeros(lenu+lenv,2N)

    if H._issymmetric
      # In symmetric case, these matrices are identical. (Interpolation is stored
      # as its transpose.)
      Hmat[1:lenu,          1:N]    = RegularizationMatrix(H,src.u,target.u)[1].M
      Hmat[lenu+1:lenu+lenv,N+1:2N] = RegularizationMatrix(H,src.v,target.v)[1].M
      return RegularizationMatrix{typeof(target),typeof(src)}(Hmat),InterpolationMatrix{typeof(target),typeof(src)}(Hmat)
    else
      Hmat[1:lenu,          1:N]    = RegularizationMatrix(H,src.u,target.u).M
      Hmat[lenu+1:lenu+lenv,N+1:2N] = RegularizationMatrix(H,src.v,target.v).M
      return RegularizationMatrix{typeof(target),typeof(src)}(Hmat)
    end
  end

  # Construct interpolation matrix
  @eval function InterpolationMatrix(H::Regularize{N,F},src::$gridtype{$(ctype...),NX,NY,T,DDT},target::$pointtype{N,S,DT}) where {N,F,NX,NY,S,T,DT,DDT}

    # note that we store interpolation matrices in the same shape as regularization matrices
    lenu = length(src.u)
    lenv = length(src.v)
    Emat = spzeros(lenu+lenv,2N)
    Emat[1:lenu,          1:N]    = InterpolationMatrix(H,src.u,target.u).M
    Emat[lenu+1:lenu+lenv,N+1:2N] = InterpolationMatrix(H,src.v,target.v).M
    InterpolationMatrix{typeof(src),typeof(target)}(Emat)
  end

end

# ======  Regularization and interpolation operators of tensor data to edge gradients ==== #
# Here, u describes both diagonal components and v the off diagonal
# We do not use the scalar-wise operations here because there is double use of
# the same ddf evaluation, so this saves us quite a bit of time.

pointtype = :TensorData
for (gridtype,ctype,dunx,duny,dvnx,dvny,shiftux,shiftuy,shiftvx,shiftvy) in @generate_collectionlist(TENSORLIST)

# Regularization
  @eval function (H::Regularize{N,F})(target::$gridtype{$(ctype...),NX,NY,T,DDT},source::$pointtype{N,S,DT}) where {N,F,NX,NY,S,T,DT,DDT}
    radius = H.ddf_radius
    fill!(target.dudx,0.0)
    fill!(target.dvdy,0.0)
    xmin, xmax = _delta_block(radius,$shiftux)
    ymin, ymax = _delta_block(radius,$shiftuy)
    @inbounds for pt in 1:N
      minx, maxx = _index_range(H.x[pt],xmin,xmax,NX,$dunx)
      miny, maxy = _index_range(H.y[pt],ymin,ymax,NY,$duny)
      prangex = _distance_list(H.x[pt],minx,maxx,$shiftux)
      prangey = _distance_list(H.y[pt],miny,maxy,$shiftuy)
      tmp = H.wgt[pt].*H.ddf(prangex,prangey)
      target.dudx[minx:maxx,miny:maxy] .+= source.dudx[pt].*tmp
      target.dvdy[minx:maxx,miny:maxy] .+= source.dvdy[pt].*tmp
    end
    fill!(target.dudy,0.0)
    fill!(target.dvdx,0.0)
    xmin, xmax = _delta_block(radius,$shiftvx)
    ymin, ymax = _delta_block(radius,$shiftvy)
    @inbounds for pt in 1:N
      minx, maxx = _index_range(H.x[pt],xmin,xmax,NX,$dvnx)
      miny, maxy = _index_range(H.y[pt],ymin,ymax,NY,$dvny)
      prangex = _distance_list(H.x[pt],minx,maxx,$shiftvx)
      prangey = _distance_list(H.y[pt],miny,maxy,$shiftvy)
      tmp = H.wgt[pt].*H.ddf(prangex,prangey)
      target.dudy[minx:maxx,miny:maxy] .+= source.dudy[pt].*tmp
      target.dvdx[minx:maxx,miny:maxy] .+= source.dvdx[pt].*tmp
    end
    target
  end

  # @eval function (H::Regularize{N,F})(target::$ctype,source::$ftype) where {N,F,NX,NY,S,T}
  #   H(target.dudx,source.dudx)
  #   H(target.dudy,source.dudy)
  #   H(target.dvdx,source.dvdx)
  #   H(target.dvdy,source.dvdy)
  #   target
  # end

# Interpolation
  @eval function (H::Regularize{N,false})(target::$pointtype{N,S,DT},source::$gridtype{$(ctype...),NX,NY,T,DDT}) where {N,NX,NY,S,T,DT,DDT}
    radius = H.ddf_radius
    fill!(target.dudx,0.0)
    fill!(target.dvdy,0.0)
    xmin, xmax = _delta_block(radius,$shiftux)
    ymin, ymax = _delta_block(radius,$shiftuy)
    @inbounds for pt in 1:N
      minx, maxx = _index_range(H.x[pt],xmin,xmax,NX,$dunx)
      miny, maxy = _index_range(H.y[pt],ymin,ymax,NY,$duny)
      prangex = _distance_list(H.x[pt],minx,maxx,$shiftux)
      prangey = _distance_list(H.y[pt],miny,maxy,$shiftuy)
      tmp = H.ddf(prangex,prangey)
      target.dudx[pt] += sum(source.dudx[minx:maxx,miny:maxy].*tmp)
      target.dvdy[pt] += sum(source.dvdy[minx:maxx,miny:maxy].*tmp)
    end
    fill!(target.dudy,0.0)
    fill!(target.dvdx,0.0)
    xmin, xmax = _delta_block(radius,$shiftvx)
    ymin, ymax = _delta_block(radius,$shiftvy)
    @inbounds for pt in 1:N
      minx, maxx = _index_range(H.x[pt],xmin,xmax,NX,$dvnx)
      miny, maxy = _index_range(H.y[pt],ymin,ymax,NY,$dvny)
      prangex = _distance_list(H.x[pt],minx,maxx,$shiftvx)
      prangey = _distance_list(H.y[pt],miny,maxy,$shiftvy)
      tmp = H.ddf(prangex,prangey)
      target.dudy[pt] += sum(source.dudy[minx:maxx,miny:maxy].*tmp)
      target.dvdx[pt] += sum(source.dvdx[minx:maxx,miny:maxy].*tmp)
    end
    target
  end

  # @eval function (H::Regularize{N,F})(target::$ftype,source::$ctype) where {N,F,NX,NY,S,T}
  #   H(target.dudx,source.dudx)
  #   H(target.dudy,source.dudy)
  #   H(target.dvdx,source.dvdx)
  #   H(target.dvdy,source.dvdy)
  #   target
  # end

# Interpolation with filtering -- need to speed up
  @eval function (H::Regularize{N,true})(target::$pointtype{N,S,DT},source::$gridtype{$(ctype...),NX,NY,T,DDT}) where {N,NX,NY,S,T,DT,DDT}
    tmp = typeof(source)()
    radius = H.ddf_radius
    fill!(target.dudx,0.0)
    fill!(target.dvdy,0.0)
    xmin, xmax = _delta_block(radius,$shiftux)
    ymin, ymax = _delta_block(radius,$shiftuy)
    @inbounds for pt in 1:N
        minx, maxx = _index_range(H.x[pt],xmin,xmax,NX,$dunx)
        miny, maxy = _index_range(H.y[pt],ymin,ymax,NY,$duny)
        prangex = _distance_list(H.x[pt],minx,maxx,$shiftux)
        prangey = _distance_list(H.y[pt],miny,maxy,$shiftuy)
        tmp.dudx[minx:maxx,miny:maxy] .+= H.wgt[pt].*H.ddf(prangex,prangey)
    end
    nzinds = findall(x -> abs(x) > eps(),tmp.dudx)
    tmp.dudx[nzinds] .= inv.(tmp.dudx[nzinds])
    tmp.dvdy .= tmp.dudx

    @inbounds for pt in 1:N
        minx, maxx = _index_range(H.x[pt],xmin,xmax,NX,$dunx)
        miny, maxy = _index_range(H.y[pt],ymin,ymax,NY,$duny)
        prangex = _distance_list(H.x[pt],minx,maxx,$shiftux)
        prangey = _distance_list(H.y[pt],miny,maxy,$shiftuy)
        d = H.ddf(prangex,prangey)
        target.dudx[pt] += sum(tmp.dudx[minx:maxx,miny:maxy].*
                               source.dudx[minx:maxx,miny:maxy].*d)
        target.dvdy[pt] += sum(tmp.dvdy[minx:maxx,miny:maxy].*
                               source.dvdy[minx:maxx,miny:maxy].*d)
    end

    fill!(target.dudy,0.0)
    fill!(target.dvdx,0.0)
    xmin, xmax = _delta_block(radius,$shiftvx)
    ymin, ymax = _delta_block(radius,$shiftvy)
    @inbounds for pt in 1:N
        minx, maxx = _index_range(H.x[pt],xmin,xmax,NX,$dvnx)
        miny, maxy = _index_range(H.y[pt],ymin,ymax,NY,$dvny)
        prangex = _distance_list(H.x[pt],minx,maxx,$shiftvx)
        prangey = _distance_list(H.y[pt],miny,maxy,$shiftvy)
        tmp.dudy[minx:maxx,miny:maxy] .+= H.wgt[pt].*H.ddf(prangex,prangey)
    end
    nzinds = findall(x -> abs(x) > eps(),tmp.dudy)
    tmp.dudy[nzinds] .= inv.(tmp.dudy[nzinds])
    tmp.dvdx .= tmp.dudy

    @inbounds for pt in 1:N
        minx, maxx = _index_range(H.x[pt],xmin,xmax,NX,$dvnx)
        miny, maxy = _index_range(H.y[pt],ymin,ymax,NY,$dvny)
        prangex = _distance_list(H.x[pt],minx,maxx,$shiftvx)
        prangey = _distance_list(H.y[pt],miny,maxy,$shiftvy)
        d = H.ddf(prangex,prangey)
        target.dudy[pt] += sum(tmp.dudy[minx:maxx,miny:maxy].*
                               source.dudy[minx:maxx,miny:maxy].*d)
        target.dvdx[pt] += sum(tmp.dvdx[minx:maxx,miny:maxy].*
                               source.dvdx[minx:maxx,miny:maxy].*d)
    end
    target
  end

  # Construct regularization matrix
  @eval function RegularizationMatrix(H::Regularize{N,F},src::$pointtype{N,S,DT},target::$gridtype{$(ctype...),NX,NY,T,DDT}) where {N,F,NX,NY,S,T,DT,DDT}

    # note that we only need to compute two distinct matrices, since there are
    # only two types of cell data in this tensor

    lenu = length(target.dudx)
    lenv = length(target.dudy)
    Hmat = spzeros(2lenu+2lenv,4N)

    if H._issymmetric
      # In symmetric case, these matrices are identical. (Interpolation is stored
      # as its transpose.)
      Hdudx, _ = RegularizationMatrix(H,src.dudx,target.dudx)
      Hdudy, _ = RegularizationMatrix(H,src.dudy,target.dudy)
      Hmat[1:lenu,          1:N]    = Hdudx.M
      Hmat[lenu+1:lenu+lenv,N+1:2N] = Hdudy.M
      Hmat[lenu+lenv+1:lenu+2lenv,2N+1:3N] = Hdudy.M
      Hmat[lenu+2lenv+1:2lenu+2lenv,3N+1:4N] = Hdudx.M
      return RegularizationMatrix{typeof(target),typeof(src)}(Hmat),InterpolationMatrix{typeof(target),typeof(src)}(Hmat)
    else
      Hdudx = RegularizationMatrix(H,src.dudx,target.dudx)
      Hdudy = RegularizationMatrix(H,src.dudy,target.dudy)
      Hmat[1:lenu,          1:N]    = Hdudx.M
      Hmat[lenu+1:lenu+lenv,N+1:2N] = Hdudy.M
      Hmat[lenu+lenv+1:lenu+2lenv,2N+1:3N] = Hdudy.M
      Hmat[lenu+2lenv+1:2lenu+2lenv,3N+1:4N] = Hdudx.M
      return RegularizationMatrix{typeof(target),typeof(src)}(Hmat)
    end
  end

  # Construct interpolation matrix
  @eval function InterpolationMatrix(H::Regularize{N,F},src::$gridtype{$(ctype...),NX,NY,T,DDT},target::$pointtype{N,S,DT}) where {N,F,NX,NY,S,T,DT,DDT}

    # note that we store interpolation matrices in the same shape as regularization matrices
    lenu = length(src.dudx)
    lenv = length(src.dudy)
    Emat = spzeros(2lenu+2lenv,4N)
    Edudx = InterpolationMatrix(H,src.dudx,target.dudx)
    Edudy = InterpolationMatrix(H,src.dudy,target.dudy)

    Emat[1:lenu,          1:N]    = Edudx.M
    Emat[lenu+1:lenu+lenv,N+1:2N] = Edudy.M
    Emat[lenu+lenv+1:lenu+2lenv,2N+1:3N] = Edudy.M
    Emat[lenu+2lenv+1:2lenu+2lenv,3N+1:4N] = Edudx.M

    InterpolationMatrix{typeof(src),typeof(target)}(Emat)
  end

end

###### Matrix multiplication extended to these Regularization/Interpolation matrices ########

# Regularize
# We need to restrict the action of the matrices to the point and grid types
# they are built for, but allow the underyling data type (DT parameter)
# to be different, since one might be, e.g., a SubArray and the other an Array.
for f in [:Nodes,:XEdges,:YEdges]
    @eval mul!(u::S1,Hmat::RegularizationMatrix{S2,F},f) where {F,S1<:$f{C,NX,NY,T},S2<:$f{C,NX,NY,T}} where {C,NX,NY,T} = _mul!(u,Hmat,f)
    @eval mul!(f,Emat::InterpolationMatrix{S2,F},u::S1) where {F,S1<:$f{C,NX,NY,T},S2<:$f{C,NX,NY,T}} where {C,NX,NY,T} = _mul!(f,Emat,u)
    @eval (*)(Emat::InterpolationMatrix{S2,F},u::S1) where {F,S1<:$f{C,NX,NY,T},S2<:$f{C,NX,NY,T}} where {C,NX,NY,T} = mul!(F(),Emat,u)
    # This is meant to generate a MethodError for non-matching point types:
    @eval (*)(Emat::InterpolationMatrix{S2,F},u::GridData) where {F,S2<:$f{C,NX,NY,T}} where {C,NX,NY,T} = mul!(F(),Emat,u)
end
# Handle all collected grid data in a stricter fashion, since underlying data
# is always stored as Vector type
mul!(u::G,Hmat::RegularizationMatrix{G,F},f) where {G <: CollectedGridData,F} = _mul!(u,Hmat,f)
mul!(f,Emat::InterpolationMatrix{G,F},u::G) where {G <: CollectedGridData,F} = _mul!(f,Emat,u)
(*)(Emat::InterpolationMatrix{G,F},u::G) where {F,G <: CollectedGridData} = mul!(F(),Emat,u)
# This is meant to generate a MethodError for non-matching point types:
(*)(Emat::InterpolationMatrix{G,F},u) where {F,G <: CollectedGridData} = mul!(F(),Emat,u)

# Now dispatch on the point data type
for f in [:ScalarData,:VectorData,:TensorData]
  @eval _mul!(u,Hmat::RegularizationMatrix{G,S1},f::S2) where {G,S1<:$f{N,T},S2<:$f{N,T}} where {N,T} = _unsafe_mul!(u,Hmat,f)
  @eval _mul!(f::S2,Emat::InterpolationMatrix{G,S1},u) where {G,S1<:$f{N,T},S2<:$f{N,T}} where {N,T} = _unsafe_mul!(f,Emat,u)
  @eval (*)(Hmat::RegularizationMatrix{G,S1},f::S2) where {G,S1<:$f{N,T},S2<:$f{N,T}} where {N,T} = mul!(G(),Hmat,f)
  # This is meant to generate a MethodError for non-matching point types:
  @eval (*)(Hmat::RegularizationMatrix{G,S1},f) where {G,S1<:$f{N,T}} where {N,T} = mul!(G(),Hmat,f)
end

# Regularization without checking types
function _unsafe_mul!(u,Hmat::RegularizationMatrix,f)
  fill!(u,0.0)
  nzv = Hmat.M.nzval
  rv = Hmat.M.rowval
  @inbounds for col = 1:Hmat.M.n
    fj = f[col]
    for j = Hmat.M.colptr[col]:(Hmat.M.colptr[col + 1] - 1)
        u[rv[j]] += nzv[j]*fj
    end
  end
  u
end

# Interpolation without checking types
function _unsafe_mul!(f,Emat::InterpolationMatrix,u)
  fill!(f,0.0)
  nzv = Emat.M.nzval
  rv = Emat.M.rowval
  @inbounds for col = 1:Emat.M.n
      tmp = zero(eltype(f))
      for j = Emat.M.colptr[col]:(Emat.M.colptr[col + 1] - 1)
          tmp += transpose(nzv[j])*u[rv[j]]
      end
      f[col] += tmp
  end
  f
end

# Interpolation of regularization, used for developing the filtering matrix
function mul!(C::Array{Float64},Emat::InterpolationMatrix{G,F},
                              Hmat::RegularizationMatrix{G,F}) where {G<:GridData,F<:PointData}
  fill!(C,0.0)
  Enzv = Emat.M.nzval
  Erv = Emat.M.rowval
  @inbounds for row = 1:Emat.M.n, col = 1:Hmat.M.n
      tmp = zero(eltype(C))
      for j = Emat.M.colptr[row]:(Emat.M.colptr[row + 1] - 1)
          tmp += transpose(Enzv[j])*Hmat[Erv[j],col]
      end
      C[row,col] += tmp
  end
  return C
end

(*)(Emat::InterpolationMatrix,Hmat::RegularizationMatrix) =
        mul!(Array{eltype(Emat),2}(undef,Emat.M.n,Hmat.M.n),Emat,Hmat)

# This fixes a method ambiguity
(*)(Emat::InterpolationMatrix{G,F},Hmat::RegularizationMatrix) where {G <: Union{VectorGridData,TensorGridData},F}=
        mul!(Array{eltype(Emat),2}(undef,Emat.M.n,Hmat.M.n),Emat,Hmat)

function Base.summary(io::IO, H::RegularizationMatrix{TU,TF}) where {TU,TF}
    print(io, "Regularization matrix acting on type $TF and returning type $TU")
end

function Base.summary(io::IO, H::InterpolationMatrix{TU,TF}) where {TU,TF}
    print(io, "Interpolation matrix acting on type $TU and returning type $TF")
end

max_stars_repo_name: JuliaIBPM/CartesianGrids

max_issues_repo_name: JuliaIBPM/CartesianGrids

max_forks_repo_name: JuliaIBPM/CartesianGrids