text:
module ExponentialLinearQuadraticTest

using Test
using ReactiveMP
using Distributions
using Random
using DomainIntegrals

@testset "ExponentialLinearQuadratic" begin

    @testset "Statistics" begin

        approximation = GaussHermiteCubature(21)
        rng = MersenneTwister(123)

        for _ in 1:3
            a, b, c, d = rand(rng, 4)
            dist = ExponentialLinearQuadratic(approximation, a, b, c, d)

            μ, v = ReactiveMP.approximate_meancov(approximation, (x) -> exp(-(a * x + b * exp(c * x + d * x ^ 2 / 2)) / 2) * exp(x^2/2), NormalMeanVariance(0.0, 1.0))

            σ = sqrt(v)
            w = inv(v)
            ξ = w * μ

            @test mean(dist)         ≈ μ
            @test var(dist)          ≈ v
            @test cov(dist)          ≈ v
            @test std(dist)          ≈ σ
            @test weightedmean(dist) ≈ ξ
            @test invcov(dist)       ≈ w
            @test precision(dist)    ≈ w

            @test all(mean_var(dist)       .≈ (μ, v))
            @test all(mean_cov(dist)       .≈ (μ, v))
            @test all(mean_invcov(dist)    .≈ (μ, w))
            @test all(mean_precision(dist) .≈ (μ, w))
            @test all(mean_std(dist)       .≈ (μ, σ))

            @test all(weightedmean_var(dist)       .≈ (ξ, v))
            @test all(weightedmean_cov(dist)       .≈ (ξ, v))
            @test all(weightedmean_invcov(dist)    .≈ (ξ, w))
            @test all(weightedmean_precision(dist) .≈ (ξ, w))
            @test all(weightedmean_std(dist)       .≈ (ξ, σ))

        end

    end

    @testset "prod" begin

        approximation = GaussHermiteCubature(51)
        rng = MersenneTwister(1234)

        for _ in 1:3

            a, b, c, d = rand(rng, 4)
            left = ExponentialLinearQuadratic(approximation, a, b, c, d)

            μ, v = 5.0 * randn(rng), 5.0 * rand(rng)
            right = NormalMeanVariance(μ, v)

            q = NormalMeanVariance(ReactiveMP.approximate_meancov(approximation, (x) -> exp(logpdf(left, x) + logpdf(right, x) + x^2/2), NormalMeanVariance(0.0, 1.0))...)

            @test all(isapprox.(mean_var(q), mean_var(prod(ProdAnalytical(), left, right)), atol = 1e-2))
            @test all(isapprox.(mean_var(q), mean_var(prod(ProdAnalytical(), right, left)), atol = 1e-2))

        end

    end

end

end


text:
@testset "Correlator" begin

    @testset "Correlator constructor" begin
        gpsl1 = GPSL1()
        correlator = @inferred EarlyPromptLateCorrelator()
        correlator_sample_shifts = SVector(-1, 0, 1)
        early_late_index_shift = 1
        @test @inferred(get_early(correlator, correlator_sample_shifts, early_late_index_shift)) == 0.0
        @test @inferred(get_prompt(correlator, correlator_sample_shifts)) == 0.0
        @test @inferred(get_late(correlator, correlator_sample_shifts, early_late_index_shift)) == 0.0

        correlator = @inferred EarlyPromptLateCorrelator(NumAnts(1))

        @test @inferred(get_early(correlator, correlator_sample_shifts, early_late_index_shift)) == 0.0
        @test @inferred(get_prompt(correlator, correlator_sample_shifts)) == 0.0
        @test @inferred(get_late(correlator, correlator_sample_shifts, early_late_index_shift)) == 0.0

        correlator = @inferred EarlyPromptLateCorrelator(NumAnts(2))

        @test @inferred(get_early(correlator, correlator_sample_shifts, early_late_index_shift)) == SVector(0.0 + 0.0im, 0.0 + 0.0im)
        @test @inferred(get_prompt(correlator, correlator_sample_shifts)) == SVector(0.0 + 0.0im, 0.0 + 0.0im)
        @test @inferred(get_late(correlator, correlator_sample_shifts, early_late_index_shift)) == SVector(0.0 + 0.0im, 0.0 + 0.0im)
    end
    @testset "Zeroing correlator" begin
        correlator = @inferred EarlyPromptLateCorrelator(
            SVector(
                SVector(1.0 + 0.0im, 1.0 + 0.0im),
                SVector(1.0 + 0.0im, 1.0 + 0.0im),
                SVector(1.0 + 0.0im, 1.0 + 0.0im)
            )
        )

        @test @inferred(zero(correlator)) == EarlyPromptLateCorrelator(
            SVector(
                SVector(0.0 + 0.0im, 0.0 + 0.0im),
                SVector(0.0 + 0.0im, 0.0 + 0.0im),
                SVector(0.0 + 0.0im, 0.0 + 0.0im)
            )
        )

        correlator = @inferred EarlyPromptLateCorrelator(
            SVector(
                1.0 + 0.0im,
                1.0 + 0.0im,
                1.0 + 0.0im
            )
        )

        @test @inferred(zero(correlator)) == EarlyPromptLateCorrelator(
            SVector(
                0.0 + 0.0im,
                0.0 + 0.0im,
                0.0 + 0.0im
            )
        )
    end
    @testset "Filter correlator" begin
        correlator = @inferred EarlyPromptLateCorrelator(
            SVector(
                SVector(1.0 + 0.0im, 1.0 + 0.0im),
                SVector(1.0 + 0.0im, 1.0 + 0.0im),
                SVector(1.0 + 0.0im, 1.0 + 0.0im)
            )
        )
        filtered_correlator = @inferred Tracking.filter(x -> x[1], correlator)
        @test filtered_correlator == EarlyPromptLateCorrelator(
            SVector(
                1.0 + 0.0im,
                1.0 + 0.0im,
                1.0 + 0.0im
            )
        )

        correlator = @inferred EarlyPromptLateCorrelator(
            SVector(
                1.0 + 0.0im,
                1.0 + 0.0im,
                1.0 + 0.0im
            )
        )
        filtered_correlator = @inferred Tracking.filter(x -> 2 * x, correlator)
        @test filtered_correlator == EarlyPromptLateCorrelator(
            SVector(
                2.0 + 0.0im,
                2.0 + 0.0im,
                2.0 + 0.0im
            )
        )
    end
    @testset "Correlator sample shifts" begin
        gpsl1 = GPSL1()
        correlator = EarlyPromptLateCorrelator()
        correlator_sample_shifts = @inferred get_correlator_sample_shifts(
            gpsl1,
            correlator,
            4e6Hz,
            0.5
        )
        @test correlator_sample_shifts ≈ round(0.5 * 4e6 / 1023e3) * SVector(-1, 0, 1)
        @test typeof(correlator_sample_shifts) <: SVector{3,<:Integer}
    end
    @testset "Early late index shift" begin
        gpsl1 = GPSL1()
        correlator = EarlyPromptLateCorrelator()
        correlator_sample_shifts = get_correlator_sample_shifts(
            gpsl1,
            correlator,
            4e6Hz,
            0.5
        )
        early_late_index_shift = @inferred get_early_late_index_shift(
           gpsl1,
           correlator_sample_shifts,
           correlator,
           4e6Hz,
           0.5
        )

        @test early_late_index_shift == 1
    end
    @testset "Early late sample spacing" begin
        @test get_early_late_sample_spacing(SVector(-2,0,2), 1) == 4
    end
    @testset "Normalize correlator" begin
        correlator = @inferred EarlyPromptLateCorrelator(
            SVector(
                1.0 + 0.0im,
                1.0 + 0.0im,
                1.0 + 0.0im
            )
        )
        normalized_correlator = @inferred Tracking.normalize(correlator, 10)
        @test normalized_correlator == EarlyPromptLateCorrelator(
            SVector(
                0.1 + 0.0im,
                0.1 + 0.0im,
                0.1 + 0.0im
            )
        )

        correlator = EarlyPromptLateCorrelator(
            SVector(
                SVector(1.0 + 0.0im, 1.0 + 0.0im),
                SVector(1.0 + 0.0im, 1.0 + 0.0im),
                SVector(1.0 + 0.0im, 1.0 + 0.0im)
            )
        )
        normalized_correlator = @inferred Tracking.normalize(correlator, 10)
        @test normalized_correlator == EarlyPromptLateCorrelator(
            SVector(
                SVector(0.1 + 0.0im, 0.1 + 0.0im),
                SVector(0.1 + 0.0im, 0.1 + 0.0im),
                SVector(0.1 + 0.0im, 0.1 + 0.0im)
            )
        )
    end
    @testset "Correlate" begin
        gpsl1 = GPSL1()
        signal = StructArray{Complex{Float32}}(
            (Float32.(get_code.(gpsl1, (1:2500) * 1023e3 / 2.5e6, 1)),
            zeros(Float32, 2500))
        )
        correlator_sample_shifts = SVector(-1, 0, 1)
        code = get_code.(
            gpsl1,
            (1 + correlator_sample_shifts[1]:2500 + correlator_sample_shifts[end]) * 1023e3 / 2.5e6,
            1
        )
        correlator = EarlyPromptLateCorrelator(SVector(0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im))

        correlator_result = Tracking.correlate(
            correlator,
            signal,
            code,
            correlator_sample_shifts,
            1,
            2500
        )
        @test get_early(correlator_result, correlator_sample_shifts, 1) ==
            get_late(correlator_result, correlator_sample_shifts, 1)
        @test get_prompt(correlator_result, correlator_sample_shifts) == 2500

        signal_mat = repeat(signal, outer = (1,3))
        correlator_sample_shifts = SVector(-1,0,1)
        correlator = EarlyPromptLateCorrelator(NumAnts(3))

        correlator_result = Tracking.correlate(
            correlator,
            signal_mat,
            code,
            correlator_sample_shifts,
            1,
            2500,
        )
        @test get_early(correlator_result, correlator_sample_shifts, 1) == get_late(correlator_result, correlator_sample_shifts, 1)
        @test all(get_early(correlator_result, correlator_sample_shifts, 1) .== 1476)
        @test all(get_prompt(correlator_result, correlator_sample_shifts) .== 2500)
    end

end

ext: jl


text:
using ClusterTrees, BEAST, CompScienceMeshes

a, h = 1.0, 0.061
Γ1 = meshsphere(a, h)
Γ2 = Γ1
#Γ2 = CompScienceMeshes.translate(Γ1, point(2.10,0,0))

X1 = raviartthomas(Γ1)
X2 = X1
#X2 = raviartthomas(Γ2)

@show numfunctions(X1)
@show numfunctions(X2)

κ = 1.0; γ = κ*im;
T = Maxwell3D.singlelayer(gamma=γ)

μ = (τ,σ) -> assemble(T, subset(X1,τ), subset(X2,σ))
A = ClusterTrees.h1compress(T,X1,X2)

sz = storedentries(A)
println("Compression: ", sz/numfunctions(X1)/numfunctions(X2))


text:
using Clustering

srand(1)

n = 200

x = vcat(randn(n, 2), randn(n, 2) .+ 10)
true_assignments = vcat(zeros(n), ones(n))

results = k_means(x, 5)


text:
using Elefridge
using ZfpCompression
using PyPlot
using PyCall
using Blosc
using ColorSchemes
using StatsBase, Statistics
using TranscodingStreams, CodecZstd
using BenchmarkTools
BenchmarkTools.DEFAULT_PARAMETERS.samples = 3
xr = pyimport("xarray")
ccrs = pyimport("cartopy.crs")

path = "/Users/milan/cams/gridded/"
filelist = filter(x->endswith(x,"_go3.grib"),readdir(path))
grib = xr.open_dataarray(joinpath(path,filelist[end]),engine="cfgrib")
X = grib.data
lat = grib.latitude.data
lon = grib.longitude.data

X = permutedims(X,[3,2,1])      # for lossless longitude first

## DEFINE COMPRESSORS
ZstdCompressors = Array{ZstdCompressor,1}(undef,22)

for i in 1:22
    ZstdCompressors[i] = ZstdCompressor(level=i)
    TranscodingStreams.initialize(ZstdCompressors[i])
end

## Compress X with round+Zstd
Xr = round(X[1:450,:,:],5)     # temperature contains 7 informationful signficant bits
Xr8 = unsafe_wrap(Array, Ptr{UInt8}(pointer(Xr)), sizeof(Xr))

# cfs = fill(0.0,22)    # compression factors
# speed_comp = fill(0.0,22)         # compression speeds and decompression
# speed_decomp = fill(0.0,22)

for i in 11:17   # loop over compression levels
    println(i)
    Xc = transcode(ZstdCompressors[i],Xr8)
    cfs[i] = sizeof(Xr)/sizeof(Xc)*2/(900*451)*(348528)    # relative to Float64 and unstructured grid
    speed_comp[i] = (@benchmark Xc = transcode(ZstdCompressors[$i],Xr8)).times[1]/1e9
    speed_decomp[i] = (@benchmark Xuc = transcode(ZstdDecompressor,$Xc)).times[1]/1e9
    # @assert Xuc == Xr8      # check whether compression is indeed reversible
end

## Compress with Zfp
BenchmarkTools.DEFAULT_PARAMETERS.samples = 10000
Xs = permutedims(X[1:450,:,:],[3,2,1])      # for zfp vert x lat x lon

Xzfp = zfp_compress(Xs,precision=10)
speed_comp_zfp = (@benchmark Xzfp = zfp_compress(Xs,precision=13)).times[1]/1e9
cfs_zfp = sizeof(Xs)/sizeof(Xzfp)*2/(900*451)*(348528)
speed_decomp_zfp = (@benchmark zfp_decompress(Xzfp)).times[1]/1e9

## Compress with LinQuant24
Xlin24 = LinQuant16Array(Xs)
Xlog16 = LogQuant16Array(Xs)

speed_comp_lin24 = (@benchmark Xlin24 = LinQuant16Array(Xs)).times[1]/1e9       # use 16 as UInt24 is slow
speed_comp_log16 = (@benchmark Xlog16 = LogQuant16Array(Xs)).times[1]/1e9

speed_decomp_lin24 = (@benchmark Array{Float32}(Xlin24)).times[1]/1e9
speed_decomp_log16 = (@benchmark Array{Float32}(Xlog16)).times[1]/1e9

## Plotting
ioff()
fig,(ax1,ax2,ax3) = subplots(1,3,sharex=true,figsize=(10,3))

# ax1.grid(axis="y")
# ax1.set_axisbelow(true)

zstd_x = 3 .+ Array(0:21)/3.0
alfa = 0.7

ax1.bar(1,64/24,color="grey",edgecolor="k")
ax1.bar(2,4,color="k")
ax1.bar(zstd_x,cfs,alpha=alfa,color="C1",width=0.25,edgecolor="k",lw=0.5)
ax1.bar(11,cfs_zfp,color="C2",edgecolor="k")

# speed in MB/s
ax2.bar(1,sizeof(Xr)/speed_comp_lin24/1e6,color="grey",edgecolor="k")
ax2.bar(2,sizeof(Xr)/speed_comp_log16/1e6,color="k")
ax2.bar(zstd_x,sizeof(Xr)./speed_comp/1e6,alpha=alfa,color="C1",width=0.25,edgecolor="k",lw=0.5)
ax2.bar(11,sizeof(Xr)/speed_comp_zfp/1e6,color="C2",edgecolor="k")

ax3.bar(1,sizeof(Xr)/speed_decomp_lin24/1e6,color="grey",edgecolor="k")
ax3.bar(2,sizeof(Xr)/speed_decomp_log16/1e6,color="k")
ax3.bar(zstd_x,sizeof(Xr)./speed_decomp/1e6,alpha=alfa,color="C1",width=0.25,edgecolor="k",lw=0.5)
ax3.bar(11,sizeof(Xr)/speed_decomp_zfp/1e6,color="C2",edgecolor="k")

ax2.set_yscale("log")
ax3.set_yscale("log")

ax1.set_ylim(1,60)
ax2.set_ylim(1,5000)
ax3.set_ylim(1,5000)
ax1.set_xlim(0,12)

ax1.set_xticks(zstd_x[[1,22]])
ax1.set_xticks(zstd_x[2:end-1],minor=true)
ax1.set_xticklabels([1,22])

ax1.text(1,8,"LinQuant24",rotation=90,ha="center")
ax1.text(2,8,"LogQuant16",rotation=90,ha="center")
ax1.text(11,52,"Zfp",rotation=90,ha="center")
ax1.text(4,26,"Zstandard,\nlevel 1-22")

ax2.text(1,2,"LinQuant24",rotation=90,ha="center",color="w")
ax2.text(2,2,"LogQuant16",rotation=90,ha="center",color="w")
ax2.text(11,600,"Zfp",rotation=90,ha="center")
ax2.text(5.6,140,"Zstandard,\nlevel 1-22")

ax3.text(1,2,"LinQuant24",rotation=90,ha="center",color="w")
ax3.text(2,2,"LogQuant16",rotation=90,ha="center",color="w")
ax3.text(11,1200,"Zfp",rotation=90,ha="center")
ax3.text(4,1200,"Zstandard,\nlevel 1-22")

ax1.set_yticks(10:10:60)
ax1.set_yticklabels(["10x","20x","30x","40x","50x","60x"])

ax2.set_yticks([1,10,100,1000,5000])
ax2.set_yticklabels([1,10,100,1000,"[MB/s]"])
ax3.set_yticks([1,10,100,1000,5000])
ax3.set_yticklabels([1,10,100,1000,"[MB/s]"])

# ax2.set_ylabel("MB/s")
# ax3.set_ylabel("MB/s")

ax1.set_title("Compression factors",loc="left")
ax2.set_title("Compression speed",loc="left")
ax3.set_title("Decompression speed",loc="left")

ax1.set_title("a",loc="right",fontweight="bold")
ax2.set_title("b",loc="right",fontweight="bold")
ax3.set_title("c",loc="right",fontweight="bold")

tight_layout()
savefig("/Users/milan/git/Elefridge.jl/plots/level_speed_o3.png",dpi=200)
close(fig)


text:
module States
export State, getstate, putstate

"""
    State(func)

    State() do state
      ....
      return value, newstate
    end

State monad, which capsulate a state within a monadic type for monadic encapsulation of the state-handling.

You can run a state by either calling it, or by using `Base.run`. If no initial state is given, `nothing` is used.
"""
struct State{T}
  "s -> (a, s)"
  func::T
end

# running a State
function (s::State)(state = nothing)
  s.func(state)
end
Base.run(state::State, initial_state = nothing) = state(initial_state)

@doc raw"""
    getstate

Standard value for returning the hidden state of the `State` Monad.

Examples
--------

```jldoctest
julia> using TypeClasses

julia> mystate = @syntax_flatmap begin
         state = getstate
         @pure println("state = $state")
       end;

julia> mystate(42)
state = 42
(nothing, 42)
```
"""
const getstate = State() do state
  state, state
end

@doc raw"""
    putstate(x)

`putstate` is a standard constructor for `State` objects which changes the underlying state to the given value.

Examples
--------

```jldoctest
julia> using TypeClasses

julia> mystate = @syntax_flatmap begin
         putstate(10)
         state = getstate
         @pure println("The state is $state, and should be 10")
       end;

julia> mystate()
The state is 10, and should be 10
(nothing, 10)
```
"""
putstate(x) = State() do state
  (), x
end

end # module


text:
function read_input(input_filename)
    lines = readlines(input_filename)
    grid = zeros(Int, length(lines), length(lines[1]))
    for iy in 1:length(lines)
        for ix in 1:length(lines[iy])
            grid[iy, ix] = parse(Int, lines[iy][ix])
        end
    end
    return grid
end

function is_local_minimum(x, y, grid)
#    local_min_x = false
    if x == 1
        local_min_x = (grid[x,y] < grid[x+1,y])
    elseif x == size(grid)[1]
        local_min_x = (grid[x,y] < grid[x-1,y])
    else
        local_min_x = (grid[x,y] < grid[x+1,y]) && (grid[x,y] < grid[x-1,y])
    end
#    local_min_y = false
    if y == 1
        local_min_y = (grid[x,y] < grid[x,y+1])
    elseif y == size(grid)[2]
        local_min_y = (grid[x,y] < grid[x,y-1])
    else
        local_min_y = (grid[x,y] < grid[x,y+1]) && (grid[x,y] < grid[x,y-1])
    end
    return (local_min_x && local_min_y)
end

function find_all_local_minima(grid)
    minima_coords = []
    for ix in 1:size(grid)[1]
        for iy in 1:size(grid)[2]
            if is_local_minimum(ix, iy, grid)
                push!(minima_coords,(ix,iy))
            end
        end
    end
    return minima_coords
end

function get_total_risk(grid)
    minima = find_all_local_minima(grid)
    total = 0
    for minimum in minima
        total += grid[minimum[1],minimum[2]] + 1
    end
    return total
end


text:
# ---
# title: 718. Maximum Length of Repeated Subarray
# id: problem718
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array, Hash Table, Binary Search, Dynamic Programming
# link: <https://leetcode.com/problems/maximum-length-of-repeated-subarray/description/>
# hidden: true
# ---
#
# Given two integer arrays `A` and `B`, return the maximum length of an subarray
# that appears in both arrays.
#
# **Example 1:**
#
#
#
#     Input:
#     A: [1,2,3,2,1]
#     B: [3,2,1,4,7]
#     Output: 3
#     Explanation:
#     The repeated subarray with maximum length is [3, 2, 1].
#
#
#
#
# **Note:**
#
#   1. 1 <= len(A), len(B) <= 1000
#   2. 0 <= A[i], B[i] < 100
#
#
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
## reader
function Base.parse(::Type{Module}, membuf::MemoryBuffer)
    out_ref = Ref{API.LLVMModuleRef}()

    status = convert(Core.Bool, API.LLVMParseBitcode2(membuf, out_ref))
    @assert !status # caught by diagnostics handler

    Module(out_ref[])
end

function Base.parse(::Type{Module}, membuf::MemoryBuffer, ctx::Context)
    out_ref = Ref{API.LLVMModuleRef}()

    status = convert(Core.Bool, API.LLVMParseBitcodeInContext2(ctx, membuf, out_ref))
    @assert !status # caught by diagnostics handler

    Module(out_ref[])
end

Base.parse(::Type{Module}, data::Vector) =
    parse(Module, MemoryBuffer(data, "", false))
Base.parse(::Type{Module}, data::Vector, ctx::Context) =
    parse(Module, MemoryBuffer(data, "", false), ctx)

## writer

Base.write(io::IOStream, mod::Module) =
    API.LLVMWriteBitcodeToFD(mod, Cint(fd(io)), convert(Bool, false), convert(Bool, true))

Base.convert(::Type{MemoryBuffer}, mod::Module) = MemoryBuffer(
    API.LLVMWriteBitcodeToMemoryBuffer(mod))

Base.convert(::Type{Vector{T}}, mod::Module) where {T<:Union{UInt8,Int8}} =
    convert(Vector{T}, convert(MemoryBuffer, mod))


text:
# Quadratic line minimization as suggested in original Beigi-Arias paper

function linmin1( pw, Vpot, α_t, psi, d, E, g )
    # compute directional derivative
    Δ = 2.0*real( sum( conj(d).*g ) )
    #
    psi_t = ortho_gram_schmidt( psi + α_t*d )
    E_trial = calc_Etot( pw, Vpot, psi_t )
    curvature = ( E_trial - ( E + α_t*Δ ) ) /α_t^2
    α = -Δ/(2*curvature)

    if abs(curvature) < eps()
        α = 0.0
    end

    return abs(α)   # need to use absolute value of α !!!
end


text: abstract type AbstractFixedPointSolver{T} <: NonlinearSolver{T} end


text:
mutable struct AdjacencyMatrix <: AbstractGraph
    cntedge::Int
    cntvertex::Int
    arr::AbstractMatrix # Should eventually be of type Vertex
    mark::Vector # 0...unvisited ; 1...visited
    function AdjacencyMatrix(n::Int)
        cntedge = 0
        cntvertex = n
        arr = zeros(Int, cntvertex, cntvertex)
        mark = zeros(Int, cntvertex)
        new(cntedge, cntvertex, arr, mark)
    end
end

countedges(g::AdjacencyMatrix) = g.cntedge
countvertices(g::AdjacencyMatrix) = g.cntvertex

function first(g::AdjacencyMatrix, v::Vertex)
    for i = 1:g.cntvertex
        (g.arr[v, i] != 0) && return i
    end
    return g.cntvertex + 1
end

function next(g::AdjacencyMatrix, v::Vertex, w::Vertex)
    for i = w+1:g.cntvertex
        (g.arr[v, i] != 0) && return i
    end
    return g.cntvertex + 1
end

function setedge!(g::AdjacencyMatrix, v::Vertex, w::Vertex, weight)
    (weight > 0) || throw(ArgumentError("Illegal weight value."))
    (g.arr[v, w] == 0) && (g.cntedge += 1)
    g.arr[v, w] = weight
end

function remove!(g::AdjacencyMatrix, v::Vertex, w::Vertex)
    (g.arr[v, w] != 0) && (g.cntedge -= 1)
    g.arr[v, w] = 0
end

isedge(g::AdjacencyMatrix, v::Vertex, w::Vertex) = g.arr[v, w] != 0
getweight(g::AdjacencyMatrix, v::Vertex, w::Vertex) = g.arr[v, w]
getmark(g::AdjacencyMatrix, v::Vertex) = g.mark[v]
setmark!(g::AdjacencyMatrix, v::Vertex, value) = (g.mark[v] = value)

==(a::AdjacencyMatrix, b::AdjacencyMatrix) = a.arr == b.arr


text:
#!/usr/bin/env julia

function load_login_passwords(filename)
    lines = readlines(filename)
    logins = Dict()

    for l in lines
        user, password = split(l, ':')
        logins[user] = password
    end

    logins
end

function login()
    login_passwords = load_login_passwords("passwords.txt")

    for i in 1:3
        println("Login:")
        user = readline(stdin)

        if haskey(login_passwords, user)
            println("Password:")
            password = readline(stdin)

            correct_password = login_passwords[user]
            if password == correct_password
                println("Correct! You are now logged in")
                break
            else
                println("Wrong password, try again...")
            end
        else
            println("There is no login account with name $user")
        end
    end

end

login()


text: include("data.jl")


text:
using Plots; gr()

# Random Walk model
M = LGSSM(diagm(1.), diagm(1.), LowerTriangular(diagm(1.)), LowerTriangular(diagm(1.)), diagm(1.))
x = cumsum(randn(50))
y = x + randn(50)
KF = DisturbanceSmoother(M, y)

plot(x)
plot!(KF[:Xf][:,1,:]')
plot!(KF[:Xs][:,1,:]')

# Local level linear model
M = LGSSM([1. 1.; 0. 1.], [1. 0.], LowerTriangular(diagm([1.,1.])), LowerTriangular(diagm(1.)) , eye(2))
S = Generate(M, 100, 1)
plot(S[:Y][:,:,1]')
plot(S[:X][:,:,1]')

# Structural Time Series Model
# Local level model
# Seasonal component with s = 4
A = blkdiag(sparse([1],[1],[1.]), sparse([1,1,1,2,3], [1,2,3,1,2],[-1.,-1.,-1.,1.,1.]))
B = [1. 1. 0. 0.]
S = LowerTriangular(diagm(0.1))
R = vcat(diagm([1., 1e-2]), zeros(2,2))
M = LGSSM(A, B, R, S, 0.5*eye(4))
sim = Generate(M, 20*4, 5)
plot(sim[:Y][1,:,:])
KF = KalmanFilter(M, sim[:Y][1,:,1])

# Structural Time Series Model
# Local level linear model
# Seasonal component with s = 4
A = blkdiag(sparse([1,1,2],[1,2,2],[1.,1.,1.]), sparse([1,1,1,2,3], [1,2,3,1,2],[-1.,-1.,-1.,1.,1.]))
B = [1. 0. 1. 0. 0.]
S = LowerTriangular(diagm(0.1))
R = vcat(diagm([1e-2, 1e-1, 1.]), zeros(2,3))
M = LGSSM(A, B, R, S, 0.5*eye(5))
sim = Generate(M, 20*4, 50)
plot(sim[:Y][1,:,:])


text:
module LeapSeconds

using Dates: DateTime, datetime2julian

export offset_tai_utc, offset_utc_tai

include(joinpath("..", "gen", "leap_seconds.jl"))

const MJD_EPOCH = 2400000.5
const SECONDS_PER_DAY = 86400

# Constants for calculating the offset between TAI and UTC for
# dates between 1960-01-01 and 1972-01-01
# See ftp://maia.usno.navy.mil/ser7/tai-utc.dat

const EPOCHS = [
    36934,
    37300,
    37512,
    37665,
    38334,
    38395,
    38486,
    38639,
    38761,
    38820,
    38942,
    39004,
    39126,
    39887,
]

const OFFSETS = [
    1.417818,
    1.422818,
    1.372818,
    1.845858,
    1.945858,
    3.240130,
    3.340130,
    3.440130,
    3.540130,
    3.640130,
    3.740130,
    3.840130,
    4.313170,
    4.213170,
]

const DRIFT_EPOCHS = [
    37300,
    37300,
    37300,
    37665,
    37665,
    38761,
    38761,
    38761,
    38761,
    38761,
    38761,
    38761,
    39126,
    39126,
]

const DRIFT_RATES = [
    0.0012960,
    0.0012960,
    0.0012960,
    0.0011232,
    0.0011232,
    0.0012960,
    0.0012960,
    0.0012960,
    0.0012960,
    0.0012960,
    0.0012960,
    0.0012960,
    0.0025920,
    0.0025920,
]

leapseconds(mjd) = LEAP_SECONDS[searchsortedlast(LS_EPOCHS, floor(Int, mjd))]

"""
    offset_tai_utc(tai1, tai2=0.0)

Returns the difference between International Atomic Time (TAI) and Coordinated
Universal Time (UTC) for a given TAI Julian day number `tai1` (optionally
split into two parts for increased precision).

``\\Delta AT = TAI - UTC``
"""
function offset_tai_utc(tai1, tai2=0.0)
    mjd = tai1 - MJD_EPOCH + tai2

    # Before 1960-01-01
    if mjd < 36934.0
        @warn "UTC is not defined for dates before 1960-01-01."
        return 0.0
    end

    # Before 1972-01-01
    if mjd < LS_EPOCHS[1]
        idx = searchsortedlast(EPOCHS, floor(Int, mjd))
        rate_utc = DRIFT_RATES[idx] / SECONDS_PER_DAY
        rate_tai = rate_utc / (1 + rate_utc) * SECONDS_PER_DAY
        offset = OFFSETS[idx]
        Δt = mjd - DRIFT_EPOCHS[idx] - offset / SECONDS_PER_DAY
        return offset + Δt * rate_tai
    end

    leapseconds(mjd)
end

"""
    offset_tai_utc(dt::DateTime)

Returns the difference between International Atomic Time (TAI) and Coordinated
Universal Time (UTC) for a given `DateTime` in TAI.

``\\Delta AT = TAI - UTC``
"""
offset_tai_utc(dt::DateTime) = offset_tai_utc(datetime2julian(dt))

"""
    offset_utc_tai(utc1, utc2=0.0)

Returns the difference between Coordinated Universal Time (UTC) and
International Atomic Time (TAI) for a given UTC pseudo-Julian day number
`utc1` (optionally split into two parts for increased precision).

``\\Delta AT = UTC - TAI``

!!! note
    This function uses the [ERFA convention](https://github.com/liberfa/erfa/blob/master/src/dtf2d.c#L49)
    for Julian day numbers representing UTC dates during leap seconds.
"""
function offset_utc_tai(utc1, utc2=0.0)
    mjd = utc1 - MJD_EPOCH + utc2

    # Before 1960-01-01
    if mjd < 36934.0
        @warn "UTC is not defined for dates before 1960-01-01."
        return 0.0
    end

    # Before 1972-01-01
    if mjd < LS_EPOCHS[1]
        idx = searchsortedlast(EPOCHS, floor(Int, mjd))
        offset = OFFSETS[idx] + (mjd - DRIFT_EPOCHS[idx]) * DRIFT_RATES[idx]
        return -offset
    end

    offset = 0.0
    for _ = 1:3
        offset = leapseconds(mjd + offset / SECONDS_PER_DAY)
    end
    return -offset
end

"""
    offset_utc_tai(dt::DateTime)

Returns the difference between Coordinated Universal Time (UTC) and
International Atomic Time (TAI) for a given `DateTime` in UTC.

``\\Delta AT = UTC - TAI``

!!! warning
    The `DateTime` type from Julia's Standard Libary cannot represent UTC dates
    during leap seconds, e.g. "2016-12-31T23:59:60.0" will not be parsed as a
    valid `DateTime` but throw an error.
    The [AstroTime.jl](https://github.com/JuliaAstro/AstroTime.jl) package
    provides a leap second-aware `Epoch` type that can be used as a replacement.
"""
offset_utc_tai(dt::DateTime) = offset_utc_tai(datetime2julian(dt))

end


text:
@testset "198.house-robber.jl" begin
    @test rob_198([1, 2, 3, 1]) == 4
    @test rob_198([2, 7, 9, 3, 1]) == 12
    @test rob_198([
        40, 28, 8, 40, 32, 49, 31, 39, 11, 24, 3, 37, 10, 16, 23, 33, 3, 48, 41, 8
    ]) == 334
    @test rob_198([
        29, 38, 35, 20, 20, 41, 38, 10, 17, 14, 32, 14, 21, 31, 15, 13, 38, 19, 50, 43
    ]) == 295
end


text:
""" HighLevelCost

Abstract base type for the two kinds of high-level MAPF costs - Makespan and Sum-of-Costs
"""
abstract type HighLevelCost end

"""
    compute_cost(::HighLevelCost, solution::Vector{PR}) where {PR <: PlanResult}

Compute the cost for the MAPF solution.
"""
function compute_cost end

"""
    accumulate_cost(::HighLevelCost, old_cost::C, new_entry_cost::C) where {C <: Number}

Given a new entry, update the high-level cost according to the metric (addition for sum-of-costs,
maximum for makespan)
"""
function accumulate_cost end

"""
    deaccumulate_cost(::HighLevelCost, old_cost::C, new_entry_cost::C) where {C <: Number}

Given an entry to be removed, update the high-level cost according to the metric (subtraction for sum-of-costs,
minimum for makespan)
"""
function deaccumulate_cost end

struct SumOfCosts <: HighLevelCost
end

# Assumes that all entries are valid
function compute_cost(::SumOfCosts, solution::Vector{PR}) where {PR <: PlanResult}
    return sum(s.cost for s in solution)
end

function accumulate_cost(::SumOfCosts, old_cost::C, new_entry_cost::C) where {C <: Number}
    return (old_cost + new_entry_cost)
end

function deaccumulate_cost(::SumOfCosts, old_cost::C, new_entry_cost::C) where {C <: Number}
    return (old_cost - new_entry_cost)
end

struct Makespan <: HighLevelCost
end

# Assumes that all entries are valid
function compute_cost(::Makespan, solution::Vector{PR}) where {PR <: PlanResult}
    return maximum(s.cost for s in solution)
end

function accumulate_cost(::Makespan, old_cost::C, new_entry_cost::C) where {C <: Number}
    return max(old_cost, new_entry_cost)
end

function deaccumulate_cost(::Makespan, old_cost::C, new_entry_cost::C) where {C <: Number}
    return min(old_cost, new_entry_cost)
end


text:
for i in 3:5, j in 1:(i-2)
    show_simulation(i, j, 10^5)
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using InternedStrings
const N = 1_000_000; const K = 100
srand(1);
@time sample_space = [string(rand(Char.(97:97+25), rand(1:8))...) for k in 1:N÷K]
@time svec = rand(sample_space, N);

using BenchmarkTools
@time codeunit.(svec,1)

include("src/three_way_radix_quick_sort.jl")

x = svec[1:10]
lo = 1;
hi = length(x);
cmppos = 1;
@time three_way_radix_qsort0!(x)
issorted(x)

lo = 1;
hi = length(svec);
cmppos = 1;
@time three_way_radix_qsort0!(svec);
issorted(svec)

const N = 100_000_000; const K = 100
srand(1);
svec = rand([string(rand(Char.(97:97+25), rand(1:8))...) for k in 1:N÷K], N);

using FastGroupBy
@time radixsort!(svec)
issorted(svec)

using BenchmarkTools, FastGroupBy
const NN = 2^9
srand(1);
svec = rand([string(rand(Char.(97:97+25), rand(1:8))...) for k in 1:NN], NN);
x = @elapsed sort!(svec, alg = InsertionSort)
srand(1);
svec = rand([string(rand(Char.(97:97+25), rand(1:8))...) for k in 1:NN], NN);
y = @elapsed radixsort!(svec)

x < y


text:
using FunctionalCollections
using JSON

import FunctionalCollections: append
export Node, node, instanceof, props

const WEBIO_NODE_MIME = MIME"application/vnd.webio.node+json"
Base.Multimedia.istextmime(::WEBIO_NODE_MIME) = true

const WEBIO_APPLICATION_MIME = MIME"application/vnd.webio.application+html"
Base.Multimedia.istextmime(::WEBIO_APPLICATION_MIME) = true

struct Node{T}
    instanceof::T # if this changes the node must be *replaced*

    children::PersistentVector{Any}
    props::Dict{Symbol, Any}

    _descendants_count::Int

end

function Node(instanceof, children::AbstractVector, props::AbstractDict)
    inst = promote_instanceof(instanceof)
    Node{typeof(inst)}(inst, _pvec(children), kwargs2props(props), descendants_count(children))
end

function node(instanceof, children...; props...)
    Node(instanceof, PersistentVector{Any}(collect(children)), Dict(props))
end

promote_instanceof(x) = x
promote_instanceof(s::Symbol) = DOM(:html, s)
promote_instanceof(s::AbstractString) = promote_instanceof(Symbol(s))

nodetype(n::Node) = typename(n.instanceof)
typename(n::T) where {T} = string(T.name.name)

function kwargs2props(propkwargs)
    props = Dict{Symbol,Any}(propkwargs)
    props # XXX IJulia/JSON bug? kernel seems to crash if this is a String not a Dict (which is obviously silly but still, it shouldn't crash the IJulia kernel)
end

struct DOM
    namespace::Symbol
    tag::Symbol
end

export children, setchildren, instanceof, setinstanceof, props, setprops

children(n::Node) = n.children
setchildren(n::Node, children) = Node(n.instanceof, PersistentVector{Any}(collect(children)), n.props)
instanceof(n::Node) = n.instanceof
setinstanceof(n::Node, instanceof) = Node(instanceof, n.children, n.props)
props(n::Node) = n.props
setprops(n::Node, props) = Node(n.instanceof, n.children, props)

######## modifying an element #######

export append, setchild, withchild, withlastchild, mergeprops

append(n::Node, cs) = setchildren(n, append(children(n), cs))
setchild(n::Node, i, c) = setchildren(n, assoc(children(n), i, c))
withchild(f, n::Node, i) = setchild(n, i, f(children(n)[i]))
withlastchild(f, n::Node) = setchild(n, length(children(n)), f(children(n)[i]))
function mergeprops(n::Node, p, ps...)
    out = recmerge(props(n), p)
    for x in ps
        recmerge!(out, x)
    end
    setprops(n, out)
end

####### Rendering to HTML ########

function JSON.lower(n::Node)
    result = Dict{String, Any}(
        "type" => "node",
        "nodeType" => nodetype(n),
        "instanceArgs" => JSON.lower(n.instanceof),
        "children" => map!(
            render,
            Vector{Any}(undef, length(children(n))),
            children(n),
        ),
        "props" => props(n),
    )
    return result
end

"""
Escape characters for a "safe" representation of JSON.

In particular, we escape '/' characters to avoid the presence of "</" (and
especially "</script>") which cause the browser to break out of the current

        """
    )
end

function Base.show(io::IO, m::WEBIO_NODE_MIME, node::Node)
    write(io, JSON.json(node))
end

Base.show(io::IO, m::MIME"text/html", x::Observable) = show(io, m, WebIO.render(x))
Base.show(io::IO, m::WEBIO_NODE_MIME, x::Union{Observable, AbstractWidget}) = show(io, m, WebIO.render(x))

### Utility

descendants_count(t::String) = 0
descendants_count(el::Node) = el._descendants_count
function descendants_count(v::AbstractArray)
    s = 0
    for i in 1:length(v)
        @inbounds s +=_count(v[i])
    end
    s
end
_count(t) = 1
_count(el::Node) = el._descendants_count + 1

## Element extension syntax

(n::Node)(x, args...) = append(n, (x, args...))
(n::Node)(props::AbstractDict...) = mergeprops(n, props...)
(n::Node)(;kwargs...) = mergeprops(n, kwargs)
(n::Node)(args...; kwargs...) = n(args...)(;kwargs...)

## Pretty printing

function showindent(io, level)
    for i=1:level
        write(io, "  ")
    end
end

function Base.show(io::IO, ::MIME"text/plain", el::Node)
    _show(io, el)
end

function showprops(io, dict)
    write(io, "{")
    write(io, ' ')
    for (k,v) in dict
        print(io, k)
        write(io, '=')
        show(io, v)
        write(io, ' ')
    end
    write(io, "}")
end

function showchildren(io, elems, indent_level)
    length(elems) == 0 && return
    write(io, "\n")
    l = length(elems)
    for i=1:l
        _show(io, elems[i], indent_level+1)
        i != l && write(io, "\n")
    end
end

function _show(io::IO, el::Node, indent_level=0)
    showindent(io, indent_level)
    write(io, "(")
    if !isa(el.instanceof, DOM)
        write(io, string(el.instanceof))
        write(io, ":")
    elseif el.instanceof.namespace != :html
        write(io, el.instanceof.namespace)
        write(io, ":")
        write(io, el.instanceof.tag)
    else
        write(io, el.instanceof.tag)
    end

    if !isempty(props(el))
        write(io, " ")
        showprops(io, props(el))
    end
    showchildren(io, children(el), indent_level)
    write(io, ")")
end

function _show(io::IO, el, indent_level)
    showindent(io, indent_level)
    show(io, el)
end


text:
# elimat(m) returns the elimination matrix Lm
# The elimination matrix Lm is for any matrix F, Vech(F)=Lm Vec(F)
# Luca Brugnolini 2014
# Based on Yun Jung Kim

function elimat(m)

A = eye(m^2);
L = A[1:m,:];

for n = 2:m
    S = A[m*(n-1)+1:n*m,:];
    S = S[n:end,:];
    L = [L;S];

end

    return L

end


text:
# question (a)
say("== (a)")
tic()
test_a()
toc()

# question (b)
say("== (b)")
m = 100000

tic()
D = [test_a(verbose=false) for i=1:m]
toc()
say("$(sum(D.>XI100)/m) percent outside $XI100.")

figure()
(range, counts) = hist(D, 0:200)
p = bar(0:199, counts)
q = plot(50:150, chi_square(50:150, 100)*m);

# question (c)
say("== (c)")
tic()
test_c()
toc()

# question (d)
say("== (d)")
m = 100000

tic()
D = [test_c(verbose=false) for i=1:m]
toc()
say("$(sum(D.>XI80)/m) percent outside $XI80.")

figure()
(range, counts) = hist(D, 0:200)
p = bar(0:199, counts)
q = plot(50:150, chi_square(50:150, 80)*m);

max_stars_repo_name: innerlee/testRand

max_issues_repo_name: innerlee/testRand

max_forks_repo_name: innerlee/testRand


text:
# usage example for package developer:
#
#     import EdgeDetectionAPI: AbstractEdgeDetectionAlgorithm,
#                             detect_edges, detect_edges!,
#                             detect_subpixel_edges, detect_subpixel_edges!

"""
    AbstractEdgeDetectionAlgorithm <: AbstractImageFilter

The root type for `ImageEdgeDetection` package.

Any concrete edge detection algorithm shall subtype it to support
[`detect_edges`](@ref), [`detect_edges!`](@ref), [`detect_subpixel_edges`](@ref)
and [`detect_subpixel_edges!`](@ref) APIs.

# Examples

All edge detection algorithms in ImageEdgeDetection are called in the
following pattern:

```julia
# first generate an algorithm instance
f = Canny()

# then pass the algorithm to `detect_edges`
img_edges = detect_edges(img, f)

# or use in-place version `detect_edges!`
img_edges = similar(img)
detect_edges!(img_edges, img, f)
```

For more examples, please check [`detect_edges`](@ref),
[`detect_edges!`](@ref) and concrete algorithms.

One can also detect edges to subpixel accuracy by specifying
`SubpixelNonmaximaSuppression` as the edge thinning algorithm and using
[`detect_subpixel_edges`](@ref) or [`detect_subpixel_edges!`](@ref). The
function returns an edge image as well as a accompanying matrix of length-2
vectors which, when added to the edge image coordinates, specify the location
of an edge to subpixel precision.

```julia
# first generate an algorithm instance
f = Canny(thinning_algorithm = SubpixelNonmaximaSuppression())

# then pass the algorithm to `detect_subpixel_edges`
img_edges, subpixel_offsets = detect_subpixel_edges(img, f)

# or use in-place version `detect_edges!`
img_edges = similar(img)
subpixel_offsets = zeros(SVector{2,Float64}, axes(img))
detect_edges!(img_edges, subpixel_offsets, img, f)
```

"""
abstract type AbstractEdgeDetectionAlgorithm <: AbstractImageFilter end

detect_edges!(out::Union{GenericGrayImage, AbstractArray{<:Color3}},
              img,
              f::AbstractEdgeDetectionAlgorithm,
              args...; kwargs...) =
        f(out, img, args...; kwargs...)

# TODO: Relax this to all color types
function detect_edges!(img::Union{GenericGrayImage, AbstractArray{<:Color3}},
                       f::AbstractEdgeDetectionAlgorithm,
                       args...; kwargs...)
    tmp = copy(img)
    f(img, tmp, args...; kwargs...)
    return img
end

function detect_edges(::Type{T},
                      img,
                      f::AbstractEdgeDetectionAlgorithm,
                      args...; kwargs...) where T
    out = similar(Array{T}, axes(img))
    detect_edges!(out, img, f, args...; kwargs...)
    return out
end

detect_edges(img::AbstractArray{T},
             f::AbstractEdgeDetectionAlgorithm,
             args...; kwargs...) where T <: Colorant =
        detect_edges(T, img, f, args...; kwargs...)

# Do not promote Number to Gray{<:Number}
detect_edges(img::AbstractArray{T},
             f::AbstractEdgeDetectionAlgorithm,
             args...; kwargs...) where T <: Number =
        detect_edges(T, img, f, args...; kwargs...)

# TODO Take care of this in a separate pull-request.
# Handle instance where the input is a sequence of images.
# detect_edges!(out_sequence::Vector{T},
#           img_sequence,
#           f::AbstractEdgeDetectionAlgorithm,
#           args...; kwargs...) where T <: Union{GenericGrayImage, AbstractArray{<:Color3}} =
#     f(out_sequence, img_sequence, args...; kwargs...)
#
# # TODO: Relax this to all color types
# function detect_edges!(img_sequence::Vector{T},
#                    f::AbstractEdgeDetectionAlgorithm,
#                    args...; kwargs...) where T <: Union{GenericGrayImage, AbstractArray{<:Color3}}
#     tmp = copy(img_sequence)
#     f(img_sequence, tmp, args...; kwargs...)
#     return img_sequence
# end
#
# function detect_edges(::Type{T},
#                   img_sequence::Vector{<:AbstractArray},
#                   f::AbstractEdgeDetectionAlgorithm,
#                   args...; kwargs...) where T
#     N  = length(img_sequence)
#     out_sequence = [similar(Array{T}, axes(img_sequence[n])) for n = 1:N]
#     detect_edges!(out_sequence, img_sequence, f, args...; kwargs...)
#     return out_sequence
# end

detect_edges(img_sequence::Vector{<:AbstractArray{T}},
             f::AbstractEdgeDetectionAlgorithm,
             args...; kwargs...) where T <: Colorant =
         detect_edges(T, img_sequence, f, args...; kwargs...)

# Do not promote Number to Gray{<:Number}
detect_edges(img_sequence::Vector{<:AbstractArray{T}},
             f::AbstractEdgeDetectionAlgorithm,
             args...; kwargs...) where T <: Number =
        detect_edges(T, img_sequence, f, args...; kwargs...)

######
detect_subpixel_edges!(out₁::Union{GenericGrayImage, AbstractArray{<:Color3}},
                       out₂::AbstractArray{<: StaticArray},
                       img,
                       f::AbstractEdgeDetectionAlgorithm,
                       args...; kwargs...) =
    f(out₁, out₂, img, args...; kwargs...)

# TODO: Relax this to all color types
function detect_subpixel_edges!(img::Union{GenericGrayImage, AbstractArray{<:Color3}},
                                f::AbstractEdgeDetectionAlgorithm,
                                args...; kwargs...)
    subpixel_offsets = zeros(SVector{2,Float64}, axes(img))
    tmp = copy(img)
    f(img, subpixel_offsets, tmp, args...; kwargs...)
    return img, subpixel_offsets
end

function detect_subpixel_edges(::Type{T},
                               img,
                               f::AbstractEdgeDetectionAlgorithm,
                               args...; kwargs...) where T
    out = similar(Array{T}, axes(img))
    subpixel_offsets = zeros(SVector{2,Float64}, axes(img))
    detect_subpixel_edges!(out, subpixel_offsets, img, f, args...; kwargs...)
    return out, subpixel_offsets
end

function detect_subpixel_edges(::Type{T₁}, ::Type{T₂},
                               img,
                               f::AbstractEdgeDetectionAlgorithm,
                               args...; kwargs...) where {T₁,T₂}
    out = similar(Array{T₁}, axes(img))
    subpixel_offsets = zeros(SVector{2,T₂}, axes(img))
    detect_subpixel_edges!(out, subpixel_offsets, img, f, args...; kwargs...)
    return out, subpixel_offsets
end

detect_subpixel_edges(img::AbstractArray{T},
                      f::AbstractEdgeDetectionAlgorithm,
                      args...; kwargs...) where T <: Colorant =
         detect_subpixel_edges(T, img, f, args...; kwargs...)

# Do not promote Number to Gray{<:Number}
detect_subpixel_edges(img::AbstractArray{T},
                 f::AbstractEdgeDetectionAlgorithm,
                 args...; kwargs...) where T <: Number =
        detect_subpixel_edges(T, img, f, args...; kwargs...)

# TODO Handle sequences properly in a separate pull-request.

# # Handle instance where the input is a sequence of images.
# detect_subpixel_edges!(out_sequence::Vector{T},
#           img_sequence,
#           f::AbstractEdgeDetectionAlgorithm,
#           args...; kwargs...) where T <: Union{GenericGrayImage, AbstractArray{<:Color3}} =
#     f(out_sequence, img_sequence, args...; kwargs...)
#
# # TODO: Relax this to all color types
# function detect_subpixel_edges!(img_sequence::Vector{T},
#                    f::AbstractEdgeDetectionAlgorithm,
#                    args...; kwargs...) where T <: Union{GenericGrayImage, AbstractArray{<:Color3}}
#     tmp = copy(img_sequence)
#     f(img_sequence, tmp, args...; kwargs...)
#     return img_sequence
# end

# function detect_subpixel_edges(::Type{T},
#                   img_sequence::Vector{<:AbstractArray},
#                   f::AbstractEdgeDetectionAlgorithm,
#                   args...; kwargs...) where T
#     N  = length(img_sequence)
#     out_sequence = [similar(Array{T}, axes(img_sequence[n])) for n = 1:N]
#     detect_subpixel_edges!(out_sequence, img_sequence, f, args...; kwargs...)
#     return out_sequence
# end
#
# detect_subpixel_edges(img_sequence::Vector{<:AbstractArray{T}},
#                  f::AbstractEdgeDetectionAlgorithm,
#                  args...; kwargs...) where T <: Colorant =
#          detect_subpixel_edges(T, img_sequence, f, args...; kwargs...)
#
# # Do not promote Number to Gray{<:Number}
# detect_subpixel_edges(img_sequence::Vector{<:AbstractArray{T}},
#                  f::AbstractEdgeDetectionAlgorithm,
#                  args...; kwargs...) where T <: Number =
#         detect_subpixel_edges(T, img_sequence, f, args...; kwargs...)

### Docstrings

"""
    detect_edges!([out,] img, f::AbstractEdgeDetectionAlgorithm, args...; kwargs...)

Detect edges of `img` using algorithm `f`;  if left unspecified, `f` is assumed
to be [`Canny`](@ref ImageEdgeDetection.Canny).

# Output

If `out` is specified, it will be changed in place. Otherwise `img` will be changed in place.

# Examples

Just simply pass an algorithm to `detect_edges!`:

```julia
img_edges = similar(img)
detect_edges!(img_edges, img, f)
```

For cases you just want to change `img` in place, you don't necessarily need to manually
allocate `img_edges`; just use the convenient method:

```julia
detect_edges!(img, f)
```

See also: [`detect_edges`](@ref)
"""
detect_edges!

"""
    out = detect_edges([T::Type,] img, f::AbstractEdgeDetectionAlgorithm, args...; kwargs...)

Detect edges of `img` using algorithm `f`;  if left unspecified, `f` is assumed
to be [`Canny`](@ref ImageEdgeDetection.Canny).

# Output

The return image `out` is an `Array{T}`. If `T` is not specified, then
it's inferred.

# Examples

Just simply pass the input image and algorithm to `detect_edges`

```julia
f = Canny(spatial_scale = 1.4)
img_edges = detect_edges(img, f)
```

This reads as "`detect_edges` of image `img` using algorithm `f`".

You can also explicitly specify the return type:

```julia
f = Canny(spatial_scale = 1.4)
img_edges_float32 = detect_edges(Gray{Float32}, img, f)
```

See also [`detect_edges!`](@ref) for in-place edge detection.
"""
detect_edges

"""
    detect_subpixel_edges!(out₁, out₂, img, f::AbstractEdgeDetectionAlgorithm, args...; kwargs...)

Detect edges of `img` to subpixel precision using algorithm `f`;  if left
unspecified, `f` is assumed to be [`Canny`](@ref ImageEdgeDetection.Canny).

# Output

The integer components of an edge correspond to non-zero row and column
entries in `out₁`. The accompanying subpixel offsets  are stored in a 2-D array
`out₂` as length-2 vectors. One can recover the  subpixel coordinates by adding
the subpixel offsets to the integer components.

# Examples

Just simply pass an algorithm to `detect_subpixel_edges!`:

```julia
f = Canny(spatial_scale = 1.4, thinning_algorithm = SubpixelNonmaximaSuppression())
img_edges = similar(img)
offsets = zeros(SVector{2,Float64}, axes(img))
detect_edges!(img_edges, offsets, img, f)
```

See also: [`detect_subpixel_edges`](@ref)
"""
detect_subpixel_edges!

"""
    out₁, out₂ = detect_subpixel_edges([T₁::Type, T₂::Type], img, f::AbstractEdgeDetectionAlgorithm, args...; kwargs...)

Detect edges of `img` to subpixel precision using algorithm `f`;  if left
unspecified, `f` is assumed to be [`Canny`](@ref ImageEdgeDetection.Canny).

# Output

The integer components of an edge correspond to non-zero row and column entries
in `out₁` which is an `Array{T₁}`. The accompanying subpixel offsets  are stored
in a 2-D array `out₂` as length-2 vectors (`Array{SVector{2, T₂}}`). One can
recover the  subpixel coordinates by adding the subpixel offsets to the integer
components.

# Examples

Just simply pass the input image and algorithm to `detect_subpixel_edges`

```julia
f = Canny(spatial_scale = 1.4, thinning_algorithm = SubpixelNonmaximaSuppression())
img_edges, offsets = detect_subpixel_edges(img, f)
```

This reads as "`detect_subpixel_edges` of image `img` using algorithm `f`".

You can also explicitly specify the return types:

```julia
f = Canny(spatial_scale = 1.4, thinning_algorithm = SubpixelNonmaximaSuppression())
img_edges, offsets = detect_subpixel_edges(Gray{Float32}, Float32, img, f)
```

See also [`detect_subpixel_edges!`](@ref) for in-place edge detection.
"""
detect_subpixel_edges


text:
include("PositionTypes.jl")
include("GreekAnalytical.jl")
include("SpreadBuilder.jl")
include("BSGreeksInverseOption.jl")

max_stars_repo_name: Stochastic-Adventure/BlackScholesGreeks

max_issues_repo_name: Stochastic-Adventure/BlackScholesGreeks

max_forks_repo_name: Stochastic-Adventure/BlackScholesGreeks


text:
module LinCombBenchs

using InplaceArrays.Arrays
using InplaceArrays.Fields
using InplaceArrays.Fields: MockField, MockBasis
using FillArrays
using InplaceArrays.TensorValues

@noinline function loop(a,cache)
  for i in eachindex(a)
    ai = getindex!(cache,a,i)
  end
end

@noinline function loop_and_evaluate(ca,cai,cx,a,x)
  for i in eachindex(a)
    ai = getindex!(ca,a,i)
    xi = getindex!(cx,x,i)
    vi = evaluate!(cai,ai,xi)
  end
end

function bench1(n)
  np = 4
  p = Point(1,2)
  x = fill(p,np)
  v = 3.0
  d = 2
  ndof = 8
  wi = 3.0
  w = fill(wi,ndof)
  l = n
  f = MockBasis{d}(v,ndof)
  af = Fill(f,l)
  ax = fill(x,l)
  aw = fill(w,l)
  az = fill(fill(1.0,np),l)
  afx = evaluate(af,ax)
  cafx = array_cache(afx)
  @time loop(afx,cafx)
  ag = lincomb(af,aw)
  cag = array_cache(ag)
  @time loop(ag,cag)
  agx = evaluate(ag,ax)
  cagx = array_cache(agx)
  @time loop(agx,cagx)
  g = testitem(ag)
  cg = field_cache(g,x)
  cax = array_cache(ax)
  @time loop_and_evaluate(cag,cg,cax,ag,ax)
end

function bench2(n)
  np = 4
  p = Point(1,2)
  x = fill(p,np)
  v = 3.0
  d = 2
  ndof = 8
  wi = 3.0
  w = fill(wi,ndof)
  l = n
  f = MockBasis{d}(v,ndof)
  af = Fill(f,l)
  ax = fill(x,l)
  aw = fill(w,l)
  az = fill(fill(1.0,np),l)
  _ag = lincomb(af,aw)
  ag = gradient(_ag)
  cag = array_cache(ag)
  @time loop(ag,cag)
  agx = evaluate(ag,ax)
  cagx = array_cache(agx)
  @time loop(agx,cagx)
  g = testitem(ag)
  cg = field_cache(g,x)
  cax = array_cache(ax)
  @time loop_and_evaluate(cag,cg,cax,ag,ax)
end

for n in (1,1,10,1000,100000)
  @eval begin
    println("+++ runing suite for n = $($n) +++")
    bench1($n)
    bench2($n)
  end
end

end # module


text:
function make_landmarks(nlandmarks::Int, boundaries, margin)
    xmin, xmax, ymin, ymax = boundaries
    bx = margin*(xmax - xmin)
    by = margin*(ymax - ymin)
    rand([xmin + bx:xmax - bx; ymin + by:ymax - by], 2, nlandmarks)
end

"""
Return indices of landmarks within sensor acceptance (r < rmax, phi < +/)
"""
function nearby_landmark_indices(pose, lm, rmax)
    dx = lm[1,:] - pose[1]
    dy = lm[2,:] - pose[2]
    phi = pose[3]

    inearby = Array{Int}(0)
    for i in 1:size(lm, 2)
        # Select landmarks within +/- 90 degrees of current heading
        if (dx[i]*cos(phi) + dy[i]*sin(phi)) > 0
            # Select landmarks inside sensor range
            if ((dx[i]^2 + dy[i]^2) < rmax^2)
                push!(inearby, i)
            end
        end
    end
    return inearby
end

"""
Add Gaussian noise to nominal control values
Noise is taken to be uncorrelated (assume diagonal control_cov)
"""
function add_control_noise!(vehicle::Vehicle, control_cov::AbstractMatrix)
    vehicle.measured_speed = vehicle.target_speed + randn(1)[1]*sqrt(control_cov[1,1])
    vehicle.measured_gamma = vehicle.target_gamma + randn(1)[1]*sqrt(control_cov[2,2])
end

function add_control_noise(speed, gamma, control_cov)
    # Add Gaussian noise to nominal control values
    # Noise is taken to be uncorrelated (assume diagonal control_cov)
    speed + randn(1)[1]*sqrt(control_cov[1,1]),
    gamma + randn(1)[1]*sqrt(control_cov[2,2])
end

"""
Return array with columns z = [range; angle] for all nearby landmark
observations as well as the corresponding list of landmark IDs.
"""
function get_observations(vehicle::Vehicle, scene::Scene, R)
    lm = scene.landmarks
    x,y,phi = vehicle.pose

    inearby = nearby_landmark_indices(vehicle.pose, lm, vehicle.sensor_range)

    dx  = lm[1, inearby] - x
    dy  = lm[2, inearby] - y

    # 2 by nz observation matrix - columns are range and bearing
    z = [sqrt(dx.^2 + dy.^2) atan2(dy, dx) - phi]'

    # Add Gaussian noise to z to simulate measurement resolution
    ncols = size(z, 2)
    if ncols > 0
        z = z + [randn(1, ncols)*sqrt(R[1,1]); randn(1, ncols)*sqrt(R[2,2])]
    end

    # Unique identifier for each landmark
    landmark_tags = 1:size(lm, 2)

    return z, landmark_tags[inearby]
end


text:
export opendb
export load_tanks, add_tank!, findalltanks, findfirsttank, tanknames
export load_engines

using SQLite
using DataFrames
using Rockets

"""
    opendb() -> SQLite.DB
Open SQLite database. Make sure you close it with `DBInterface.close!(db)`
"""
function opendb()
    path = joinpath(datadir, "rockets.db")
    return SQLite.DB(path)
end

# path = joinpath(datadir, "rockets.db")
# db = SQLite.DB(path)
# DataFrame(DBInterface.execute(db, "SELECT * FROM Tanks"))

"""
    load_tanks() -> DataFrame
Load tank data as a data frame from data/rockets.db SQLite database
"""
function load_engines()
    path = joinpath(datadir, "rockets.db")
    db = SQLite.DB(path)

    query = DBInterface.prepare(db, "SELECT * FROM Engines")
    df = DBInterface.execute(query) |> DataFrame

    DBInterface.close!(db)  # close database for further interaction
    return df
end

"""
    load_tanks() -> DataFrame
Load tank data as a data frame from data/rockets.db SQLite database
"""
function load_tanks()
    path = joinpath(datadir, "rockets.db")
    db = SQLite.DB(path)

    query = DBInterface.prepare(db, "SELECT * FROM Tanks")
    df = DBInterface.execute(query) |> DataFrame

    DBInterface.close!(db)  # close database for further interaction
    return df
end

"""
    add_tank(db::Connection, name, tank::Tank)
Adds a tank with given `name` to the `db` database.
"""
function add_tank!(db::DBInterface.Connection, name::AbstractString, tank::Tank)
    query = DBInterface.prepare(db, "INSERT INTO Tanks VALUES(?, ?, ?)")
    DBInterface.execute(query, [name, tank.total_mass, tank.dry_mass])
end

"Find all tanks with given tank name. Should only be one"
function findalltanks(db::DBInterface.Connection, tankname::AbstractString)
    query  = DBInterface.prepare(db, "SELECT * FROM Tanks WHERE Tank LIKE ?")
    result = DBInterface.execute(query, [tankname])
    rows   = Tables.rows(result)        # Get row iterator
    tanks  = map(rows) do row           # Iterate over each row in table
        Tank(row.DryMass, row.TotalMass)
    end

    tanks
end

"""
    findfirsttank(db, name)
Finds first tank with `name`. If no such tank is found return `nothing`
"""
function findfirsttank(db::DBInterface.Connection, tankname::AbstractString)
    tanks = findalltanks(db, tankname)
    if isempty(tanks)
        nothing
    else
        tanks[1]
    end
end

function tanknames(db::DBInterface.Connection)
    result = DBInterface.execute(db, "SELECT Tank FROM Tanks")
    rows   = Tables.rows(result)        # Get row iterator
    map(rows) do row
        row.Tank
    end
end


text:
# ## Preliminary steps
#
# Let's start by loading the relevant packages and generating some dummy data.

using MLJ
import DataFrames
import Statistics
using PrettyPrinting
using StableRNGs

MLJ.color_off() # hide

rng = StableRNG(512)
Xraw = rand(rng, 300, 3)
y = exp.(Xraw[:,1] - Xraw[:,2] - 2Xraw[:,3] + 0.1*rand(rng, 300))
X = DataFrames.DataFrame(Xraw)

train, test = partition(eachindex(y), 0.7);

# Let's also load a simple model:

@load KNNRegressor
knn_model = KNNRegressor(K=10)

# As before, let's instantiate a machine that wraps the model and data:

knn = machine(knn_model, X, y)

# and fit it

fit!(knn, rows=train)
ŷ = predict(knn, X[test, :]) # or use rows=test
rms(ŷ, y[test])

# The few steps above are equivalent to just calling `evaluate!`:

evaluate!(knn, resampling=Holdout(fraction_train=0.7, rng=StableRNG(666)),
          measure=rms) |> pprint

# ## Homogenous ensembles
#
# MLJ offers basic support for ensembling such as [_bagging_](https://en.wikipedia.org/wiki/Bootstrap_aggregating).
# Defining such an ensemble of simple "atomic" models is done via the `EnsembleModel` constructor:

ensemble_model = EnsembleModel(atom=knn_model, n=20);

# where the `n=20` indicates how many models are present in the ensemble.
#
# ### Training and testing an ensemble
#
# Now that we've instantiated an ensemble, it can be trained and tested the same as any other model:

ensemble = machine(ensemble_model, X, y)
estimates = evaluate!(ensemble, resampling=CV())
estimates |> pprint

# here the implicit measure is the `rms` (default for regressions). The `measurement` is the mean taken over the folds:

@show estimates.measurement[1]
@show mean(estimates.per_fold[1])

# Note that multiple measurements can be specified jointly. Here only on measurement is (implicitly) specified but we still have to select the corresponding results (whence the `[1]` for both  the `measurement` and `per_fold`).
#
# ### Systematic tuning
#
# Let's simultaneously tune the ensemble's `bagging_fraction` and the K-Nearest neighbour hyperparameter `K`. Since one of our models is  a field of the  other, we have nested hyperparameters:

params(ensemble_model) |> pprint

# To define a tuning grid, we construct ranges for the two parameters and collate these ranges:

B_range = range(ensemble_model, :bagging_fraction,
                lower=0.5, upper=1.0)
K_range = range(ensemble_model, :(atom.K),
                lower=1, upper=20);

# the scale for a tuning grid is linear by default but can be specified to `:log10` for logarithmic ranges.
# Now we have to define a `TunedModel` and fit it:

tm = TunedModel(model=ensemble_model,
                tuning=Grid(resolution=10), # 10x10 grid
                resampling=Holdout(fraction_train=0.8, rng=StableRNG(42)),
                ranges=[B_range, K_range])

tuned_ensemble = machine(tm, X, y)
fit!(tuned_ensemble, rows=train);

# Note the `rng=42` seeds the random number generator for reproducibility of this example.
#
# ### Reporting results
#
# The best model can be accessed like so:

best_ensemble = fitted_params(tuned_ensemble).best_model
@show best_ensemble.atom.K
@show best_ensemble.bagging_fraction

# The `report` method gives more detailed information on the tuning process:

r = report(tuned_ensemble);

# For instance, `r.measurements` are the measurements for all pairs of hyperparameters which you could visualise nicely:

using PyPlot

figure(figsize=(8,6))

res = r.plotting
vals_b = res.parameter_values[:, 1]
vals_k = res.parameter_values[:, 2]

tricontourf(vals_b, vals_k, res.measurements)
xticks(0.5:0.1:1, fontsize=12)
xlabel("Bagging fraction", fontsize=14)
yticks([1, 5, 10, 15, 20], fontsize=12)
ylabel("Number of neighbors - K", fontsize=14)

savefig(joinpath(@OUTPUT, "A-ensembles-heatmap.svg")) # hide

# \figalt{Hyperparameter heatmap}{A-ensembles-heatmap.svg}
#
# Finally you can always just evaluate the model by reporting `rms` on the test set:

ŷ = predict(tuned_ensemble, rows=test)
rms(ŷ, y[test])

PyPlot.close_figs() # hide


text:
# Show how to do joins. First joins on two distinct columns, then joins on a range of columns.
using PyPlot

# Set max elements in iterator
MaxElem = 1000

# Pick two columns to join
col1 = "1,"
col2 = "100,"

# Find all rows with these columns
Ajoin = Tadj[getrow(sum(logical(Tadj[:,col1*col2]),2) == 2),:]

# Display
figure()
spy(Ajoin)
axis("auto")
xlabel("end vertex")
ylabel("start vertex")

# Set column ranges
colRange1 = StartsWith("111,")
colRange2 = StartsWith("222,")

# Set up query iterators
TadjIt1 = getiterator(Tadj,MaxElem)
TadjIt2 = getiterator(Tadj,MaxElem)

# Start first query iterator
A1 = logical(TadjIt1[:,colRange1])
A1outDeg = Assoc("","","");

while nnz(A1) > 0
    global A1, A1outDeg
    # Combine
    A1outDeg = A1outDeg + sum(A1,2)
    # Run next query iterator
    A1 = logical(TadjIt1[])
end

# Start second query iterator
A2 = logical(TadjIt2[:,colRange2])
A2outDeg = Assoc("","","")

while nnz(A2) > 0
    global A2, A2outDeg
    # Combine
    A2outDeg = A2outDeg + sum(A2,2)
    # Run next query iterator
    A2 = logical(TadjIt2[])
end

# Join columns
AjoinRange = Tadj[getrow(A1outDeg[getrow(A2outDeg),:]),:]


text:
let
    roadway = get_test_roadway()
    trajdata = get_test_trajdata(roadway)
    veh = get(trajdata, 1, 1)

    let
        s = VehicleState()
        @test VehicleState() == propagate(veh, s, roadway, NaN)
    end

    let
        a = AccelTurnrate(0.1,0.2)
        @test a == convert(AccelTurnrate, [0.1,0.2])
        @test copy!([NaN, NaN], AccelTurnrate(0.1,0.2)) == [0.1,0.2]

        s = propagate(veh, AccelTurnrate(0.0,0.0), roadway, 1.0)
        @test isapprox(s.posG.x, veh.state.v)
        @test isapprox(s.posG.y, 0.0)
        @test isapprox(s.posG.θ, 0.0)
    end

    let
        a = AccelDesang(0.1,0.2)
        @test a == convert(AccelDesang, [0.1,0.2])
        @test copy!([NaN, NaN], AccelDesang(0.1,0.2)) == [0.1,0.2]

        s = propagate(veh, AccelDesang(0.0,0.0), roadway, 1.0)
        @test isapprox(s.posG.x, veh.state.v*1.0)
        @test isapprox(s.posG.y, 0.0)
        @test isapprox(s.posG.θ, 0.0)
        @test isapprox(s.v, veh.state.v)
    end

    let
        a = LatLonAccel(0.1,0.2)
        @test a == convert(LatLonAccel, [0.1,0.2])
        @test copy!([NaN, NaN], LatLonAccel(0.1,0.2)) == [0.1,0.2]

        Δt = 0.1
        s = propagate(veh, LatLonAccel(0.0,0.0), roadway, Δt)
        @test isapprox(s.posG.x, veh.state.v*Δt)
        @test isapprox(s.posG.y, 0.0)
        @test isapprox(s.posG.θ, 0.0)
        @test isapprox(s.v, veh.state.v)
    end
end


text:
for N in [Float64, Float32, Rational{Int}]
    # unit vector
    for n in 1:3
        for i in 1:n
            e = SingleEntryVector(i, n, one(N))
            for j in 1:n
                @test e[j] == (i == j ? one(N) : zero(N))
            end
        end
    end

    # check that #1359 is fixed
    e = SingleEntryVector(2, 3, N(1)) # vector [0, 1, 0]
    D = Diagonal(N[1, 2, 3])
    De = D * e
    @test De isa SingleEntryVector && Vector(De) == N[0, 2, 0]

    # template direction approximation
    for n in 1:3
        B = BallInf(zeros(N, n), N(2))
        A = Matrix{N}(2I, n, n) + ones(N, n, n)
        X = A * B

        # box directions
        dir = BoxDirections{N}(n)
        @test isbounding(dir)
        @test dim(dir) == n
        box = overapproximate(X, dir)
        @test length(dir) == length(box.constraints) == 2*n
        box = overapproximate(X, BoxDirections)
        @test length(dir) == length(box.constraints) == 2*n

        # octagon directions
        dir = OctDirections{N}(n)
        @test isbounding(dir)
        @test dim(dir) == n
        oct = overapproximate(X, dir)
        @test length(dir) == length(oct.constraints) == 2 * n^2
        oct = overapproximate(X, OctDirections)
        @test length(dir) == length(oct.constraints) == 2 * n^2

        # box-diagonal directions
        dir = BoxDiagDirections{N}(n)
        @test isbounding(dir)
        @test dim(dir) == n
        boxdiag = overapproximate(X, dir)
        @test length(dir) == length(boxdiag.constraints) ==
              (n == 1 ? 2 : 2^n + 2*n)
        boxdiag = overapproximate(X, BoxDiagDirections)
        @test length(dir) == length(boxdiag.constraints) ==
              (n == 1 ? 2 : 2^n + 2*n)

        # spherical directions approximation
        if n == 2 && N in [Float32, Float64]
            dir = PolarDirections{N}(2)
            @test !isbounding(dir)
            dir = PolarDirections{N}(5)
            @test isbounding(dir)
            @test dim(dir) == 2
            polar = overapproximate(X, dir)
        end

        # spherical directions approximation
        if n == 3 && N in [Float32, Float64]
            dir = SphericalDirections{N}(2, 2)
            @test !isbounding(dir)
            dir = SphericalDirections{N}(5, 5)
            @test isbounding(dir)
            @test dim(dir) == 3
            spherical = overapproximate(X, dir)
        end

        # overapproximate lazy polyhedral intersections
        if N in [Float64]
            Y = B ∩ Ball1(zeros(N, n), N(1))
            Z = B ∩ HalfSpace(ones(N, n), N(1))
            for dir in [BoxDirections, OctDirections, BoxDiagDirections]
                overapproximate(Y, dir)
                overapproximate(Z, dir)
            end
        end

        # custom directions
        # empty list of directions
        dir = CustomDirections(Vector{N}[]; n=n)
        @test dim(dir) == n
        P = overapproximate(X, dir)
        @test isempty(constraints_list(P))
        # minimal number of bounded directions
        if n == 1
            dirs = [N[-1], N[1]]
        elseif n == 2
            dirs = [N[-1, 0], N[0, -1], N[1, 1]]
        elseif n == 3
            dirs = [N[-1, 0, 0], N[0, -1, 0], N[0, 0, -1], N[1, 1, 1]]
        end
        dir = CustomDirections(dirs)
        @test isbounding(dir)
        @test !isbounding(CustomDirections(dirs[1:end-1]))
        P = overapproximate(X, dir)
        @test P isa HPolytope && length(constraints_list(P)) == length(dirs)
    end
end

# default Float64 constructors
BoxDirections(3)
OctDirections(3)
BoxDiagDirections(3)
SphericalDirections(3)


text:
"""
Linear evolution of the NP equations: this includes the
Teukolsky equation and the Metric reconstruction equations.

For more details see
   Loutrel et. al. Phys.Rev.D 103 (2021) 10, 104017, arXiv:2008.11770

   Ripley et. al. Phys.Rev.D 103 (2021) 1040180, arXiv:2010.00162
"""
module LinearEvolution

include("Radial.jl")
include("Sphere.jl")
include("GHP.jl")
include("BackgroundNP.jl")
include("Fields.jl")
include("Evolution.jl")

using .Radial: set_d1!, set_d2!, filter!
using .Sphere: angular_matrix_mult!
using .GHP: GHP_ops
using .Fields: Field
using .BackgroundNP: NP_0
using .Evolution: Evo_lin_f

export Linear_evolution!

const half  = 1.0/2.0
const third = 1.0/3.0
const sixth = 1.0/6.0

function set_psi4_p_k!(;
      kp::Array{ComplexF64,2},
      f_rd1::Array{ComplexF64,2},
      f_rd2::Array{ComplexF64,2},
      sph_lap::Array{ComplexF64,2},
      p_rd1::Array{ComplexF64,2},
      f::Array{ComplexF64,2},
      p::Array{ComplexF64,2},
      A_pp::Array{Float64,2},
      A_pq::Array{Float64,2},
      B_pp::Array{ComplexF64,2},
      B_pq::Array{ComplexF64,2},
      B_pf::Array{ComplexF64,2}
   )::Nothing
   nx, ny = size(kp)

   for j=1:ny
      for i=1:nx
         kp[i,j] = (
               A_pp[i,j] * p_rd1[i,j]
             + A_pq[i,j] * f_rd2[i,j]
             + B_pp[i,j] * p[i,j]
             + B_pq[i,j] * f_rd1[i,j]
             + B_pf[i,j] * f[i,j]

             + sph_lap[i,j]
      )
      end
   end
   return nothing
end

function set_psi4_f_k!(;
      k::Array{ComplexF64,2},
      p::Array{ComplexF64,2}
     )::Nothing
   nx, ny = size(k)
   for j=1:ny
      for i=1:nx
         k[i,j] = p[i,j]
      end
   end
   return nothing
end

function set_psi3_k!(
      k::Array{ComplexF64},
      psi3::Array{ComplexF64},
      psi4::Array{ComplexF64},
      psi4_edth::Array{ComplexF64},
      mu_0::Array{ComplexF64},
      tau_0::Array{ComplexF64},
      R::Vector{Float64},
   )::Nothing
   nx, ny = size(k)
   for j=1:ny
      for i=1:nx
         k[i,j] = (-
                   4.0*R[i]*mu_0[i,j]*psi3[i,j]
                   -
                   R[i]*tau_0[i,j]*psi4[i,j]
                   +
                   psi4_edth[i,j]
                  )
      end
   end
   return nothing
end

function set_lam_k!(
      k::Array{ComplexF64,2},
      lam::Array{ComplexF64,2},
      psi4::Array{ComplexF64,2},
      mu_0::Array{ComplexF64},
      R::Vector{Float64}
   )::Nothing
   nx, ny = size(k)
   for j=1:ny
      for i=1:nx
         k[i,j] = (-
                   R[i]*(mu_0[i,j] + conj(mu_0[i,j]))*lam[i,j]
                   -
                   psi4[i,j]
                  )
      end
   end
   return nothing
end

function set_psi2_k!(
      k::Array{ComplexF64,2},
      psi2::Array{ComplexF64,2},
      psi3::Array{ComplexF64,2},
      psi3_edth::Array{ComplexF64,2},
      mu_0::Array{ComplexF64,2},
      tau_0::Array{ComplexF64,2},
      R::Vector{Float64}
   )::Nothing
   nx, ny = size(k)
   for j=1:ny
      for i=1:nx
         k[i,j] = (-
                   3.0*R[i]*mu_0[i,j] *psi2[i,j]
                   -
                   2.0*R[i]*tau_0[i,j]*psi3[i,j]
                   +
                   psi3_edth[i,j]
                  )
      end
   end
   return nothing
end

function set_hmbmb_k!(
      k::Array{ComplexF64,2},
      hmbmb::Array{ComplexF64,2},
      lam::Array{ComplexF64,2},
      mu_0::Array{ComplexF64,2},
      R::Vector{Float64}
   )::Nothing
   nx, ny = size(k)
   for j=1:ny
      for i=1:nx
         k[i,j] = (
                   R[i]*(mu_0[i,j] - conj(mu_0[i,j]))*hmbmb[i,j]
                   -
                   2.0*lam[i,j]
                  )
      end
   end
   return nothing
end

function set_pi_k!(
      k::Array{ComplexF64,2},
      lam::Array{ComplexF64,2},
      hmbmb::Array{ComplexF64,2},
      psi3::Array{ComplexF64,2},
      mu_0::Array{ComplexF64,2},
      tau_0::Array{ComplexF64,2},
      pi_0::Array{ComplexF64,2},
      R::Vector{Float64}
   )::Nothing
   nx, ny = size(k)
   for j=1:ny
      for i=1:nx
         k[i,j] = (-
                   R[i]*(conj(pi_0[i,j]) + tau_0[i,j])*lam[i,j]
                   +
                   (R[i]^2)*0.5*mu_0[i,j]*(
                                           conj(pi_0[i,j])
                                           +
                                           tau_0[i,j]
                                          )*hmbmb[i,j]
                   -
                   psi3[i,j]
                  )
      end
   end
   return nothing
end

function set_hlmb_k!(
      k::Array{ComplexF64,2},
      lam::Array{ComplexF64,2},
      hlmb::Array{ComplexF64,2},
      hmbmb::Array{ComplexF64,2},
      psi3::Array{ComplexF64,2},
      pi_l::Array{ComplexF64,2},
      mu_0::Array{ComplexF64,2},
      tau_0::Array{ComplexF64,2},
      R::Vector{Float64}
   )::Nothing
   nx, ny = size(k)
   for j=1:ny
      for i=1:nx
         k[i,j] = (-
                   R[i]*conj(mu_0[i,j])*hlmb[i,j]
                   -
                   2.0*pi_l[i,j]
                   -
                   R[i]*tau_0[i,j]*hmbmb[i,j]
                  )
      end
   end
   return nothing
end

function set_muhll_k!(
      k::Array{ComplexF64,2},
      muhll::Array{ComplexF64,2},
      hlmb_edth::Array{ComplexF64,2},
      hlmb::Array{ComplexF64,2},
      pi_edth::Array{ComplexF64,2},
      pi_l::Array{ComplexF64,2},
      psi2::Array{ComplexF64,2},
      pi_nm::Array{ComplexF64,2},
      hmbmb_edth_nm::Array{ComplexF64,2},
      hmbmb_nm::Array{ComplexF64,2},
      hlmb_edth_nm::Array{ComplexF64,2},
      hlmb_nm::Array{ComplexF64,2},
      mu_0::Array{ComplexF64,2},
      tau_0::Array{ComplexF64,2},
      pi_0::Array{ComplexF64,2},
      R::Vector{Float64}
   )::Nothing
   nx, ny = size(k)
   for j=1:ny
      for i=1:nx
         k[i,j] = (-
                   R[i]*conj(mu_0[i,j])*muhll[i,j]
                   -
                   R[i]*mu_0[i,j]*hlmb_edth[i,j]
                   -
                   (R[i]^2)*mu_0[i,j]*(
                                       conj(pi_0[i,j])
                                       +
                                       2.0*tau_0[i,j]
                                      )*hlmb[i,j]
                   -
                   2.0*pi_edth[i,j]
                   -
                   2.0*R[i]*conj(pi_0[i,j])*pi_l[i,j]
                   -
                   2.0*psi2[i,j]
                   -
                   2.0*R[i]*pi_0[i,j]*conj(pi_nm[i,j])
                   -
                   R[i]*pi_0[i,j]*conj(hmbmb_edth_nm[i,j])
                   +
                   (R[i]^2)*(pi_0[i,j]^2)*conj(hmbmb_nm[i,j])
                   +
                   R[i]*mu_0[i,j]*conj(hlmb_edth_nm[i,j])
                   +
                   (R[i]^2)*(-
                             3.0*mu_0[i,j]*pi_0[i,j]
                             +
                             2.0*conj(mu_0[i,j])*pi_0[i,j]
                             -
                             2.0*mu_0[i,j]*conj(tau_0[i,j])
                            )*conj(hlmb_nm[i,j])
                  )
      end
   end
   return nothing
end

function finish_linear_reconstruction_k!(;
      k::Array{ComplexF64,2},
      f::Array{ComplexF64,2},
      DR::Array{ComplexF64,2},
      R::Vector{Float64},
      falloff::Int64,
      bhm::Float64,
      cl::Float64,
      dr::Float64
     )::Nothing
   nx, ny = size(k)
   set_d1!(DR,f,dr)
   for j=1:ny
      for i=1:nx
         k[i,j] += (-
                    ((R[i]/cl)^2)*DR[i,j]
                    -
                    (falloff*R[i]/(cl^2))*f[i,j]
                   )/(
                      2.0+(4.0*bhm*R[i]/(cl^2))
                     )
      end
   end
   return nothing
end

function set_linear_k!(;
      psi4_f,  psi4_f_l::Array{ComplexF64,2},
      psi4_p,  psi4_p_l::Array{ComplexF64,2},
      psi3_f,  psi3_l  ::Array{ComplexF64,2},
      psi2_f,  psi2_l  ::Array{ComplexF64,2},
      lam_f,   lam_l   ::Array{ComplexF64,2},
      pi_f,    pi_l    ::Array{ComplexF64,2}, pi_nm_l::Array{ComplexF64,2},
      hmbmb_f,    hmbmb_l   ::Array{ComplexF64,2},
      hmbmb_f_nm, hmbmb_nm_l::Array{ComplexF64,2},
      hlmb_f,     hlmb_l    ::Array{ComplexF64,2},
      hlmb_f_nm,  hlmb_nm_l ::Array{ComplexF64,2},
      muhll_f, muhll_l::Array{ComplexF64,2},
      Evo, Op, NP,
      R::Vector{Float64},
      m_ang::Int64,
      bhm::Float64,
      cl::Float64,
      dr::Float64
     )::Nothing
   psi4_f_k       = psi4_f.k
   psi4_f_rd1     = psi4_f.rad_d1
   psi4_f_rd2     = psi4_f.rad_d2
   psi4_f_sph_lap = psi4_f.sph_lap
   psi4_f_raised  = psi4_f.raised
   psi4_f_edth    = psi4_f.edth

   psi4_p_k   = psi4_p.k
   psi4_p_rd1 = psi4_p.rad_d1

   psi3_k      = psi3_f.k
   psi3_raised = psi3_f.raised
   psi3_edth   = psi3_f.edth
   psi3_DR     = psi3_f.rad_d1

   psi2_k      = psi2_f.k
   psi2_DR     = psi2_f.rad_d1

   lam_k      = lam_f.k
   lam_DR     = lam_f.rad_d1

   pi_k      = pi_f.k
   pi_DR     = pi_f.rad_d1
   pi_raised = pi_f.raised
   pi_edth   = pi_f.edth

   hmbmb_k      = hmbmb_f.k
   hmbmb_DR     = hmbmb_f.rad_d1

   hmbmb_k_nm      = hmbmb_f_nm.k
   hmbmb_raised_nm = hmbmb_f_nm.raised
   hmbmb_edth_nm   = hmbmb_f_nm.edth

   hlmb_k      = hlmb_f.k
   hlmb_DR     = hlmb_f.rad_d1
   hlmb_raised = hlmb_f.raised
   hlmb_edth   = hlmb_f.edth

   hlmb_k_nm      = hlmb_f_nm.k
   hlmb_raised_nm = hlmb_f_nm.raised
   hlmb_edth_nm   = hlmb_f_nm.edth

   muhll_k      = muhll_f.k
   muhll_DR     = muhll_f.rad_d1

   A_pp = Evo.A_pp
   A_pq = Evo.A_pq
   B_pp = Evo.B_pp
   B_pq = Evo.B_pq
   B_pf = Evo.B_pf

   laplM = Evo.S_lapl

   mu_0  = NP.mu_0
   tau_0 = NP.tau_0
   pi_0  = NP.pi_0

   set_d1!(psi4_f_rd1, psi4_f_l, dr)
   set_d1!(psi4_p_rd1, psi4_p_l, dr)
   set_d2!(psi4_f_rd2, psi4_f_l, dr)

   angular_matrix_mult!(psi4_f_sph_lap, psi4_f_l, laplM)

   set_psi4_p_k!(kp=psi4_p_k,
      f_rd1=psi4_f_rd1, f_rd2=psi4_f_rd2,
      sph_lap=psi4_f_sph_lap, p_rd1=psi4_p_rd1,
      f=psi4_f_l, p=psi4_p_l,
      A_pp=A_pp, A_pq=A_pq, B_pp=B_pp, B_pq=B_pq, B_pf=B_pf)

   set_psi4_f_k!(k=psi4_f_k, p=psi4_p_l)
   ##===============
   GHP.set_edth!(edth=psi4_f_edth,spin=psi4_f.spin,boost=psi4_f.boost,m_ang=m_ang,f=psi4_f_l,DT=psi4_f_k,raised=psi4_f_raised,Op=Op)

   finish_linear_reconstruction_k!(k=psi3_k, f=psi3_l, DR=psi3_DR, R=R, falloff=psi3_f.falloff, bhm=bhm, cl=cl, dr=dr)
   ##===============
   set_lam_k!(lam_k, lam_l, psi4_f_l, mu_0, R)
   finish_linear_reconstruction_k!(k=lam_k, f=lam_l, DR=lam_DR, R=R, falloff=lam_f.falloff, bhm=bhm, cl=cl, dr=dr)
   ##===============
   GHP.set_edth!(edth=psi3_edth,spin=psi3_f.spin,boost=psi3_f.boost,m_ang=m_ang,f=psi3_l,DT=psi3_k,raised=psi3_raised,Op=Op)

   set_psi2_k!(psi2_k, psi2_l, psi3_l, psi3_edth, mu_0, tau_0, R)
   finish_linear_reconstruction_k!(k=psi2_k, f=psi2_l, DR=psi2_DR, R=R, falloff=psi2_f.falloff, bhm=bhm, cl=cl, dr=dr)
   ##===============
   set_hmbmb_k!(hmbmb_k, hmbmb_l, lam_l, mu_0, R)
   finish_linear_reconstruction_k!(k=hmbmb_k, f=hmbmb_l, DR=hmbmb_DR, R=R, falloff=hmbmb_f.falloff, bhm=bhm, cl=cl, dr=dr)
   ##===============
   set_pi_k!(pi_k, lam_l, hmbmb_l, psi3_l, mu_0, tau_0, pi_0, R)
   finish_linear_reconstruction_k!(k=pi_k, f=pi_l, DR=pi_DR, R=R, falloff=pi_f.falloff, bhm=bhm, cl=cl, dr=dr)
   ##===============
   set_hlmb_k!(hlmb_k, lam_l, hlmb_l, hmbmb_l, psi3_l, pi_l, mu_0, tau_0, R)
   finish_linear_reconstruction_k!(k=hlmb_k, f=hlmb_l, DR=hlmb_DR, R=R, falloff=hlmb_f.falloff, bhm=bhm, cl=cl, dr=dr)
   ##===============
   GHP.set_edth!(edth=hlmb_edth,    spin=hlmb_f.spin, boost=hlmb_f.boost, m_ang= m_ang,f=hlmb_l,    DT=hlmb_k,    raised=hlmb_raised,    Op=Op)
   GHP.set_edth!(edth=pi_edth,      spin=pi_f.spin,   boost=pi_f.boost,   m_ang= m_ang,f=pi_l,      DT=pi_k,      raised=pi_raised,      Op=Op)
   GHP.set_edth!(edth=hmbmb_edth_nm,spin=hmbmb_f.spin,boost=hmbmb_f.boost,m_ang=-m_ang,f=hmbmb_nm_l,DT=hmbmb_k_nm,raised=hmbmb_raised_nm,Op=Op)
   GHP.set_edth!(edth=hlmb_edth_nm, spin=hlmb_f.spin, boost=hlmb_f.boost, m_ang=-m_ang,f=hlmb_nm_l, DT=hlmb_k_nm, raised=hlmb_raised_nm, Op=Op)

   set_muhll_k!(muhll_k, muhll_l, hlmb_edth, hlmb_l, pi_edth, pi_l, psi2_l,
                pi_nm_l, hmbmb_edth_nm, hmbmb_nm_l, hlmb_edth_nm, hlmb_nm_l,
                mu_0, tau_0, pi_0,
                R
               )
   finish_linear_reconstruction_k!(k=muhll_k, f=muhll_l, DR=muhll_DR, R=R, falloff=muhll_f.falloff, bhm=bhm, cl=cl, dr=dr)

   return nothing
end

"""
One giant RK4 evolution for all of the linear variables
"""
function Linear_evolution!(;
      psi4_f_pm, psi4_f_nm,
      psi4_p_pm, psi4_p_nm,
      psi3_pm,   psi3_nm,
      psi2_pm,   psi2_nm,
      lam_pm,    lam_nm,
      pi_pm,     pi_nm,
      hmbmb_pm,  hmbmb_nm,
      hlmb_pm,   hlmb_nm,
      muhll_pm,  muhll_nm,
      Evo_pm,    Evo_nm,
      Op_pm,     Op_nm,
      NP,
      R::Vector{Float64},
      m_ang::Int64,
      bhm::Float64,
      cl::Float64,
      dr::Float64,
      dt::Float64
   )::Nothing
   psi4_f_pm_n   = psi4_f_pm.n  ; psi4_f_nm_n   = psi4_f_nm.n
   psi4_f_pm_tmp = psi4_f_pm.tmp; psi4_f_nm_tmp = psi4_f_nm.tmp
   psi4_f_pm_np1 = psi4_f_pm.np1; psi4_f_nm_np1 = psi4_f_nm.np1
   psi4_f_pm_k   = psi4_f_pm.k;   psi4_f_nm_k   = psi4_f_nm.k

   psi4_p_pm_n   = psi4_p_pm.n;   psi4_p_nm_n   = psi4_p_nm.n
   psi4_p_pm_tmp = psi4_p_pm.tmp; psi4_p_nm_tmp = psi4_p_nm.tmp
   psi4_p_pm_np1 = psi4_p_pm.np1; psi4_p_nm_np1 = psi4_p_nm.np1
   psi4_p_pm_k   = psi4_p_pm.k;   psi4_p_nm_k   = psi4_p_nm.k

   psi3_pm_n   = psi3_pm.n;   psi3_nm_n   = psi3_nm.n
   psi3_pm_tmp = psi3_pm.tmp; psi3_nm_tmp = psi3_nm.tmp
   psi3_pm_np1 = psi3_pm.np1; psi3_nm_np1 = psi3_nm.np1
   psi3_pm_k   = psi3_pm.k;   psi3_nm_k   = psi3_nm.k

   psi2_pm_n   = psi2_pm.n;   psi2_nm_n   = psi2_nm.n
   psi2_pm_tmp = psi2_pm.tmp; psi2_nm_tmp = psi2_nm.tmp
   psi2_pm_np1 = psi2_pm.np1; psi2_nm_np1 = psi2_nm.np1
   psi2_pm_k   = psi2_pm.k;   psi2_nm_k   = psi2_nm.k

   lam_pm_n   = lam_pm.n;   lam_nm_n   = lam_nm.n
   lam_pm_tmp = lam_pm.tmp; lam_nm_tmp = lam_nm.tmp
   lam_pm_np1 = lam_pm.np1; lam_nm_np1 = lam_nm.np1
   lam_pm_k   = lam_pm.k;   lam_nm_k   = lam_nm.k

   pi_pm_n   = pi_pm.n;   pi_nm_n   = pi_nm.n
   pi_pm_tmp = pi_pm.tmp; pi_nm_tmp = pi_nm.tmp
   pi_pm_np1 = pi_pm.np1; pi_nm_np1 = pi_nm.np1
   pi_pm_k   = pi_pm.k;   pi_nm_k   = pi_nm.k

   hmbmb_pm_n   = hmbmb_pm.n;   hmbmb_nm_n   = hmbmb_nm.n
   hmbmb_pm_tmp = hmbmb_pm.tmp; hmbmb_nm_tmp = hmbmb_nm.tmp
   hmbmb_pm_np1 = hmbmb_pm.np1; hmbmb_nm_np1 = hmbmb_nm.np1
   hmbmb_pm_k   = hmbmb_pm.k;   hmbmb_nm_k   = hmbmb_nm.k

   hlmb_pm_n   = hlmb_pm.n;   hlmb_nm_n   = hlmb_nm.n
   hlmb_pm_tmp = hlmb_pm.tmp; hlmb_nm_tmp = hlmb_nm.tmp
   hlmb_pm_np1 = hlmb_pm.np1; hlmb_nm_np1 = hlmb_nm.np1
   hlmb_pm_k   = hlmb_pm.k;   hlmb_nm_k   = hlmb_nm.k

   muhll_pm_n   = muhll_pm.n;   muhll_nm_n   = muhll_nm.n
   muhll_pm_tmp = muhll_pm.tmp; muhll_nm_tmp = muhll_nm.tmp
   muhll_pm_np1 = muhll_pm.np1; muhll_nm_np1 = muhll_nm.np1
   muhll_pm_k   = muhll_pm.k;   muhll_nm_k   = muhll_nm.k

   nx, ny = size(psi4_f_pm_n)

   fltrM_pm = Evo_pm.S_fltr
   fltrM_nm = Evo_nm.S_fltr

   set_linear_k!(
      psi4_f=psi4_f_pm,    psi4_f_l=psi4_f_pm_n,
      psi4_p=psi4_p_pm,    psi4_p_l=psi4_p_pm_n,
      psi3_f=psi3_pm,      psi3_l=psi3_pm_n,
      psi2_f=psi2_pm,      psi2_l=psi2_pm_n,
      lam_f=lam_pm,        lam_l=lam_pm_n,
      pi_f=pi_pm,          pi_l=pi_pm_n,                            pi_nm_l=pi_nm_n,
      hmbmb_f=hmbmb_pm,    hmbmb_l=hmbmb_pm_n, hmbmb_f_nm=hmbmb_nm, hmbmb_nm_l=hmbmb_nm_n,
      hlmb_f=hlmb_pm,      hlmb_l=hlmb_pm_n,   hlmb_f_nm=hlmb_nm,   hlmb_nm_l=hlmb_nm_n,
      muhll_f=muhll_pm,    muhll_l=muhll_pm_n,
      Evo=Evo_pm, Op=Op_pm, NP=NP,
      R=R,
      m_ang=m_ang,
      bhm=bhm, cl=cl, dr=dr
     )

   set_linear_k!(
      psi4_f=psi4_f_nm,    psi4_f_l=psi4_f_nm_n,
      psi4_p=psi4_p_nm,    psi4_p_l=psi4_p_nm_n,
      psi3_f=psi3_nm,      psi3_l=psi3_nm_n,
      psi2_f=psi2_nm,      psi2_l=psi2_nm_n,
      lam_f=lam_nm,        lam_l=lam_nm_n,
      pi_f=pi_nm,          pi_l=pi_nm_n,                            pi_nm_l=pi_pm_n,
      hmbmb_f=hmbmb_nm,    hmbmb_l=hmbmb_nm_n, hmbmb_f_nm=hmbmb_pm, hmbmb_nm_l=hmbmb_pm_n,
      hlmb_f=hlmb_nm,      hlmb_l=hlmb_nm_n,   hlmb_f_nm=hlmb_pm,   hlmb_nm_l=hlmb_pm_n,
      muhll_f=muhll_nm,    muhll_l=muhll_nm_n,
      Evo=Evo_pm, Op=Op_nm, NP=NP,
      R=R,
      m_ang=-m_ang,
      bhm=bhm, cl=cl, dr=dr
     )

   for j=1:ny
      for i=1:nx
         psi4_f_pm_tmp[i,j] = psi4_f_pm_n[i,j] + half*dt*psi4_f_pm_k[i,j]
         psi4_p_pm_tmp[i,j] = psi4_p_pm_n[i,j] + half*dt*psi4_p_pm_k[i,j]
         psi3_pm_tmp[i,j]   = psi3_pm_n[i,j]   + half*dt*psi3_pm_k[i,j]
         psi2_pm_tmp[i,j]   = psi2_pm_n[i,j]   + half*dt*psi2_pm_k[i,j]
         lam_pm_tmp[i,j]    = lam_pm_n[i,j]    + half*dt*lam_pm_k[i,j]
         pi_pm_tmp[i,j]     = pi_pm_n[i,j]     + half*dt*pi_pm_k[i,j]
         hmbmb_pm_tmp[i,j]  = hmbmb_pm_n[i,j]  + half*dt*hmbmb_pm_k[i,j]
         hlmb_pm_tmp[i,j]   = hlmb_pm_n[i,j]   + half*dt*hlmb_pm_k[i,j]
         muhll_pm_tmp[i,j]  = muhll_pm_n[i,j]  + half*dt*muhll_pm_k[i,j]

         psi4_f_pm_np1[i,j] = psi4_f_pm_n[i,j] + sixth*dt*psi4_f_pm_k[i,j]
         psi4_p_pm_np1[i,j] = psi4_p_pm_n[i,j] + sixth*dt*psi4_p_pm_k[i,j]
         psi3_pm_np1[i,j]   = psi3_pm_n[i,j]   + sixth*dt*psi3_pm_k[i,j]
         psi2_pm_np1[i,j]   = psi2_pm_n[i,j]   + sixth*dt*psi2_pm_k[i,j]
         lam_pm_np1[i,j]    = lam_pm_n[i,j]    + sixth*dt*lam_pm_k[i,j]
         pi_pm_np1[i,j]     = pi_pm_n[i,j]     + sixth*dt*pi_pm_k[i,j]
         hmbmb_pm_np1[i,j]  = hmbmb_pm_n[i,j]  + sixth*dt*hmbmb_pm_k[i,j]
         hlmb_pm_np1[i,j]   = hlmb_pm_n[i,j]   + sixth*dt*hlmb_pm_k[i,j]
         muhll_pm_np1[i,j]  = muhll_pm_n[i,j]  + sixth*dt*muhll_pm_k[i,j]

         psi4_f_nm_tmp[i,j] = psi4_f_nm_n[i,j] + half*dt*psi4_f_nm_k[i,j]
         psi4_p_nm_tmp[i,j] = psi4_p_nm_n[i,j] + half*dt*psi4_p_nm_k[i,j]
         psi3_nm_tmp[i,j]   = psi3_nm_n[i,j]   + half*dt*psi3_nm_k[i,j]
         psi2_nm_tmp[i,j]   = psi2_nm_n[i,j]   + half*dt*psi2_nm_k[i,j]
         lam_nm_tmp[i,j]    = lam_nm_n[i,j]    + half*dt*lam_nm_k[i,j]
         pi_nm_tmp[i,j]     = pi_nm_n[i,j]     + half*dt*pi_nm_k[i,j]
         hmbmb_nm_tmp[i,j]  = hmbmb_nm_n[i,j]  + half*dt*hmbmb_nm_k[i,j]
         hlmb_nm_tmp[i,j]   = hlmb_nm_n[i,j]   + half*dt*hlmb_nm_k[i,j]
         muhll_nm_tmp[i,j]  = muhll_nm_n[i,j]  + half*dt*muhll_nm_k[i,j]

         psi4_f_nm_np1[i,j] = psi4_f_nm_n[i,j] + sixth*dt*psi4_f_nm_k[i,j]
         psi4_p_nm_np1[i,j] = psi4_p_nm_n[i,j] + sixth*dt*psi4_p_nm_k[i,j]
         psi3_nm_np1[i,j]   = psi3_nm_n[i,j]   + sixth*dt*psi3_nm_k[i,j]
         psi2_nm_np1[i,j]   = psi2_nm_n[i,j]   + sixth*dt*psi2_nm_k[i,j]
         lam_nm_np1[i,j]    = lam_nm_n[i,j]    + sixth*dt*lam_nm_k[i,j]
         pi_nm_np1[i,j]     = pi_nm_n[i,j]     + sixth*dt*pi_nm_k[i,j]
         hmbmb_nm_np1[i,j]  = hmbmb_nm_n[i,j]  + sixth*dt*hmbmb_nm_k[i,j]
         hlmb_nm_np1[i,j]   = hlmb_nm_n[i,j]   + sixth*dt*hlmb_nm_k[i,j]
         muhll_nm_np1[i,j]  = muhll_nm_n[i,j]  + sixth*dt*muhll_nm_k[i,j]
      end
   end

   set_linear_k!(
      psi4_f=psi4_f_pm,    psi4_f_l=psi4_f_pm_tmp,
      psi4_p=psi4_p_pm,    psi4_p_l=psi4_p_pm_tmp,
      psi3_f=psi3_pm,      psi3_l=psi3_pm_tmp,
      psi2_f=psi2_pm,      psi2_l=psi2_pm_tmp,
      lam_f=lam_pm,        lam_l=lam_pm_tmp,
      pi_f=pi_pm,          pi_l=pi_pm_tmp,                            pi_nm_l=pi_nm_tmp,
      hmbmb_f=hmbmb_pm,    hmbmb_l=hmbmb_pm_tmp, hmbmb_f_nm=hmbmb_nm, hmbmb_nm_l=hmbmb_nm_tmp,
      hlmb_f=hlmb_pm,      hlmb_l=hlmb_pm_tmp,   hlmb_f_nm=hlmb_nm,   hlmb_nm_l=hlmb_nm_tmp,
      muhll_f=muhll_pm,    muhll_l=muhll_pm_tmp,
      Evo=Evo_pm, Op=Op_pm, NP=NP,
      R=R,
      m_ang=m_ang,
      bhm=bhm, cl=cl, dr=dr
     )

   set_linear_k!(
      psi4_f=psi4_f_nm,    psi4_f_l=psi4_f_nm_tmp,
      psi4_p=psi4_p_nm,    psi4_p_l=psi4_p_nm_tmp,
      psi3_f=psi3_nm,      psi3_l=psi3_nm_tmp,
      psi2_f=psi2_nm,      psi2_l=psi2_nm_tmp,
      lam_f=lam_nm,        lam_l=lam_nm_tmp,
      pi_f=pi_nm,          pi_l=pi_nm_tmp,                            pi_nm_l=pi_pm_tmp,
      hmbmb_f=hmbmb_nm,    hmbmb_l=hmbmb_nm_tmp, hmbmb_f_nm=hmbmb_pm, hmbmb_nm_l=hmbmb_pm_tmp,
      hlmb_f=hlmb_nm,      hlmb_l=hlmb_nm_tmp,   hlmb_f_nm=hlmb_pm,   hlmb_nm_l=hlmb_pm_tmp,
      muhll_f=muhll_nm,    muhll_l=muhll_nm_tmp,
      Evo=Evo_pm, Op=Op_nm, NP=NP,
      R=R,
      m_ang=-m_ang,
      bhm=bhm, cl=cl, dr=dr
     )

   for j=1:ny
      for i=1:nx
         psi4_f_pm_tmp[i,j] = psi4_f_pm_n[i,j] + half*dt*psi4_f_pm_k[i,j]
         psi4_p_pm_tmp[i,j] = psi4_p_pm_n[i,j] + half*dt*psi4_p_pm_k[i,j]
         psi3_pm_tmp[i,j]   = psi3_pm_n[i,j]   + half*dt*psi3_pm_k[i,j]
         psi2_pm_tmp[i,j]   = psi2_pm_n[i,j]   + half*dt*psi2_pm_k[i,j]
         lam_pm_tmp[i,j]    = lam_pm_n[i,j]    + half*dt*lam_pm_k[i,j]
         pi_pm_tmp[i,j]     = pi_pm_n[i,j]     + half*dt*pi_pm_k[i,j]
         hmbmb_pm_tmp[i,j]  = hmbmb_pm_n[i,j]  + half*dt*hmbmb_pm_k[i,j]
         hlmb_pm_tmp[i,j]   = hlmb_pm_n[i,j]   + half*dt*hlmb_pm_k[i,j]
         muhll_pm_tmp[i,j]  = muhll_pm_n[i,j]  + half*dt*muhll_pm_k[i,j]

         psi4_f_pm_np1[i,j] += third*dt*psi4_f_pm_k[i,j]
         psi4_p_pm_np1[i,j] += third*dt*psi4_p_pm_k[i,j]
         psi3_pm_np1[i,j]   += third*dt*psi3_pm_k[i,j]
         psi2_pm_np1[i,j]   += third*dt*psi2_pm_k[i,j]
         lam_pm_np1[i,j]    += third*dt*lam_pm_k[i,j]
         pi_pm_np1[i,j]     += third*dt*pi_pm_k[i,j]
         hmbmb_pm_np1[i,j]  += third*dt*hmbmb_pm_k[i,j]
         hlmb_pm_np1[i,j]   += third*dt*hlmb_pm_k[i,j]
         muhll_pm_np1[i,j]  += third*dt*muhll_pm_k[i,j]

         psi4_f_nm_tmp[i,j] = psi4_f_nm_n[i,j] + half*dt*psi4_f_nm_k[i,j]
         psi4_p_nm_tmp[i,j] = psi4_p_nm_n[i,j] + half*dt*psi4_p_nm_k[i,j]
         psi3_nm_tmp[i,j]   = psi3_nm_n[i,j]   + half*dt*psi3_nm_k[i,j]
         psi2_nm_tmp[i,j]   = psi2_nm_n[i,j]   + half*dt*psi2_nm_k[i,j]
         lam_nm_tmp[i,j]    = lam_nm_n[i,j]    + half*dt*lam_nm_k[i,j]
         pi_nm_tmp[i,j]     = pi_nm_n[i,j]     + half*dt*pi_nm_k[i,j]
         hmbmb_nm_tmp[i,j]  = hmbmb_nm_n[i,j]  + half*dt*hmbmb_nm_k[i,j]
         hlmb_nm_tmp[i,j]   = hlmb_nm_n[i,j]   + half*dt*hlmb_nm_k[i,j]
         muhll_nm_tmp[i,j]  = muhll_nm_n[i,j]  + half*dt*muhll_nm_k[i,j]

         psi4_f_nm_np1[i,j] += third*dt*psi4_f_nm_k[i,j]
         psi4_p_nm_np1[i,j] += third*dt*psi4_p_nm_k[i,j]
         psi3_nm_np1[i,j]   += third*dt*psi3_nm_k[i,j]
         psi2_nm_np1[i,j]   += third*dt*psi2_nm_k[i,j]
         lam_nm_np1[i,j]    += third*dt*lam_nm_k[i,j]
         pi_nm_np1[i,j]     += third*dt*pi_nm_k[i,j]
         hmbmb_nm_np1[i,j]  += third*dt*hmbmb_nm_k[i,j]
         hlmb_nm_np1[i,j]   += third*dt*hlmb_nm_k[i,j]
         muhll_nm_np1[i,j]  += third*dt*muhll_nm_k[i,j]
      end
   end

   set_linear_k!(
      psi4_f=psi4_f_pm,    psi4_f_l=psi4_f_pm_tmp,
      psi4_p=psi4_p_pm,    psi4_p_l=psi4_p_pm_tmp,
      psi3_f=psi3_pm,      psi3_l=psi3_pm_tmp,
      psi2_f=psi2_pm,      psi2_l=psi2_pm_tmp,
      lam_f=lam_pm,        lam_l=lam_pm_tmp,
      pi_f=pi_pm,          pi_l=pi_pm_tmp,                            pi_nm_l=pi_nm_tmp,
      hmbmb_f=hmbmb_pm,    hmbmb_l=hmbmb_pm_tmp, hmbmb_f_nm=hmbmb_nm, hmbmb_nm_l=hmbmb_nm_tmp,
      hlmb_f=hlmb_pm,      hlmb_l=hlmb_pm_tmp,   hlmb_f_nm=hlmb_nm,   hlmb_nm_l=hlmb_nm_tmp,
      muhll_f=muhll_pm,    muhll_l=muhll_pm_tmp,
      Evo=Evo_pm, Op=Op_pm, NP=NP,
      R=R,
      m_ang=m_ang,
      bhm=bhm, cl=cl, dr=dr
     )

   set_linear_k!(
      psi4_f=psi4_f_nm,    psi4_f_l=psi4_f_nm_tmp,
      psi4_p=psi4_p_nm,    psi4_p_l=psi4_p_nm_tmp,
      psi3_f=psi3_nm,      psi3_l=psi3_nm_tmp,
      psi2_f=psi2_nm,      psi2_l=psi2_nm_tmp,
      lam_f=lam_nm,        lam_l=lam_nm_tmp,
      pi_f=pi_nm,          pi_l=pi_nm_tmp,                            pi_nm_l=pi_pm_tmp,
      hmbmb_f=hmbmb_nm,    hmbmb_l=hmbmb_nm_tmp, hmbmb_f_nm=hmbmb_pm, hmbmb_nm_l=hmbmb_pm_tmp,
      hlmb_f=hlmb_nm,      hlmb_l=hlmb_nm_tmp,   hlmb_f_nm=hlmb_pm,   hlmb_nm_l=hlmb_pm_tmp,
      muhll_f=muhll_nm,    muhll_l=muhll_nm_tmp,
      Evo=Evo_pm, Op=Op_nm, NP=NP,
      R=R,
      m_ang=-m_ang,
      bhm=bhm, cl=cl, dr=dr
     )

   for j=1:ny
      for i=1:nx
         psi4_f_pm_tmp[i,j] = psi4_f_pm_n[i,j] + dt*psi4_f_pm_k[i,j]
         psi4_p_pm_tmp[i,j] = psi4_p_pm_n[i,j] + dt*psi4_p_pm_k[i,j]
         psi3_pm_tmp[i,j]   = psi3_pm_n[i,j]   + dt*psi3_pm_k[i,j]
         psi2_pm_tmp[i,j]   = psi2_pm_n[i,j]   + dt*psi2_pm_k[i,j]
         lam_pm_tmp[i,j]    = lam_pm_n[i,j]    + dt*lam_pm_k[i,j]
         pi_pm_tmp[i,j]     = pi_pm_n[i,j]     + dt*pi_pm_k[i,j]
         hmbmb_pm_tmp[i,j]  = hmbmb_pm_n[i,j]  + dt*hmbmb_pm_k[i,j]
         hlmb_pm_tmp[i,j]   = hlmb_pm_n[i,j]   + dt*hlmb_pm_k[i,j]
         muhll_pm_tmp[i,j]  = muhll_pm_n[i,j]  + dt*muhll_pm_k[i,j]

         psi4_f_pm_np1[i,j] += third*dt*psi4_f_pm_k[i,j]
         psi4_p_pm_np1[i,j] += third*dt*psi4_p_pm_k[i,j]
         psi3_pm_np1[i,j]   += third*dt*psi3_pm_k[i,j]
         psi2_pm_np1[i,j]   += third*dt*psi2_pm_k[i,j]
         lam_pm_np1[i,j]    += third*dt*lam_pm_k[i,j]
         pi_pm_np1[i,j]     += third*dt*pi_pm_k[i,j]
         hmbmb_pm_np1[i,j]  += third*dt*hmbmb_pm_k[i,j]
         hlmb_pm_np1[i,j]   += third*dt*hlmb_pm_k[i,j]
         muhll_nm_np1[i,j]  += third*dt*muhll_nm_k[i,j]

         psi4_f_nm_tmp[i,j] = psi4_f_nm_n[i,j] + dt*psi4_f_nm_k[i,j]
         psi4_p_nm_tmp[i,j] = psi4_p_nm_n[i,j] + dt*psi4_p_nm_k[i,j]
         psi3_nm_tmp[i,j]   = psi3_nm_n[i,j]   + dt*psi3_nm_k[i,j]
         psi2_nm_tmp[i,j]   = psi2_nm_n[i,j]   + dt*psi2_nm_k[i,j]
         lam_nm_tmp[i,j]    = lam_nm_n[i,j]    + dt*lam_nm_k[i,j]
         pi_nm_tmp[i,j]     = pi_nm_n[i,j]     + dt*pi_nm_k[i,j]
         hmbmb_nm_tmp[i,j]  = hmbmb_nm_n[i,j]  + dt*hmbmb_nm_k[i,j]
         hlmb_nm_tmp[i,j]   = hlmb_nm_n[i,j]   + dt*hlmb_nm_k[i,j]
         muhll_nm_tmp[i,j]  = muhll_nm_n[i,j]  + dt*muhll_nm_k[i,j]

         psi4_f_nm_np1[i,j] += third*dt*psi4_f_nm_k[i,j]
         psi4_p_nm_np1[i,j] += third*dt*psi4_p_nm_k[i,j]
         psi3_nm_np1[i,j]   += third*dt*psi3_nm_k[i,j]
         psi2_nm_np1[i,j]   += third*dt*psi2_nm_k[i,j]
         lam_nm_np1[i,j]    += third*dt*lam_nm_k[i,j]
         pi_nm_np1[i,j]     += third*dt*pi_nm_k[i,j]
         hmbmb_nm_np1[i,j]  += third*dt*hmbmb_nm_k[i,j]
         hlmb_nm_np1[i,j]   += third*dt*hlmb_nm_k[i,j]
         muhll_nm_np1[i,j]  += third*dt*muhll_nm_k[i,j]
      end
   end

   set_linear_k!(
      psi4_f=psi4_f_pm,    psi4_f_l=psi4_f_pm_tmp,
      psi4_p=psi4_p_pm,    psi4_p_l=psi4_p_pm_tmp,
      psi3_f=psi3_pm,      psi3_l=psi3_pm_tmp,
      psi2_f=psi2_pm,      psi2_l=psi2_pm_tmp,
      lam_f=lam_pm,        lam_l=lam_pm_tmp,
      pi_f=pi_pm,          pi_l=pi_pm_tmp,                            pi_nm_l=pi_nm_tmp,
      hmbmb_f=hmbmb_pm,    hmbmb_l=hmbmb_pm_tmp, hmbmb_f_nm=hmbmb_nm, hmbmb_nm_l=hmbmb_nm_tmp,
      hlmb_f=hlmb_pm,      hlmb_l=hlmb_pm_tmp,   hlmb_f_nm=hlmb_nm,   hlmb_nm_l=hlmb_nm_tmp,
      muhll_f=muhll_pm,    muhll_l=muhll_pm_tmp,
      Evo=Evo_pm, Op=Op_pm, NP=NP,
      R=R,
      m_ang=m_ang,
      bhm=bhm, cl=cl, dr=dr
     )

   set_linear_k!(
      psi4_f=psi4_f_nm,    psi4_f_l=psi4_f_nm_tmp,
      psi4_p=psi4_p_nm,    psi4_p_l=psi4_p_nm_tmp,
      psi3_f=psi3_nm,      psi3_l=psi3_nm_tmp,
      psi2_f=psi2_nm,      psi2_l=psi2_nm_tmp,
      lam_f=lam_nm,        lam_l=lam_nm_tmp,
      pi_f=pi_nm,          pi_l=pi_nm_tmp,                            pi_nm_l=pi_pm_tmp,
      hmbmb_f=hmbmb_nm,    hmbmb_l=hmbmb_nm_tmp, hmbmb_f_nm=hmbmb_pm, hmbmb_nm_l=hmbmb_pm_tmp,
      hlmb_f=hlmb_nm,      hlmb_l=hlmb_nm_tmp,   hlmb_f_nm=hlmb_pm,   hlmb_nm_l=hlmb_pm_tmp,
      muhll_f=muhll_nm,    muhll_l=muhll_nm_tmp,
      Evo=Evo_pm, Op=Op_nm, NP=NP,
      R=R,
      m_ang=-m_ang,
      bhm=bhm, cl=cl, dr=dr
     )

   for j=1:ny
      for i=1:nx
         psi4_f_pm_np1[i,j] += sixth*dt*psi4_f_pm_k[i,j]
         psi4_p_pm_np1[i,j] += sixth*dt*psi4_p_pm_k[i,j]
         psi3_pm_np1[i,j]   += third*dt*psi3_pm_k[i,j]
         psi2_pm_np1[i,j]   += third*dt*psi2_pm_k[i,j]
         lam_pm_np1[i,j]    += third*dt*lam_pm_k[i,j]
         pi_pm_np1[i,j]     += third*dt*pi_pm_k[i,j]
         hmbmb_pm_np1[i,j]  += third*dt*hmbmb_pm_k[i,j]
         hlmb_pm_np1[i,j]   += third*dt*hlmb_pm_k[i,j]
         muhll_pm_np1[i,j]  += third*dt*muhll_pm_k[i,j]

         psi4_f_nm_np1[i,j] += sixth*dt*psi4_f_nm_k[i,j]
         psi4_p_nm_np1[i,j] += sixth*dt*psi4_p_nm_k[i,j]
         psi3_nm_np1[i,j]   += third*dt*psi3_nm_k[i,j]
         psi2_nm_np1[i,j]   += third*dt*psi2_nm_k[i,j]
         lam_nm_np1[i,j]    += third*dt*lam_nm_k[i,j]
         pi_nm_np1[i,j]     += third*dt*pi_nm_k[i,j]
         hmbmb_nm_np1[i,j]  += third*dt*hmbmb_nm_k[i,j]
         hlmb_nm_np1[i,j]   += third*dt*hlmb_nm_k[i,j]
         muhll_nm_np1[i,j]  += third*dt*muhll_nm_k[i,j]
      end
   end

   ## Only filter psi4 for now

   filter!(psi4_f_pm_np1,psi4_f_pm_tmp,0.5)
   filter!(psi4_p_pm_np1,psi4_p_pm_tmp,0.5)

   filter!(psi4_f_nm_np1,psi4_f_nm_tmp,0.5)
   filter!(psi4_p_nm_np1,psi4_p_nm_tmp,0.5)
   for j=1:ny
      for i=1:nx
         psi4_f_pm_tmp[i,j] = psi4_f_pm_np1[i,j]
         psi4_p_pm_tmp[i,j] = psi4_p_pm_np1[i,j]

         psi4_f_nm_tmp[i,j] = psi4_f_nm_np1[i,j]
         psi4_p_nm_tmp[i,j] = psi4_p_nm_np1[i,j]
      end
   end
   angular_matrix_mult!(psi4_f_pm_np1,psi4_f_pm_tmp,fltrM_pm)
   angular_matrix_mult!(psi4_p_pm_np1,psi4_p_pm_tmp,fltrM_pm)

   angular_matrix_mult!(psi4_f_nm_np1,psi4_f_nm_tmp,fltrM_nm)
   angular_matrix_mult!(psi4_p_nm_np1,psi4_p_nm_tmp,fltrM_nm)

   return nothing
end

end