text:
using Oceananigans
using Oceananigans.Units

using Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities:
    CATKEVerticalDiffusivity,
    SurfaceTKEFlux,
    MixingLength

using OceanTurbulenceParameterEstimation.Observations: FieldTimeSeriesCollector
using JLD2
using GLMakie
using Printf
using Statistics

dir = "six_day_suite"
name = "weak_wind_strong_cooling"
xz_filepath = joinpath(dir, name * "_xz_slice.jld2")
yz_filepath = joinpath(dir, name * "_yz_slice.jld2")
xy_filepath = joinpath(dir, name * "_xy_slice.jld2")
statistics_filepath = joinpath(dir, name * "_instantaneous_statistics.jld2")

#####
##### Metadata
#####

xy_file = jldopen(xy_filepath)

Nx = xy_file["grid/Nx"]
Ny = xy_file["grid/Ny"]
Nz = xy_file["grid/Nz"]
Hx = xy_file["grid/Hx"]
Hy = xy_file["grid/Hy"]
Hz = xy_file["grid/Hz"]
Lx = xy_file["grid/Lx"]
Ly = xy_file["grid/Ly"]
Lz = xy_file["grid/Lz"]

iterations = parse.(Int, keys(xy_file["timeseries/t"]))
times = [xy_file["timeseries/t/$i"] for i in iterations]
Nt = length(times)

Qᵀ = xy_file["parameters/temperature_flux"]
Qᵘ = xy_file["parameters/momentum_flux"]
f  = xy_file["parameters/coriolis_parameter"]
α  = xy_file["parameters/thermal_expansion_coefficient"]
f  = xy_file["parameters/coriolis_parameter"]
dTdz = xy_file["parameters/dθdz_deep"]

close(xy_file)

H = Hz
chop(a::AbstractArray{<:Any, 1}) = a[1+H:end-H]
chop(a::AbstractArray{<:Any, 2}) = a[1+H:end-H, 1+H:end-H]

function extract_slices(filepath; dims, name="w")
    file = jldopen(filepath)
    slices = [chop(dropdims(file["timeseries/$name/$i"]; dims)) for i in iterations]
    close(file)
    return slices
end

#####
##### Construct 1D model
#####

# Data
u_avg_series = extract_slices(statistics_filepath, dims=(1, 2), name="u")
v_avg_series = extract_slices(statistics_filepath, dims=(1, 2), name="v")
T_avg_series = extract_slices(statistics_filepath, dims=(1, 2), name="T")

observations = SyntheticObservations(statistics_filepath, field_names=(:u, :v, :T, :e))

catke_mixing_length = MixingLength(Cᴬu=0.0, Cᴬc=0.0, Cᴬe=0.0,
                                   Cᴷuʳ=0.0, Cᴷcʳ=0.0, Cᴷeʳ=0.0)

catke = CATKEVerticalDiffusivity(mixing_length=catke_mixing_length)

simulation = ensemble_column_model_simulation(observations;
                                              Nensemble = 100,
                                              architecture = CPU(),
                                              equation_of_state = LinearEquationOfState(; α),
                                              buoyancy = SeawaterBuoyancy(; equation_of_state, constant_salinity=true),
                                              tracers = (:T, :e),
                                              closure = catke)

simulation.Δt = 10.0

Qᵘ = simulation.model.velocities.u.boundary_conditions.top.condition
Qᵀ = simulation.model.tracers.T.boundary_conditions.top.condition
dTdz = simulation.model.tracers.T.boundary_conditions.bottom.condition

observations = [observations]
for (case, obs) in enumerate(observations)
    @show case cases[case]
    @show obs.metadata.parameters.momentum_flux
    @show obs.metadata.parameters.temperature_flux
    @show obs.metadata.parameters.buoyancy_flux
    @show f = obs.metadata.parameters.coriolis_parameter

    view(Qᵘ, :, case) .= obs.metadata.parameters.momentum_flux
    view(Qᵀ, :, case) .= obs.metadata.parameters.temperature_flux
    view(dTdz, :, case) .= obs.metadata.parameters.dθdz_deep
    view(simulation.model.coriolis, :, case) .= Ref(FPlane(f=f))
end

set!(simulation.model, observations)

times = round.(times)
collector = FieldTimeSeriesCollector(merge(model.velocities, model.tracers), times)
simulation.callbacks[:collector] = Callback(collector, SpecifiedTimes(times...))

progress(sim) = @info "Sim time: " * prettytime(sim) * ", wall time: " * prettytime(sim.run_wall_time)
simulation.callbacks[:progress] = Callback(progress, IterationInterval(100))

@info "Running 1D simulation..."
start = time_ns()
run!(simulation)
elapsed = 1e-9 * (time_ns() - start)
@info "    ... done (" * prettytime(elapsed) * ")."
=#

#####
##### 3D visualization
#####

grid = RectilinearGrid(size=(Nx, Ny, Nz), halo=(Hx, Hy, Hz), x=(0, Lx), y=(0, Ly), z=(-Lz, 0))

u_yz_series = extract_slices(yz_filepath, dims=1, name="u")
u_xz_series = extract_slices(xz_filepath, dims=2, name="u")
u_xy_series = extract_slices(xy_filepath, dims=3, name="u")

w_yz_series = extract_slices(yz_filepath, dims=1, name="w")
w_xz_series = extract_slices(xz_filepath, dims=2, name="w")
w_xy_series = extract_slices(xy_filepath, dims=3, name="w")

T_yz_series = extract_slices(yz_filepath, dims=1, name="T")
T_xz_series = extract_slices(xz_filepath, dims=2, name="T")
T_xy_series = extract_slices(xy_filepath, dims=3, name="T")

x, y, z = nodes((Center, Center, Center), grid)

x_xz = repeat(x, 1, Nz)
y_xz = 0.995 * Ly * ones(Nx, Nz)
z_xz = repeat(reshape(z, 1, Nz), Nx, 1)

x_yz = 0.995 * Lx * ones(Ny, Nz)
y_yz = repeat(y, 1, Nz)
z_yz = repeat(reshape(z, 1, Nz), grid.Ny, 1)

# Slight displacements to "stitch" the cube together
x_xy = x
y_xy = y
z_xy = - 0.001 * Lz * ones(grid.Nx, grid.Ny)

#####
##### Animate!
#####

fig = Figure(resolution=(1400, 1200))

n = Observable(1)

###
### Vertical velocity
###

w_yz = @lift w_yz_series[$n]
w_xz = @lift w_xz_series[$n]
w_xy = @lift w_xy_series[$n]

u_yz = @lift u_yz_series[$n]
u_xz = @lift u_xz_series[$n]
u_xy = @lift u_xy_series[$n]

T_yz = @lift T_yz_series[$n]
T_xz = @lift T_xz_series[$n]
T_xy = @lift T_xy_series[$n]

u_catke = @lift interior(collector.field_time_serieses.u[$n])[1, 1, :]
v_catke = @lift interior(collector.field_time_serieses.v[$n])[1, 1, :]
T_catke = @lift interior(collector.field_time_serieses.T[$n])[1, 1, :]

u_LES = @lift u_avg_series[$n]
v_LES = @lift v_avg_series[$n]
T_LES = @lift T_avg_series[$n]

colormap_T = :oslo
colormap_u = :balance
colormap_w = :balance

w_max = maximum(maximum.(abs, w_xz_series[end]))
colorrange_w = (-w_max/2, w_max/2)

u_max = maximum(maximum.(abs, u_xz_series[end]))
colorrange_u = (-u_max/2, u_max/2)

# Compute maximum temperature for colorrange with a moving average
ΔT = 0.04 # Absolute range
Δn = 2 # Half-width of moving average
Tf_max = maximum(T_xy_series[end])
T0_max = maximum(T_xy_series[1])

colorrange_T = @lift begin
    if $n < Δn + 1
        T_max = T0_max
    elseif $n < Nt - Δn
        T_max = mean([maximum(T_xy_series[nn]) for nn in $n-Δn:$n+Δn])
    else
        T_max = Tf_max
    end

    T_min = T_max - ΔT

    (T_min, T_max)
end

ax_T = Axis3(fig[1:8,  2], aspect=:data, xlabel="x (m)", ylabel="y (m)", zlabel="z (m)", azimuth=0.86, elevation=0.38)
ax_u = Axis3(fig[9:16, 2], aspect=:data, xlabel="x (m)", ylabel="y (m)", zlabel="z (m)", azimuth=0.86, elevation=0.38)

#sfc_w = surface!(ax_w, x_xz, y_xz, z_xz; color=w_xz, colormap=colormap_w, colorrange=colorrange_w)
#        surface!(ax_w, x_yz, y_yz, z_yz; color=w_yz, colormap=colormap_w, colorrange=colorrange_w)
#        surface!(ax_w, x_xy, y_xy, z_xy; color=w_xy, colormap=colormap_w, colorrange=colorrange_w)

sfc_u = surface!(ax_u, x_xz, y_xz, z_xz; color=u_xz, colormap=colormap_u, colorrange=colorrange_u)
        surface!(ax_u, x_yz, y_yz, z_yz; color=u_yz, colormap=colormap_u, colorrange=colorrange_u)
        surface!(ax_u, x_xy, y_xy, z_xy; color=u_xy, colormap=colormap_u, colorrange=colorrange_u)

sfc_T = surface!(ax_T, x_xz, y_xz, z_xz; color=T_xz, colormap=colormap_T, colorrange=colorrange_T)
        surface!(ax_T, x_yz, y_yz, z_yz; color=T_yz, colormap=colormap_T, colorrange=colorrange_T)
        surface!(ax_T, x_xy, y_xy, z_xy; color=T_xy, colormap=colormap_T, colorrange=colorrange_T)

cp_T = fig[3:6,   1] = Colorbar(fig, sfc_T, flipaxis=false, label="Temperature (ᵒC)")
cp_u = fig[11:14, 1] = Colorbar(fig, sfc_u, flipaxis=false, label="u (x-velocity component, m s⁻¹)")

ax_T_avg = Axis(fig[3:7, 3],
                xlabel = "Horizontally-averaged \n temperature (ᵒC)",
                ylabel = "z (m)",
                yaxisposition = :right)

ax_u_avg = Axis(fig[11:15, 3],
                xlabel = "Horizontally-averaged velocity \n components (m s⁻¹)",
                ylabel = "z (m)",
                yaxisposition = :right)

T_catke_lines = lines!(ax_T_avg, T_catke, z, label="CATKE turbulence model")
T_LES_lines = lines!(ax_T_avg, T_LES, z, label="Large eddy simulation data")

model_color = T_catke_lines.attributes.color.val
LES_color = T_LES_lines.attributes.color.val

lines!(ax_u_avg, u_catke, z, color=model_color, label="u")
lines!(ax_u_avg, v_catke, z, color=model_color, linestyle=:dash, label="v")

lines!(ax_u_avg, u_LES, z, color=LES_color, label="u")
lines!(ax_u_avg, v_LES, z, color=LES_color, linestyle=:dash, label="v")

xlims!(ax_u_avg, -0.2, 0.2)

Legend(fig[8:9, 3], ax_T_avg)
axislegend(ax_u_avg, position=:rb, merge=true)

title = @lift "Ocean surface boundary layer at t = " * prettytime(times[$n])
lab = Label(fig[0, :], title, textsize=18)

record(fig, joinpath(dir, "$name.mp4"), 1:Nt; framerate=12) do nn
    @info "Drawing frame $nn of $Nt..."
    n[] = nn
end

max_stars_repo_name: glwagner/LocalOceanClosureCalibration


text:
using Test, MPI
using ClimateMachine.Mesh.Topologies
using ClimateMachine.Mesh.Grids
using ClimateMachine.Mesh.Geometry
using StaticArrays

MPI.Initialized() || MPI.Init()

@testset "LocalGeometry" begin
    FT = Float64
    ArrayType = Array

    xmin = 0
    ymin = 0
    zmin = 0
    xmax = 2000
    ymax = 400
    zmax = 2000

    Ne = (20, 2, 20)

    polynomialorder = 4

    brickrange = (
        range(FT(xmin); length = Ne[1] + 1, stop = xmax),
        range(FT(ymin); length = Ne[2] + 1, stop = ymax),
        range(FT(zmin); length = Ne[3] + 1, stop = zmax),
    )
    topl = StackedBrickTopology(
        MPI.COMM_SELF,
        brickrange,
        periodicity = (true, true, false),
    )

    grid = DiscontinuousSpectralElementGrid(
        topl,
        FloatType = FT,
        DeviceArray = ArrayType,
        polynomialorder = polynomialorder,
    )

    S = (
        (xmax - xmin) / (polynomialorder * Ne[1]),
        (ymax - ymin) / (polynomialorder * Ne[2]),
        (zmax - zmin) / (polynomialorder * Ne[3]),
    )
    Savg = cbrt(prod(S))
    M = SDiagonal(S .^ -2)

    for e in 1:size(grid.vgeo, 3)
        for n in 1:size(grid.vgeo, 1)
            g = LocalGeometry(Val(polynomialorder), grid.vgeo, n, e)
            @test lengthscale(g) ≈ Savg
            @test Geometry.resolutionmetric(g) ≈ M
        end
    end
end

ext: jl


text:
abstract type ClapeyronParam end

"""
    SingleParam{T}

Struct designed to contain single parameters. Basically a vector with some extra info.

## Creation:
```julia-repl
julia> mw = SingleParam("molecular weight",["water","ammonia"],[18.01,17.03])
SingleParam{Float64}("molecular weight") with 2 components:
 "water" => 18.01
 "ammonia" => 17.03

julia> mw.values
2-element Vector{Float64}:
 18.01
 17.03

julia> mw.components
2-element Vector{String}:
 "water"
 "ammonia"

julia> mw2 = SingleParam(mw,"new name")
SingleParam{Float64}("new name") with 2 components:
 "water" => 18.01
 "ammonia" => 17.03

julia> has_oxigen = [true,false]; has_o = SingleParam(mw2,has_oxigen)
SingleParam{Bool}("new name") with 2 components:
 "water" => true
 "ammonia" => false

```

## Example usage in models:

```
function molecular_weight(model,molar_frac)
    mw = model.params.mw.values
    res = zero(eltype(molarfrac))
    for i in @comps #iterating through all components
        res += molar_frac[i]*mw[i]
    end
    return res
end
```
"""
struct SingleParam{T} <: ClapeyronParam
    name::String
    components::Array{String,1}
    values::Array{T,1}
    ismissingvalues::Array{Bool,1}
    sourcecsvs::Array{String,1}
    sources::Array{String,1}
end

function Base.show(io::IO, param::SingleParam)
    print(io, typeof(param), "(\"", param.name, "\")[")
    for component in param.components
        component != first(param.components) && print(io, ",")
        print(io, "\"", component, "\"")
    end
    print(io, "]")
end

function Base.show(io::IO, ::MIME"text/plain", param::SingleParam)
    len = length(param.values)
    print(io, typeof(param), "(\"", param.name)
    println(io, "\") with ", len, " component", ifelse(len==1, ":", "s:"))
    i = 0
    for (name, val, miss) in zip(param.components, param.values, param.ismissingvalues)
        i += 1
        if i > 1
            println(io)
        end
        if miss == false
            if typeof(val) <: AbstractString
                print(io, " \"", name, "\" => \"", val, "\"")
            else
                print(io, " \"", name, "\" => ", val)
            end
        else
            print(io, " \"", name, " => -")
        end
    end
end

function SingleParam(x::SingleParam{T},name=x.name) where T
    return SingleParam(name, x.components,deepcopy(x.values), deepcopy(x.ismissingvalues), x.sourcecsvs, x.sources)
end

#a barebones constructor, in case we dont build from csv
function SingleParam(
    name::String,
    components::Vector{String},
    values::Vector{T},
    sourcecsvs = String[],
    sources = String[]
    ;default = _zero(T)) where T
    _values,_ismissingvalues = defaultmissing(values,default)
    TT = eltype(_values)
    return  SingleParam{TT}(name,components, _values, _ismissingvalues, sourcecsvs, sources)
end

function SingleParam(x::SingleParam, v::Vector)
    _values,_ismissingvalues = defaultmissing(v)
    return SingleParam(x.name, x.components,_values, _ismissingvalues , x.sourcecsvs, x.sources)
end
"""
    PairParam{T}

Struct designed to contain pair data. used a matrix as underlying data storage.

## Creation:
```julia-repl
julia> kij = PairParam("interaction params",["water","ammonia"],[0.1 0.0;0.1 0.0])
PairParam{Float64}["water", "ammonia"]) with values:
2×2 Matrix{Float64}:
 0.1  0.0
 0.1  0.0

julia> kij.values
2×2 Matrix{Float64}:
 0.1  0.0
 0.1  0.0

julia> kij.diagvalues
2-element view(::Vector{Float64},
1:3:4) with eltype Float64:
 0.1
 0.0
```

## Example usage in models:

```julia
#lets compute ∑xᵢxⱼkᵢⱼ
function alpha(model,x)
    kij = model.params.kij.values
    ki = model.params.kij.diagvalues
    res = zero(eltype(molarfrac))
    for i in @comps
        @show ki[i] #diagonal values
        for j in @comps
            res += x[i]*x[j]*kij[i,j]
        end
    end
    return res
end
```
"""
struct PairParam{T} <: ClapeyronParam
    name::String
    components::Array{String,1}
    values::Array{T,2}
    diagvalues::SubArray{T, 1, Vector{T}, Tuple{StepRange{Int64, Int64}}, true}
    ismissingvalues::Array{Bool,2}
    sourcecsvs::Array{String,1}
    sources::Array{String,1}
end

function PairParam(name::String,
                    components::Array{String,1},
                    values::Array{T,2},
                    sourcecsvs::Array{String,1} = String[],
                    sources::Array{String,1} = String[]) where T

    _values,_ismissingvalues = defaultmissing(values)
    diagvalues = view(_values, diagind(_values))

    return PairParam{T}(name, components,_values, diagvalues, _ismissingvalues, sourcecsvs, sources)
end

function PairParam(x::PairParam,name::String=x.name)
    values = deepcopy(x.values)
    diagvalues = view(values,diagind(values))
    return PairParam(name, x.components,values ,diagvalues, deepcopy(x.ismissingvalues), x.sourcecsvs, x.sources)
end

function PairParam(x::SingleParam,name::String=x.name)
    pairvalues = singletopair(x.values,missing)
    for i in 1:length(x.values)
        if x.ismissingvalues[i]
            pairvalues[i,i] = missing
        end
    end
    _values,_ismissingvalues = defaultmissing(pairvalues)
    diagvalues = view(_values, diagind(_values))
    return PairParam(name, x.components, _values,diagvalues,_ismissingvalues,x.sourcecsvs, x.sources)
end

function PairParam(x::PairParam, v::Matrix,name::String=x.name)
    return PairParam(name, x.components,deepcopy(v), x.sourcecsvs, x.sources)
end
function PairParam(x::SingleParam, v::Vector,name::String=x.name)
    pairvalues = singletopair(v,missing)
    return PairParam(x.name, x.components, pairvalues,x.sourcecsvs, x.sources)
end
function PairParam(x::SingleParam, v::Matrix,name::String=x.name)
    return PairParam(x.name, x.components, deepcopy(v),x.sourcecsvs, x.sources)
end

#barebones constructor by list of pairs.

function Base.show(io::IO,mime::MIME"text/plain",param::PairParam)
    print(io,"PairParam{",string(typeof(param)),"}")
    show(io,param.components)
    println(io,") with values:")
    show(io,mime,param.values)
end

function Base.show(io::IO,param::PairParam)
    print(io,"PairParam(",param.name)
    print(io,"\"",param.name,"\"",")[")
    for (name,val,miss,i) in zip(param.components,param.values,param.ismissingvalues,1:length(param.values))
        i != 1 && print(io,",")
        if miss == false
            print(io,name,"=",val)
        else
            print(io,name,"=","-")
        end
    end
    print(io,"]")
end
"""
    AssocParam{T}

Struct holding association parameters.
"""
struct AssocParam{T} <: ClapeyronParam
    name::String
    components::Array{String,1}
    values::Array{Array{T,2},2}
    ismissingvalues::Array{Array{Bool,2},2}
    allcomponentsites::Array{Array{String,1},1}
    sourcecsvs::Array{String,1}
    sources::Array{String,1}
end

function AssocParam(x::AssocParam{T}) where T
    return PairParam{T}(x.name,x.components, deepcopy(x.values), deepcopy(x.ismissingvalues), x.allcomponentsites, x.sourcecsvs, x.sources)
end

function AssocParam{T}(x::AssocParam, v::Matrix{Matrix{T}}) where T
    return AssocParam{T}(x.name, x.components,deepcopy(v), deepcopy(x.ismissingvalues), x.allcomponentsites, x.sourcecsvs, x.sources)
end

function Base.show(io::IO,mime::MIME"text/plain",param::AssocParam{T}) where T
    print(io,"AssocParam{",string(T),"}")
    show(io,param.components)
    println(io,") with values:")
    show(io,mime,param.values)
end

function Base.show(io::IO,param::AssocParam)
    print(io,"AssocParam(",param.name)
    print(io,"\"",param.name,"\"",")[")
    for (name,val,miss,i) in zip(param.components,param.values,param.ismissingvalues,1:length(param.values))
        i != 1 && print(io,",")
        if miss == false
            print(io,name,"=",val)
        else
            print(io,name,"=","-")
        end
    end
    print(io,"]")
end
const PARSED_GROUP_VECTOR_TYPE =  Vector{Tuple{String, Vector{Pair{String, Int64}}}}

"""
    GroupParam

Struct holding group parameters.contains:
* `components`: a list of all components
* `groups`: a list of groups names for each component
* `i_groups`: a list containing the number of groups for each component
* `n_groups`: a list of the group multiplicity of each group corresponding to each group in `i_groups`
* `flattenedgroups`: a list of all unique groups--the parameters correspond to this list
* `n_flattenedgroups`: the group multiplicities corresponding to each group in `flattenedgroups`
* `i_flattenedgroups`: an iterator that goes through the indices for each flattened group

You can create a group param by passing a `Vector{Tuple{String, Vector{Pair{String, Int64}}}}.
For example:
```julia-repl
julia> grouplist = [
           ("ethanol", ["CH3"=>1, "CH2"=>1, "OH"=>1]),
           ("nonadecanol", ["CH3"=>1, "CH2"=>18, "OH"=>1]),
           ("ibuprofen", ["CH3"=>3, "COOH"=>1, "aCCH"=>1, "aCCH2"=>1, "aCH"=>4])];

julia> groups = GroupParam(grouplist)
GroupParam with 3 components:
 "ethanol": "CH3" => 1, "CH2" => 1, "OH" => 1
 "nonadecanol": "CH3" => 1, "CH2" => 18, "OH" => 1
 "ibuprofen": "CH3" => 3, "COOH" => 1, "aCCH" => 1, "aCCH2" => 1, "aCH" => 4

julia> groups.flattenedgroups
7-element Vector{String}:
 "CH3"
 "CH2"
 "OH"
 "COOH"
 "aCCH"
 "aCCH2"
 "aCH"

julia> groups.i_groups
3-element Vector{Vector{Int64}}:
 [1, 2, 3]
 [1, 2, 3]
 [1, 4, 5, 6, 7]

julia> groups.n_groups
3-element Vector{Vector{Int64}}:
 [1, 1, 1]
 [1, 18, 1]
 [3, 1, 1, 1, 4]

julia> groups.n_flattenedgroups
 3-element Vector{Vector{Int64}}:
 [1, 1, 1, 0, 0, 0, 0]
 [1, 18, 1, 0, 0, 0, 0]
 [3, 0, 0, 1, 1, 1, 4]
```

if you have CSV with group data, you can also pass those, to automatically query the missing groups in your input vector:

```julia-repl
julia> grouplist = [
           "ethanol",
           ("nonadecanol", ["CH3"=>1, "CH2"=>18, "OH"=>1]),
           ("ibuprofen", ["CH3"=>3, "COOH"=>1, "aCCH"=>1, "aCCH2"=>1, "aCH"=>4])];

           julia> groups = GroupParam(grouplist, ["SAFT/SAFTgammaMie/SAFTgammaMie_groups.csv"])
           GroupParam with 3 components:
            "ethanol": "CH2OH" => 1, "CH3" => 1
            "nonadecanol": "CH3" => 1, "CH2" => 18, "OH" => 1
            "ibuprofen": "CH3" => 3, "COOH" => 1, "aCCH" => 1, "aCCH2" => 1, "aCH" => 4
```
In this case, `SAFTGammaMie` files support the second order group `CH2OH`.

"""
struct GroupParam <: ClapeyronParam
    components::Array{String,1}
    groups::Array{Array{String,1},1}
    n_groups::Array{Array{Int,1},1}
    i_groups::Array{Array{Int,1},1}
    flattenedgroups::Array{String,1}
    n_flattenedgroups::Array{Array{Int,1},1}
    i_flattenedgroups::UnitRange{Int}
    sourcecsvs::Array{String,1}
end

function GroupParam(input::PARSED_GROUP_VECTOR_TYPE,sourcecsvs::Vector{String}=String[],options::ParamOptions = ParamOptions())
    components = [first(i) for i ∈ input]
    raw_groups =  [last(i) for i ∈ input]
    groups = [first.(grouppairs) for grouppairs ∈ raw_groups]
    n_groups = [last.(grouppairs) for grouppairs ∈ raw_groups]
    flattenedgroups = unique!(reduce(vcat,groups))
    i_groups = [[findfirst(isequal(group), flattenedgroups) for group ∈ componentgroups] for componentgroups ∈ groups]
    len_flattenedgroups = length(flattenedgroups)
    i_flattenedgroups = 1:len_flattenedgroups
    n_flattenedgroups = [zeros(Int,len_flattenedgroups) for _ ∈ 1:length(input)]
    for i in length(input)
        setindex!.(n_flattenedgroups,n_groups,i_groups)
    end

    return GroupParam(components,
    groups,
    n_groups,
    i_groups,
    flattenedgroups,
    n_flattenedgroups,
    i_flattenedgroups,
    sourcecsvs)
end

function Base.show(io::IO, mime::MIME"text/plain", param::GroupParam)
    print(io,"GroupParam ")
    len = length(param.components)
    println(io,"with ", len, " component", ifelse(len==1, ":", "s:"))

    for i in 1:length(param.components)

        print(io, " \"", param.components[i], "\": ")
        firstloop = true
        for j in 1:length(param.n_groups[i])
            firstloop == false && print(io, ", ")
            print(io, "\"", param.groups[i][j], "\" => ", param.n_groups[i][j])
            firstloop = false
        end
        i != length(param.components) && println(io)
    end
end

function Base.show(io::IO, param::GroupParam)
    print(io,"GroupParam[")
    len = length(param.components)

    for i in 1:length(param.components)

        print(io, "\"", param.components[i], "\" => [")
        firstloop = true
        for j in 1:length(param.n_groups[i])
            firstloop == false && print(io, ", ")
            print(io, "\"", param.groups[i][j], "\" => ", param.n_groups[i][j])
            firstloop = false
        end
        print(io,']')
        i != length(param.components) && print(io,", ")
    end
    print(io,"]")
end

"""
    SiteParam

Struct holding site parameters.
Is built by parsing all association parameters in the input CSV files.
It has the following fields:
* `components`: a list of all components (or groups in Group Contribution models)
* `sites`: a list containing a list of all sites corresponding to each component (or group) in the components field
* `n_sites`: a list of the site multiplicities corresponding to each site in `flattenedsites`
* `flattenedsites`: a list of all unique sites
* `i_sites`: an iterator that goes through the indices corresponding  to each site in `flattenedsites`
* `n_flattenedsites`: the site multiplicities corresponding to each site in `flattenedsites`
* `i_flattenedsites`: an iterator that goes through the indices for each flattened site

Let's explore the sites in a 3-component `SAFTGammaMie` model:

```julia

julia> model3 = SAFTgammaMie([
                "ethanol",
                ("nonadecanol", ["CH3"=>1, "CH2"=>18, "OH"=>1]),
                ("ibuprofen", ["CH3"=>3, "COOH"=>1, "aCCH"=>1, "aCCH2"=>1, "aCH"=>4])
                               ])

SAFTgammaMie{BasicIdeal} with 3 components:
 "ethanol"
 "nonadecanol"
 "ibuprofen"
Contains parameters: segment, shapefactor, lambda_a, lambda_r, sigma, epsilon, epsilon_assoc, bondvol

julia> model3.sites
SiteParam with 8 sites:
 "CH2OH": "H" => 1, "e1" => 2
 "CH3": (no sites)
 "CH2": (no sites)
 "OH": "H" => 1, "e1" => 2
 "COOH": "e2" => 2, "H" => 1, "e1" => 2
 "aCCH": (no sites)
 "aCCH2": (no sites)
 "aCH": (no sites)

julia> model3.sites.flattenedsites
3-element Vector{String}:
 "H"
 "e1"
 "e2"

julia> model3.sites.i_sites
8-element Vector{Vector{Int64}}:
 [1, 2]
 []
 []
 [1, 2]
 [1, 2, 3]
 []
 []
 []

julia> model3.sites.n_sites
8-element Vector{Vector{Int64}}:
 [1, 2]
 []
 []
 [1, 2]
 [2, 1, 2]
 []
 []
 []
```
"""
struct SiteParam <: ClapeyronParam
    components::Array{String,1}
    sites::Array{Array{String,1},1}
    n_sites::Array{Array{Int,1},1}
    i_sites::Array{Array{Int,1},1}
    flattenedsites::Array{String,1}
    n_flattenedsites::Array{Array{Int,1},1}
    i_flattenedsites::UnitRange{Int}
    sourcecsvs::Array{String,1}
end

function Base.show(io::IO, mime::MIME"text/plain", param::SiteParam)
    print(io,"SiteParam ")
    len = length(param.components)
    println(io,"with ", len, " site", ifelse(len==1, ":", "s:"))

    for i in 1:length(param.components)

        print(io, " \"", param.components[i], "\": ")
        firstloop = true
        if length(param.n_sites[i]) == 0
            print(io,"(no sites)")
        end
        for j in 1:length(param.n_sites[i])
            firstloop == false && print(io, ", ")
            print(io, "\"", param.sites[i][j], "\" => ", param.n_sites[i][j])
            firstloop = false
        end
        i != length(param.components) && println(io)
    end
end

function Base.show(io::IO, param::SiteParam)
    print(io,"SiteParam[")
    len = length(param.components)

    for i in 1:length(param.components)

        print(io, "\"", param.components[i], "\" => [")
        firstloop = true

        for j in 1:length(param.n_sites[i])
            firstloop == false && print(io, ", ")
            print(io, "\"", param.sites[i][j], "\" => ", param.n_sites[i][j])
            firstloop = false
        end
        print(io,']')
        i != length(param.components) && print(io,", ")
    end
    print(io,"]")
end

function SiteParam(pairs::Dict{String,SingleParam{Int}},allcomponentsites)
    arbitraryparam = first(values(pairs))
    components = arbitraryparam.components
    sites = allcomponentsites

    sourcecsvs = String[]
    for x in values(pairs)
        vcat(sourcecsvs,x.sourcecsvs)
    end
    if length(sourcecsvs) >0
        unique!(sourcecsvs)
    end
    n_sites = [[pairs[sites[i][j]].values[i] for j ∈ 1:length(sites[i])] for i ∈ 1:length(components)]  # or groupsites
    length_sites = [length(componentsites) for componentsites ∈ sites]
    i_sites = [1:length_sites[i] for i ∈ 1:length(components)]
    flattenedsites = unique!(reduce(vcat,sites,init = String[]))
    len_flattenedsites = length(flattenedsites)
    i_flattenedsites = 1:len_flattenedsites
    n_flattenedsites = [zeros(Int,len_flattenedsites) for _ ∈ 1:length(components)]
    for i in length(components)
        setindex!.(n_flattenedsites,n_sites,i_sites)
    end
    return SiteParam(components,
    sites,
    n_sites,
    i_sites,
    flattenedsites,
    n_flattenedsites,
    i_flattenedsites,
    sourcecsvs)
end

function SiteParam(input::PARSED_GROUP_VECTOR_TYPE,sourcecsvs::Vector{String}=String[])
    components = [first(i) for i ∈ input]
    raw_sites =  [last(i) for i ∈ input]
    sites = [first.(sitepairs) for sitepairs ∈ raw_sites]
    n_sites = [last.(sitepairs) for sitepairs ∈ raw_sites]
    flattenedsites = unique!(reduce(vcat,sites,init = String[]))
    i_sites = [[findfirst(isequal(site), flattenedsites) for site ∈ componentsites] for componentsites ∈ sites]
    len_flattenedsites = length(flattenedsites)
    i_flattenedsites = 1:len_flattenedsites
    n_flattenedsites = [zeros(Int,len_flattenedsites) for _ ∈ 1:length(input)]
    for i in length(input)
        setindex!.(n_flattenedsites,n_sites,i_sites)
    end

    return SiteParam(components,
    sites,
    n_sites,
    i_sites,
    flattenedsites,
    n_flattenedsites,
    i_flattenedsites,
    sourcecsvs)
end

function SiteParam(components::Vector{String})
    n = length(components)
    return SiteParam(
    components,
    [String[] for _ ∈ 1:n],
    [Int[] for _ ∈ 1:n],
    [Int[] for _ ∈ 1:n],
    String[],
    [Int[] for _ ∈ 1:n],
    1:0,
    String[])
end

paramvals(param::ClapeyronParam) = param.values
paramvals(x) = x


text:
include(GLVisualize.dir("test", "ExampleRunner.jl"))
using ExampleRunner
using Base.Test

config = ExampleRunner.RunnerConfig(
    number_of_frames = 1060,
    record=false,
    interactive_time = 5.0,
    resolution = (500, 500),
    directory = GLVisualize.dir("examples", "sprites", "image_texture_atlas.jl")
)
#config.directory = GLVisualize.dir("examples", "interactive", "image_processing.jl")
ExampleRunner.run(config)
x = first(config.attributes)[2][:timings]
using Plots;glvisualize()
plot(x)


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

using NeuralPDE, Flux, ModelingToolkit, GalacticOptim, Optim, DiffEqFlux

@parameters x y θ
@variables u1(..)
@derivatives Dxx''~x
@derivatives Dyy''~y

# 2D PDE
eq  = Dxx(u1(x,y,θ)) + Dyy(u1(x,y,θ)) ~ -sin(pi*x)*sin(pi*y)

# Boundary conditions
bcs = [u1(0,y,θ) ~ 0.f0, u1(1,y,θ) ~ -sin(pi*1)*sin(pi*y),
       u1(x,0,θ) ~ 0.f0, u1(x,1,θ) ~ -sin(pi*x)*sin(pi*1)]
# Space and time domains
domains = [x ∈ IntervalDomain(0.0,1.0),
           y ∈ IntervalDomain(0.0,1.0)]
# Discretization
dx = 0.1

# Neural network
dim = 2 # number of dimensions
chain = FastChain(FastDense(dim,16,Flux.σ),FastDense(16,16,Flux.σ),FastDense(16,1))

discretization = PhysicsInformedNN(dx,chain)

pde_system = PDESystem(eq,bcs,domains,[x,y],[u1])
prob = discretize(pde_system,discretization)

cb = function (p,l)
    println("Current loss is: $l")
    return false
end

res1 = GalacticOptim.solve(prob, Optim.BFGS(); cb = cb, maxiters=1000)
phi = discretization.phi

##

xs,ys = [domain.domain.lower:dx/10:domain.domain.upper for domain in domains]
analytic_sol_func(x,y) = (sin(pi*x)*sin(pi*y))/(2pi^2)

u_predict = reshape([first(phi([x,y],res1.minimizer)) for x in xs for y in ys],(length(xs),length(ys)))
u_real = reshape([analytic_sol_func(x,y) for x in xs for y in ys], (length(xs),length(ys)))
diff_u = abs.(u_predict .- u_real)

using Plots
p1 = plot(xs, ys, u_real, linetype=:contourf,title = "analytic");
p2 = plot(xs, ys, u_predict, linetype=:contourf,title = "predict");
p3 = plot(xs, ys, diff_u,linetype=:contourf,title = "error");
plot(p1,p2,p3)

max_stars_repo_name: stevengogogo/ProgrammableAging

max_issues_repo_name: stevengogogo/ProgrammableAging

max_forks_repo_name: stevengogogo/ProgrammableAging


text:
module TestDoctrines

using Test
using Catlab.Doctrines
using Catlab.Syntax

sexpr(expr::GATExpr) = sprint(show_sexpr, expr)
unicode(expr::GATExpr) = sprint(show_unicode, expr)
latex(expr::GATExpr) = sprint(show_latex, expr)

@testset "Categories" begin
  include("Category.jl")
end

@testset "MonoidalCategories" begin
  include("Monoidal.jl")
end

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

end


text:
import Distributions
using LoopVectorization
using Random

export MvNormalFast
"""
    MvNormalFast
This is a specialized MvNormal that assumes a **diagonal** covariance. This is substantially
faster than the usual Distributions.jl version that for some reason has large overhead.

# Notes
Only logpdf is defined so no calling rand!
"""
struct MvNormalFast{T<:AbstractVector, S<:AbstractVector, N} <: Distributions.AbstractMvNormal
    μ::T
    Σ::S
    lnorm::N
    function MvNormalFast(μ::T, Σ::S) where {T,S}
        lnorm = -0.5*sum(log.(Σ)) - length(μ)/2*log(2π)
        return new{T,S,typeof(lnorm)}(μ, Σ, lnorm)
    end
end

Base.length(d::MvNormalFast) = length(d.μ)
function Distributions.sampler(d::MvNormalFast)
    return sampler(Distributions.MvNormal(d.μ, sqrt.(d.Σ)))
end

Base.eltype(::MvNormalFast{T,N}) where {T,N} = eltype(T)

@inline function Distributions._logpdf(d::MvNormalFast, x::AbstractVector{T}) where {T}
    μ,Σ = d.μ, d.Σ
    acc = zero(eltype(x))
    @inbounds @fastmath for i in eachindex(x)
        acc += -(x[i]-μ[i])^2/Σ[i]/2
    end
    return acc + d.lnorm
end

function Distributions._rand!(rng::Random.AbstractRNG, d::MvNormalFast, x::AbstractArray)
    return Distributions._rand!(rng, Distributions.MvNormal(d.μ, sqrt.(d.Σ)), x)
end

# These are some soss stuff to make Dirichlet and other stuff work.
#Soss.xform(d::Distributions.Dirichlet, _data::NamedTuple=NamedTuple()) = Soss.TV.UnitSimplex(length(d.alpha))
#Soss.xform(d::Distributions.LogNormal, _data::NamedTuple=NamedTuple()) = Soss.TV.as_positive_real
#Soss.TV.as(d::Distributions.LogNormal) = Soss.TV.as_positive_real


text:
### A Pluto.jl notebook ###
# v0.14.9

using Markdown
using InteractiveUtils

# ╔═╡ e23acf2d-d4b1-4e5c-bbff-1d0e954aee09
using CSV

# ╔═╡ 7d7537ab-d969-4a33-918e-66fcd3d90984
using DataFrames

# ╔═╡ 20ccd69d-e6ea-4d8b-bc2f-ffd4ea8ff4d4
md"""
## Covid with Julia
Soourced from
"""

# ╔═╡ 1464bdcd-3f82-459e-ba2f-6d03f97e2a40
md"""
---------------------------------------------------------------------------------
"""

# ╔═╡ 4407b13c-6893-4910-ab0a-e6042bb2d7c6
md"""
https://github.com/mitmath/6S083

https://www.youtube.com/watch?v=Yr8Og1Rre1I&t=7s

https://github.com/CSSEGISandData/COVID-19

COVID-19/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_confirmed_global.csv
"""

# ╔═╡ 81917f5d-ff6f-4fe0-8448-71239a44f7e7
md"""
https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_confirmed_global.csv
"""

# ╔═╡ 1e3a99cc-d200-4988-abc9-7fac03bf7a69
md"""
---------------------------------------------------------------------------------
"""

# ╔═╡ a1709fd4-3729-4193-af68-15ea361f8d96
md"""
store above URL CSV in a varaible
"""

# ╔═╡ 7f8f8bac-6bfc-4eca-9dc4-dd3296e1f825
url = "https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_confirmed_global.csv"

# ╔═╡ 7c696ef1-34f3-4a83-a4fa-cef998dde43c
download(url, "confirmed_global.csv")

# ╔═╡ ac55e60d-ccfa-4712-8e72-d5269f90f979
data = CSV.read("confirmed_global.csv", DataFrame);

# ╔═╡ f04bc82d-65f6-41e9-b75e-b86addf98e80
md"""
Change Header titles
"""

# ╔═╡ 948e2842-2d99-4d39-b9ad-eadc4d96517c
data2 = rename(data, 1 => "Provice", 2 => "Country" );

# ╔═╡ 6ce29b14-bdf9-41e3-8178-ab38c68745b6
md"""
extract unique countries name
"""

# ╔═╡ e63f364e-4690-4222-a02f-4b6391746774
cntr_list = unique(collect(data2[:,2]));

# ╔═╡ 44900fd7-fa9e-4acb-bae5-92feb117db7c
for i in 1:length(cntr_list)
	#@show cntr_list[4] #output goes on REPL Screen
end

# ╔═╡ dfb08a44-4316-4175-9180-71f462061874
cntr_list[4]

# ╔═╡ 25d9d261-afb2-45e5-8174-4074967f840f
@show cntr_list[7]

# ╔═╡ 0bdf37a8-f935-4293-9243-00a57a8eb752

# ╔═╡ b02a566d-0618-498c-b4ef-49afc52cc6b5

# ╔═╡ 600f7e94-3cfe-4872-b04a-aa376b7edb94

# ╔═╡ 1d084057-2b63-426b-87f6-d94dd9a79e17

# ╔═╡ Cell order:
# ╟─20ccd69d-e6ea-4d8b-bc2f-ffd4ea8ff4d4
# ╟─1464bdcd-3f82-459e-ba2f-6d03f97e2a40
# ╟─4407b13c-6893-4910-ab0a-e6042bb2d7c6
# ╟─81917f5d-ff6f-4fe0-8448-71239a44f7e7
# ╟─1e3a99cc-d200-4988-abc9-7fac03bf7a69
# ╟─a1709fd4-3729-4193-af68-15ea361f8d96
# ╟─7f8f8bac-6bfc-4eca-9dc4-dd3296e1f825
# ╟─7c696ef1-34f3-4a83-a4fa-cef998dde43c
# ╠═e23acf2d-d4b1-4e5c-bbff-1d0e954aee09
# ╠═7d7537ab-d969-4a33-918e-66fcd3d90984
# ╠═ac55e60d-ccfa-4712-8e72-d5269f90f979
# ╟─f04bc82d-65f6-41e9-b75e-b86addf98e80
# ╠═948e2842-2d99-4d39-b9ad-eadc4d96517c
# ╟─6ce29b14-bdf9-41e3-8178-ab38c68745b6
# ╠═e63f364e-4690-4222-a02f-4b6391746774
# ╠═44900fd7-fa9e-4acb-bae5-92feb117db7c
# ╠═dfb08a44-4316-4175-9180-71f462061874
# ╠═25d9d261-afb2-45e5-8174-4074967f840f
# ╠═0bdf37a8-f935-4293-9243-00a57a8eb752
# ╠═b02a566d-0618-498c-b4ef-49afc52cc6b5
# ╠═600f7e94-3cfe-4872-b04a-aa376b7edb94
# ╠═1d084057-2b63-426b-87f6-d94dd9a79e17


text:
export implot, plot_encoding, plot_centroids
export plot_buffer, plot_active_buffer, plot_pipeline

using Plots
using Plots.PlotMeasures
# using ImageView

"""
    function implot(img::AbstractArray; kwargs...)

Plot input image using heatmap.
Magnitude parameter may be precised using the `clim` keyword. The default value
is set to `clim=(0,255)`. Using `clim="auto"` amounts to take the maximum of
the input image as maximum magnitude.

Examples:

    implot(m)
    implot(m, clim="auto")
    implot(m, clim=(1, 10))
"""
function implot(img::AbstractArray; kwargs...)
    kwargs_dict = Dict(kwargs)
    if haskey(kwargs_dict, :clim)
        if kwargs_dict[:clim] == "auto"
            clim = (0, maximum(img))
        else
            clim = kwargs_dict[:clim]
        end
    else
        clim = (0, 255)
    end
    if ndims(img) == 3
        img = img[1,:,:]
    end
    heatmap(transpose(img), yflip=true, color=:grays, clim=clim,
            ratio=:equal)
end

"""
    function display_buffer(
        ind::CGPInd,
        enlargement::E=1;
        indexes=eachindex(ind.buffer)
    ) where {E <: Union{Int64, Float64}}

Display the images contained in each node in the input IPCGP individual.

Examples:

    IICGP.display_buffer(ind)
    IICGP.display_buffer(ind, 2)
    IICGP.display_buffer(ind, indexes=1:3)
    IICGP.display_buffer(ind, 2, indexes=1:3)
"""
function plot_buffer(
    ind::CGPInd, # enlargement::E=1;
    indexes=eachindex(ind.buffer)
) where {E <: Union{Int64, Float64}}
    for i in indexes
        implot(ind.buffer[i]) # , enlargement)
    end
end

"""
    plot_pipeline

Generic pipeline plotting function.
"""
function plot_pipeline(s::Array{UInt8,3}) #, enco::CGPInd)
    # display(implot(s[1]))
    imshow(s, canvassize=(500, 500))
end

"""
    plot_active_buffer(
        ind::CGPInd,
        enlargement::E=1
    ) where {E <: Union{Int64, Float64}}

Display the images contained in each active node in the given IPCGP individual.
"""
function plot_active_buffer(ind::CGPInd)
    for i in eachindex(ind.buffer)
        if ind.nodes[i].active
            #plt = implot(ind.buffer[i]) # , enlargement)
            #=plt = plot(
                ind.buffer[i], seriestype=:heatmap, flip=true,
                ratio=:equal, color=:grays, leg=false, framestyle=:none
            )=#
            plt = heatmap(
                ind.buffer[i], yflip=true, color=:grays, clim=(0,255),
                ratio=:equal, leg=false, framestyle=:none,
                padding = (0.0, 0.0) # , margin=-5mm
            )
            # println(ind.buffer[i][1:3])
            display(plt)
            savefig(plt, "buffer.png")
        end
    end
end

"""
    plot_encoding(
        n_in::Int64,
        buffer::Array{Array{UInt8,2},1},
        features::AbstractArray
    )

Plot the complete encoding pipeline from input to projection on feature space.
"""
function plot_encoding(
        n_in::Int64,
        buffer::Array{Array{UInt8,2},1},
        features::AbstractArray
    )
    n_cols = max(n_in, length(features), length(buffer)-n_in)
    mrg = -4mm
    pdg = (0.0, 0.0)
    p = plot(layout=grid(3, n_cols), leg=false, framestyle=:none, margin=mrg) #axis=nothing)
    for i in 1:n_in
        plot!(p[i], buffer[i], seriestype=:heatmap, flip=true, ratio=:equal,
              color=:grays)
    end
    for i in n_in+1:length(buffer)
        plot!(p[2,i-n_in], buffer[i], seriestype=:heatmap, flip=true, ratio=:equal, color=:grays)
    end
    for i in eachindex(features)
        plot!(p[3,i], features[i], seriestype=:heatmap, flip=true, ratio=:equal, color=:grays)
    end
    p
end

"""
    plot_centroids(
        images::Array{Array{UInt8,2},1},
        centroids::Array{Array{Tuple{Float64,Float64},1},1}
    )

Given an image and a set of centroids, plot the image as a heatmap along with
the centroids.
"""
function plot_centroids(
        images::Array{Array{UInt8,2},1},
        centroids::Array{Array{Tuple{Float64,Float64},1},1}
    )
    x = images[1]
    img_size = size(x)
    mlt = 1
    mrg = -100mm
    pdg = (0.0, 0.0)
    final_size = (mlt*img_size[2], mlt*img_size[1])
    centro = IICGP.scaled_centroids(centroids[1], img_size)

    xs = [c[1] for c in centro]
    ys = [c[2] for c in centro]
    pal = palette([:blue, :red, :orange, :yellow, :green], length(xs))
    plt = heatmap(x, color=:grays, ratio=:equal, yflip=true, leg=false,
                  framestyle=:none, padding=pdg, margin=mrg,
                  size=final_size)
    for i in eachindex(xs)
        scatter!(plt, [ys[i]], [xs[i]], padding=pdg,
                 margin=mrg, color=pal[i])
    end
    plt
end


text:
@testset "wiener" begin
    k_1 = WienerKernel(; i=-1)
    @test typeof(k_1) <: WhiteKernel

    k0 = WienerKernel()
    @test k0 isa WienerKernel{0}

    k1 = WienerKernel(; i=1)
    @test k1 isa WienerKernel{1}

    k2 = WienerKernel(; i=2)
    @test k2 isa WienerKernel{2}

    k3 = WienerKernel(; i=3)
    @test k3 isa WienerKernel{3}

    @test_throws ArgumentError WienerKernel(; i=4)
    @test_throws ArgumentError WienerKernel(; i=-2)

    v1 = rand(4)
    v2 = rand(4)

    X = sqrt(sum(abs2, v1))
    Y = sqrt(sum(abs2, v2))
    minXY = min(X, Y)

    @test k0(v1, v2) ≈ minXY
    @test k1(v1, v2) ≈ 1 / 3 * minXY^3 + 1 / 2 * minXY^2 * euclidean(v1, v2)
    @test k2(v1, v2) ≈
          1 / 20 * minXY^5 + 1 / 12 * minXY^3 * euclidean(v1, v2) * (X + Y - 1 / 2 * minXY)
    @test k3(v1, v2) ≈
          1 / 252 * minXY^7 +
          1 / 720 *
          minXY^4 *
          euclidean(v1, v2) *
          (5 * max(X, Y)^2 + 2 * X * Y + 3 * minXY^2)

    # Standardised tests. Requires careful input choice.
    x0 = rand(3)
    x1 = rand(3)
    x2 = rand(2)
    TestUtils.test_interface(k0, x0, x1, x2)
    TestUtils.test_interface(k1, x0, x1, x2)
    TestUtils.test_interface(k2, x0, x1, x2)
    TestUtils.test_interface(k3, x0, x1, x2)
    # test_ADs(()->WienerKernel(i=1))
    @test_broken "No tests passing"
end


text:
# Calculate boundary "forces" in advection

import PDESolver._evalFunctional

@doc """
### AdvectionEquationMod.evalFunctional

Hight level function that evaluates functionals specified in the options
dictionary. The user must call this function for functional evaluation.This
function is agnostic which type of a functional is being computed and calls a
mid level type specific function for the actual functional evaluation.

**Arguments**

*  `mesh` :  Abstract mesh object
*  `sbp`  : Summation-By-Parts operator
*  `eqn`  : Advection equation object
*  `opts` : Options dictionary
*  `functionalData` : Object of the functional being computed.
"""->
function _evalFunctional(mesh::AbstractMesh{Tmsh},
            sbp::AbstractOperator, eqn::AdvectionData{Tsol}, opts,
            functionalData::AbstractBoundaryFunctional) where {Tmsh, Tsol}
  if mesh.isDG
    boundaryinterpolate!(mesh.sbpface, mesh.bndryfaces, eqn.q, eqn.q_bndry)
  end

  calcBndryFunctional(mesh, sbp, eqn, opts, functionalData)

  return functionalData.val
end

@doc """
AdvectionEquationMod.calcBndryfunctional

This function calculates the functional on a geometric boundary of a the
computational space. This is a mid level function that should not be called
from outside the module. Depending on the functional being computed, it
may be necessary to define another method for this function based on a
different boundary functional type.

**Arguments**

*  `mesh` :  Abstract mesh object
*  `sbp`  : Summation-By-Parts operator
*  `eqn`  : Advection equation object
*  `opts` : Options dictionary
*  `functionalData` : Object of the functional being computed

"""->
#=  TODO: uncomment and run test when mesh.bndry_geo_nums gets added to CG meshes

function calcBndryFunctional(mesh::AbstractCGMesh{Tmsh},sbp::AbstractOperator,
                         eqn::AdvectionData{Tsol}, opts, functor, functional_edges) where {Tmsh, Tsol}

  # Specify the boundary conditions for the edge on which the force needs to be
  # computed separately. Use that boundary number to access the boundary
  # offset array. Then proceed the same as bndryflux to get the forces using
  # boundaryintegrate!

  functional_val = zero(Tsol)
  alpha_x = eqn.params.alpha_x
  alpha_y = eqn.params.alpha_y

  for itr = 1:length(functional_edges)
    g_edge_number = functional_edges[itr] # Extract geometric edge number

    # get the boundary array associated with the geometric edge
    itr2 = 0
    for itr2 = 1:mesh.numBC
      if findfirst(mesh.bndry_geo_nums[itr2],g_edge_number) > 0
        break
      end
    end
    start_index = mesh.bndry_offsets[itr2]
    end_index = mesh.bndry_offsets[itr2 + 1]
    idx_range = start_index:(end_index-1)  # Index range
    bndry_facenums = sview(mesh.bndryfaces, idx_range) # faces on geometric edge i

    nfaces = length(bndry_facenums)
    boundary_integrand = zeros(Tsol, 1, mesh.numNodesPerFace, nfaces)
    boundary_functional = zeros(Tsol, 1, sbp.numnodes, mesh.numEl)

    for i = 1:nfaces
      bndry_i = bndry_facenums[i]
      for j = 1:mesh.numNodesPerFace
        k = mesh.facenodes[j, bndry_i.face]
        q = eqn.q[1,k,bndry_i.element]
        x = ro_sview(mesh.coords, :, k, bndry_i.element)
        dxidx = ro_sview(mesh.dxidx, :, :, k, bndry_i.element)
        nrm = ro_sview(mesh.sbpface.normal, :, bndry_i.face)
        nx = dxidx[1,1]*nrm[1] + dxidx[2,1]*nrm[2]
        ny = dxidx[1,2]*nrm[1] + dxidx[2,2]*nrm[2]
        boundary_integrand[1,j,i] = functor(eqn.params, nx, ny, q) # Boundary Flux
      end
    end

    boundaryintegrate!(mesh.sbpface, mesh.bndryfaces[idx_range], boundary_integrand, boundary_functional)
    # Add all boundary_force nodal values along the edge to get the nodal force value
    edge_functional_val = zero(Tsol) # functional value over a geometric edge
    for (bindex, bndry) in enumerate(mesh.bndryfaces[idx_range])
      for i = 1:mesh.numNodesPerFace
        k = mesh.facenodes[i, bndry.face]
        edge_functional_val += boundary_functional[1,k,bndry.element]
      end  # end for i = 1:mesh.numNodesPerFace
    end    # end enumerate

    functional_val += edge_functional_val
  end      # for itr = 1:length(functional_edges)

  return functional_val
end
=#

function calcBndryFunctional(mesh::AbstractDGMesh{Tmsh},
          sbp::AbstractOperator,
          eqn::AdvectionData{Tsol}, opts,
          functionalData::AbstractBoundaryFunctional{Topt}) where {Tmsh, Tsol, Topt}

  # Specify the boundary conditions for the edge on which the force needs to be
  # computed separately. Use that boundary number to access the boundary
  # offset array. Then proceed the same as bndryflux to get the integrand. Finally
  # use integratefunctional! to get the solution.

  functionalData.val = zero(Topt)
  local_functional_val = zero(Tsol)
  alpha_x = eqn.params.alpha_x
  alpha_y = eqn.params.alpha_y

  for itr = 1:length(functionalData.bcnums)
    bcnum = functionalData.bcnums[itr]

    start_index = mesh.bndry_offsets[bcnum]
    end_index = mesh.bndry_offsets[bcnum+1]
    idx_range = start_index:(end_index-1)
    bndry_facenums = sview(mesh.bndryfaces, idx_range) # faces on geometric edge i

    nfaces = length(bndry_facenums)
    boundary_integrand = zeros(Tsol, 1, mesh.sbpface.numnodes, nfaces)

    for i = 1:nfaces
      global_facenum = idx_range[i]

      for j = 1:mesh.sbpface.numnodes
        q = eqn.q_bndry[ 1, j, global_facenum]
        coords = ro_sview(mesh.coords_bndry, :, j, global_facenum)
        nrm = ro_sview(mesh.nrm_bndry, :, j, global_facenum)
        nx = nrm[1]
        ny = nrm[2]
        boundary_integrand[1,j,i] = calcBoundaryFunctionalIntegrand(eqn.params, nx, ny, q,
                                    functionalData) # Boundary Flux
      end
    end

    val_per_geom_edge = zeros(Tsol, 1)
    integratefunctional!(mesh.sbpface, mesh.bndryfaces[idx_range],
                         boundary_integrand, val_per_geom_edge)

    # Add contributions of multiple geometric edges to the final functional value
    local_functional_val += val_per_geom_edge[1]

  end  # end loop itr

  functionalData.val = MPI.Allreduce(local_functional_val, MPI.SUM, eqn.comm)

  return nothing
end

#-----    Computing Functional Integrand    -----#
@doc """
### AdvectionEquationMod.calcBoundaryFunctionalIntegrand

Computes the integrand for boundary functional at a surface SBP node. Every
functional needs to have its own method and the functional type determines
which method is called.

**Inputs**

*  `params` : eqn.params object
*  `nx` : X component of face normal vector
*  `ny` : Y component of face normal vector
*  `q`  : Nodal solution variable
*  `functionalData` : Object of the functional being computed

**Outputs**

* `functional_integrand` : Computed integrand at the surface node

"""->

function calcBoundaryFunctionalIntegrand(params::ParamType2, nx, ny, q,
                                         functionalData::QfluxData)

  alpha_x = params.alpha_x
  alpha_y = params.alpha_y

  return functional_integrand = (alpha_x*nx + alpha_y*ny)*q
end

"""
  Method for [`IntegralQData`](@ref) functional.

  Note that q should be scaled by the length of the normal vector so the
  integration works correctly
"""
function calcBoundaryFunctionalIntegrand(params::ParamType2, nx, ny, q,
                                         functionalData::IntegralQData)

  fac = sqrt(nx*nx + ny*ny)
  return q*fac
end


text:
function SP3Parser()
    return SP3Parser(())
end

function SP3Parser(arg0::jdouble, arg1::jint, arg2::Function)
    return SP3Parser((jdouble, jint, Function), arg0, arg1, arg2)
end

function SP3Parser(arg0::jdouble, arg1::jint, arg2::Function, arg3::TimeScales)
    return SP3Parser((jdouble, jint, Function, TimeScales), arg0, arg1, arg2, arg3)
end

function parse(obj::SP3Parser, arg0::BufferedReader, arg1::JString)
    return jcall(obj, "parse", SP3File, (BufferedReader, JString), arg0, arg1)
end

function parse(obj::SP3Parser, arg0::InputStream)
    return jcall(obj, "parse", SP3File, (InputStream,), arg0)
end

function parse(obj::SP3Parser, arg0::JString)
    return jcall(obj, "parse", SP3File, (JString,), arg0)
end


text:
const FILE_FORMAT = ("2.2", "0", "8")

struct GeomType
    gmsh_code :: Int
    dimen     :: Int
    nonodes   :: Int
end

const LINE        = GeomType(1, 1, 2)
const TRIANGLE    = GeomType(2, 2, 3)
const TETRAHEDRON = GeomType(4, 3, 4)

const GETGEOMTYPE = Dict(1 => LINE,
                         2 => TRIANGLE,
                         4 => TETRAHEDRON)

# Mesh data type with fields
#
# coord[:,n]          = x, y, z coordinates of node n
# physname[m]         = string label of physical group with
#                       the numerical label 'm'
# physnum[s]          = numerical label of physical group
#                       the string label 's'
# elmtype[physname]   = element type for elements in physname
# elms_of[physname]   = array whose columns give the node numbers
#                       of the elements in physname
# elmnbrs_of[physnam] = Gmsh element numbers
#                       of the elements in physname
# nodes_of[physname]  = nodes in physname
#

struct Mesh
    coord      :: Array{Float64, 2}
    physdim    :: Dict{String, Int}
    physnum    :: Dict{String, Int}
    physname   :: Dict{Int, String}
    elmtype    :: Dict{String, GeomType}
    elms_of    :: Dict{String, Matrix{Int}}
    elmnbrs_of :: Dict{String, Vector{Int}}
    nodes_of   :: Dict{String, Vector{Int}}
end

function verify_next_line(s, f)
    line = readline(f)
    @assert strip(line) == s
end

"""
mesh = read_msh_file("filename.msh")

Read a .msh file created by Gmsh and return the corresponding
Mesh object.
"""
function read_msh_file(fname::String)
    coord = Float64[]
    nonodes = 0
    line = ""
    elms = String[]
    elmcount = Dict{Int,Int}()
    physdim  = Dict{String,Int}()
    physnum  = Dict{String,Int}()
    physname = Dict{Int,String}()
    f = open(fname,"r")
    while !eof(f)
        header = strip( readline(f) )
        if header == "\$MeshFormat"
            fmt = split( readline(f) )
            for j = 1:3
                @assert fmt[j] == FILE_FORMAT[j]
            end
            verify_next_line("\$EndMeshFormat", f)
        elseif header == "\$Nodes"
            line = readline(f)
            nonodes = parse(Int, line)
            for n = 1:nonodes
                line = readline(f)
                s = split(line)
                @assert n == parse(Int, s[1])
                for j = 1:3
                   push!(coord, parse(Float64, s[j+1]))
                end
            end
            verify_next_line("\$EndNodes", f)
        elseif header == "\$Elements"
            line = readline(f)
            noelms = parse(Int, line)
            for n = 1:noelms
                line = readline(f)
                push!(elms, line)
                row = split(line)
                physno = parse(Int, row[4])
                elmcount[physno] = get(elmcount, physno, 0) + 1
            end
            verify_next_line("\$EndElements", f)
        elseif header == "\$PhysicalNames"
            line = readline(f)
            nonames = parse(Int, line)
            for n = 1:nonames
                line = readline(f)
                s = split(line)
                dimen = parse(Int, s[1])
                num   = parse(Int, s[2])
                name  = strip(s[3], '"')
                physdim[name] = dimen
                physnum[name] = num
                physname[num] = name
            end
            verify_next_line("\$EndPhysicalNames", f)
        elseif header == "\$Comment"
            while true
                line = readline(f)
                if strip(line) == "\$EndComment"
                    break
                end
            end
        else
            msg = string("Unknown section header: ", line)
            error(msg)
        end
    end
    close(f)
    coord = reshape(coord, (3,nonodes))

    elmtype    = Dict{String, GeomType}()
    elms_of    = Dict{String,Matrix{Int}}()
    elmnbrs_of = Dict{String,Vector{Int}}()
    nodes_of   = Dict{String,Vector{Int}}()
    k = Dict{String,Int}()
    for elm in elms
        s = split(elm)
        name = physname[parse(Int, s[4])]
        gmsh_code = parse(Int, s[2])
        if haskey(elmtype, name)
            @assert elmtype[name] == GETGEOMTYPE[gmsh_code]
        else
            elmtype[name] = GETGEOMTYPE[gmsh_code]
        end
    end
    for (name, element_type) in elmtype
        noelms = elmcount[physnum[name]]
        elms_of[name]    = zeros(Int64, element_type.nonodes, noelms)
        elmnbrs_of[name] = zeros(Int64, noelms)
        k[name] = 0
    end
    for elm in elms
        line = [parse(Int, s) for s in split(elm)]
        # line uses the format
        # elm-number elm-type number-of-tags <tags> node-number-list
        # where the first tag is the numerical label of the physical
        # name.
        notags = line[3]
        name = physname[line[4]]
        k[name] += 1
        elms_of[name][:,k[name]]  = line[notags+4:end]
        elmnbrs_of[name][k[name]] = line[1]
    end
    for (name, elm) in elms_of
        nodes_of[name] = noduplicates(elm[:])
    end
    return Mesh(coord, physdim, physnum, physname, elmtype,
                elms_of, elmnbrs_of, nodes_of)
end

"""
b = noduplicates(a)

Return a vector b containing the elements of a with no duplicates,
and sorted in increasing order.
"""
function noduplicates(a::Vector{Int})
    sort!(a)
    b = Array{Int}(undef, length(a))
    j = 1
    b[j] = a[1]
    for k = 2:length(a)
        if a[k] != a[k-1]
            j += 1
            b[j] = a[k]
        end
    end
    return b[1:j]
end

function save_nodal_scalar_field(u::Vector{Float64},
                                 name::String,
                                 fid::IOStream,
                                 time=0.0, timeidx=0)
    write(fid, "\$NodeData\n")
    write(fid, "1\n")             # number of string tags
    write(fid, "\"$name\"\n")     # label for scalar field u
    write(fid, "1\n")             # number of real tags
    s = @sprintf("%.15e\n", time)
    write(fid, s)
    write(fid, "3\n")             # number of integer tags
    s = @sprintf("%d\n", timeidx)
    write(fid, s)
    write(fid, "1\n")             # number of field components
    nonodes = length(u)
    s = @sprintf("%d\n", nonodes)
    write(fid, s)
    for nd = 1:nonodes
        s = @sprintf("%d  %.15e\n", nd, u[nd])
        write(fid, s)
    end
    write(fid, "\$EndNodeData\n")
end

"""
Gmsh kludge to facilitate surface plotting of 2D scalar field via
the warp plugin.
"""
function save_warp_nodal_scalar_field(u::Array{Float64},
                                      name::String,
                                      fid::IOStream,
                                      time=0.0, timeidx=0)
    n = length(u)
    uwarp = zeros(3, length(u))
    uwarp[3,:] = u
    save_nodal_scalar_field(u, name, fid, time, timeidx)
    save_nodal_vector_field(uwarp, "$(name)_warp", fid, time, timeidx)
end

function write_format_version(fid::IOStream)
    write(fid, "\$MeshFormat\n")
    write(fid, "$(FILE_FORMAT[1]) $(FILE_FORMAT[2]) $(FILE_FORMAT[3])\n")
    write(fid, "\$EndMeshFormat\n")
end

function save_nodal_vector_field(u::Array{Float64,2}, name::String,
                                 fid::IOStream,
                                 time=0.0, timeidx=0)
    write(fid, "\$NodeData\n")
    write(fid, "1\n")             # number of string tags
    write(fid, "\"$name\"\n")     # label for scalar field u
    write(fid, "1\n")             # number of real tags
    s = @sprintf("%.15e\n", time)
    write(fid, s)
    write(fid, "3\n")             # number of integer tags
    s = @sprintf("%d\n", timeidx)
    write(fid, s)
    write(fid, "3\n")             # number of field components
    nonodes = size(u, 2)
    s = @sprintf("%d\n", nonodes)
    write(fid, s)
    @assert size(u,1) == 3
    for nd = 1:nonodes
        s = @sprintf("%d  %.15e  %.15e  %.15e\n",
                     nd, u[1,nd], u[2,nd], u[3,nd])
        write(fid, s)
    end
    write(fid, "\$EndNodeData\n")
end

function save_elm_scalar_field(u::Dict{String,Vector{Float64}},
                               name::String,
                               mesh::Mesh, elmtype::GeomType,
                               fid::IOStream,
                               time=0.0, timeidx=0)
    write(fid, "\$ElementData\n")
    write(fid, "1\n")             # number of string tags
    write(fid, "\"$name\"\n")     # label for scalar field u
    write(fid, "1\n")             # number of real tags
    s = @sprintf("%.15e\n", time)
    write(fid, s)
    write(fid, "3\n")             # number of integer tags
    s = @sprintf("%d\n", timeidx)
    write(fid, s)
    write(fid, "1\n")             # number of field components
    noelms = count_elms(mesh, elmtype)
    s = @sprintf("%d\n", noelms)
    write(fid, s)
    for (name, tp) in mesh.elmtype
        if tp == elmtype
            elms    = mesh.elms_of[name]
            elmnbrs = mesh.elmnbrs_of[name]
            for k = 1:size(elms,2)
                nbr = elmnbrs[k]
                s = @sprintf("%d  %.15e\n", nbr, u[name][k])
                write(fid, s)
            end
        end
    end
    write(fid, "\$EndElementData\n")
end

"""
Convenience function for use with do-syntax.
"""
function write_pos_file(f::Function, fname::String)
    fid = open(fname, "w")
    try
        write_format_version(fid)
        f(fid)
    finally
        close(fid)
    end
end

"""
Convenience function for use with do-syntax.
Include the mesh data.
"""
function write_pos_file(f::Function, posfname::String,
                        meshfname::String)
    cp(meshfname, posfname)
    fid = open(posfname, "a")
    try
        f(fid)
    finally
        close(fid)
    end
end

function get_node_coords(mesh::Mesh)
   nonodes = size(mesh.coord, 2)
   x = zeros(nonodes)
   y = zeros(nonodes)
   z = zeros(nonodes)
   for nd = 1:nonodes
       x[nd] = mesh.coord[1,nd]
       y[nd] = mesh.coord[2,nd]
       z[nd] = mesh.coord[3,nd]
   end
   return x, y, z
end

"""
Convenience function returns mesh parameter h.
"""
function get_mesh_h(mesh::Mesh)
    h = 0.0
    for name in keys(mesh.physnum)
        next_h = get_mesh_h(mesh, name)
	h = h > next_h ? h : next_h
    end
    return h
end

function get_mesh_h(mesh::Mesh, name)
    x, y, z = get_node_coords(mesh)
    h = 0.0
    elm = mesh.elms_of[name]
    etype = mesh.elmtype[name]
    if etype == LINE
        pt = zeros(3,2)
        for k = 1:size(elm,2)
            for nd = 1:2
                pt[1,nd] = x[elm[nd,k]]
                pt[2,nd] = y[elm[nd,k]]
                pt[3,nd] = z[elm[nd,k]]
            end
            h = max(h, norm(pt[:,2]-pt[:,1]))
        end
    elseif etype == TRIANGLE
        pt = zeros(3,3)
        for k = 1:size(elm,2)
            for nd = 1:3
                pt[1,nd] = x[elm[nd,k]]
                pt[2,nd] = y[elm[nd,k]]
                pt[3,nd] = z[elm[nd,k]]
            end
            s1 = norm(pt[:,2]-pt[:,1])
            s2 = norm(pt[:,3]-pt[:,2])
            s3 = norm(pt[:,1]-pt[:,3])
            diam = max(s1, s2, s3)
            h = max(h, diam)
        end
    else
        error("Not implemented for elements of type $etype")
    end
    return h
end

function count_elms(mesh::Mesh, elmtype::GeomType)
    n = 0
    for name in keys(mesh.physnum)
        if mesh.elmtype[name] == elmtype
	    n += count_elms(mesh, name)
	end
    end
    return n
end

function count_elms(mesh::Mesh, name)
    return size(mesh.elms_of[name], 2)
end

"""
Returns the vector of nodal values of f
"""
function get_nodal_vals(f::Function, mesh::Mesh)
    nonodes = size(mesh.coord, 2)
    x = zeros(3)
    fvec = zeros(nonodes)
    for nd = 1:nonodes
        for j = 1:3
            x[j] = mesh.coord[j,nd]
        end
        fvec[nd] = f(x)
    end
    return fvec
end

"""
Read the file shape.geo and writes a sequence of files
shape0.msh, shape1.msh, ...,
by generating an initial mesh (shape0.msh) with maximum
element size h0, and performing the specified number of
refinements.  Thus, refinements+1 files are generated
altogether.
"""
function successive_refine(shape::String, dimen::Integer,
                           refinements::Integer, h0::Float64)
    infile = "$shape.geo"
    outfile = "$(shape)0.msh"
    run(`gmsh -$dimen -clmax $h0 -o $outfile -format msh22 $infile`)
    for k = 1:refinements
        infile = "$shape$(k-1).msh"
        outfile = "$shape$k.msh"
        run(`gmsh -v 0 -refine -o $outfile -format msh22 $infile`)
    end
    return
end


text:
### Fit a full rank model with sketched Frank Wolfe

# todo: make it work for mpca

import FirstOrderOptimization: frank_wolfe_sketched, FrankWolfeParams, DecreasingStepSize,
                               AbstractSketch, AsymmetricSketch,
                               IndexingOperator, LowRankOperator,
                               IndexedLowRankOperator, mysvds
import LowRankModels: fit!, ConvergenceHistory, get_yidxs, grad, evaluate
import Base: axpy!, scale!

export fit_sketch!, FrankWolfeParams

### FITTING
function fit_sketch!(gfrm::GFRM, params::FrankWolfeParams = FrankWolfeParams();
			  ch::ConvergenceHistory=ConvergenceHistory("FrankWolfeGFRM"),
        z::AbstractVector = zeros(sum(map(length, gfrm.observed_examples))),
        sketch::AbstractSketch = AsymmetricSketch(size(gfrm.A)..., gfrm.k),
        fenchel_tol::Float64 = 1e-5,
			  verbose=true,
			  kwargs...)

    if !isa(gfrm.r, TraceNormConstraint)
        error("Frank Wolfe fitting is only implemented for trace norm constrained problems")
    end

    # the functions below close over all this problem data
    yidxs = get_yidxs(gfrm.losses)
    d = maximum(yidxs[end])
    m,n = size(gfrm.A)
    alpha = gfrm.r.scale

    nobsj = map(length, gfrm.observed_examples)
    startobsj = cumsum(vcat(1, nobsj))
    nobs = sum(nobsj)
    obs = Array(Int, nobs)

    Threads.@threads for j=1:n
        obs[startobsj[j]:(startobsj[j+1]-1)] = m*(j-1) + gfrm.observed_examples[j]
    end
    indexing_operator = IndexingOperator(m, n, obs)
    @assert size(indexing_operator, 1) == nobs
    @assert length(z) == nobs

    objs = zeros(Threads.nthreads())
    function f(z::Array{Float64,1})
        scale!(objs, 0)
        # println("obj")
        # @time
        Threads.@threads for j=1:n
            objs[Threads.threadid()] += evaluate(gfrm.losses[j],
                         z[startobsj[j]:(startobsj[j+1]-1)],
                         gfrm.A[gfrm.observed_examples[j],j])
        end
        return sum(objs)
    end

    ## Grad of f
    g = Array(Float64,nobs) # working variable for computing gradient; grad_f mutates this
    function grad_f(z::Array{Float64,1})
        # println("grad")
        # @time
        Threads.@threads for j=1:n
            ii = startobsj[j]:(startobsj[j+1]-1)
            g[ii] = grad(gfrm.losses[j], z[ii], gfrm.A[gfrm.observed_examples[j],j])
        end
        # return G = A'*g
        # LowRankOperator(indexing_operator, g, transpose = Symbol[:T, :N])
        # @show g
        return IndexedLowRankOperator(indexing_operator, g)
    end

    const_nucnorm(z) = alpha # we'll always saturate the constraint, don't bother computing it
    # returns solution, optval of min <G, Delta> st ||Delta||_* \leq alpha
    function min_lin_st_nucnorm_sketched(G, alpha, tol::Float64 = fenchel_tol)
        # this looks inefficient, but ga is a sparse matrix, so it's not bad
        ga = Array(G)
        # println("svds")
        # @time
        u,s,v = mysvds(ga, nsv=1, tol=tol)
        # u,s,v = mysvds(G, nsv=1) # right now there's no method to directly do svds efficiently on an IndexedLowRankOperator
        return LowRankOperator(-alpha*u, v'), -alpha*s[1]
    end

    # recover
    t = time()
    X_sketched = frank_wolfe_sketched(
        z,
        f, grad_f,
        const_nucnorm,
        alpha,
        min_lin_st_nucnorm_sketched,
        sketch,
        params,
        ch,
        verbose=verbose
        )

    t = time() - t
    update_ch!(ch, t, f(z))

    gfrm.W = X_sketched

    return X_sketched, ch
end


text:
@reexport module TSMLTypes

using DataFrames

export 	Transformer,
		TSLearner,
		fit!,
		transform!

abstract type Transformer end
abstract type TSLearner <: Transformer end

"""
    fit!(tr::Transformer, instances::T, labels::Vector) where {T<:Union{Vector,Matrix,DataFrame}}

Generic `fit!` function to be redefined using multidispatch in  different subtypes of `Transformer`.
"""
function fit!(tr::Transformer, instances::T, labels::Vector) where {T<:Union{Vector,Matrix,DataFrame}}
	error(typeof(tr)," not implemented yet: fit!")
end

"""
    transform!(tr::Transformer, instances::T) where {T<:Union{Vector,Matrix,DataFrame}}

Generic `transform!` function to be redefined using multidispatch in  different subtypes of `Transformer`.
"""
function transform!(tr::Transformer, instances::T) where {T<:Union{Vector,Matrix,DataFrame}}
	error(typeof(tr)," not implemented yet: transform!")
end

end


text:
if basename(pwd()) == "aoc"
    cd("2017/13")
end

using OffsetArrays

function loadfirewall(filename::AbstractString)
    data = collect(parse.(Int, split(line, ": ")) for line in eachline(filename))
    firewall = OffsetArray(zeros(Int, first(last(data)) + 1), -1)
    for (i, v) in data
        firewall[i] = v
    end
    firewall
end

function part1(firewall, start = 0)
    total = 0
    for (depth, range) in zip(Iterators.countfrom(start), firewall)
        if !iszero(range) && iszero(scannerposition(range, depth))
            total += depth * range
        end
    end
    total
end

function scannerposition(range, step)
    if range < 2
        return 0
    end
    roundrange = 2 * range - 2
    roundstep = step % roundrange
    if roundstep < range
        roundstep
    else
        roundrange - roundstep
    end
end

@assert collect(scannerposition(5, x) for x in 0:10) == [0:4..., 3:-1:0..., 1:2...]
@assert part1(loadfirewall("example.txt")) == 24
part1(loadfirewall("input.txt"))

function part2(filename::AbstractString)
    firewall = loadfirewall(filename)
    for start in Iterators.countfrom(1)
        if all(
            iszero(range) || !iszero(scannerposition(range, depth))
            for (depth, range) in zip(Iterators.countfrom(start), firewall)
        )
            return start
        end
    end
end

@assert part2("example.txt") == 10
part2("input.txt")

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
isdefined(:__precompile__) && __precompile__(true)
module InplaceOps

include("common.jl")
include("inplace_math.jl")

export @in1!, @in2!, @into!

inplace_sym(s::Symbol) = inplace_sym(Operator{s}())

inplace_sym{S}(::Operator{S}) = error("invalid operator")

inplace_sym(::Operator{:+}) = :add!
inplace_sym(::Operator{:-}) = :sub!
inplace_sym(::Operator{:*}) = :mul!
inplace_sym(::Operator{:\}) = :ldiv!
inplace_sym(::Operator{:/}) = :rdiv!

inplace_sym(::Operator{:.+}) = :badd!
inplace_sym(::Operator{:.-}) = :bsub!
inplace_sym(::Operator{:.*}) = :bmul!
inplace_sym(::Operator{:.\}) = :bldiv!
inplace_sym(::Operator{:./}) = :brdiv!

op_ctranspose{T}(x::AbstractArray{T}) = CTranspose(x)
op_ctranspose{T<:Real}(x::AbstractArray{T}) = Transpose(x)

op_transpose{T}(x::AbstractArray{T}) = Transpose(x)

typealias AbstractVMF Union{AbstractVecOrMat,Factorization}

#TODO: Most of the 2-argument A_foo_B methods overwrite B, though there are some exceptions (e.g. QRPackedQ)
mul!(::Type{Inplace{2}}, A::AbstractVMF, B::AbstractVMF) = A_mul_B!(A,B)
mul!(::Type{Inplace{2}}, A::Transpose  , B::AbstractVMF) = At_mul_B!(A.obj,B)
mul!(::Type{Inplace{2}}, A::CTranspose , B::AbstractVMF) = Ac_mul_B!(A.obj,B)

ldiv!(::Type{Inplace{2}}, A::AbstractVMF, B::AbstractVMF) = A_ldiv_B!(A,B)
ldiv!(::Type{Inplace{2}}, A::Transpose  , B::AbstractVMF) = At_ldiv_B!(A.obj,B)
ldiv!(::Type{Inplace{2}}, A::CTranspose , B::AbstractVMF) = Ac_ldiv_B!(A.obj,B)

rdiv!(::Type{Inplace{2}}, A::AbstractVMF, B::AbstractVMF) = A_rdiv_B!(A,B)
rdiv!(::Type{Inplace{2}}, A::Transpose  , B::AbstractVMF) = At_rdiv_B!(A.obj,B)
rdiv!(::Type{Inplace{2}}, A::CTranspose , B::AbstractVMF) = Ac_rdiv_B!(A.obj,B)

mul!(O::AbstractVMF, A::AbstractVMF, B::AbstractVMF) = A_mul_B!(O,A,B)
mul!(O::AbstractVMF, A::Transpose  , B::AbstractVMF) = At_mul_B!(O,A.obj,B)
mul!(O::AbstractVMF, A::CTranspose , B::AbstractVMF) = Ac_mul_B!(O,A.obj,B)
mul!(O::AbstractVMF, A::AbstractVMF, B::Transpose  ) = A_mul_Bt!(O,A,B.obj)
mul!(O::AbstractVMF, A::Transpose  , B::Transpose  ) = At_mul_Bt!(O,A.obj,B.obj)
mul!(O::AbstractVMF, A::CTranspose , B::Transpose  ) = Ac_mul_Bt!(O,A.obj,B.obj)
mul!(O::AbstractVMF, A::AbstractVMF, B::CTranspose ) = A_mul_Bc!(O,A,B.obj)
mul!(O::AbstractVMF, A::Transpose  , B::CTranspose ) = At_mul_Bc!(O,A.obj,B.obj)
mul!(O::AbstractVMF, A::CTranspose , B::CTranspose ) = Ac_mul_Bc!(O,A.obj,B.obj)

ldiv!(O::AbstractVMF, A::AbstractVMF, B::AbstractVMF) = A_ldiv_B!(O,A,B)
ldiv!(O::AbstractVMF, A::Transpose  , B::AbstractVMF) = At_ldiv_B!(O,A.obj,B)
ldiv!(O::AbstractVMF, A::CTranspose , B::AbstractVMF) = Ac_ldiv_B!(O,A.obj,B)
ldiv!(O::AbstractVMF, A::AbstractVMF, B::Transpose  ) = A_ldiv_Bt!(O,A,B.obj)
ldiv!(O::AbstractVMF, A::Transpose  , B::Transpose  ) = At_ldiv_Bt!(O,A.obj,B.obj)
ldiv!(O::AbstractVMF, A::CTranspose , B::Transpose  ) = Ac_ldiv_Bt!(O,A.obj,B.obj)
ldiv!(O::AbstractVMF, A::AbstractVMF, B::CTranspose ) = A_ldiv_Bc!(O,A,B.obj)
ldiv!(O::AbstractVMF, A::Transpose  , B::CTranspose ) = At_ldiv_Bc!(O,A.obj,B.obj)
ldiv!(O::AbstractVMF, A::CTranspose , B::CTranspose ) = Ac_ldiv_Bc!(O,A.obj,B.obj)

rdiv!(O::AbstractVMF, A::AbstractVMF, B::AbstractVMF) = A_rdiv_B!(O,A,B)
rdiv!(O::AbstractVMF, A::Transpose  , B::AbstractVMF) = At_rdiv_B!(O,A.obj,B)
rdiv!(O::AbstractVMF, A::CTranspose , B::AbstractVMF) = Ac_rdiv_B!(O,A.obj,B)
rdiv!(O::AbstractVMF, A::AbstractVMF, B::Transpose  ) = A_rdiv_Bt!(O,A,B.obj)
rdiv!(O::AbstractVMF, A::Transpose  , B::Transpose  ) = At_rdiv_Bt!(O,A.obj,B.obj)
rdiv!(O::AbstractVMF, A::CTranspose , B::Transpose  ) = Ac_rdiv_Bt!(O,A.obj,B.obj)
rdiv!(O::AbstractVMF, A::AbstractVMF, B::CTranspose ) = A_rdiv_Bc!(O,A,B.obj)
rdiv!(O::AbstractVMF, A::Transpose  , B::CTranspose ) = At_rdiv_Bc!(O,A.obj,B.obj)
rdiv!(O::AbstractVMF, A::CTranspose , B::CTranspose ) = Ac_rdiv_Bc!(O,A.obj,B.obj)

badd!{N}(::Type{Inplace{N}}, As...) = broadcast!(+,As[N],As...)
bsub!{N}(::Type{Inplace{N}}, As...) = broadcast!(-,As[N],As...)
bmul!{N}(::Type{Inplace{N}}, As...) = broadcast!(*,As[N],As...)
bldiv!{N}(::Type{Inplace{N}}, As...) = broadcast!(\,As[N],As...)
brdiv!{N}(::Type{Inplace{N}}, As...) = broadcast!(/,As[N],As...)
add!{N}(t::Type{Inplace{N}}, As...) = _add!(t, As...)
sub!{N}(t::Type{Inplace{N}}, As...) = _sub!(t, As...)

badd!(O::AbstractArray, As...) = broadcast!(+,O,As...)
bsub!(O::AbstractArray, As...) = broadcast!(-,O,As...)
bmul!(O::AbstractArray, As...) = broadcast!(*,O,As...)
bldiv!(O::AbstractArray, As...) = broadcast!(\,O,As...)
brdiv!(O::AbstractArray, As...) = broadcast!(/,O,As...)
add!(O::AbstractArray, As...) = _add!(O, As...)
sub!(O::AbstractArray, As...) = _sub!(O, As...)

replace_t(ex) = esc(ex)
function replace_t(ex::Expr)
    if ex.head == Symbol("'")
        :(op_ctranspose($(esc(ex.args[1]))))
    elseif ex.head == Symbol(".'")
        :(op_transpose($(esc(ex.args[1]))))
    else
        esc(ex)
    end
end

macro in1!(ex)
    ex.head == :call || error("incorrect macro usage")
    Expr(:call,inplace_sym(ex.args[1]),:(Inplace{1}),[replace_t(a) for a in ex.args[2:end]]...)
end

macro in2!(ex)
    ex.head == :call || error("incorrect macro usage")
    Expr(:call,inplace_sym(ex.args[1]),:(Inplace{2}),[replace_t(a) for a in ex.args[2:end]]...)
end

macro into!(ex)
    ex.head == :(=) || error("incorrect macro usage")
    out, ex = ex.args
    Expr(:call,inplace_sym(ex.args[1]),esc(out),[replace_t(a) for a in ex.args[2:end]]...)
end

end  # module


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    UniformSampling(size, replace=false)

Sample elements uniformly from a given domain/data. Produce a
sample of given `size` with or without replacement depending on
the `replace` option.
"""
struct UniformSampling <: DiscreteSamplingMethod
  size::Int
  replace::Bool
end

UniformSampling(size::Int) = UniformSampling(size, false)

function sample(rng::AbstractRNG, Ω::DomainOrData, method::UniformSampling)
  n = nelements(Ω)
  s = method.size
  r = method.replace
  if s > n && r == false
    @error "invalid sample size for sampling without replacement"
  end
  view(Ω, sample(rng, 1:n, s, replace=r))
end


text:
include("helpers/listing.jl")
using Colors

function main(window)
    push!(window.assets, "widgets")
    push!(window.assets, "date")
    push!(window.assets, "codemirror")
    push!(window.assets, "tex")
    colors = ["violet", "blue", "green", "yellow", "orange", "red"]
    colortile(c, w, h) = fillcolor(c, size(w, h, empty))

    vbox(
        title(3, "Widgets"),
        md"""Escher provides a set of input widgets that can be used to create
        dynamic UIs that take inputs from users and are interactive. Widgets
        require the `widgets` asset to be pushed to the window to render and
        behave correctly. This is done with `push!(window.assets, "widgets")`.

        Escher follows the concept of reactive programming and uses Reactive.jl
        to implement this. The Escher workflow for creating dynamic UIs can be
        broken down into the following workflow.

        - Create a `Signal`
        - Create a `Widget` and subscribe the `Signal` to the `Widget`. Updates on the widgets will now reflect
        to the `Signal`.
        - Use `map` to update a UI based on the value of the `Signal`

        """,
        h3(md"Sliders"),
        md"""
        Sliders are very useful input widgets for taking inputs in a certain range.
        They are well suited to be used for users to input numerical inputs, such
        as the number of iterations for a simulation to run. Sliders can be created
        using the `slider` function and providing a range as the input. The starting
        value can also be provided using the `value` keyword argument.
        """,
        listing(
            """
            # Create a Signal
            inp = Signal(5.0)
            # Create a slider and subscribe it to the signal
            s = subscribe(
                    inp,
                    slider(1.0:9.0, value=5.0)
                )
            # Map the signal and create a dynamic UI
            map(inp) do val
                vbox(
                    "Adjust the slider to see the value update",
                    s, #the slider!
                    string(val) #the value
                )
            end
            """
        ),
        md"""The above pattern is followed for the rest of the input widgets in
        Escher too. We use the slider to present a more involved example below.
        """,
        listing(
            """
            # Colortile set up
            colors = ["violet", "blue", "green", "yellow", "orange", "red"]
            colortile(c, w, h) = fillcolor(c, size(w, h, empty))
            colortiles = [colortile(c, 4em, 4em) for c in colors]
            # Create a Signal
            n = Signal(1)
            # Shortcut to subscribe
            nslider = slider(1:length(colors)) >>> n
            # Mapping the signal and creating a dynamic UI
            map(n) do val
                vbox(
                    "Adjust the slider to display that many colortiles",
                    nslider,
                    vskip(1em),
                    hbox(
                    colortiles[1:val]...
                    )
                )
            end
            """
        ),
        md"""Notice how we use `>>>` as syntactic sugar instead of `subscribe`.""",
        h3("Buttons"),
        md"""Buttons can be created using `button` and providing a label. We can
        also set the button to be disabled, or raised using keywords. Buttons can
        be used in UI's to allow users to update the UI by clicking it, like run
        a simulation again.
        """,
        listing(
        """
            colors = ["violet", "blue", "green", "yellow", "orange", "red"]
            colortile(c, w, h) = fillcolor(c, size(w, h, empty))
            inp = Signal{Any}(nothing)
            b = button("Shuffle tiles") >>> inp
            map(inp) do _
                vbox(
                    b,
                    vskip(1em),
                    hbox(
                    [
                    colortile(color, 4em, 4em)
                    for color in shuffle(colors)
                    ]
                    )
                )
            end
        """
        ),
        h3(md"Checkbox"),
        md"""
        Checkboxes are a great way to take boolean inputs like allowing users to specify what content to display.
        They act like `Bool`s and can have values of either `true` and `false`.
        """,
        listing(
        """
            colors = ["violet", "blue", "green", "yellow", "orange", "red"]
            colortile(c, w, h) = fillcolor(c, size(w, h, empty))
            signals = [Signal(true) for i=1:length(colors)]
            checkboxes = [checkbox(true, colors[i]) >>> signals[i] for i=1:length(colors)]
            function checkboxtile(signal, checkbox, color)
                map(signal) do val
                    vbox(
                    checkbox,
                    if val
                        colortile(color, 4em, 4em)
                    else
                        colortile("white", 4em, 4em)
                    end
                    ) |> width(5em)
                end
            end

            hbox(
                [checkboxtile(signals[i], checkboxes[i], colors[i]) for i in 1:length(colors)]
            )
        """
        ),
        h3(md"Togglebuttons"),
        md"""
        Togglebuttons can be created using `togglebutton`. They are also like checkboxes,
        but do not have labels. They also act like `Bool`s and can have values of either `true` and `false`.
        """,
        listing(
        """
            colors = ["violet", "blue", "green", "yellow", "orange", "red"]
            colortile(c, w, h) = fillcolor(c, size(w, h, empty))
            colortiles = [colortile(c, 4em, 4em) for c in colors]
            disp = Signal(false)
            t = togglebutton(false) >>> disp
            map(disp) do disp
                vbox(
                    hbox("Hide",t),
                    vskip(1em),
                    if disp
                        "Tiles Hidden!"
                    else
                        hbox(colortiles)
                    end
                )
            end
        """),
        h3(md"Radio Groups"),
        md"""
        Radio buttons can be created using `radio`. We combine several of these
        `radio` buttons in a `radiogroup`. In one group, only one button can
        be selected at any one time. This can be used to provide a list of choices
        to users from which only one can be selected.
        """,
        listing(
        """
            colors = ["violet", "blue", "green", "yellow", "orange", "red"]
            colortile(c, w, h) = fillcolor(c, size(w, h, empty))

            selected = Signal("violet")
            r = radiogroup([radio(c, c) for c in colors], selected="violet") >>> selected
            map(selected) do c
                vbox(
                    r,
                    colortile(c, 4em, 4em)
                )
            end
        """
        ),
        h3(md"Text input"),
        md"""
        The text input widget can be used to take text input from users. It can
        be created using `textinput`. We can
        set pattern validation and display error messages in the case of invalid
        input. Character counters and multiline inputs can also set!
        """,
        listing(
            """
            text = Signal("")
            textwidget = textinput(
                label="Type here", multiline=true, charcounter=true
                ) >>> text
            map(text) do val
                vbox(
                    textwidget,
                    val
                )
            end
            """
        ),
        #=
        h3("Date Picker"),
        md"""
        Date pickers are a great way of allowing users to choose dates. `datepicker`
        can be used to create such a date picker in Escher.
        """,
        listing(
            """
            datechosen = Signal(Dates.today())
            datewidget = datepicker(Dates.today())
            map(datechosen) do d
                vbox(
                    datewidget,
                    d
                )
            end
            """
        ),=#
        h3(md"Codemirror"),
        md"""
        We can set up a codemirror input widget using `codemirror`. Syntax highlighting,
        tab indents, themes, line indents etc can be set using the keywork arguments. This requires
        us to push the `codemirror` asset first by doing `push!(window.assets, "codemirror")`.
        """,
        listing(
            """
            code = Signal("1+1")
            codewidget = codemirror("1+1") >>> code
            map(code) do code
                vbox(
                    codewidget,
                    eval(parse(code))
                )
            end
            """
        )
    ) |> pad(2em)
end


text:
import Oceananigans.Fields: set!

using Oceananigans.TimeSteppers: update_state!

function set!(model::ShallowWaterModel; kwargs...)
    for (fldname, value) in kwargs
        if fldname ∈ propertynames(model.solution)
            ϕ = getproperty(model.solution, fldname)
        elseif fldname ∈ propertynames(model.tracers)
            ϕ = getproperty(model.tracers, fldname)
        else
            throw(ArgumentError("name $fldname not found in model.solution or model.tracers."))
        end
        set!(ϕ, value)
    end

    update_state!(model)

    return nothing
end


text:
########## deconvolve ##########

function AbstractNFFTs.deconvolve!(p::NFFTPlan{T,D,1}, f::AbstractArray{U,D}, g::StridedArray{Complex{T},D}) where {D,T,U}
  if !p.params.storeApodizationIdx
    deconvolve_alloc_free!(p, f, g)
  else
    p.tmpVecHat[:] .= vec(f) .* p.windowHatInvLUT[1]
    g[p.deconvolveIdx] = p.tmpVecHat
  end
  return
end

function deconvolve_alloc_free!(p::NFFTPlan{T,D,1}, f::AbstractArray{U,D}, g::StridedArray{Complex{T},D}) where {D,T,U}
  if D == 1
    _deconvolve_alloc_free!(p, f, g, 1)
  else
    @cthreads for o = 1:p.N[end]
        _deconvolve_alloc_free!(p, f, g, o)
    end
  end
end

@generated function _deconvolve_alloc_free!(p::NFFTPlan{T,D,1}, f::AbstractArray{U,D}, g::StridedArray{Complex{T},D}, o) where {D,T,U}
  quote
    @nexprs 1 d -> gidx_{$D-1} = rem(o+p.n[$D] - p.N[$D]÷2 - 1, p.n[$D]) + 1
    @nexprs 1 d -> l_{$D-1} = o
    @nloops $(D-2) l d->(1:size(f,d+1)) d-> begin
        gidx_d = rem(l_d+p.n[d+1] - p.N[d+1]÷2 - 1, p.n[d+1]) + 1
      end begin
      N2 = p.N[1]÷2
      @inbounds @simd for i = 1:N2
        v = f[i, CartesianIndex(@ntuple $(D-1) l)]
        v *= p.windowHatInvLUT[1][i]
        @nexprs $(D-1) d -> v *= p.windowHatInvLUT[d+1][l_d]
        g[i-N2+p.n[1], CartesianIndex(@ntuple $(D-1) gidx)] = v

        v = f[i+N2, CartesianIndex(@ntuple $(D-1) l)]
        v *= p.windowHatInvLUT[1][i+N2]
        @nexprs $(D-1) d -> v *= p.windowHatInvLUT[d+1][l_d]
        g[i, CartesianIndex(@ntuple $(D-1) gidx)] = v
      end
    end
  end
end

########## deconvolve adjoint ##########

function AbstractNFFTs.deconvolve_transpose!(p::NFFTPlan{T,D,1}, g::AbstractArray{Complex{T},D}, f::StridedArray{U,D}) where {D,T,U}
  if !p.params.storeApodizationIdx
    deconvolve_transpose_alloc_free!(p, g, f)
  else
    p.tmpVecHat[:] = g[p.deconvolveIdx]
    f[:] .= vec(p.tmpVecHat) .* p.windowHatInvLUT[1]
  end
  return
end

function deconvolve_transpose_alloc_free!(p::NFFTPlan{T,D,1}, g::AbstractArray{Complex{T},D}, f::StridedArray{U,D}) where {D,T,U}
  if D == 1
    _deconvolve_transpose_alloc_free!(p, g, f, 1)
  else
    @cthreads for o = 1:p.N[end]
      _deconvolve_transpose_alloc_free!(p, g, f, o)
    end
  end
end

@generated function _deconvolve_transpose_alloc_free!(p::NFFTPlan{T,D,1}, g::AbstractArray{Complex{T},D}, f::StridedArray{U,D}, o) where {D,T,U}
  quote
    @nexprs 1 d -> gidx_{$D-1} = rem(o+p.n[$D] - p.N[$D]÷2 - 1, p.n[$D]) + 1
    @nexprs 1 d -> l_{$D-1} = o
    @nloops $(D-2) l d->(1:size(f,d+1)) d-> begin
        gidx_d = rem(l_d+p.n[d+1] - p.N[d+1]÷2 - 1, p.n[d+1]) + 1
      end begin
      N2 = p.N[1]÷2
      @inbounds @simd for i = 1:N2
        v = g[i-N2+p.n[1], CartesianIndex(@ntuple $(D-1) gidx)]
        v *= p.windowHatInvLUT[1][i]
        @nexprs $(D-1) d -> v *= p.windowHatInvLUT[d+1][l_d]
        f[i, CartesianIndex(@ntuple $(D-1) l)] = v

        v = g[i, CartesianIndex(@ntuple $(D-1) gidx)]
        v *= p.windowHatInvLUT[1][i+N2]
        @nexprs $(D-1) d -> v *= p.windowHatInvLUT[d+1][l_d]
        f[i+N2, CartesianIndex(@ntuple $(D-1) l)] = v
      end
    end
  end
end

########## convolve ##########

function AbstractNFFTs.convolve!(p::NFFTPlan{T,D,1}, g::AbstractArray{Complex{T},D}, fHat::StridedVector{U}) where {D,T,U}
  if isempty(p.B)
    if p.params.blocking
      convolve_blocking!(p, g, fHat)
    else
      convolve_nonblocking!(p, g, fHat)
    end
  else
    convolve_sparse_matrix!(p, g, fHat)
  end
  return
end

function convolve_nonblocking!(p::NFFTPlan{T,D,1}, g::AbstractArray{Complex{T},D}, fHat::StridedVector{U}) where {D,T,U}
  L = Val(2*p.params.m)
  scale = Int(p.params.LUTSize/(p.params.m+2))

  @cthreads for k in 1:p.M
      fHat[k] = _convolve_nonblocking(p, g, L, scale, k)
  end
end

function precomputeOneNode(p::NFFTPlan{T,D,1}, scale, k, d, L::Val{Z}) where {T,D,Z}
    return precomputeOneNode(p.windowLinInterp, p.x, p.n,
                    p.params.m, p.params.σ, scale, k, d, L, p.params.LUTSize)
end

@generated function _convolve_nonblocking(p::NFFTPlan{T,D,1}, g::AbstractArray{Complex{T},D}, L::Val{Z}, scale, k) where {D,T,Z}
  quote
    @nexprs $(D) d -> ((tmpIdx_d, tmpWin_d) = precomputeOneNode(p, scale, k, d, L) )

    fHat = zero(Complex{T})

    @nexprs 1 d -> prodWin_{$D} = one(T)
    @nloops_ $D l d -> 1:$Z d->begin
      # preexpr
      prodWin_{d-1} = prodWin_d * tmpWin_d[l_d]
      gidx_d = tmpIdx_d[l_d]
    end begin
      # bodyexpr
      fHat += prodWin_0 * @nref $D g gidx
    end
    return fHat
  end
end

function convolve_sparse_matrix!(p::NFFTPlan{T,D,1}, g::AbstractArray{Complex{T},D}, fHat::StridedVector{U}) where {D,T,U}
  threaded_mul!(fHat, transpose(p.B), vec(g))
  #mul!(fHat, transpose(p.B), vec(g))
end

########## convolve adjoint ##########

function AbstractNFFTs.convolve_transpose!(p::NFFTPlan{T,D,1}, fHat::AbstractVector{U}, g::StridedArray{Complex{T},D}) where {D,T,U}
  if isempty(p.B)
    if p.params.blocking
      convolve_transpose_blocking!(p, fHat, g)
    else
      convolve_transpose_nonblocking!(p, fHat, g)
    end
  else
    convolve_transpose_sparse_matrix!(p, fHat, g)
  end
end

function convolve_transpose_nonblocking!(p::NFFTPlan{T,D,1}, fHat::AbstractVector{U}, g::StridedArray{Complex{T},D}) where {D,T,U}
  fill!(g, zero(T))
  L = Val(2*p.params.m)
  scale = Int(p.params.LUTSize/(p.params.m+2))

  @inbounds @simd for k in 1:p.M
    _convolve_transpose_nonblocking!(p, fHat, g, L, scale, k)
  end
end

@generated function _convolve_transpose_nonblocking!(p::NFFTPlan{T,D,1}, fHat::AbstractVector{U}, g::StridedArray{Complex{T},D}, L::Val{Z}, scale, k) where {D,T,U,Z}
  quote
    @nexprs $(D) d -> ((tmpIdx_d, tmpWin_d) = precomputeOneNode(p, scale, k, d, L) )

    @nexprs 1 d -> prodWin_{$D} = one(T)
    @nloops_ $D l d -> 1:$Z d->begin
      # preexpr
      prodWin_{d-1} = prodWin_d * tmpWin_d[l_d]
      gidx_d = tmpIdx_d[l_d]
    end begin
      # bodyexpr
      (@nref $D g gidx) += prodWin_0 * fHat[k]
    end
  end
end

function convolve_transpose_sparse_matrix!(p::NFFTPlan{T,D,1},
                        fHat::AbstractVector{U}, g::StridedArray{Complex{T},D}) where {D,T,U}
  #threaded_mul!(vec(g), p.B, fHat)
  mul!(vec(g), p.B, fHat)
end

######## blocked multi-threading #########

function makePolyArrayStatic(p::NFFTPlan)
  if p.params.precompute == POLYNOMIAL
    P = p.windowPolyInterp
    return ntuple(d-> ntuple(g-> P[g,d], size(P,1)), size(P,2))
    #return SArray{Tuple{12,12}}(P) # Static Array Version
  else
    return nothing
  end
end

function convolve_blocking!(p::NFFTPlan, g, fHat)
  L = Val(2*p.params.m)
  winPoly = makePolyArrayStatic(p)
  _convolve_blocking!(p, g, fHat, L, winPoly)
end

function _convolve_blocking!(p::NFFTPlan{T,D,1}, g, fHat, L, winPoly) where {D,T}
  scale = Int(p.params.LUTSize/(p.params.m+2))

  @cthreads for l in CartesianIndices(size(p.blocks))
    if !isempty(p.nodesInBlock[l])
      toBlock!(p, g, p.blocks[l], p.blockOffsets[l])
      winTensor = (p.params.precompute == TENSOR) ? p.windowTensor[l] : nothing
      calcOneNode!(p, fHat, p.nodesInBlock[l], p.blocks[l], p.blockOffsets[l], L,
                   scale, p.idxInBlock[l], winTensor, winPoly)
    end
  end
end

@generated function toBlock!(p::NFFTPlan{T,D,1}, g, block, off) where {T,D}
quote
  if off[1] + 2 < 1
    LA = -(off[1] + 2) + 1
    offA = size(g,1) - LA
    offB = -LA
  else
    LA = 0
    offB = off[1] + 1
  end

  if offB + size(block,1) > size(g,1)
    LB = size(g,1) - (offB )
    offC = -LB
  else
    LB = size(block,1)
  end

  @nloops_ $(D-1)  (d->l_{d+1})  (d -> 1:size(block,d+1)) d->begin
    # preexpr
    idx_{d+1} = ( rem(l_{d+1}  + off[d+1] + p.n[d+1], p.n[d+1]) + 1)
  end begin
    # bodyexpr
    @inbounds for l_1 = 1:LA
      idx_1 = l_1 + offA
      (@nref $D block l) = (@nref $D g idx)
    end
    @inbounds for l_1 = (LA+1):LB
      idx_1 = l_1 + offB
      (@nref $D block l) = (@nref $D g idx)
    end
    @inbounds for l_1 = (LB+1):size(block,1)
      idx_1 = l_1 + offC
      (@nref $D block l) = (@nref $D g idx)
    end
  end
  return
 end
end

@noinline function calcOneNode!(p::NFFTPlan{T,D,1}, fHat, nodesInBlock, block,
         off, L, scale, idxInBlock, winTensor, winPoly) where {D,T}
  for (kLocal,k) in enumerate(nodesInBlock)
    fHat[k] = _convolve_nonblocking(p, block, off, L, scale, k, kLocal,
                               idxInBlock, winTensor, winPoly)
  end
  return
end

@generated function _convolve_nonblocking(p::NFFTPlan{T,D,1}, block, off, L::Val{Z}, scale,
                   k, kLocal, idxInBlock, winTensor, winPoly) where {D,T,Z}
  quote
    @nexprs $(D) d -> ((off_d, tmpWin_d) =  precomputeOneNodeBlocking(p.windowLinInterp, winTensor, winPoly, scale,
                      kLocal, d, L, idxInBlock) )

    fHat = zero(Complex{T})

    @nexprs 1 d -> prodWin_{$D} = one(T)
    @nloops_ $(D-1)  (d->l_{d+1})  (d -> 1:$Z) d->begin
      # preexpr
      prodWin_{d} = prodWin_{d+1} * tmpWin_{d+1}[l_{d+1}]
      block_idx_{d+1} = off_{d+1} + l_{d+1}
    end begin
      # bodyexpr
      @inbounds for l_1 = 1:$Z
        block_idx_1 = off_1 + l_1
        prodWin_0 = prodWin_1 * tmpWin_1[l_1]
        fHat += prodWin_0 * (@nref $D block block_idx)
      end
    end
    return fHat
  end
end

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

function convolve_transpose_blocking!(p::NFFTPlan, fHat, g)
  L = Val(2*p.params.m)
  winPoly = makePolyArrayStatic(p)
  _convolve_transpose_blocking!(p, fHat, g, L, winPoly)
end

function _convolve_transpose_blocking!(p::NFFTPlan{T,D,1}, fHat, g, L, winPoly) where {D,T}
  #g .= zero(T)
  ccall(:memset, Ptr{Cvoid}, (Ptr{Cvoid}, Cint, Csize_t), g, 0, sizeof(g))
  scale = Int(p.params.LUTSize/(p.params.m+2))

  lk = ReentrantLock()
  @cthreads for l in CartesianIndices(size(p.blocks))
    if !isempty(p.nodesInBlock[l])
      # p.blocks[l] .= zero(T)
      ccall(:memset, Ptr{Cvoid}, (Ptr{Cvoid}, Cint, Csize_t), p.blocks[l], 0, sizeof(p.blocks[l]))

      winTensor = (p.params.precompute == TENSOR) ? p.windowTensor[l] : nothing
      fillBlock!(p, fHat, p.blocks[l], p.nodesInBlock[l], p.blockOffsets[l], L, scale,
                 p.idxInBlock[l], winTensor, winPoly)

      lock(lk) do
        addBlock!(p, g, p.blocks[l], p.blockOffsets[l])
      end
    end
  end
end

@generated function addBlock!(p::NFFTPlan{T,D,1}, g, block, off) where {T,D}
  quote
    # addBlock! needs to wrap indices. The easies is to apply a modulo operation (rem)
    # but doing so for each index is slow. We thus only do this for the outer dimensions
    # For the inner dimension we determine the indices where a wrap occurs. This can happen
    # 0 time, 1 time, or 2 times. The last is an extreme case where the block size matches
    # the dimension of g.
    #
    # The following code divides the range 1:size(block,1) into up to three parts and then
    # applies three dedicated for loops. The wrap indices are LA and LB.

    if off[1] + 2 < 1
      # negative indices: first wrapping
      LA = -(off[1] + 2) + 1
      offA = size(g,1) - LA
      offB = -LA
    else
      # no negative indices: skip first sum
      LA = 0
      offB = off[1] + 1
    end

    if offB + size(block,1) > size(g,1)
      # out of bounds indices: second (or first) wrapping
      LB = size(g,1) - (offB )
      offC = -LB
    else
      # no out of bounds indices
      LB = size(block,1)
    end

    @nloops_ $(D-1)  (d->l_{d+1})  (d -> 1:size(block,d+1)) d->begin
      # preexpr
      idx_{d+1} = ( rem(l_{d+1}  + off[d+1] + p.n[d+1], p.n[d+1]) + 1)
    end begin
      # bodyexpr
      @inbounds for l_1 = 1:LA
        idx_1 = l_1 + offA
        (@nref $D g idx) += (@nref $D block l)
      end
      @inbounds for l_1 = (LA+1):LB
        idx_1 = l_1 + offB
        (@nref $D g idx) += (@nref $D block l)
      end
      @inbounds for l_1 = (LB+1):size(block,1)
        idx_1 = l_1 + offC
        (@nref $D g idx) += (@nref $D block l)
      end
    end
    return
  end
end

@noinline function fillBlock!(p::NFFTPlan{T,D,1}, fHat, block, nodesInBlock, off, L::Val{Z}, scale,
                              idxInBlock, winTensor, winPoly) where {T,D,Z}
  if (D >= 3 && Z >= 10) || (D == 2 && Z >= 16) # magic
    for (kLocal,k) in enumerate(nodesInBlock)
      fillOneNodeReal!(p, fHat, block, off, L, scale, k, kLocal, idxInBlock, winTensor, winPoly)
    end
  else
    for (kLocal,k) in enumerate(nodesInBlock)
      fillOneNode!(p, fHat, block, off, L, scale, k, kLocal, idxInBlock, winTensor, winPoly)
    end
  end
  return
end

@generated function fillOneNode!(p::NFFTPlan{T,D,1}, fHat, block, off, L::Val{Z}, scale,
              k, kLocal, idxInBlock, winTensor, winPoly) where {D,T,Z}
  quote
    fHat_ = fHat[k]

    @nexprs $(D) d -> ((off_d, tmpWin_d) = precomputeOneNodeBlocking(p.windowLinInterp, winTensor, winPoly, scale, kLocal, d, L,
                                             idxInBlock) )

    innerWin = @ntuple $(Z) l -> tmpWin_1[l] * fHat_

    @nexprs 1 d -> prodWin_{$D} = one(T)
    @nloops_ $(D-1)  (d->l_{d+1})  (d -> 1:$Z) d->begin
      # preexpr
      prodWin_{d} = prodWin_{d+1} * tmpWin_{d+1}[l_{d+1}]
      block_idx_{d+1} = off_{d+1} + l_{d+1}
    end begin
      # bodyexpr
      @inbounds @simd for l_1 = 1:$(Z)
        block_idx_1 = off_1 + l_1
        (@nref $D block block_idx) += innerWin[l_1] * prodWin_1
      end
    end
    return
  end
end

@generated function fillOneNodeReal!(p::NFFTPlan{T,D,1}, fHat, block, off, L::Val{Z}, scale, k,
                kLocal, idxInBlock, winTensor, winPoly) where {D,T,Z}
  quote
    fHat_ = fHat[k]

    @nexprs $(D) d -> ((off_d, tmpWin_d) =  precomputeOneNodeBlocking(p.windowLinInterp, winTensor, winPoly, scale, kLocal,
                                              d, L, idxInBlock) )

    innerWin = Vector{Float64}(undef,$(2*Z)) # Probably cache me somewhere
    for l=1:$Z
      innerWin[2*(l-1)+1] = tmpWin_1[l] * real(fHat_)
      innerWin[2*(l-1)+2] = tmpWin_1[l] * imag(fHat_)
    end
    #innerWin = @ntuple $(2*Z) l -> isodd(l) ? tmpWin_1[(l-1)÷2+1] * real(fHat_) : tmpWin_1[(l-1)÷2+1] * imag(fHat_)

    # Convert to real
    # We use unsafe_wrap here because reinterpret(T, block) is much slower
    # and reinterpret(T, reshape, block) is slightly slower than unsafe_wrap
    # In fact the conversion to real only pays off when using unsafe_wrap,
    # otherwise sticking with complex is faster
    U = @ntuple $(D) d -> d==1 ? 2*size(block,1) : size(block,d)
    blockReal = unsafe_wrap(Array,reinterpret(Ptr{T},pointer(block)), U)
    off_1 *= 2

    @nexprs 1 d -> prodWin_{$D} = one(T)
    @nloops_ $(D-1)  (d->l_{d+1})  (d -> 1:$Z) d->begin
      # preexpr
      prodWin_{d} = prodWin_{d+1} * tmpWin_{d+1}[l_{d+1}]
      block_idx_{d+1} = off_{d+1} + l_{d+1}
    end begin
      # bodyexpr
      @inbounds @simd for l_1 = 1:$(2*Z)
        block_idx_1 = off_1 + l_1
        (@nref $(D) blockReal block_idx) += innerWin[l_1] * prodWin_1
      end
    end
    return
  end
end