text:
# MIT License
# Copyright (c) 2017: Xavier Gandibleux, Anthony Przybylski, Gauthier Soleilhac, and contributors.
mutable struct Triangle
    xr::solution
    xs::solution
    XΔ::Vector{solution}
    λ::Tuple{Int,Int}
    ub::Int
    lb::Int
    pending::Bool
end

Triangle(xr::solution,xs::solution) = begin
    λ1,λ2 = obj_2(xr) - obj_2(xs), obj_1(xs) - obj_1(xr)
    Triangle(xr,xs,solution[xr,xs],(λ1,λ2),λ1*obj_1(xr) + λ2*obj_2(xr),λ1*obj_1(xr) + λ2*obj_2(xs), true)
end

#const PartitionHeap = BinaryHeap{Partition, DataStructures.GreaterThan} # update 22-07-2021
const PartitionHeap = BinaryMaxHeap{Partition}

#Returns the partition p in τ maximizing p.zλ
function parent_partition!(τ::PartitionHeap)
    pop!(τ)
end

function clean!(τ::PartitionHeap, lb::Int)
    deleteat!(τ.valtree, findall(elt -> elt.zλ < lb, τ.valtree))
    heapify!(τ.valtree, Base.Order.Reverse)
    # @assert all(top(τ).zλ .>= zλ.(τ.valtree[2:end]))
end

#constructs the optimal path for a partition p=<v, U>
#new partitions are added to τ if their value for zλ is >= lb
function const_Path_New_Partition(p::Partition, τ::PartitionHeap, lb::Int, mono_pb::mono_problem)::BitArray
    t = p.v #The vertex from which we're constructing the path
    U = p.arcs #The imposed arcs
    U_prime = reverse(U)
    ϕ = falses(size(mono_pb)) #The path that will be returned

    while t.layer != 1
        if !isempty(U_prime) && last(head(U_prime)) == t #if an arc (s=>t) is imposed
            s = first(head(U_prime))
            U_prime = tail(U_prime)
        elseif inner_degree(t) == 1 #elseif the vertex only has one parent
            s = parent(t)
        else
            s1,s2 = parents(t)

            #(s => t) is the optimal arc, (s' => t) is the secondary arc, not imposed yet
            if s2.zλ + zλ(s2, t, mono_pb) == t.zλ
                s = s2
                s_prime = s1
            else
                s = s1
                s_prime = s2
            end

            if isempty(U_prime) #Create a partition imposing the secondary arc (s' => t)
                ϕ_zλ = zλ(p) - zλ(t) + zλ(s_prime) + zλ(s_prime, t, mono_pb)
                if ϕ_zλ >= lb #Add it only if it's interesting
                    ϕ_z1 = z1(p) - z1(t) + z1(s_prime) + z1(s_prime, t, mono_pb)
                    ϕ_z2 = z2(p) - z2(t) + z2(s_prime) + z2(s_prime, t, mono_pb)
                    #@assert ϕ_zλ <= p.zλ
                    ♇ = Partition(p.v, List(s_prime=>t, U), p.nbArcs+1, ϕ_zλ, ϕ_z1, ϕ_z2 )
                    push!(τ,♇)
                end
            end
        end
        #If the weight of the arc isn't 0 (-> if we picked the item)
        if weight(s, t) != 0
            ϕ[s.layer] = true #We add s.i to the list of items picked
        end
        t = s
    end
    ϕ
end

in_triangle(y::Tuple{Int,Int}, yr::Tuple{Int,Int}, ys::Tuple{Int,Int}) = yr[1] <= y[1] <= ys[1] && ys[2] <= y[2] <= yr[2]
in_triangle(y::Tuple{Int,Int}, xr::solution, xs::solution) = in_triangle(y, obj(xr), obj(xs))
in_triangle(xt::solution, xr::solution, xs::solution) = in_triangle(obj(xt), obj(xr), obj(xs))
in_triangle(y::solution, Δ::Triangle) = in_triangle(y, Δ.XΔ[1], Δ.XΔ[end])
in_triangle(y::Tuple{Int,Int}, Δ::Triangle) = in_triangle(y, Δ.XΔ[1], Δ.XΔ[end])

function update_lb!(Δ::Triangle)
    # @assert issorted(Δ.XΔ, by = x -> obj_1(x))

    XΔ = Δ.XΔ
    λ1, λ2 = Δ.λ

    lb = λ1*obj_1(XΔ[1]) + λ2*obj_2(XΔ[1])

    for i = 2:length(XΔ)
        zλ = λ1*obj_1(XΔ[i]) + λ2*obj_2(XΔ[i])
        lb = min(lb,zλ)
    end

    for i = 1:length(XΔ)-1
        xr = XΔ[i]
        xs = XΔ[i+1]
        zλ = λ1*(obj_1(xr)+1) + λ2*(obj_2(xs)+1)
        lb = min(lb,zλ)
    end
    Δ.lb = lb
end

function explore_triangle(Δ::Triangle, output::Bool)
    Δ.pending = false
    XΔ = Δ.XΔ
    OΔ = solution[] #Nondominated points located outside of Δ(yr, ys)

    GKP = build_graph(Δ, output)

    isempty(GKP) && return OΔ
    mono_pb = GKP.mono_pb

#    τ = binary_maxheap([Partition(v) for v in GKP.layer]) # updated 23-07-2021
    τ = BinaryMaxHeap([Partition(v) for v in GKP.layer])
    Tk = parent_partition!(τ)
    while Tk.zλ >= Δ.lb

        ϕk = const_Path_New_Partition(Tk, τ, Δ.lb, mono_pb)

        let yk = (Tk.z1, Tk.z2)

            if in_triangle(yk, Δ)#If the solution is in the triangle
                if !any(x -> dominates(obj(x), yk), XΔ) #and no other solution in the triangle dominates it
                    s = solution(ϕk, GKP.pb, GKP.mono_pb) #Create the solution( O(n) )
                    deleteat!(XΔ, findall(x -> s>=x, XΔ))#Delete solutions dominated by this one

                    #Add the solution to the list
                    indinsert = searchsortedfirst(XΔ, s, by = obj_1)
                    insert!(XΔ, indinsert, s)

                    if Δ.lb != update_lb!(Δ)#And update the lower bound
                        clean!(τ, Δ.lb) #If the lower bound has been updated, we can remove some partitions from τ
                    end
                end
            else
                if !any(y -> dominates(obj(y), yk), OΔ)#If the solution isn't in the triangle and isn't dominated by any other in OΔ
                    deleteat!(OΔ, findall(y -> dominates(yk, obj(y)), OΔ))#Delete solutions dominated by this one
                    push!(OΔ, solution(ϕk, GKP.pb, GKP.mono_pb))#push it in OΔ
                end
            end
        end

        length(τ)==0 && break #If there are no more partitions to explore, stop
        Tk = parent_partition!(τ)#Else, pick the one with the better zλ
    end
    OΔ
end

function second_phase(XSE::Vector{solution}, output::Bool)

    res = solution[]

    Δlist = [Triangle(XSE[i], XSE[i+1]) for i = 1:length(XSE)-1]

    for i = 1:length(Δlist)

        Δ = select_triangle(Δlist)
        output && println("exploring triangle $i/$(length(Δlist)) : Δ$(obj(Δ.xr))$(obj(Δ.xs)) ($((Δ.ub-Δ.lb)/2))")
        # output && plot_triangle(Δ)

        OΔ = explore_triangle(Δ, output)
        append!(res, Δ.XΔ)

        for sol in OΔ
            for j = 1:length(Δlist)
                Δj = Δlist[j]
                if in_triangle(sol,Δj)
                    let sol=sol #15938
                        !(Δj.pending) && break
                        any(x -> x>=sol, Δj.XΔ) && break

                        deleteat!(Δj.XΔ, findall(x -> sol>=x, Δj.XΔ))

                        indinsert = searchsortedfirst(Δj.XΔ, sol, by = obj_1)
                        insert!(Δj.XΔ, indinsert, sol)
                    end

                    update_lb!(Δj)
                    break
                end
            end
        end

    end
    sort!(res, by = obj_1, alg=QuickSort)
    return unique(elt -> elt.x, res)
end

function select_triangle(Δlist::Vector{Triangle})
    res = Δlist[findfirst(x->x.pending, Δlist)]
    val_min = res.ub - res.lb
    for i in findall(x -> x.pending, Δlist)
        val = Δlist[i].ub - Δlist[i].lb
        if val < val_min
            val_min = val
            res = Δlist[i]
        end
    end
    return res
end

max_stars_repo_name: vOptSolver/vOptSpecific

max_issues_repo_name: vOptSolver/vOptSpecific

max_forks_repo_name: vOptSolver/vOptSpecific


text:
# module PushInterface.DisplayInterface

import ColorTypes
import FixedPointNumbers
import Cairo: CairoSurface, CairoContext, CairoARGBSurface, read_from_png, set_source_surface, set_source_rgb, rectangle, fill, save, restore, width, height, paint
import Cairo: _jl_libcairo

mutable struct Canvas
    surface::CairoSurface
    ctx::CairoContext
    function Canvas()
        surface = CairoARGBSurface(DisplayPixelWidth, DisplayPixelHeight)
        ctx = CairoContext(surface)
        new(surface, ctx)
    end
end

function draw(f::Function, canvas::Canvas)
    save(canvas.ctx)
    f(canvas.ctx)
    restore(canvas.ctx)
end

import Base: fill!
function fill!(canvas::Canvas, color::ColorTypes.RGB{FixedPointNumbers.Normed{UInt8,8}})
    draw(canvas) do ctx
        set_source_rgb(ctx, color.r, color.g, color.b)
        rectangle(ctx, 0, 0, width(canvas.surface), height(canvas.surface))
        fill(ctx)
    end
end

function draw!(canvas::Canvas, path::String)
    image = read_from_png(path)
    draw(canvas) do ctx
        set_source_surface(ctx, image, 0, 0)
        paint(ctx)
    end
end

# module PushInterface.DisplayInterface

ext: jl


text:
function movesCounter!(r, list)
    push!(list, goToWall!(r, Sud))
    push!(list, goToWall!(r, West))
    while list[end] != 0 && list[end-1] != 0
        push!(list, goToWall!(r, Sud))
        push!(list, goToWall!(r, West))
    end
end

function goToWall!(r, side)
    k = 0
    while !isborder(r, side)
        k+=1
        move!(r, side)
    end
    return k
end

function countXY(list)
    return sum(list[2:2:end]), sum(list[1:2:end])
end

function rightSide!(r, x, y)
    sides = [Ost, Nord]
    xy = [x, y]
    result = [0, 0]

    for el = 1:2
        makeMove(r, sides[el], xy[el])
        putmarker!(r)
        while !isborder(r, sides[el])
            result[el]+=1
            move!(r, sides[el])
        end
    end

    return result[1], result[2]
end

function leftSide!(r, kx, ky)
    sides = [West, Sud]
    kxy = [kx, ky]

    for el = 1:2
        makeMove(r, sides[el], kxy[el])
        putmarker!(r)
        while !isborder(r, sides[el])
            move!(r, sides[el])
        end
    end
end

function goToHome!(r, list)
    reverse!(list)
    for n = 1:length(list)
        if n%2==1
            makeMove(r, Ost, list[n])
        else
            makeMove(r, Nord, list[n])
        end
    end
end

function makeMove(r, side, mov)
    for _ = 1:mov
        move!(r, side)
    end
end

function draw!(r)
    list = []

    movesCounter!(r, list)

    x, y = countXY(list)

    kx, ky = rightSide!(r, x, y)

    leftSide!(r, kx, ky)

    goToHome!(r, list)
end

max_stars_repo_name: Trilonka/MireaAl

max_issues_repo_name: Trilonka/MireaAl

max_forks_repo_name: Trilonka/MireaAl


text:
using LabelledArrays, Test, StaticArrays

@testset "Basic interface" begin
    ABC = @SLSliced (3,2) (:a,:b,:c), (:x, :y)
    b = ABC(1,2,3,4,5,6)

    @test b.a.x == 1
    @test b.b.x == 2
    @test b.c.y == 6
    @test b[1,1] == b.a.x
    @test b[2,2] == b.b.y
    @test b[3,1] == b.c.x
    @test_throws UndefVarError fill!(a,1)
    @test typeof(b.__x) == SArray{Tuple{3,2},Int,2,6}

    # Type stability tests
    ABC_fl = @SLSliced Float64 (3,2) (:a, :b, :c), (:x, :y)
    ABC_int = @SLSliced Int (3,2) (:a, :b, :c), (:x, :y)

    @test @inferred(similar_type(b, Float64)) === ABC_fl
    @test @inferred(similar_type(b, Float64, Size(3,2))) === @SLSliced Float64 (3,2) (:a, :b, :c), (:x, :y)
    @test @inferred(similar_type(b, Float64, Size(3,3))) === SArray{Tuple{3,3},Float64,2,9}

    @test typeof(@inferred(similar(b))) === LArray{Int,2,Array{Int,2},Tuple{(:a,:b,:c),(:x,:y)}}
    @test typeof(@inferred(similar(b, Float64))) === LArray{Float64,2,Array{Float64,2},Tuple{(:a,:b,:c),(:x,:y)}}
    @test typeof(@inferred(similar(b, Size(3,2)))) === LArray{Int,2,Array{Int,2},Tuple{(:a, :b, :c),(:x,:y)}}
    @test typeof(@inferred(similar(b, Size(3,3)))) === MArray{Tuple{3,3},Int,2,9}

    @test @inferred(copy(b)) === ABC_int(Tuple(b))

    @test @inferred(broadcast(Float64, b)) === ABC_fl(Tuple(b))
    @test @inferred(broadcast(+, b, b)) === ABC_int(Tuple(b.__x .+ b.__x))
    @test @inferred(broadcast(+, b, 1.0)) === ABC_fl(Tuple(b.__x .+ 1.0))

    @test @inferred(zero(b)) === ABC_int(zero(b))
end

@testset "accessing dimensions symbols of SLSliced" begin
    ABC = @SLSliced (3,2) (:a,:b,:c), (:x, :y)
    A = ABC([1 2; 3 4; 5 6])
    ret = symbols(A)
    @test ret == ((:a,:b,:c),(:x,:y))
end


text:
# ---
# title: 273. Integer to English Words
# id: problem273
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Math, String
# link: <https://leetcode.com/problems/integer-to-english-words/description/>
# hidden: true
# ---
#
# Convert a non-negative integer `num` to its English words representation.
#
#
#
# **Example 1:**
#
#
#
#     Input: num = 123
#     Output: "One Hundred Twenty Three"
#
#
# **Example 2:**
#
#
#
#     Input: num = 12345
#     Output: "Twelve Thousand Three Hundred Forty Five"
#
#
# **Example 3:**
#
#
#
#     Input: num = 1234567
#     Output: "One Million Two Hundred Thirty Four Thousand Five Hundred Sixty Seven"
#
#
# **Example 4:**
#
#
#
#     Input: num = 1234567891
#     Output: "One Billion Two Hundred Thirty Four Million Five Hundred Sixty Seven Thousand Eight Hundred Ninety One"
#
#
#
#
# **Constraints:**
#
#   * `0 <= num <= 231 - 1`
#
#
## @lc code=start
using LeetCode

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


text:
# This takes too long on Travis:
# Pkg.add("ArgParse")
# load_only=true
# include(Knet.dir("examples","charlm.jl"))
# CharLM.main("--gcheck 3")

using Knet
srand(42)

predict(w,x)=(w*x)

loss(w,x,y)=(sum((y-predict(w,x)).^2) / size(x,2))

lossgradient = grad(loss)

function train(w, data; lr=.02, epochs=10)
    for epoch=1:epochs
        for (x,y) in data
            g = lossgradient(w, x, y)
            w -= lr * g
        end
    end
    return w
end

function test(w, data)
    sumloss = numloss = 0
    for (x,y) in data
        sumloss += loss(w,x,y)
        numloss += 1
    end
    return sumloss/numloss
end

# Data generator:
import Base: start, next, done

type Data; w; batchsize; epochsize; noise; rng; atype; end

function Data(outputdims,inputdims; batchsize=20, epochsize=10000, noise=.01, rng=Base.GLOBAL_RNG, atype=Array)
    Data(convert(atype, randn(rng,outputdims,inputdims)),batchsize,epochsize,noise,rng,atype)
end

function next(l::Data, n)
    (outputdims, inputdims) = size(l.w)
    x = convert(l.atype, rand(l.rng, inputdims, l.batchsize))
    y = l.w * x + convert(l.atype, l.noise * randn(l.rng, outputdims, l.batchsize))
    return ((x,y), n+l.batchsize)
end

start(l::Data)=0
done(l::Data,n)=(n >= l.epochsize)

# Run:

data = Data(10,100)
w = 0.1*randn(10,100)
println((:epoch,0,:loss,test(w,data)))
@time for epoch=1:10
    w = train(w, data; epochs=1, lr=0.02)
    println((:epoch,epoch,:loss,test(w,data)))
end


text:
include("header.jl")

@testset "loss" begin
    for f in (logp, logsumexp)
        a = rand(10,10)
        @test gradcheck(f,a)
        @test gradcheck(f,a,kw=(:dims=>1,))
        @test gradcheck(f,a,kw=(:dims=>2,))
        if gpu() >= 0
            k = KnetArray(a)
            @test gradcheck(f,k)
            @test gradcheck(f,k,kw=(:dims=>1,))
            @test gradcheck(f,k,kw=(:dims=>2,))
            @test isapprox(f(a),f(k))
            @test isapprox(f(a,dims=1),f(k,dims=1))
            @test isapprox(f(a,dims=2),f(k,dims=2))
        end

        a = rand(10,10,10)
        @test gradcheck(f,a)
        @test gradcheck(f,a,kw=(:dims=>1,))
        @test gradcheck(f,a,kw=(:dims=>2,))
        @test gradcheck(f,a,kw=(:dims=>3,))
        @test gradcheck(f,a,kw=(:dims=>(1,2),))
        @test gradcheck(f,a,kw=(:dims=>(3,2),))
        @test gradcheck(f,a,kw=(:dims=>(1,3),))

        if gpu() >= 0
            k = KnetArray(a)
            @test gradcheck(f,k)
            @test isapprox(f(a),f(k))
            for d in [1,2,3]
                @test gradcheck(f,k,kw=(:dims=>d,))
                @test isapprox(f(a,dims=d),f(k,dims=d))
            end
            for dims in [(1,2), (1,3), (3,1), [1,3,2]]
                @test isapprox(f(a,dims=dims),f(k,dims=dims))
            end
        end

        a = rand(10,10, 10)
        for d in [1, 2, 3, (1,2), (1,3), [2,3], [1,2,3]]
            @test softmax(a, dims=d) ≈ exp.(logsoftmax(a, dims=d))
            @test all(sum(softmax(a, dims=d), dims=d) .≈ 1)
            if gpu() > 0
                k = KnetArray(a)
                @test softmax(k, dims=d) ≈ exp.(logsoftmax(k, dims=d))
                @test all(Array(sum(softmax(k, dims=d), dims=d)) .≈ 1)
            end
        end
    end

    a = rand(10,10)
    indices = rand(1:10,10)
    @test gradcheck(nll, a, indices, kw=(:dims=>1,), args=1)
    @test gradcheck(nll, a, indices, kw=(:dims=>2,), args=1)
    if gpu() >= 0
        k = KnetArray(a)
        @test gradcheck(nll, k, indices, kw=(:dims=>1,), args=1)
        @test gradcheck(nll, k, indices, kw=(:dims=>2,), args=1)
        @test isapprox(nll(k, indices, dims=1), nll(a, indices, dims=1))
        @test isapprox(nll(k, indices, dims=2), nll(a, indices, dims=2))
    end
    @test gradcheck(logistic,a[:],a[:])
    @test gradcheck(bce,a[:],a[:])

    # Issue 439: highorder derivatives
    using Knet: _softmax, _softback, _logp, _logpback, cudnnSoftmaxForward, cudnnSoftmaxBackward
    x = randn(3,4); y1 = _softmax(x,dims=1); y2 = _logp(x,dims=1); dy = randn(3,4)
    @test @gcheck _softmax(Param(x),dims=1)
    @test @gcheck _softback(Param(x),Param(y1),Param(dy),dims=1)
    @test @gcheck _logp(Param(x),dims=1)
    @test @gcheck _logpback(Param(x),Param(y2),Param(dy),dims=1)
    if gpu() >= 0
        x = KnetArray(x); y1 = KnetArray(y1); y2 = KnetArray(y2); dy = KnetArray(dy)
        @test isapprox(_softmax(x,dims=1), cudnnSoftmaxForward(x,algo=0))
        @test isapprox(_softmax(x,dims=1), cudnnSoftmaxForward(x,algo=1))
        @test isapprox(_logp(x,dims=1), cudnnSoftmaxForward(x,algo=2))
        @test isapprox(_softback(x,y1,dy,dims=1), cudnnSoftmaxBackward(y1,dy,algo=0))
        @test isapprox(_softback(x,y1,dy,dims=1), cudnnSoftmaxBackward(y1,dy,algo=1))
        @test isapprox(_logpback(x,y2,dy,dims=1), cudnnSoftmaxBackward(y2,dy,algo=2))
        @test @gcheck cudnnSoftmaxForward(Param(x),algo=0)
        @test @gcheck cudnnSoftmaxForward(Param(x),algo=1)
        @test @gcheck cudnnSoftmaxForward(Param(x),algo=2)
        @test @gcheck cudnnSoftmaxBackward(Param(y1),Param(dy),algo=0)
        @test @gcheck cudnnSoftmaxBackward(Param(y1),Param(dy),algo=1)
        @test @gcheck cudnnSoftmaxBackward(Param(y2),Param(dy),algo=2)

        # Broken example from Alkan's notebook:
        f(w,x,y) = nll(w*x,y)
        ∇f = grad(f)
        ∇fj(w,x,y,j) = ∇f(w,x,y)[j]
        ∇∇fj = grad(∇fj)
        a = rand(10,10); b = rand(10,10); c = rand(1:10,10)
        A = KnetArray(a); B = KnetArray(b); C = c
        @test isapprox(f(a,b,c), f(A,B,C))
        @test isapprox(∇f(a,b,c), ∇f(A,B,C))
        i = 10; j = 20; d = 1e-4
        @test isapprox(∇∇fj(a,b,c,i), ∇∇fj(A,B,C,i))
    end
end


text:
"""
    function dinic_impl(residual_graph, source, target, capacity_matrix)

Compute the maximum flow between the `source` and `target` for `residual_graph`
with edge flow capacities in `capacity_matrix` using
[Dinic\'s Algorithm](https://en.wikipedia.org/wiki/Dinic%27s_algorithm).
Return the value of the maximum flow as well as the final flow matrix.
"""
function dinic_impl end
@traitfn function dinic_impl(
        residual_graph::::lg.IsDirected,               # the input graph
        source::Integer,                       # the source vertex
        target::Integer,                       # the target vertex
        capacity_matrix::AbstractMatrix{T}    # edge flow capacities
    ) where {T}
    n = lg.nv(residual_graph)                     # number of vertexes
    flow_matrix = zeros(T, n, n)           # initialize flow matrix
    P = zeros(Int, n)                      # Sharable parent vector

    flow = 0

    while true
        augment = blocking_flow!(residual_graph, source, target, capacity_matrix, flow_matrix, P)
        augment == 0 && break
        flow += augment
    end
    return flow, flow_matrix
end

"""
    blocking_flow!(residual_graph, source, target, capacity_matrix, flow-matrix, P)

Like `blocking_flow`, but requires a preallocated parent vector `P`.
"""
function blocking_flow! end
@traitfn function blocking_flow!(
        residual_graph::::lg.IsDirected,     # the input graph
        source::Integer,                     # the source vertex
        target::Integer,                     # the target vertex
        capacity_matrix::AbstractMatrix{T},  # edge flow capacities
        flow_matrix::AbstractMatrix,         # the current flow matrix
        P::AbstractVector{Int}               # Parent vector to store Level Graph
    ) where {T}
    n = lg.nv(residual_graph)                # number of vertexes
    fill!(P, -1)
    P[source] = -2

    Q = [source]
    sizehint!(Q, n)

    while length(Q) > 0                   # Construct the Level Graph using BFS
        u = pop!(Q)
        for v in lg.outneighbors(residual_graph, u)
            if P[v] == -1 && capacity_matrix[u, v] > flow_matrix[u, v]
                P[v] = u
                pushfirst!(Q, v)
            end
        end
    end

    P[target] == -1 && return 0                    # BFS couldn't reach the target

    total_flow = 0

    for bv in lg.inneighbors(residual_graph, target)    # Trace all possible routes to source
        flow = typemax(T)
        v = target
        u = bv
        while v != source
            if u == -1                    # Vertex unreachable from source
                flow = 0
                break
            else
                flow = min(flow, capacity_matrix[u, v] - flow_matrix[u, v])
                v = u
                u = P[u]
            end
        end

        flow == 0 && continue                      # Flow cannot be augmented along path

        v = target
        u = bv
        while v != source             # Augment flow along path
            flow_matrix[u, v] += flow
            flow_matrix[v, u] -= flow
            v = u
            u = P[u]
        end

        total_flow += flow
    end
    return total_flow
end

"""
    blocking_flow(residual_graph, source, target, capacity_matrix, flow-matrix)

Use BFS to identify a blocking flow in the `residual_graph` with current flow
matrix `flow_matrix`and then backtrack from `target` to `source`,
augmenting flow along all possible paths.
"""
blocking_flow(
    residual_graph::lg.AbstractGraph,  # the input graph
    source::Integer,                   # the source vertex
    target::Integer,                   # the target vertex
    capacity_matrix::AbstractMatrix,   # edge flow capacities
    flow_matrix::AbstractMatrix,       # the current flow matrix
    ) = blocking_flow!(
            residual_graph,
            source,
            target,
            capacity_matrix,
            flow_matrix,
            zeros(Int, lg.nv(residual_graph)))


text:
#################################################################################################
# AIM 4: signal-to-noise ratio study of nanopolish model file
#################################################################################################
## Deps
using Pkg
Pkg.activate("./") # <=
using CpelNano
using StatsBase
using StatsPlots
using Distributed
using LinearAlgebra
using Distributions
using DelimitedFiles

## Constants

# IO
const aim_dir = "/Users/jordiabante/OneDrive - Johns Hopkins/CpelNano/Data/Simulations/Aim-4/" # <=

## Default attributes
default(titlefont=(14, "arial"),guidefont=(16, "arial"),tickfont=(12, "arial"))

#################################################################################################
# Calls
#################################################################################################

# IO
dataDir = "/Users/jordiabante/OneDrive - Johns Hopkins/CpelNano/Data/Real-Data/Nanopolish-Pore-Model/"
nano = "$(dataDir)/r9.4_450bps.cpg.6mer.template.model"

# Get nanopolish model dictionary
nano_dic = CpelNano.get_nano_pore_mod(nano)

# Plot histogram of Δμ
Δμs = [nano_dic[key][1] - nano_dic[key][3] for key in keys(nano_dic)]
Δμ_range = -5.0:0.2:5.0
Δμ_histo = fit(Histogram, Δμs, Δμ_range)
Δμ_histo = normalize(Δμ_histo, mode=:probability)
sts_Δμ = collect(Δμ_range)
cent_Δμ = sts_Δμ .+ step(Δμ_range) / 2.0
p_Δμ = plot(cent_Δμ, Δμ_histo.weights * 100, seriestype=:bar, alpha=0.25, xlabel="\\Delta\\mu", ylabel="Percentage (%)", label="");

# Plot histogram of σ̄
σ̄s = [0.5 * (nano_dic[key][2] + nano_dic[key][4]) for key in keys(nano_dic)]
σ̄_range = 0.0:0.1:5.0
σ̄_histo = fit(Histogram, σ̄s, σ̄_range)
σ̄_histo = normalize(σ̄_histo, mode=:probability)
sts_σ̄ = collect(σ̄_range)
cent_σ̄ = sts_σ̄ .+ step(σ̄_range) / 2.0
p_σ̄ = plot(cent_σ̄, σ̄_histo.weights * 100, seriestype=:bar, alpha=0.25, xlabel="Average \\sigma", ylabel="Percentage (%)", label="");

# Plot histogram of SNR
snrs = abs.(Δμs ./ σ̄s)
snr_range = 0.0:0.1:5.0
snr_histo = fit(Histogram, snrs, snr_range)
snr_histo = normalize(snr_histo, mode=:probability)
sts_snr = collect(snr_range)
cent_snr = sts_snr .+ step(snr_range) / 2.0
p_snr = plot(cent_snr, snr_histo.weights * 100, seriestype=:bar, alpha=0.25, xlabel="SNR", ylabel="Percentage (%)", label="");

# Plot all of them and save
plot(p_Δμ,p_σ̄,p_snr,layout=(3, 1),size=(600, 800))
savefig("$(aim_dir)/Histogram-Nanopolish-r9-model-file.pdf")

# Median SNR
snr_med = median(snrs)

# Δμ for simulations
Δμ_sims_1 = snr_med * 3.0 / sqrt(4000 / 450)
Δμ_sims_2 = median(Δμs)


text:
module MinimalRides

using Dates
using Geodesy
using TimeZones
using ZipFile
using ProgressMeter
using CSV
using JSON3
using StructTypes

include("animals.jl")

########################################
# JSON patch
########################################

struct JSONInfo
    id::String
    car_id::String
    category::String
    car_passengers::Int
    cat_carry_weight::Int
    start_dt::String
    end_dt::String
end

struct JSONRide
    info::JSONInfo
    data::Vector{Vector{Float64}}
end
StructTypes.StructType(::Type{JSONInfo}) = StructTypes.Struct()
StructTypes.StructType(::Type{JSONRide}) = StructTypes.Struct()

########################################
# Structures
########################################

struct Location{T}
    ts::Int
    coord::T
end

Location{LatLon{T}}(v::AbstractVector) where T = Location{LatLon{Float64}}(v[1], LatLon{Float64}(v[3], v[2]))
Location(v::AbstractVector{T}) where T = Location{LatLon{T}}(v[1], LatLon{Float64}(v[3], v[2]))
Location{ECEF}(v::AbstractVector{T}, datum = wgs84) where T = Location{ECEF{T}}(v[1], ECEF(LLA(LatLon(v[3], v[2])), datum))
Location{ECEF{T}}(v::AbstractVector, datum = wgs84) where T = Location{ECEF{T}}(v[1], ECEF(LLA(LatLon(v[3], v[2])), datum))
Location{ECEF}(l::Location{LatLon{T}}, datum = wgs84) where T = Location{ECEF{T}}(l.ts, ECEF(LLA(l.coord), datum))

const GEO_RADIUS = 6371.009

function great_circle(p1::LatLon, p2::LatLon)
    lat1, lng1 = deg2rad(p1.lat), deg2rad(p1.lon)
    lat2, lng2 = deg2rad(p2.lat), deg2rad(p2.lon)

    sin_lat1, cos_lat1 = sin(lat1), cos(lat1)
    sin_lat2, cos_lat2 = sin(lat2), cos(lat2)

    delta_lng = lng2 - lng1
    cos_delta_lng, sin_delta_lng = cos(delta_lng), sin(delta_lng)

    d = atan(sqrt((cos_lat2 * sin_delta_lng) ^ 2 + (cos_lat1 * sin_lat2 - sin_lat1 * cos_lat2 * cos_delta_lng) ^ 2), sin_lat1 * sin_lat2 + cos_lat1 * cos_lat2 * cos_delta_lng)

    return GEO_RADIUS * d
end

# this definition twice as fast as native Geodesy version `distance(l1.coord, l2.coord)`
distance(l1::Location{T}, l2::Location{T}) where {T <: LatLon} = great_circle(l1.coord, l2.coord)
distance(l1::Location{T}, l2::Location{T}) where {T <: ECEF} = @inbounds (l1.coord[1] - l2.coord[1])^2 + (l1.coord[2] - l2.coord[2])^2 + (l1.coord[3] - l2.coord[3])^2
# distance(l1::Location, l2::Location) = distance(l1.coord, l2.coord)

struct Info
    id::String
    car_id::String
    category::String
    car_passengers::Int
    cat_carry_weight::Int
    car_type::Symbol
    start_ts::ZonedDateTime
    end_ts::ZonedDateTime
end

function car_type(car)
    res = :other
    res = car.car_passengers >= 8 ? :bus : res

    if (res == :other) & ((car.category == "Специальный\\Автобус ") | (occursin("Легковой", car.category)))
        res = :passenger
    end

    res = occursin("Грузовой", car.category) ? :freight : res
    res = car.category == "Строительный\\Автокран" ? :special : res
    # Weird logic, Специальный\\Вахтовая а/м accounted as :bus, not :special
    if (res == :other) & (occursin("Специальный", car.category))
        res = :special
    end

    return res
end

dateformatter() = Dates.DateFormat("yyyy-mm-dd HH:MM:SSzz")

function Info(car)
    Info(
        car.id,
        car.car_id,
        car.category,
        car.car_passengers,
        car.cat_carry_weight,
        car_type(car),
        ZonedDateTime(car.start_dt, dateformatter()),
        ZonedDateTime(car.end_dt, dateformatter())
    )
end

function Info(info::JSONInfo)
    Info(
        info.id,
        info.car_id,
        info.category,
        info.car_passengers,
        info.cat_carry_weight,
        car_type(info),
        ZonedDateTime(info.start_dt, dateformatter()),
        ZonedDateTime(info.end_dt, dateformatter())
    )
end

struct Ride{T}
    info::Info
    route::Vector{Location{T}}
end
Base.getindex(ride::Ride, i) = @inbounds ride.route[i]
Base.length(ride::Ride) = length(ride.route)
Base.lastindex(ride::Ride) = length(ride)

function Ride(ride)
    Ride{LatLon{Float64}}(
        Info(ride.info),
        sort!(Location{LatLon{Float64}}.(ride.data), by = x -> x.ts)
    )
end

function Ride{ECEF}(ride, datum = wgs84)
    Ride{ECEF{Float64}}(
        Info(ride.info),
        sort!(Location{ECEF{Float64}}.(ride.data, datum), by = x -> x.ts)
    )
end

function Ride{ECEF}(ride::Ride{LatLon{T}}, datum = wgs84) where T
    Ride{ECEF{T}}(ride.info, Location{ECEF}.(ride.route))
end

########################################
# Download utility functions
########################################

function load_raw(url, cache, force = false)
    if !isfile(cache) | force
        @info "Downloading rides from $(url)"
        download(url, cache)
    end

    res = JSONRide[]

    r = ZipFile.Reader(cache)
    @info "Processing $(length(r.files)) rides from $(cache)"
    @showprogress for f in r.files
        v = Vector{UInt8}(undef, f.uncompressedsize)
        read!(f, v)
        push!(res, JSON3.read(v, JSONRide))
    end
    close(r)

    return res
end

function load(url, cache, force = false)
    if !isfile(cache) | force
        @info "Downloading rides from $(url)"
        download(url, cache)
    end

    res = Ride[]
    r = ZipFile.Reader(cache)
    @info "Processing $(length(r.files)) rides from $(cache)"
    @showprogress for f in r.files
        v = Vector{UInt8}(undef, f.uncompressedsize)
        read!(f, v)
        push!(res, Ride(JSON3.read(v)))
    end
    close(r)
    return identity.(res)
end

########################################
# Simple Ride processing functions
########################################

function make_subset(rides, days, types)
    filter(x -> (Date(x.info.start_ts) in days) & (x.info.car_type in types), rides)
end

mileage(ride::Ride) = mileage(ride.route)
function mileage(route::AbstractVector{Location{LatLon{T}}}) where T
    res = zero(T)
    @inbounds @simd for i in 2:length(route)
        res += distance(route[i - 1], route[i])
    end

    return res
end

function mileage(route::AbstractVector{Location{ECEF{T}}}) where T
    res = zero(T)
    @inbounds @simd for i in 2:length(route)
        res += sqrt(distance(route[i - 1], route[i]))
    end

    return res/1000
end

########################################
# Smoothing
########################################

abstract type SmoothAlgorithm end

struct DistanceSmoothing <: SmoothAlgorithm
    r::Float64
end

function smooth(alg::DistanceSmoothing, ride::Ride{T}) where {T <: LatLon}
    r = alg.r
    segment = 0.
    route = [ride[1]]
    for i in 2:length(ride)
        segment += distance(ride[i - 1], ride[i])
        if segment >= r
            segment = 0.
            push!(route, ride[i])
        end
    end

    if ride[end] != route[end]
        push!(route, ride[end])
    end

    return Ride(ride.info, route)
end

function smooth(alg::DistanceSmoothing, ride::Ride{T}) where {T <: ECEF}
    r = alg.r*1000
    segment = 0.
    route = [ride[1]]
    for i in 2:length(ride)
        segment += sqrt(distance(ride[i - 1], ride[i]))
        if segment >= r
            segment = 0.
            push!(route, ride[i])
        end
    end

    if ride[end] != route[end]
        push!(route, ride[end])
    end

    return Ride(ride.info, route)
end

struct SegmentSmoothing <: SmoothAlgorithm
    n::Int
end

function smooth(alg::SegmentSmoothing, ride)
    route = ride[1:alg.n:length(ride)]

    if ride[end] != route[end]
        push!(route, ride[end])
    end

    return Ride(ride.info, route)
end

########################################
# Coverage and Results
########################################

function coverage(ride1::Ride{T}, ride2::Ride{T}, radius) where {T <: LatLon}
    dist1 = fill(Inf, length(ride1))
    dist2 = fill(Inf, length(ride2))

    @inbounds for i in 1:length(ride1)
        for j in 1:length(ride2)
            d = distance(ride1[i], ride2[j])
            dist1[i] = d < dist1[i] ? d : dist1[i]
            dist2[j] = d < dist2[j] ? d : dist2[j]
        end
    end

    cov1 = sum(<(radius), dist1)/length(ride1)
    cov2 = sum(<(radius), dist2)/length(ride2)

    return (;cov1, cov2)
end

function coverage(ride1::Ride{T}, ride2::Ride{T}, radius) where {T <: ECEF}
    dist1 = fill(Inf, length(ride1))
    dist2 = fill(Inf, length(ride2))
    radius = (radius * 1000.0)^2

    @inbounds for i in 1:length(ride1)
        for j in 1:length(ride2)
            d = distance(ride1[i], ride2[j])
            dist1[i] = d < dist1[i] ? d : dist1[i]
            dist2[j] = d < dist2[j] ? d : dist2[j]
        end
    end

    cov1 = sum(<(radius), dist1)/length(ride1)
    cov2 = sum(<(radius), dist2)/length(ride2)

    return (;cov1, cov2)
end

function results(rides; smoothalg = DistanceSmoothing(0.5), search_radius = 1.0)
    smooth_rides = smooth.(Ref(smoothalg), rides)
    mileages = mileage.(rides)
    # Too lazy to define new structure, NamedTuple is fine at this stage
    res = NamedTuple{(:i, :j, :cov1, :cov2, :cov, :len1, :len2, :len, :op),Tuple{Int64,Int64,Float64,Float64,Float64,Float64,Float64,Float64,Float64}}[]
    Threads.@threads for i in 1:length(rides)
        for j in i+1:length(rides)
            cov1, cov2 = coverage(smooth_rides[i], smooth_rides[j], search_radius)
            len1 = mileages[i]
            len2 = mileages[j]
            cov = cov1 + cov2
            len = len1 + len2
            op = (cov1 * len1 + cov2 * len2)/(len1 + len2)
            push!(res, (; i, j, cov1, cov2, cov, len1, len2, len, op))
        end
    end

    return res
end

end # module


text:
using Plots; pyplot(fmt = :png)

include("ReadGeant4Hits.jl")
include("SetupSSD.jl")
include("EventSimulation.jl")
include("WaveformAnalysis.jl")
include("DigitalFilters.jl")

@info "Start"

# Setup the simulation
simFilename = "save/hex_pixel_Si_ring_circle.h5"
geometryFile = "config_files/simple_si_pixel_ring_circle.json"
chargeDriftConfigFile = "config_files/drift_velocity_config_linear_si.json"
max_ref = 3

# Get either a saved simulation or a new one
sim = ReadSimulation(simFilename)
if isnothing(sim)
    sim = SetupDefaultSimulation(geometryFile)
    CalculateDefaultDetectorFields!(sim, chargeDriftConfigFile, max_ref)
    SaveSimulation(simFilename, sim)
end

SetChargeDriftModel!(sim, chargeDriftConfigFile)

geantFilename = "../data/e_500keV.root"
gdf = GetHitInformation(geantFilename)
#
offset0 = CartesianPoint{T}(0, 0, 0)
events0 = DriftGeant4Events(gdf, sim, offset0)
riseTimes0 = CollectRiseTimes(events0)

#
# offsetEdge = CartesianPoint{T}(3e-3, 0, 0)
# eventsEdge = DriftGeant4Events(gdf, sim, offsetEdge)
# riseTimesEdge = CollectRiseTimes(eventsEdge)
#
histogram(riseTimes0 / 1u"ns", bins = 50, alpha = 0.5, label="Center")
# histogram!(riseTimesEdge / 1u"ns", bins = 50, alpha= 0.5, label="Edge")
# xlabel!("Rise time [ns]")

@info "Stop"

max_stars_repo_name: rjtayl/NDP

max_issues_repo_name: rjtayl/NDP

max_forks_repo_name: rjtayl/NDP


text:
# export LD_LIBRARY_PATH=/usr/local/cuda/lib64:$LD_LIBRARY_PATH
using TensorFlowDiffEq
using Base.Test

using Plots; plotly()
using DiffEqBase, ParameterizedFunctions
using DiffEqProblemLibrary, DiffEqDevTools

# Toy problem 1

f = (t,u) -> -u/5 + exp(-t/5).*cos(t)
(::typeof(f))(::Type{Val{:analytic}},t,u0) =  exp(-t/5)*(u0 + sin(t))
prob = ODEProblem(f,Float32(0.0),(Float32(0.0),Float32(2.0)))
sol = solve(prob,odetf(),dt=0.02)

sol(0.5)
sol(0.5,idxs=1)
plot(sol,plot_analytic=true)

prob = ODEProblem(f,0.0,(0.0,2.0))
sol = solve(prob,odetf(),dt=0.02)

sol = solve(prob,odetf(hl_width=[128,64,32,64,128]),dt=0.02,maxiters=Int(1e5))

dts = 1./2.^(8:-1:2) #14->7 good plot
prob = ODEProblem(f,0.0,(0.0,1.0))
sim  = test_convergence(dts,prob,odetf(hl_width=[128,64,32,64,128]))
@test abs(sim.𝒪est[:l2]-1) < 0.2

# Standard Convergence Problem

dts = 1./2.^(14:-1:7) #14->7 good plot
sim_linear  = test_convergence(dts,prob_ode_linear,odetf(),maxiters=Int(1e5))
@test abs(sim_linear.𝒪est[:l2]-1) < 0.2

# Problem 2

function lotka_voltera_tf(t,u)
  du1 = 1.5 .* u[:,1] - 1.0 .* u[:,1].*u[:,2]
  du2 = -3 .* u[:,2] + u[:,1].*u[:,2]
  [du1 du2]
end
prob = ODEProblem(lotka_voltera_tf,Float32[1.0,1.0],(Float32(0.0),Float32(10.0)))
sol = solve(prob,odetf(hl_width=[256]),dt=0.01,maxiters=Int(1e4),progress_steps=500)

sol(0.5,idxs=1)
sol([0.5,0.52])
plot(sol)

sol = solve(prob,odetf(hl_width=[128,128]),dt=0.01,maxiters=Int(1e4),progress_steps=500)
plot(sol)

sol = solve(prob,odetf(hl_width=[128,64, 64]),dt=0.01,maxiters=Int(1e4),progress_steps=500)
plot(sol)

# Toy problem 3

function lorenz_tf(t,u)
    du1 = 10.0(u[:, 2]-u[:,1])
    du2 = u[:,1].*(28.0-u[:, 3]) - u[:, 2]
    du3 = u[:,1].*u[:,2] - (8/3)*u[:,3]

    [du1 du2 du3]
end
prob = ODEProblem(lorenz_tf,Float32[1.0,0.0,0.0],(Float32(0.0),Float32(100.0)))
sol = solve(prob,odetf(hl_width=256),dt=0.1,maxiters=Int(1e5),progress_steps=500)

plot(sol,vars=(1,2,3))

# Macro Usage

g = @ode_def LorenzExample begin
  dx = σ*(y-x)
  dy = x*(ρ-z) - y
  dz = x*y - β*z
end σ=10.0 ρ=28.0 β=(8/3)

# g.vector_ex and g.vector_ex_return are vectorized forms
# g(Val{:vec},t,u,du) and g(Val{:vec},t,u)


text:
"""
测试上VHS的kagome
"""

using Plots
using MARY_fRG.Fermi.Kagome
using MARY_fRG.FlowEquation

#阻止画图
ENV["GKSwstype"] = "100"

function run_tf()
    model = upperband_kagome_lattice(0.)
    Γ4 = TFGamma4(
        model, 6.0, 24, 223
    )
    println(length(Γ4.ltris))
    Γ4.V[1, 1, 1, 1, :, :, :] .+= get_kagome_U_mt(3.0, Γ4.patches[1])
    println(
        all(isfinite.(Γ4.V[1, 1, 1, 1, :, :, :]))
    )
    #plt = heatmap(Γ4.V[1, 1, 1, 1, :, :, 1])
    #png(plt, "Gamma41")
    #return
    lval = 0.
    lstep = 0.01
    for idx in 1:1:3501
        bubb_pp, bubb_fs, bubb_ex = all_bubble_tf_mt(Γ4, lval)
        dl = dl_tf_mt(Γ4, bubb_pp, bubb_fs, bubb_ex)
        Γ4.V .+= dl .* lstep
        if idx % 50 == 1
            plt = heatmap(Γ4.V[1, 1, 1, 1, :, :, 1])
            png(plt, "Gamma4"*string(idx))
        end
        lval += lstep
    end
end

run_tf()


text:
import Base: *, copyto!
import LinearAlgebra: mul!
using SparseArrays

abstract type AbstractDerivOperator{T,N} end

abstract type AbstractFiniteDiff{T,N} <: AbstractDerivOperator{T,N} end

struct PeriodicFD{T<:AbstractFloat,N,S} <: AbstractFiniteDiff{T,N}
    derivative_order        :: Int
    approximation_order     :: Int
    dx                      :: T
    len                     :: Int
    stencil_length          :: Int
    stencil_offset          :: Int
    stencil_coefs           :: S
end

struct FiniteDiffDeriv{T<:AbstractFloat,N,S1,S2,S3} <: AbstractFiniteDiff{T,N}
    derivative_order        :: Int
    approximation_order     :: Int
    dx                      :: T
    len                     :: Int
    stencil_length          :: Int
    stencil_offset          :: Int
    stencil_coefs           :: S1
    boundary_stencil_length :: Int
    # these are assumed to be ordered from the boundary inwards
    low_boundary_offsets    :: Vector{Int}
    low_boundary_coefs      :: S2
    high_boundary_offsets   :: Vector{Int}
    high_boundary_coefs     :: S3
end

struct SpectralDeriv{T<:Real,N,S} <: AbstractDerivOperator{T,N}
    derivative_order        :: Int
    len                     :: Int
    D                       :: S
end

struct CenteredDiff{N} end

function CenteredDiff{N}(derivative_order::Int,
                         approximation_order::Int, dx::T,
                         len::Int) where {T<:AbstractFloat,N}
    @assert approximation_order > 1 "approximation_order must be greater than 1."

    stencil_length = derivative_order + approximation_order - 1 +
        (derivative_order+approximation_order)%2

    stencil_offset = div(stencil_length-1,2)
    weights = calculate_weights(derivative_order, stencil_length-1, stencil_offset)

    stencil_coefs = (1/dx^derivative_order) .* weights

    PeriodicFD{T,N,typeof(stencil_coefs)}(derivative_order, approximation_order,
                                          dx, len, stencil_length,
                                          stencil_offset, stencil_coefs)
end

CenteredDiff(args...) = CenteredDiff{1}(args...)

struct EqualSizeStencilFD{N} end

"""
    EqualSizeStencilFD{N}(derivative_order, approximation_order, dx, len)

Finite difference operator where interior and boundary stencils have the same
length and same approximation order. Low (high) boundary operators will use as
many points to the left (right) as possible. For approximation_order > 2 there
is more than one boundary point, and each boundary point uses a different
stencil. Note that these are *not* truncation matched.
"""
function EqualSizeStencilFD{N}(derivative_order::Int,
                               approximation_order::Int, dx::T,
                               len::Int) where {T<:AbstractFloat,N}
    @assert approximation_order > 1 "approximation_order must be greater than 1."

    stencil_length = derivative_order + approximation_order - 1 +
        (derivative_order+approximation_order)%2
    stencil_offset = div(stencil_length-1,2)
    weights = calculate_weights(derivative_order, stencil_length-1, stencil_offset)

    stencil_coefs = (1/dx^derivative_order) .* weights

    boundary_stencil_length = stencil_length

    low_boundary_size  = div(approximation_order,2)
    high_boundary_size = div(approximation_order,2)

    # remember to order from the boundary inwards, which means the offset array
    # is increasing in the low boundary and decreasing in the high boundary
    low_boundary_offsets  = collect(0:low_boundary_size-1)
    high_boundary_offsets = collect(boundary_stencil_length-1:-1:high_boundary_size+1)

    low_weights = [calculate_weights(derivative_order, stencil_length-1, offset)
                   for offset in low_boundary_offsets]
    low_boundary_coefs = (1/dx^derivative_order) .* low_weights

    high_weights = [calculate_weights(derivative_order, stencil_length-1, offset)
                    for offset in high_boundary_offsets]
    high_boundary_coefs = (1/dx^derivative_order) .* high_weights

    FiniteDiffDeriv{T,N,typeof(stencil_coefs),typeof(low_boundary_coefs),
                    typeof(high_boundary_coefs)}(
                        derivative_order, approximation_order,
                        dx, len, stencil_length,
                        stencil_offset, stencil_coefs,
                        boundary_stencil_length, low_boundary_offsets,
                        low_boundary_coefs, high_boundary_offsets, high_boundary_coefs
                    )
end

EqualSizeStencilFD(args...) = EqualSizeStencilFD{1}(args...)

struct ChebDeriv{N} end

function ChebDeriv{N}(derivative_order::Int, xmin::T, xmax::T, len::Int) where {T<:Real,N}
    x, Dx, Dxx = cheb(xmin, xmax, len)

    if derivative_order == 1
        D = Dx
    elseif derivative_order == 2
        D = Dxx
    else
        error("derivative_order not implemented")
    end

    SpectralDeriv{T,N,typeof(D)}(derivative_order, len, D)
end

ChebDeriv(args...) = ChebDeriv{1}(args...)

struct FourierDeriv{N} end

function FourierDeriv{N}(derivative_order::Int, xsize::T, len::Int) where {T<:Real,N}
    Dx, Dxx = fourier(xsize, len)

    if derivative_order == 1
        D = Dx
    elseif derivative_order == 2
        D = Dxx
    else
        error("derivative_order not implemented")
    end

    SpectralDeriv{T,N,typeof(D)}(derivative_order, len, D)
end

FourierDeriv(args...) = FourierDeriv{1}(args...)

function Base.getindex(A::FiniteDiffDeriv{T}, i::Int, j::Int) where {T}
    s = A.stencil_offset + 1

    # interior points
    if s <= i <= (A.len-s+1)
        idx  = j - i + s
        coeffs = A.stencil_coefs
        ss     = s
        stencil_length = A.stencil_length
    elseif i < s
        coeffs = A.low_boundary_coefs[i]
        ss     = A.low_boundary_offsets[i] + 1
        idx    = j - i + ss
        stencil_length = A.boundary_stencil_length
    else # i > (A.len-s+1)
        ii     = A.len-i+1
        coeffs = A.high_boundary_coefs[ii]
        ss     = A.high_boundary_offsets[ii] + 1
        idx    = j - i + ss
        stencil_length = A.boundary_stencil_length
    end

    if idx < 1 || idx > stencil_length
        return zero(T)
    else
        return coeffs[idx]
    end
end

function Base.getindex(A::PeriodicFD{T}, i::Int, j::Int) where {T}
    N      = A.len
    coeffs = A.stencil_coefs
    s      = A.stencil_offset + 1

    if 1 <= j - i + s <= N
        idx  = j - i + s
    else
        # imposing periodicity
        idx  = mod1(j - i + s, N)
    end

    if idx < 1 || idx > A.stencil_length
        return zero(T)
    else
        return coeffs[idx]
    end
end

@inline Base.getindex(A::AbstractFiniteDiff, i::Int, ::Colon) =
    [A[i,j] for j in 1:A.len]

@inline Base.getindex(A::SpectralDeriv, i, j) = A.D[i,j]

# make it a callable struct, to compute derivatives at a given point
function (A::AbstractFiniteDiff{T,N})(f::AbstractArray{T,M},
                                      idx::Vararg{Int,M}) where {T<:Real,N,M}
    # make sure axis of differentiation is contained in the dimensions of f
    @assert N <= M

    s = A.stencil_offset + 1
    i = idx[N] # point where derivative will be taken (with respect to the N-axis)

    if s <= i <= (A.len-s+1)
        return _D_interior(A, f, idx)
    elseif i < s
        return _D_low_bdr(A, f, idx)
    else # i > (A.len-s+1)
        return _D_high_bdr(A, f, idx)
    end
end

# interior points
function _D_interior(A::AbstractFiniteDiff{T,N}, f::AbstractArray, idx) where {T<:Real,N}
    coeffs = A.stencil_coefs
    s      = A.stencil_offset + 1
    i      = idx[N] # point where derivative will be taken (with respect to the N-axis)

    sum_i  = zero(T)
    @fastmath @inbounds for aa in 1:A.stencil_length
        i_circ = i - (s - aa)
        I      = Base.setindex(idx, i_circ, N)
        sum_i += coeffs[aa] * f[I...]
    end

    sum_i
end

# boundary points
function _D_low_bdr(A::PeriodicFD{T,N}, f::AbstractArray, idx) where {T<:Real,N}
    coeffs = A.stencil_coefs
    s      = A.stencil_offset + 1
    i      = idx[N]

    sum_i  = zero(T)
    @fastmath @inbounds for aa in 1:A.stencil_length
        # imposing periodicity
        i_circ = mod1(i - (s-aa), A.len)
        I      = Base.setindex(idx, i_circ, N)
        sum_i += coeffs[aa] * f[I...]
    end

    sum_i
end
function _D_high_bdr(A::PeriodicFD{T,N}, f::AbstractArray, idx) where {T<:Real,N}
    coeffs = A.stencil_coefs
    s      = A.stencil_offset + 1
    i      = idx[N]

    sum_i  = zero(T)
    @fastmath @inbounds for aa in 1:A.stencil_length
        # imposing periodicity
        i_circ = mod1(i - (s-aa), A.len)
        I      = Base.setindex(idx, i_circ, N)
        sum_i += coeffs[aa] * f[I...]
    end

    sum_i
end

# low boundary points
function _D_low_bdr(A::FiniteDiffDeriv{T,N}, f::AbstractArray, idx) where {T<:Real,N}
    i      = idx[N]

    sum_i  = zero(T)
    @fastmath @inbounds for aa in 1:A.boundary_stencil_length
        coeffs = A.low_boundary_coefs[i]
        s      = A.low_boundary_offsets[i] + 1
        i_circ = i - (s - aa)
        I      = Base.setindex(idx, i_circ, N)
        sum_i += coeffs[aa] * f[I...]
    end

    sum_i
end

# high boundary points
function _D_high_bdr(A::FiniteDiffDeriv{T,N}, f::AbstractArray, idx) where {T<:Real,N}
    i      = idx[N]
    ii     = A.len-i+1

    sum_i  = zero(T)
    @fastmath @inbounds for aa in 1:A.boundary_stencil_length
        coeffs = A.high_boundary_coefs[ii]
        s      = A.high_boundary_offsets[ii] + 1
        i_circ = i - (s - aa)
        I      = Base.setindex(idx, i_circ, N)
        sum_i += coeffs[aa] * f[I...]
    end

    sum_i
end

function LinearAlgebra.mul!(df::AbstractVector, A::AbstractFiniteDiff, f::AbstractVector)
    @fastmath @inbounds for idx in eachindex(f)
        df[idx] = A(f,idx)
    end
    nothing
end

function LinearAlgebra.mul!(df::AbstractArray, A::AbstractFiniteDiff, f::AbstractArray)
    @fastmath @inbounds for idx in CartesianIndices(f)
        df[idx] = A(f,idx.I...)
    end
    nothing
end

# mul! done by standard matrix multiplication for Chebyshev differentiation
# matrices. This can also be done through FFT, but after some testing the FFT
# route actually seemed to be performing slower, so let's stick with this
function LinearAlgebra.mul!(df::AbstractVector, A::SpectralDeriv, f::AbstractVector)
    mul!(df, A.D, f)
end

# and now for Arrays
function LinearAlgebra.mul!(df::AbstractArray{T}, A::SpectralDeriv{T,N,S},
                            f::AbstractArray{T}) where {T,N,S}
    Rpre  = CartesianIndices(axes(f)[1:N-1])
    Rpost = CartesianIndices(axes(f)[N+1:end])

    _mul_loop!(df, A, f, Rpre, size(f,N), Rpost)
end

# this works as follows. suppose we are taking the derivative of a 4-dimensional
# array g, with coordinates w,x,y,z, and size
#
#   size(g) = (nw, nx, ny, nz)
#
# let's further suppose that we want the derivative along the x-direction (here,
# the second entry). this fixes
#
#   N = 2
#
# and then
#
#   Rpre  = CartesianIndices((nw,))
#   Rpost = CartesianIndices((ny,nz))
#
# now, the entry [Ipre,:,Ipost]
#
# slices *only* along the x-direction, for fixed Ipre and Ipost. Ipre and Ipost
# loop along Rpre and Rpost (respectively) without touching the x-direction. we
# further take care to loop in memory-order, since julia's arrays are faster in
# the first dimension.
#
# adapted from http://julialang.org/blog/2016/02/iteration
#
@noinline function _mul_loop!(df::AbstractArray{T}, A, f, Rpre, n, Rpost) where {T}
    @fastmath @inbounds for Ipost in Rpost
        @inbounds for i in 1:n
            @inbounds for Ipre in Rpre
                sum_i = zero(T)
                @inbounds @simd for j in 1:n
                    sum_i += A.D[i,j] * f[Ipre,j,Ipost]
                end
                df[Ipre,i,Ipost] = sum_i
            end
        end
    end
    nothing
end

function *(A::AbstractDerivOperator, x::AbstractArray)
    y = similar(x)
    mul!(y, A, x)
    y
end

# make SpectralDeriv a callable struct, to compute derivatives only at a given point.
function (A::SpectralDeriv{T,N,S})(f::AbstractArray{T,M},
                                   idx::Vararg{Int,M}) where {T<:Real,N,M,S}
    # make sure axis of differentiation is contained in the dimensions of f
    @assert N <= M

    i  = idx[N] # point where derivative will be taken (with respect to the N-axis)

    sum_i = zero(T)
    @fastmath @inbounds for ii in 1:A.len
        I = Base.setindex(idx, ii, N)
        sum_i += A.D[i,ii] * f[I...]
    end
    sum_i
end

# now for cross-derivatives. we assume that A acts on the first and B on the
# second axis of the x Matrix.
function (A::AbstractFiniteDiff{T,N1})(B::AbstractFiniteDiff{T,N2},
                                       f::AbstractArray{T,M},
                                       idx::Vararg{Int,M}) where {T<:Real,N1,N2,M}
    NA   = A.len
    NB   = B.len
    sA   = A.stencil_offset + 1
    sB   = B.stencil_offset + 1

    # make sure axes of differentiation are contained in the dimensions of f
    @assert N1 < N2 <= M

    # points where derivative will be taken (along their respective axes)
    i  = idx[N1]
    j  = idx[N2]

    if sA <= i <= NA-sA+1 && sB <= j <= NB-sB+1
        return _D_interior(A, B, f, idx)

    elseif sA <= i <= NA-sA+1 && j < sB
        return _D_intA_lowB(A, B, f, idx)

    elseif sA <= i <= NA-sA+1 && j > NB-sB+1
        return _D_intA_highB(A, B, f, idx)

    elseif i < sA && sB <= j <= NB-sB+1
        return _D_lowA_intB(A, B, f, idx)

    elseif i > NA-sA+1 && sB <= j <= NB-sB+1
        return _D_highA_intB(A, B, f, idx)

    elseif i < sA && j < sB
        return _D_lowA_lowB(A, B, f, idx)

    elseif i < sA && j > NB-sB+1
        return _D_lowA_highB(A, B, f, idx)

    elseif i > NA-sA+1 && j < sB
        return _D_highA_lowB(A, B, f, idx)

    elseif i > NA-sA+1 && j > NB-sB+1
        return _D_highA_highB(A, B, f, idx)

    else
        error("I shouldn't be here...")
    end

end

# interior points
function _D_interior(A::AbstractFiniteDiff{T,N1}, B::AbstractFiniteDiff{T,N2},
                     f::AbstractArray, idx) where {T<:Real,N1,N2}
    NA   = A.len
    NB   = B.len
    qA   = A.stencil_coefs
    qB   = B.stencil_coefs
    sA   = A.stencil_offset + 1
    sB   = B.stencil_offset + 1

    # points where derivative will be taken (along their respective axes)
    i  = idx[N1]
    j  = idx[N2]

    sum_ij = zero(T)
    @fastmath @inbounds for jj in 1:B.stencil_length
        j_circ = j - (sB-jj)
        Itmp   = Base.setindex(idx, j_circ, N2)
        sum_i  = zero(T)
        @inbounds for ii in 1:A.stencil_length
            i_circ = i - (sA-ii)
            I      = Base.setindex(Itmp, i_circ, N1)
            sum_i += qA[ii] * qB[jj] * f[I...]
        end
        sum_ij += sum_i
    end

    sum_ij
end

# boundary points
function _D_intA_lowB(A::AbstractFiniteDiff{T,N1}, B::PeriodicFD{T,N2},
                      f::AbstractArray, idx) where {T<:Real,N1,N2}
    NA   = A.len
    NB   = B.len
    qA   = A.stencil_coefs
    qB   = B.stencil_coefs
    sA   = A.stencil_offset + 1
    sB   = B.stencil_offset + 1

    # points where derivative will be taken (along their respective axes)
    i  = idx[N1]
    j  = idx[N2]

    sum_ij = zero(T)
    @fastmath @inbounds for jj in 1:B.stencil_length
        # imposing periodicity
        j_circ = mod1(j - (sB-jj), NB)
        Itmp   = Base.setindex(idx, j_circ, N2)
        sum_i  = zero(T)
        @inbounds for ii in 1:A.stencil_length
            i_circ = i - (sA-ii)
            I      = Base.setindex(Itmp, i_circ, N1)
            sum_i += qA[ii] * qB[jj] * f[I...]
        end
        sum_ij += sum_i
    end

    sum_ij
end
function _D_intA_highB(A::AbstractFiniteDiff{T,N1}, B::PeriodicFD{T,N2},
                       f::AbstractArray, idx) where {T<:Real,N1,N2}
    NA   = A.len
    NB   = B.len
    qA   = A.stencil_coefs
    qB   = B.stencil_coefs
    sA   = A.stencil_offset + 1
    sB   = B.stencil_offset + 1

    # points where derivative will be taken (along their respective axes)
    i  = idx[N1]
    j  = idx[N2]

    sum_ij = zero(T)
    @fastmath @inbounds for jj in 1:B.stencil_length
        # imposing periodicity
        j_circ = mod1(j - (sB-jj), NB)
        Itmp   = Base.setindex(idx, j_circ, N2)
        sum_i  = zero(T)
        @inbounds for ii in 1:A.stencil_length
            i_circ = i - (sA-ii)
            I      = Base.setindex(Itmp, i_circ, N1)
            sum_i += qA[ii] * qB[jj] * f[I...]
        end
        sum_ij += sum_i
    end

    sum_ij
end
function _D_lowA_intB(A::PeriodicFD{T,N1}, B::AbstractFiniteDiff{T,N2},
                      f::AbstractArray, idx) where {T<:Real,N1,N2}
    NA   = A.len
    NB   = B.len
    qA   = A.stencil_coefs
    qB   = B.stencil_coefs
    sA   = A.stencil_offset + 1
    sB   = B.stencil_offset + 1

    # points where derivative will be taken (along their respective axes)
    i  = idx[N1]
    j  = idx[N2]

    sum_ij = zero(T)
    @fastmath @inbounds for jj in 1:B.stencil_length
        j_circ = j - (sB-jj)
        Itmp   = Base.setindex(idx, j_circ, N2)
        sum_i  = zero(T)
        @inbounds for ii in 1:A.stencil_length
            # imposing periodicity
            i_circ = mod1(i - (sA-ii), NA)
            I      = Base.setindex(Itmp, i_circ, N1)
            sum_i += qA[ii] * qB[jj] * f[I...]
        end
        sum_ij += sum_i
    end

    sum_ij
end
function _D_highA_intB(A::PeriodicFD{T,N1}, B::AbstractFiniteDiff{T,N2},
                       f::AbstractArray, idx) where {T<:Real,N1,N2}
    NA   = A.len
    NB   = B.len
    qA   = A.stencil_coefs
    qB   = B.stencil_coefs
    sA   = A.stencil_offset + 1
    sB   = B.stencil_offset + 1

    # points where derivative will be taken (along their respective axes)
    i  = idx[N1]
    j  = idx[N2]

    sum_ij = zero(T)
    @fastmath @inbounds for jj in 1:B.stencil_length
        j_circ = j - (sB-jj)
        Itmp   = Base.setindex(idx, j_circ, N2)
        sum_i  = zero(T)
        @inbounds for ii in 1:A.stencil_length
            # imposing periodicity
            i_circ = mod1(i - (sA-ii), NA)
            I      = Base.setindex(Itmp, i_circ, N1)
            sum_i += qA[ii] * qB[jj] * f[I...]
        end
        sum_ij += sum_i
    end

    sum_ij
end
function _D_lowA_lowB(A::PeriodicFD{T,N1}, B::PeriodicFD{T,N2},
                      f::AbstractArray, idx) where {T<:Real,N1,N2}
    NA   = A.len
    NB   = B.len
    qA   = A.stencil_coefs
    qB   = B.stencil_coefs
    sA   = A.stencil_offset + 1
    sB   = B.stencil_offset + 1

    # points where derivative will be taken (along their respective axes)
    i  = idx[N1]
    j  = idx[N2]

    sum_ij = zero(T)
    @fastmath @inbounds for jj in 1:B.stencil_length
        # imposing periodicity
        j_circ = mod1(j - (sB-jj), NB)
        Itmp   = Base.setindex(idx, j_circ, N2)
        sum_i  = zero(T)
        @inbounds for ii in 1:A.stencil_length
            # imposing periodicity
            i_circ = mod1(i - (sA-ii), NA)
            I      = Base.setindex(Itmp, i_circ, N1)
            sum_i += qA[ii] * qB[jj] * f[I...]
        end
        sum_ij += sum_i
    end

    sum_ij
end
function _D_lowA_highB(A::PeriodicFD{T,N1}, B::PeriodicFD{T,N2},
                       f::AbstractArray, idx) where {T<:Real,N1,N2}
    NA   = A.len
    NB   = B.len
    qA   = A.stencil_coefs
    qB   = B.stencil_coefs
    sA   = A.stencil_offset + 1
    sB   = B.stencil_offset + 1

    # points where derivative will be taken (along their respective axes)
    i  = idx[N1]
    j  = idx[N2]

    sum_ij = zero(T)
    @fastmath @inbounds for jj in 1:B.stencil_length
        # imposing periodicity
        j_circ = mod1(j - (sB-jj), NB)
        Itmp   = Base.setindex(idx, j_circ, N2)
        sum_i  = zero(T)
        @inbounds for ii in 1:A.stencil_length
            # imposing periodicity
            i_circ = mod1(i - (sA-ii), NA)
            I      = Base.setindex(Itmp, i_circ, N1)
            sum_i += qA[ii] * qB[jj] * f[I...]
        end
        sum_ij += sum_i
    end

    sum_ij
end
function _D_highA_lowB(A::PeriodicFD{T,N1}, B::PeriodicFD{T,N2},
                       f::AbstractArray, idx) where {T<:Real,N1,N2}
    NA   = A.len
    NB   = B.len
    qA   = A.stencil_coefs
    qB   = B.stencil_coefs
    sA   = A.stencil_offset + 1
    sB   = B.stencil_offset + 1

    # points where derivative will be taken (along their respective axes)
    i  = idx[N1]
    j  = idx[N2]

    sum_ij = zero(T)
    @fastmath @inbounds for jj in 1:B.stencil_length
        # imposing periodicity
        j_circ = mod1(j - (sB-jj), NB)
        Itmp   = Base.setindex(idx, j_circ, N2)
        sum_i  = zero(T)
        @inbounds for ii in 1:A.stencil_length
            # imposing periodicity
            i_circ = mod1(i - (sA-ii), NA)
            I      = Base.setindex(Itmp, i_circ, N1)
            sum_i += qA[ii] * qB[jj] * f[I...]
        end
        sum_ij += sum_i
    end

    sum_ij
end
function _D_highA_highB(A::PeriodicFD{T,N1}, B::PeriodicFD{T,N2},
                        f::AbstractArray, idx) where {T<:Real,N1,N2}
    NA   = A.len
    NB   = B.len
    qA   = A.stencil_coefs
    qB   = B.stencil_coefs
    sA   = A.stencil_offset + 1
    sB   = B.stencil_offset + 1

    # points where derivative will be taken (along their respective axes)
    i  = idx[N1]
    j  = idx[N2]

    sum_ij = zero(T)
    @fastmath @inbounds for jj in 1:B.stencil_length
        # imposing periodicity
        j_circ = mod1(j - (sB-jj), NB)
        Itmp   = Base.setindex(idx, j_circ, N2)
        sum_i  = zero(T)
        @inbounds for ii in 1:A.stencil_length
            # imposing periodicity
            i_circ = mod1(i - (sA-ii), NA)
            I      = Base.setindex(Itmp, i_circ, N1)
            sum_i += qA[ii] * qB[jj] * f[I...]
        end
        sum_ij += sum_i
    end

    sum_ij
end

# now for FiniteDiffDeriv structs

function _D_intA_lowB(A::AbstractFiniteDiff{T,N1}, B::FiniteDiffDeriv{T,N2},
                      f::AbstractArray, idx) where {T<:Real,N1,N2}
    NA   = A.len
    NB   = B.len
    qA   = A.stencil_coefs
    sA   = A.stencil_offset + 1

    # points where derivative will be taken (along their respective axes)
    i  = idx[N1]
    j  = idx[N2]

    qB  = B.low_boundary_coefs[j]
    sB  = B.low_boundary_offsets[j] + 1

    sum_ij = zero(T)
    @fastmath @inbounds for jj in 1:B.boundary_stencil_length
        j_circ = j - (sB-jj)
        Itmp   = Base.setindex(idx, j_circ, N2)
        sum_i  = zero(T)
        @inbounds for ii in 1:A.stencil_length
            i_circ = i - (sA-ii)
            I      = Base.setindex(Itmp, i_circ, N1)
            sum_i += qA[ii] * qB[jj] * f[I...]
        end
        sum_ij += sum_i
    end

    sum_ij
end

function _D_intA_highB(A::AbstractFiniteDiff{T,N1}, B::FiniteDiffDeriv{T,N2},
                       f::AbstractArray, idx) where {T<:Real,N1,N2}
    NB   = B.len
    qA   = A.stencil_coefs
    sA   = A.stencil_offset + 1

    i     = idx[N1]
    j     = idx[N2]
    jhigh = NB-j+1

    qB    = B.high_boundary_coefs[jhigh]
    sB    = B.high_boundary_offsets[jhigh] + 1

    sum_ij = zero(T)
    @fastmath @inbounds for jj in 1:B.boundary_stencil_length
        j_circ = j - (sB-jj)
        Itmp   = Base.setindex(idx, j_circ, N2)
        sum_i  = zero(T)
        @inbounds for ii in 1:A.stencil_length
            i_circ = i - (sA-ii)
            I      = Base.setindex(Itmp, i_circ, N1)
            sum_i += qA[ii] * qB[jj] * f[I...]
        end
        sum_ij += sum_i
    end

    sum_ij
end

function _D_lowA_intB(A::FiniteDiffDeriv{T,N1}, B::AbstractFiniteDiff{T,N2},
                      f::AbstractArray, idx) where {T<:Real,N1,N2}
    qB = B.stencil_coefs
    sB = B.stencil_offset + 1

    i  = idx[N1]
    j  = idx[N2]

    qA = A.low_boundary_coefs[i]
    sA = A.low_boundary_offsets[i] + 1

    sum_ij = zero(T)
    @fastmath @inbounds for jj in 1:B.stencil_length
        j_circ = j - (sB-jj)
        Itmp   = Base.setindex(idx, j_circ, N2)
        sum_i  = zero(T)
        @inbounds for ii in 1:A.boundary_stencil_length
            i_circ = i - (sA-ii)
            I      = Base.setindex(Itmp, i_circ, N1)
            sum_i += qA[ii] * qB[jj] * f[I...]
        end
        sum_ij += sum_i
    end

    sum_ij
end

function _D_highA_intB(A::FiniteDiffDeriv{T,N1}, B::AbstractFiniteDiff{T,N2},
                       f::AbstractArray, idx) where {T<:Real,N1,N2}
    NA   = A.len
    qB   = B.stencil_coefs
    sB   = B.stencil_offset + 1

    i     = idx[N1]
    j     = idx[N2]
    ihigh = NA-i+1

    qA    = A.high_boundary_coefs[ihigh]
    sA    = A.high_boundary_offsets[ihigh] + 1

    sum_ij = zero(T)
    @fastmath @inbounds for jj in 1:B.stencil_length
        j_circ = j - (sB-jj)
        Itmp   = Base.setindex(idx, j_circ, N2)
        sum_i  = zero(T)
        @inbounds for ii in 1:A.boundary_stencil_length
            i_circ = i - (sA-ii)
            I      = Base.setindex(Itmp, i_circ, N1)
            sum_i += qA[ii] * qB[jj] * f[I...]
        end
        sum_ij += sum_i
    end

    sum_ij
end

function _D_lowA_lowB(A::FiniteDiffDeriv{T,N1}, B::FiniteDiffDeriv{T,N2},
                      f::AbstractArray, idx) where {T<:Real,N1,N2}
    i  = idx[N1]
    j  = idx[N2]

    qA = A.low_boundary_coefs[i]
    sA = A.low_boundary_offsets[i] + 1
    qB = B.low_boundary_coefs[j]
    sB = B.low_boundary_offsets[j] + 1

    sum_ij = zero(T)
    @fastmath @inbounds for jj in 1:B.boundary_stencil_length
        j_circ = j - (sB-jj)
        Itmp   = Base.setindex(idx, j_circ, N2)
        sum_i  = zero(T)
        @inbounds for ii in 1:A.boundary_stencil_length
            i_circ = i - (sA-ii)
            I      = Base.setindex(Itmp, i_circ, N1)
            sum_i += qA[ii] * qB[jj] * f[I...]
        end
        sum_ij += sum_i
    end

    sum_ij
end

function _D_lowA_highB(A::FiniteDiffDeriv{T,N1}, B::FiniteDiffDeriv{T,N2},
                       f::AbstractArray, idx) where {T<:Real,N1,N2}
    NB = B.len

    i     = idx[N1]
    j     = idx[N2]
    jhigh = NB-j+1

    qA    = A.low_boundary_coefs[i]
    sA    = A.low_boundary_offsets[i] + 1
    qB    = B.high_boundary_coefs[jhigh]
    sB    = B.high_boundary_offsets[jhigh] + 1

    sum_ij = zero(T)
    @fastmath @inbounds for jj in 1:B.boundary_stencil_length
        j_circ = j - (sB-jj)
        Itmp   = Base.setindex(idx, j_circ, N2)
        sum_i  = zero(T)
        @inbounds for ii in 1:A.boundary_stencil_length
            i_circ = i - (sA-ii)
            I      = Base.setindex(Itmp, i_circ, N1)
            sum_i += qA[ii] * qB[jj] * f[I...]
        end
        sum_ij += sum_i
    end

    sum_ij
end

function _D_highA_lowB(A::FiniteDiffDeriv{T,N1}, B::FiniteDiffDeriv{T,N2},
                       f::AbstractArray, idx) where {T<:Real,N1,N2}
    NA = A.len

    i     = idx[N1]
    j     = idx[N2]
    ihigh = NA-i+1

    qA    = A.high_boundary_coefs[ihigh]
    sA    = A.high_boundary_offsets[ihigh] + 1
    qB    = B.low_boundary_coefs[j]
    sB    = B.low_boundary_offsets[j] + 1

    sum_ij = zero(T)
    @fastmath @inbounds for jj in 1:B.boundary_stencil_length
        j_circ = j - (sB-jj)
        Itmp   = Base.setindex(idx, j_circ, N2)
        sum_i  = zero(T)
        @inbounds for ii in 1:A.boundary_stencil_length
            i_circ = i - (sA-ii)
            I      = Base.setindex(Itmp, i_circ, N1)
            sum_i += qA[ii] * qB[jj] * f[I...]
        end
        sum_ij += sum_i
    end

    sum_ij
end

function _D_highA_highB(A::FiniteDiffDeriv{T,N1}, B::FiniteDiffDeriv{T,N2},
                        f::AbstractArray, idx) where {T<:Real,N1,N2}
    NA = A.len
    NB = B.len

    i     = idx[N1]
    j     = idx[N2]
    ihigh = NA-i+1
    jhigh = NB-j+1

    qA    = A.high_boundary_coefs[ihigh]
    sA    = A.high_boundary_offsets[ihigh] + 1
    qB    = B.high_boundary_coefs[jhigh]
    sB    = B.high_boundary_offsets[jhigh] + 1

    sum_ij = zero(T)
    @fastmath @inbounds for jj in 1:B.boundary_stencil_length
        j_circ = j - (sB-jj)
        Itmp   = Base.setindex(idx, j_circ, N2)
        sum_i  = zero(T)
        @inbounds for ii in 1:A.boundary_stencil_length
            i_circ = i - (sA-ii)
            I      = Base.setindex(Itmp, i_circ, N1)
            sum_i += qA[ii] * qB[jj] * f[I...]
        end
        sum_ij += sum_i
    end

    sum_ij
end

# cross derivatives with FiniteDiffDeriv and PeriodicFD

function _D_lowA_lowB(A::PeriodicFD{T,N1}, B::FiniteDiffDeriv{T,N2},
                      f::AbstractArray, idx) where {T<:Real,N1,N2}
    NA = A.len
    qA = A.stencil_coefs
    sA = A.stencil_offset + 1

    i  = idx[N1]
    j  = idx[N2]

    qB = B.low_boundary_coefs[j]
    sB = B.low_boundary_offsets[j] + 1

    sum_ij = zero(T)
    @fastmath @inbounds for jj in 1:B.boundary_stencil_length
        j_circ = j - (sB-jj)
        Itmp   = Base.setindex(idx, j_circ, N2)
        sum_i  = zero(T)
        @inbounds for ii in 1:A.stencil_length
            # imposing periodicity
            i_circ = mod1(i - (sA-ii), NA)
            I      = Base.setindex(Itmp, i_circ, N1)
            sum_i += qA[ii] * qB[jj] * f[I...]
        end
        sum_ij += sum_i
    end

    sum_ij
end
function _D_lowA_lowB(A::FiniteDiffDeriv{T,N1}, B::PeriodicFD{T,N2},
                      f::AbstractArray, idx) where {T<:Real,N1,N2}
    NB = B.len
    qB = B.stencil_coefs
    sB = B.stencil_offset + 1

    i  = idx[N1]
    j  = idx[N2]

    qA = A.low_boundary_coefs[i]
    sA = A.low_boundary_offsets[i] + 1

    sum_ij = zero(T)
    @fastmath @inbounds for jj in 1:B.stencil_length
        # imposing periodicity
        j_circ = mod1(j - (sB-jj), NB)
        Itmp   = Base.setindex(idx, j_circ, N2)
        sum_i  = zero(T)
        @inbounds for ii in 1:A.boundary_stencil_length
            i_circ = i - (sA-ii)
            I      = Base.setindex(Itmp, i_circ, N1)
            sum_i += qA[ii] * qB[jj] * f[I...]
        end
        sum_ij += sum_i
    end

    sum_ij
end

function _D_lowA_highB(A::PeriodicFD{T,N1}, B::FiniteDiffDeriv{T,N2},
                       f::AbstractArray, idx) where {T<:Real,N1,N2}
    NA    = A.len
    qA    = A.stencil_coefs
    sA    = A.stencil_offset + 1

    i     = idx[N1]
    j     = idx[N2]
    jhigh = NB-j+1

    qB    = B.high_boundary_coefs[jhigh]
    sB    = B.high_boundary_offsets[jhigh] + 1

    sum_ij = zero(T)
    @fastmath @inbounds for jj in 1:B.boundary_stencil_length
        j_circ = j - (sB-jj)
        Itmp   = Base.setindex(idx, j_circ, N2)
        sum_i  = zero(T)
        @inbounds for ii in 1:A.stencil_length
            # imposing periodicity
            i_circ = mod1(i - (sA-ii), NA)
            I      = Base.setindex(Itmp, i_circ, N1)
            sum_i += qA[ii] * qB[jj] * f[I...]
        end
        sum_ij += sum_i
    end

    sum_ij
end
function _D_lowA_highB(A::FiniteDiffDeriv{T,N1}, B::PeriodicFD{T,N2},
                       f::AbstractArray, idx) where {T<:Real,N1,N2}
    NB    = B.len
    qB    = B.stencil_coefs
    sB    = B.stencil_offset + 1

    i     = idx[N1]
    j     = idx[N2]
    jhigh = NB-j+1

    qA    = A.low_boundary_coefs[i]
    sA    = A.low_boundary_offsets[i] + 1

    sum_ij = zero(T)
    @fastmath @inbounds for jj in 1:B.stencil_length
        # imposing periodicity
        j_circ = mod1(j - (sB-jj), NB)
        Itmp   = Base.setindex(idx, j_circ, N2)
        sum_i  = zero(T)
        @inbounds for ii in 1:A.boundary_stencil_length
            i_circ = i - (sA-ii)
            I      = Base.setindex(Itmp, i_circ, N1)
            sum_i += qA[ii] * qB[jj] * f[I...]
        end
        sum_ij += sum_i
    end

    sum_ij
end

function _D_highA_lowB(A::PeriodicFD{T,N1}, B::FiniteDiffDeriv{T,N2},
                       f::AbstractArray, idx) where {T<:Real,N1,N2}
    NA    = A.len
    qA    = A.stencil_coefs
    sA    = A.stencil_offset + 1

    i     = idx[N1]
    j     = idx[N2]
    ihigh = NA-i+1

    qB    = B.low_boundary_coefs[j]
    sB    = B.low_boundary_offsets[j] + 1

    sum_ij = zero(T)
    @fastmath @inbounds for jj in 1:B.boundary_stencil_length
        j_circ = j - (sB-jj)
        Itmp   = Base.setindex(idx, j_circ, N2)
        sum_i  = zero(T)
        @inbounds for ii in 1:A.stencil_length
            # imposing periodicity
            i_circ = mod1(i - (sA-ii), NA)
            I      = Base.setindex(Itmp, i_circ, N1)
            sum_i += qA[ii] * qB[jj] * f[I...]
        end
        sum_ij += sum_i
    end

    sum_ij
end
function _D_highA_lowB(A::FiniteDiffDeriv{T,N1}, B::PeriodicFD{T,N2},
                       f::AbstractArray, idx) where {T<:Real,N1,N2}
    NA    = A.len
    NB    = B.len
    qB    = B.stencil_coefs
    sB    = B.stencil_offset + 1

    i     = idx[N1]
    j     = idx[N2]
    ihigh = NA-i+1

    qA    = A.high_boundary_coefs[ihigh]
    sA    = A.high_boundary_offsets[ihigh] + 1

    sum_ij = zero(T)
    @fastmath @inbounds for jj in 1:B.stencil_length
        # imposing periodicity
        j_circ = mod1(j - (sB-jj), NB)
        Itmp   = Base.setindex(idx, j_circ, N2)
        sum_i  = zero(T)
        @inbounds for ii in 1:A.boundary_stencil_length
            i_circ = i - (sA-ii)
            I      = Base.setindex(Itmp, i_circ, N1)
            sum_i += qA[ii] * qB[jj] * f[I...]
        end
        sum_ij += sum_i
    end

    sum_ij
end

function _D_highA_highB(A::PeriodicFD{T,N1}, B::FiniteDiffDeriv{T,N2},
                        f::AbstractArray, idx) where {T<:Real,N1,N2}
    NA    = A.len
    NB    = B.len

    qA    = A.stencil_coefs
    sA    = A.stencil_offset + 1

    i     = idx[N1]
    j     = idx[N2]
    ihigh = NA-i+1
    jhigh = NB-j+1

    qB    = B.high_boundary_coefs[jhigh]
    sB    = B.high_boundary_offsets[jhigh] + 1

    sum_ij = zero(T)
    @fastmath @inbounds for jj in 1:B.boundary_stencil_length
        j_circ = j - (sB-jj)
        Itmp   = Base.setindex(idx, j_circ, N2)
        sum_i  = zero(T)
        @inbounds for ii in 1:A.stencil_length
            # imposing periodicity
            i_circ = mod1(i - (sA-ii), NA)
            I      = Base.setindex(Itmp, i_circ, N1)
            sum_i += qA[ii] * qB[jj] * f[I...]
        end
        sum_ij += sum_i
    end

    sum_ij
end
function _D_highA_highB(A::FiniteDiffDeriv{T,N1}, B::PeriodicFD{T,N2},
                        f::AbstractArray, idx) where {T<:Real,N1,N2}
    NA    = A.len
    NB    = B.len
    qB    = B.stencil_coefs
    sB    = B.stencil_offset + 1

    i     = idx[N1]
    j     = idx[N2]
    ihigh = NA-i+1
    jhigh = NB-j+1

    qA    = A.high_boundary_coefs[ihigh]
    sA    = A.high_boundary_offsets[ihigh] + 1

    sum_ij = zero(T)
    @fastmath @inbounds for jj in 1:B.stencil_length
        # imposing periodicity
        j_circ = mod1(j - (sB-jj), NB)
        Itmp   = Base.setindex(idx, j_circ, N2)
        sum_i  = zero(T)
        @inbounds for ii in 1:A.boundary_stencil_length
            i_circ = i - (sA-ii)
            I      = Base.setindex(Itmp, i_circ, N1)
            sum_i += qA[ii] * qB[jj] * f[I...]
        end
        sum_ij += sum_i
    end

    sum_ij
end

# Casting to matrix types

copyto!(M::AbstractMatrix{T}, A::SpectralDeriv) where {T<:Real} =
    copyto!(M, A.D)

function copyto!(M::AbstractMatrix{T}, A::AbstractFiniteDiff) where {T<:Real}
    for idx in CartesianIndices(M)
        M[idx] = A[idx.I...]
    end
    M
end

LinearAlgebra.Array(A::AbstractDerivOperator{T}) where {T} =
    copyto!(zeros(T, A.len, A.len), A)

SparseArrays.SparseMatrixCSC(A::AbstractFiniteDiff{T}) where {T} =
    copyto!(spzeros(T, A.len, A.len), A)


text:
set_theme!(
    rowgap = 5,
    colgap = 5,
    fontsize = 16,
    font =:sans,
    Axis = (
        xgridcolor = :black,
        ygridcolor = :black,
        backgroundcolor = :white,
        xgridstyle=:dash,
        ygridstyle=:dash,
        xgridwidth=0.15,
        ygridwidth=0.15,
        xtickalign=1,
        ytickalign=1, #xtickcolor = :red, ytickcolor = :white,
        xgridvisible = true,
        ygridvisible = true),
    Colorbar = (
        topspinevisible = false,
        rightspinevisible = false,
        bottomspinevisible = false,
        leftspinevisible = false,
        width = 12,
        height = Relative(4/4),
        tickalign = 1,
        labelpadding = -5),
    Legend = (
        tellheight = false,
        tellwidth = false,
        halign = :left,
        valign = :bottom,
        labelsize = 14,
        linewidth = 2,
        margin = (10, 10, 10, 10),
        framevisible = true)
)

max_stars_repo_name: aivodji/MakieNotebooks

max_issues_repo_name: aivodji/MakieNotebooks

max_forks_repo_name: aivodji/MakieNotebooks


text:
@testset "Forward Models: Convolution" begin

    N = 5
    psf = zeros((N, N))
    psf[1, 1] = 1
    img = randn((N, N))

    c(img, psf) = conv(img, psf, [1, 2])
    conv_temp = c
    @test conv_temp(img, psf) ≈ img
    s(img, psf) = sum(conv_temp(img, psf))

    @test all(1 .≈ Zygote.gradient(s, img, psf)[1])

end


text:
#!/usr/bin/env julia --project

# This file is the entry point for a generalized GCM driver that is used
# to run a variety of GCM experiments.
# Currently available default experiments are:
#   1. `baroclinicwave_problem`: initial value problem (no sources to
#      maintain equilibration)
#   2. `heldsuarez_problem`: runs to equilibration
#
# Experiment name is required to be specified in the command line:
#   e.g.: julia --project experiments/AtmosGCM/GCMDriver/jl --experiment=heldsuarez_problem
#
# The initial / boundary conditions of each experiment (defaults of which are
# defined in the experiment`_problem.jl` files) can be mixed and matched, as long
# as they are specified in:
#   - gcm_bcs.jl (use `--surface-flux`)
#   - gcm_perturbations.jl (use `--init-perturbation`)
#   - gcm_base_states.jl (use `--init-base-state`)
#   - gcm_moisture_profiles.jl (use `--init-moisture-profile`)
# Default sources cannot be currently overriden from command line and are defined in `gcm_sources.jl`

using ArgParse
using LinearAlgebra
using StaticArrays
using Test

using ClimateMachine
using ClimateMachine.Atmos
using ClimateMachine.Atmos: recover_thermo_state
using ClimateMachine.ConfigTypes
using ClimateMachine.Diagnostics
using ClimateMachine.GenericCallbacks
using ClimateMachine.Mesh.Filters
using ClimateMachine.Mesh.Grids
using ClimateMachine.ODESolvers
using ClimateMachine.Orientations
using ClimateMachine.SystemSolvers
using ClimateMachine.TemperatureProfiles
using ClimateMachine.Thermodynamics
using ClimateMachine.TurbulenceClosures
using ClimateMachine.VariableTemplates

using CLIMAParameters
using CLIMAParameters.Planet
using CLIMAParameters.Atmos.SubgridScale

struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

# Menu for initial conditions, boundary conditions and sources
include("gcm_bcs.jl")                # Boundary conditions
include("gcm_perturbations.jl")      # Initial perturbation
include("gcm_base_states.jl")        # Initial base state
include("gcm_moisture_profiles.jl")  # Initial moisture profile
include("gcm_sources.jl")            # GCM-specific sources and parametrizations

# Set default initial conditions, boundary conditions and sources for each experiment type
include("baroclinicwave_problem.jl") # initial value problem (no sources to maintain equilibration)
include("heldsuarez_problem.jl")     # runs to equilibration

# Initial conditions (common to all GCM experiments)
function init_gcm_experiment!(problem, bl, state, aux, coords, t)
    FT = eltype(state)

    # General parameters
    M_v::FT = molmass_ratio(bl.param_set) - 1 # constant for virtual temperature conversion - FIX: this assumes no initial liq/ice

    # Select initial perturbation
    u′, v′, w′, rand_pert =
        init_perturbation(problem.perturbation, bl, state, aux, coords, t)

    # Select initial base state
    T_v, p, u_ref, v_ref, w_ref =
        init_base_state(problem.base_state, bl, state, aux, coords, t)

    # Select initial moisture profile
    q_tot = init_moisture_profile(
        problem.moisture_profile,
        bl,
        state,
        aux,
        coords,
        t,
        p,
    )

    # Calculate initial total winds
    u_sphere = SVector{3, FT}(u_ref + u′, v_ref + v′, w_ref + w′)
    u_cart = sphr_to_cart_vec(bl.orientation, u_sphere, aux)

    # Calculate initial temperature and density
    phase_partition = PhasePartition(q_tot)
    T::FT = T_v / (1 + M_v * q_tot) # this needs to be adapted for ice and liq

    ρ::FT = air_density(bl.param_set, T, p, phase_partition)

    ## potential & kinetic energy
    e_pot::FT = gravitational_potential(bl.orientation, aux)
    e_kin::FT = 0.5 * u_cart' * u_cart
    e_tot::FT = total_energy(bl.param_set, e_kin, e_pot, T, phase_partition)

    ## Assign state variables
    state.ρ = ρ
    state.ρu = ρ * u_cart
    state.ρe = ρ * e_tot * rand_pert

    if bl.moisture isa EquilMoist
        state.moisture.ρq_tot = ρ * q_tot
    end

    return nothing
end

lowercasearg(arg::Nothing) = nothing
lowercasearg(arg) = lowercase(arg)

# Helper for parsing `--experiment` command line argument to initialize defaults for the specified problem
function parse_experiment_arg(arg)
    if arg == "baroclinic_wave"
        problem_type = BaroclinicWaveProblem
    elseif arg == "heldsuarez"
        problem_type = HeldSuarezProblem
    else
        error("unknown experiment: " * arg)
    end

    return problem_type
end

# General GCM configuration setup
function config_gcm_experiment(
    ::Type{FT},
    poly_order,
    resolution,
    experiment_arg,
    perturbation_arg,
    base_state_arg,
    moisture_profile_arg,
    surface_flux_arg,
) where {FT}
    # Set up a reference state for linearization of equations
    temp_profile_ref =
        DecayingTemperatureProfile{FT}(param_set, FT(290), FT(220), FT(8e3))
    ref_state = HydrostaticState(temp_profile_ref)

    # Distance between surface and top of atmosphere (m)
    domain_height::FT = 30e3

    # make the orientation explicit (for surface fluxes)
    orientation = SphericalOrientation()

    # Determine the problem type
    problem_type = parse_experiment_arg(experiment_arg)

    # Set up problem components
    perturbation = parse_perturbation_arg(perturbation_arg)
    base_state = parse_base_state_arg(base_state_arg)
    moisture_profile = parse_moisture_profile_arg(moisture_profile_arg)

    # Choose the moisture model
    if moisture_profile isa NoMoistureProfile
        hyperdiffusion = DryBiharmonic(FT(8 * 3600))
        moisture = DryModel()
    else
        hyperdiffusion = EquilMoistBiharmonic(FT(8 * 3600))
        moisture = EquilMoist{FT}()
    end

    # Set up the boundary conditions
    boundarycondition = parse_surface_flux_arg(
        surface_flux_arg,
        FT,
        param_set,
        orientation,
        moisture,
    )

    # Set up the problem
    problem = problem_type(
        boundarycondition = boundarycondition,
        perturbation = perturbation,
        base_state = base_state,
        moisture_profile = moisture_profile,
    )

    # Create the Atmosphere model
    model = AtmosModel{FT}(
        AtmosGCMConfigType,
        param_set;
        problem = problem,
        orientation = orientation,
        ref_state = ref_state,
        turbulence = ConstantKinematicViscosity(FT(0)),
        hyperdiffusion = hyperdiffusion,
        moisture = moisture,
        source = setup_source(problem),
    )

    # Create the GCM driver configuration
    config = ClimateMachine.AtmosGCMConfiguration(
        problem_name(problem),
        poly_order,
        resolution,
        domain_height,
        param_set,
        init_gcm_experiment!;
        model = model,
    )

    return config
end

# Diagnostics configuration setup
function config_diagnostics(::Type{FT}, driver_config) where {FT}
    interval = "40000steps"

    _planet_radius = planet_radius(param_set)::FT

    info = driver_config.config_info
    boundaries = [
        FT(-90.0) FT(-180.0) _planet_radius
        FT(90.0) FT(180.0) FT(_planet_radius + info.domain_height)
    ]
    resolution = (FT(1), FT(1), FT(1000)) # in (deg, deg, m)
    interpol = ClimateMachine.InterpolationConfiguration(
        driver_config,
        boundaries,
        resolution,
    )

    dgngrp = setup_atmos_default_diagnostics(
        AtmosGCMConfigType(),
        interval,
        driver_config.name,
        interpol = interpol,
    )

    return ClimateMachine.DiagnosticsConfiguration([dgngrp])
end

# Entry point
function main()
    # Helper for command-line arguments to specify the experiment and
    # the override experiment defaults
    # various options
    # TODO: some of this must move to a future namelist functionality
    exp_args = ArgParseSettings(autofix_names = true) # hyphens replaced with underscores in cl_args
    add_arg_group!(exp_args, "GCMDriver")
    @add_arg_table! exp_args begin
        "--experiment"
        help = """
            baroclinic_wave defaults:
                init-perturbation: deterministic,
                init-base-state: bc_wave,
                init-moisture-profile: moist_low_tropics;
            heldsuarez defaults:
                init-perturbation: deterministic,
                init-base-state: heldsuarez,
                init-moisture-profile: moist_low_tropics
            """
        metavar = "baroclinic_wave|heldsuarez"
        arg_type = String
        required = true
        "--init-perturbation"
        help = "select the perturbation to use"
        metavar = "deterministic|random|zero"
        arg_type = String
        "--init-base-state"
        help = "select the initial state to use"
        metavar = "bc_wave|heldsuarez|zero"
        arg_type = String
        "--init-moisture-profile"
        help = "specify the moisture profile"
        metavar = "moist_low_tropics|zero|dry"
        arg_type = String
        "--surface-flux"
        help = "specify surface flux for energy and moisture"
        metavar = "default|bulk"
        arg_type = String
        default = "default"
    end

    cl_args = ClimateMachine.init(parse_clargs = true, custom_clargs = exp_args)
    experiment_arg = lowercasearg(cl_args["experiment"])
    perturbation_arg = lowercasearg(cl_args["init_perturbation"]) # use "_" not "-" in cl_args
    base_state_arg = lowercasearg(cl_args["init_base_state"])
    moisture_profile_arg = lowercasearg(cl_args["init_moisture_profile"])
    surface_flux_arg = lowercasearg(cl_args["surface_flux"])

    # Driver configuration parameters
    FT = Float64                             # floating type precision
    poly_order = 3                           # discontinuous Galerkin polynomial order
    n_horz = 12                              # horizontal element number
    n_vert = 6                               # vertical element number
    n_days::FT = 0.1
    timestart::FT = 0                        # start time (s)
    timeend::FT = n_days * day(param_set)    # end time (s)

    # Set up driver configuration
    driver_config = config_gcm_experiment(
        FT,
        poly_order,
        (n_horz, n_vert),
        experiment_arg,
        perturbation_arg,
        base_state_arg,
        moisture_profile_arg,
        surface_flux_arg,
    )

    # Choose time stepper
    ode_solver_type = ClimateMachine.IMEXSolverType(
        implicit_model = AtmosAcousticGravityLinearModel,
        implicit_solver = ManyColumnLU,
        solver_method = ARK2GiraldoKellyConstantinescu,
        split_explicit_implicit = true,
        discrete_splitting = false,
    )

    # The target acoustic CFL number; time step is computed such that the
    # horizontal acoustic Courant number is CFL.
    CFL = FT(0.1)

    # Set up the solver
    solver_config = ClimateMachine.SolverConfiguration(
        timestart,
        timeend,
        driver_config,
        Courant_number = CFL,
        ode_solver_type = ode_solver_type,
        CFL_direction = HorizontalDirection(),
        diffdir = HorizontalDirection(),
    )

    # Set up diagnostics
    dgn_config = config_diagnostics(FT, driver_config)

    # Set up filters as user-defined callbacks
    filterorder = 20
    filter = ExponentialFilter(solver_config.dg.grid, 0, filterorder)
    cbfilter = GenericCallbacks.EveryXSimulationSteps(1) do
        Filters.apply!(
            solver_config.Q,
            AtmosFilterPerturbations(driver_config.bl),
            solver_config.dg.grid,
            filter,
            state_auxiliary = solver_config.dg.state_auxiliary,
        )
        nothing
    end

    cbtmarfilter = GenericCallbacks.EveryXSimulationSteps(1) do
        Filters.apply!(
            solver_config.Q,
            ("moisture.ρq_tot",),
            solver_config.dg.grid,
            TMARFilter(),
        )
        nothing
    end

    # Run the model
    result = ClimateMachine.invoke!(
        solver_config;
        diagnostics_config = dgn_config,
        user_callbacks = (cbtmarfilter, cbfilter),
        check_euclidean_distance = true,
    )

    return result, solver_config
end

result, solver_config = main()

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
__precompile__(true)

module PyConversions

using Reexport

using DataTables
@reexport using PyCall

const PyPickle = PyNULL()
const PyPandas = PyNULL()

function __init__()
    copy!(PyPickle, pyimport("pickle"))
    copy!(PyPandas, pyimport("pandas"))
end

include("dfutils.jl")
include("utils.jl")
include("macros.jl")

end


text:
# 2D SLAM with velocity states

"""
$(TYPEDEF)
"""
mutable struct DynPoint2VelocityPrior{T <: SamplableBelief} <: AbstractPrior
  z::T
  DynPoint2VelocityPrior{T}() where {T <: SamplableBelief} = new{T}()
  DynPoint2VelocityPrior{T}(z1::T) where {T <: SamplableBelief} = new{T}(z1)
end
DynPoint2VelocityPrior(z1::T) where {T <: SamplableBelief} = DynPoint2VelocityPrior{T}(z1)

getSample(cfo::CalcFactor{<:DynPoint2VelocityPrior}) = rand(cfo.factor.z)

"""
$(TYPEDEF)
"""
mutable struct DynPoint2DynPoint2{T <: SamplableBelief} <: AbstractRelativeRoots
  z::T
  DynPoint2DynPoint2{T}() where {T <: SamplableBelief} = new{T}()
  DynPoint2DynPoint2(z1::T) where {T <: SamplableBelief} = new{T}(z1)
end

getSample(cfo::CalcFactor{<:DynPoint2DynPoint2}) = rand(cfo.factor.z)

function (cfo::CalcFactor{<:DynPoint2DynPoint2})(z, xi, xj)
  #
  dt = Dates.value(cfo.metadata.fullvariables[2].nstime - cfo.metadata.fullvariables[1].nstime)*1e-9   # roughly the intended use of userdata
  res12 = z[1:2] - (xj[1:2] - (xi[1:2]+dt*xi[3:4]))
  res34 = z[3:4] - (xj[3:4] - xi[3:4])
  return [res12; res34]
end

"""
$(TYPEDEF)
"""
mutable struct Point2Point2Velocity{T <: IIF.SamplableBelief} <: IIF.AbstractRelativeMinimize
  z::T
  Point2Point2Velocity{T}() where {T <: IIF.SamplableBelief} = new{T}()
  Point2Point2Velocity(z1::T) where {T <: IIF.SamplableBelief} = new{T}(z1)
end

getSample(cfo::CalcFactor{<:Point2Point2Velocity}) = rand(cfo.factor.z)
function (cfo::CalcFactor{<:Point2Point2Velocity})( z,
                                                    xi,
                                                    xj  )
  #
  dt = (cfo.metadata.fullvariables[2].nstime - cfo.metadata.fullvariables[1].nstime)*1e-9     # roughly the intended use of userdata
  dp = (xj[1:2] .- xi[1:2])
  dv = (xj[3:4] .- xi[3:4])

  res12 = z[1:2] .- dp
  res34 =  dp/dt .- 0.5*(xj[3:4] .+ xi[3:4])  # (dp/dt - 0.5*(xj[3:4]+xi[3:4])) # midpoint integration

  return [res12; res34]
end

## Packing Types================================================================

"""
$(TYPEDEF)
"""
mutable struct PackedDynPoint2VelocityPrior <: IIF.PackedInferenceType
  str::String
  PackedDynPoint2VelocityPrior() = new()
  PackedDynPoint2VelocityPrior(z1::String) = new(z1)
end

function convert(::Type{PackedDynPoint2VelocityPrior}, d::DynPoint2VelocityPrior)
  return PackedDynPoint2VelocityPrior(convert(PackedSamplableBelief, d.z))
end
function convert(::Type{DynPoint2VelocityPrior}, d::PackedDynPoint2VelocityPrior)
  distr = convert(SamplableBelief, d.str)
  return DynPoint2VelocityPrior(distr)
end

#