text:
using Test

using Seis
using Beamforming, Beamforming.Synth

@testset "Beamforming" begin
    @testset "Conventional" begin
        let sxi = 2rand()-1, syi = 2rand()-1
            s = Synth.synthetic_arrival(sxi, syi)
            nsta = length(s)
            s_max = 2
            ds = 0.1
            @testset "$method" for (method, n) in zip(
                    (:linear, :nthroot, :phaseweight), (nothing, 3, 3))
                bf = Beamforming.beamform(s, 40, 60, -s_max, s_max, -s_max, s_max, ds)
                @test bf isa Beamforming.BeamformGrid
                @test length(bf.x) == length(bf.y) == nsta
                @test bf.sx == -s_max:ds:s_max
                @test bf.sy == -s_max:ds:s_max
                @test size(bf.power) == (length(bf.sx), length(bf.sy))
                @test bf.nmax == 1
                @test length(bf.bazmax) == length(bf.sx_max) ==
                    length(bf.sy_max) == length(bf.pmax)
                @test bf.sx_max[1] ≈ sxi atol=ds
                @test bf.sy_max[1] ≈ syi atol=ds
            end
        end
    end

    @testset "Cross-correlation BF" begin
        let sx = 20*(rand() - 0.5), sy = 20*(rand() - 0.5), nsta = 8
            s = Synth.synthetic_arrival(sx, sy, nsta=nsta,
                lon=2rand(nsta), lat=2rand(nsta), distance=90)
            t1, t2 = maximum(starttime.(s)), minimum(endtime.(s))
            s_max = 11
            ds = 0.5
            t = 25
            bf = crosscorrelation_beamform(s, t1, t2,
                -s_max, s_max, -s_max, s_max, ds, t=t)
            @test bf isa Beamforming.BeamformGrid
            @test length(bf.x) == length(bf.y) == nsta
            @test bf.sx == -s_max:ds:s_max
            @test bf.sy == -s_max:ds:s_max
            @test size(bf.power) == (length(bf.sx), length(bf.sy))
            @test bf.nmax == 1
            @test length(bf.bazmax) == length(bf.sx_max) ==
                length(bf.sy_max) == length(bf.pmax)
            @test bf.sx_max[1] ≈ sx atol=ds
            @test bf.sy_max[1] ≈ sy atol=ds
            # Check the max slowness version is the same
            bf′ = crosscorrelation_beamform(s, t1, t2, s_max, ds, t=t)
            @testset "$f" for f in fieldnames(Beamforming.BeamformGrid)
                @test getfield(bf, f) == getfield(bf′, f)
            end
        end
    end

    @testset "Cross-correlation BF with correlations" begin
        let sx = 5*(rand() - 0.5), sy = 5*(rand() - 0.5), nsta = 8
            s = Synth.synthetic_arrival(sx, sy, nsta=nsta,
                lon=rand(nsta), lat=rand(nsta), distance=30)
            t1, t2 = 40, 60
            corrs = Beamforming.compute_crosscorrelations(s, t1, t2)
            win = 1
            smax = 3
            ds = 0.1
            bf = crosscorrelation_beamform(s, t1, t2, smax, ds, t=win)
            bf′ = crosscorrelation_beamform_corrs(corrs, win, smax, ds)
            # Skip some fields in comparison
            @testset "$f" for f in fieldnames(Beamforming.BeamformGrid)
                f in (:traces, :t1, :t2) && continue
                if f in (:x, :y)
                    # Order of unique stations will be different,
                    # but values should be the same.
                    @test sort(getfield(bf, f)) ≈ sort(getfield(bf′, f))
                else
                    @test getfield(bf, f) == getfield(bf′, f) ||
                        # Catch floating-point roundoff errors in different order
                        getfield(bf, f) ≈ getfield(bf′, f)
                end
            end
        end
    end
end

max_stars_repo_name: anowacki/Beamforming

max_issues_repo_name: anowacki/Beamforming

max_forks_repo_name: anowacki/Beamforming


text:
using PyPlot

include("FPDE_problem.jl")

Nx = 20
Nt = 30
Nquad = 15
x = range(0, L, length=Nx)
t = Float64[ max_t * (j/Nt)^1.5 for j = 1:Nt ]
X = Float64[ x[i] for i = 1:Nx, j = 1:Nt ]
T = Float64[ t[j] for i = 1:Nx, j = 1:Nt ]
u = Array{Float64}(undef, Nx, Nt)
u[:,1] = u0.(x)
for j = 2:Nt, i = 1:Nx
    u[i,j] = refsoln(x[i], t[j], α, Nquad)
end

figure(1)
mesh(X, T, u)
xlabel(L"$x$")
ylabel(L"$t$")
ax = gca()
ax.view_init(elev=30, azim=35)
savefig("fig5.eps")

ext: jl


text:
# Basic surface operators

"""
    regularize!(s::Nodes{Primal},f::ScalarData,cache::BasicILMCache)
    regularize!(s::Nodes{Primal},f::ScalarData,sys::ILMSystem)

The operation ``s = R_c f``, which regularizes scalar surface data `f`
onto the grid in the form of scalar grid data `s`. This is the adjoint
to [`interpolate!`](@ref)
"""
@inline regularize!(s::Nodes{Primal},f::ScalarData,cache::BasicILMCache) = regularize!(s,f,cache.R)

function regularize!(s::Nodes{Primal},f::ScalarData,Rc::RegularizationMatrix)
  return s .= Rc*f
end

"""
    regularize!(v::Edges{Primal},vb::VectorData,cache::BasicILMCache)
    regularize!(v::Edges{Primal},vb::VectorData,sys::ILMSystem)

The operation ``\\mathbf{v} = R_f \\mathbf{v}_b``, which regularizes vector surface data `vb`
onto the grid in the form of scalar grid data `v`. This is the adjoint
to [`interpolate!`](@ref)
"""
@inline regularize!(v::Edges{Primal},vb::VectorData,cache::BasicILMCache) = regularize!(v,vb,cache.R)

function regularize!(v::Edges{Primal},vb::VectorData,Rf::RegularizationMatrix)
  return v .= Rf*vb
end

"""
    interpolate!(f::ScalarData,s::Nodes{Primal},cache::BasicILMCache)
    interpolate!(f::ScalarData,s::Nodes{Primal},sys::ILMSystem)

The operation ``f = R_c^T s``, which interpolates scalar grid data `s`
onto the surface points in the form of scalar point data `f`. This is the adjoint
to [`regularize!`](@ref)
"""
@inline interpolate!(f::ScalarData,s::Nodes{Primal},cache::BasicILMCache) = interpolate!(f,s,cache.E)

function interpolate!(f::ScalarData,s::Nodes{Primal},Ec::InterpolationMatrix)
  return f .= Ec*s
end

"""
    interpolate!(vb::VectorData,v::Edges{Primal},cache::BasicILMCache)
    interpolate!(vb::VectorData,v::Edges{Primal},sys::ILMSystem)

The operation ``\\mathbf{v}_b = R_c^T \\mathbf{v}``, which interpolates vector grid data `v`
onto the surface points in the form of scalar point data `vb`. This is the adjoint
to [`regularize!`](@ref)
"""
@inline interpolate!(vb::VectorData,v::Edges{Primal},cache::BasicILMCache) = interpolate!(vb,v,cache.E)

function interpolate!(vb::VectorData,v::Edges{Primal},Ef::InterpolationMatrix)
  return vb .= Ef*v
end

"""
    regularize_normal!(v::Edges{Primal},f::ScalarData,cache::BasicILMCache)
    regularize_normal!(v::Edges{Primal},f::ScalarData,sys::ILMSystem)

The operation ``\\mathbf{v} = R_f \\mathbf{n}\\circ f``, which maps scalar surface data `f` (like
a jump in scalar potential) to grid data `v` (like velocity). This is the adjoint
to [`normal_interpolate!`](@ref).
"""
@inline regularize_normal!(q::Edges{Primal},f::ScalarData,cache::BasicILMCache) = regularize_normal!(q,f,cache.nrm,cache.Rsn,cache.snorm_cache)

function regularize_normal!(q::Edges{Primal,NX,NY},f::ScalarData{N},nrm::VectorData{N},Rf::RegularizationMatrix,snorm_cache::VectorData{N}) where {NX,NY,N}
    product!(snorm_cache,nrm,f)
    q .= Rf*snorm_cache
end

"""
    regularize_normal!(qt::EdgeGradient{Primal},v::VectorData,cache::BasicILMCache)
    regularize_normal!(qt::EdgeGradient{Primal},v::VectorData,sys::ILMSystem)

The operation ``\\mathbf{q}_t = R_t \\mathbf{n}\\circ \\mathbf{v}``, which maps scalar vector data `v` (like
a jump in velocity) to grid data `qt` (like velocity-normal tensor). This is the adjoint
to [`normal_interpolate!`](@ref).
"""
@inline regularize_normal!(q::EdgeGradient{Primal},f::VectorData,cache::BasicILMCache) = regularize_normal!(q,f,cache.nrm,cache.Rsn,cache.snorm_cache,cache.snorm2_cache)

function regularize_normal!(q::EdgeGradient{Primal,Dual,NX,NY},f::VectorData{N},nrm::VectorData{N},Rf::RegularizationMatrix,snorm_cache::TensorData{N},snorm2_cache::TensorData{N}) where {NX,NY,N}
    tensorproduct!(snorm_cache,nrm,f)
    transpose!(snorm2_cache,snorm_cache)
    snorm_cache .+= snorm2_cache
    q .= Rf*snorm_cache
end

"""
    regularize_normal_cross!(v::Edges{Primal},f::ScalarData,cache::BasicILMCache)
    regularize_normal_cross!(v::Edges{Primal},f::ScalarData,sys::ILMSystem)

The operation ``\\mathbf{v} = R_f \\mathbf{n}\\times f \\mathbf{e}_z``, which maps scalar surface data `f` (like
a jump in streamfunction, endowed with the out-of-plane unit vector ``\\mathbf{e}_z``) to grid
data `v` (like velocity). This is the adjoint to [`normal_cross_interpolate!`](@ref).
"""
@inline regularize_normal_cross!(q::Edges{Primal},f::ScalarData,cache::BasicILMCache) =
           regularize_normal_cross!(q,f,cache.nrm,cache.Rsn,cache.snorm_cache)

function regularize_normal_cross!(q::Edges{Primal,NX,NY},f::ScalarData{N},nrm::VectorData{N},Rf::RegularizationMatrix,snorm_cache::VectorData{N}) where {NX,NY,N}
    cross!(snorm_cache,nrm,f)
    q .= Rf*snorm_cache
end

"""
    normal_interpolate!(vn::ScalarData,v::Edges{Primal},cache::BasicILMCache)
    normal_interpolate!(vn::ScalarData,v::Edges{Primal},sys::ILMSystem)

The operation ``v_n = \\mathbf{n} \\cdot R_f^T \\mathbf{v}``, which maps grid data `v` (like velocity) to scalar
surface data `vn` (like normal component of surface velocity). This is the
adjoint to [`regularize_normal!`](@ref).
"""
@inline normal_interpolate!(vn::ScalarData,q::Edges{Primal},cache::BasicILMCache) = normal_interpolate!(vn,q,cache.nrm,cache.Esn,cache.snorm_cache)

function normal_interpolate!(vn::ScalarData{N},q::Edges{Primal,NX,NY},nrm::VectorData{N},Ef::InterpolationMatrix,snorm_cache::VectorData{N}) where {NX,NY,N}
    snorm_cache .= Ef*q
    pointwise_dot!(vn,nrm,snorm_cache)
end

"""
    normal_interpolate!(τ::VectorData,A::EdgeGradient{Primal},cache::BasicILMCache)
    normal_interpolate!(τ::VectorData,A::EdgeGradient{Primal},sys::ILMSystem)

The operation ``\\mathbf{\\tau} = \\mathbf{n} \\cdot R_t^T (\\mathbf{A} +\\mathbf{A}^T)``, which maps grid tensor data `A` (like velocity gradient tensor) to vector
surface data `τ` (like traction). This is the adjoint to [`regularize_normal!`](@ref).
"""
@inline normal_interpolate!(vn::VectorData,q::EdgeGradient{Primal},cache::BasicILMCache) = normal_interpolate!(vn,q,cache.nrm,cache.Esn,cache.gsnorm2_cache,cache.snorm_cache)

function normal_interpolate!(vn::VectorData{N},q::EdgeGradient{Primal,Dual,NX,NY},nrm::VectorData{N},Ef::InterpolationMatrix,gsnorm_cache::EdgeGradient{Primal,Dual,NX,NY},snorm_cache::TensorData{N}) where {NX,NY,N}
    transpose!(gsnorm_cache,q)
    gsnorm_cache .+= q
    snorm_cache .= Ef*gsnorm_cache
    pointwise_dot!(vn,nrm,snorm_cache)
end

"""
    normal_cross_interpolate!(wn::ScalarData,v::Edges{Primal},cache::BasicILMCache)
    normal_cross_interpolate!(wn::ScalarData,v::Edges{Primal},sys::ILMSystem)

The operation ``w_n = e_z\\cdot (\\mathbf{n} \\times R_f^T \\mathbf{v})``, which maps grid data `v` (like velocity) to scalar
surface data `wn` (like vorticity in the surface). This is the
adjoint to [`regularize_normal!`](@ref).
"""
@inline normal_cross_interpolate!(vn::ScalarData,q::Edges{Primal},cache::BasicILMCache) =
          normal_cross_interpolate!(vn,q,cache.nrm,cache.Esn,cache.snorm_cache)

function normal_cross_interpolate!(vn::ScalarData{N},q::Edges{Primal,NX,NY},nrm::VectorData{N},Ef::InterpolationMatrix,snorm_cache::VectorData{N}) where {NX,NY,N}
    snorm_cache .= Ef*q
    cross!(vn,nrm,snorm_cache)
end

"""
    surface_curl!(w::Nodes{Dual},f::ScalarData,cache::BasicILMCache)
    surface_curl!(w::Nodes{Dual},f::ScalarData,sys::ILMSystem)

The operation ``w = C_s^T f = C^T R_f \\mathbf{n}\\circ f``, which maps scalar surface data `f` (like
a jump in scalar potential) to grid data `w` (like vorticity). This is the adjoint
to ``C_s``, also given by `surface_curl!` (but with arguments switched). Note that
the differential operations are divided either by 1 or by the grid cell size,
depending on whether `sys` has been designated with `IndexScaling` or `GridScaling`,
respectively.
"""
function surface_curl!(w::Nodes{Dual},f::ScalarData,cache::BasicILMCache)
  _unscaled_surface_curl!(w,f,cache.nrm,cache.Rsn,cache.gsnorm_cache,cache.snorm_cache)
  _scale_derivative!(w,cache)
end

function _unscaled_surface_curl!(w::Nodes{Dual,NX,NY},f::ScalarData{N},nrm::VectorData{N},Rf::RegularizationMatrix,q_cache::Edges{Primal,NX,NY},snorm_cache::VectorData{N}) where {NX,NY,N}
    fill!(w,0.0)
    fill!(q_cache,0.0)
    regularize_normal!(q_cache,f,nrm,Rf,snorm_cache)
    curl!(w,q_cache)
end

"""
    surface_curl!(w::Nodes{Dual},v::VectorData,cache::BasicILMCache)
    surface_curl!(w::Nodes{Dual},v::VectorData,sys::ILMSystem)

The operation ``w = C_s^T \\mathbf{v} = C^T R_f \\mathbf{v}``, which maps vector surface data `v` (like
velocity) to grid data `w` (like vorticity). This is the adjoint
to ``C_s``, also given by `surface_curl!` (but with arguments switched). Note that
the differential operations are divided either by 1 or by the grid cell size,
depending on whether `sys` has been designated with `IndexScaling` or `GridScaling`,
respectively.
"""
function surface_curl!(w::Nodes{Dual},v::VectorData,cache::BasicILMCache)
  _unscaled_surface_curl!(w,v,cache.R,cache.gdata_cache)
  _scale_derivative!(w,cache)
end

function _unscaled_surface_curl!(w::Nodes{Dual,NX,NY},v::VectorData{N},Rf::RegularizationMatrix,q_cache::Edges{Primal,NX,NY}) where {NX,NY,N}
    fill!(w,0.0)
    fill!(q_cache,0.0)
    regularize!(q_cache,v,Rf)
    curl!(w,q_cache)
end

"""
    surface_curl!(vn::ScalarData,s::Nodes{Dual},cache::BasicILMCache)
    surface_curl!(vn::ScalarData,s::Nodes{Dual},sys::ILMSystem)

The operation ``v_n = C_s s = \\mathbf{n} \\cdot R_f^T C s``, which maps grid data `s` (like
streamfunction) to scalar surface data `vn` (like normal component of velocity).
This is the adjoint to ``C_s^T``, also given by `surface_curl!`, but with
arguments switched.  Note that
the differential operations are divided either by 1 or by the grid cell size,
depending on whether `sys` has been designated with `IndexScaling` or `GridScaling`,
respectively.
"""
function surface_curl!(vn::ScalarData,s::Nodes{Dual},cache::BasicILMCache)
  _unscaled_surface_curl!(vn,s,cache.nrm,cache.Esn,cache.gsnorm_cache,cache.snorm_cache)
  _scale_derivative!(vn,cache)
end

function _unscaled_surface_curl!(vn::ScalarData{N},s::Nodes{Dual,NX,NY},nrm::VectorData{N},Ef::InterpolationMatrix,q_cache::Edges{Primal,NX,NY},snorm_cache::VectorData{N}) where {NX,NY,N}
    fill!(vn,0.0)
    fill!(q_cache,0.0)
    curl!(q_cache,s)
    normal_interpolate!(vn,q_cache,nrm,Ef,snorm_cache)
end

"""
    surface_curl!(v::VectorData,s::Nodes{Dual},cache::BasicILMCache)
    surface_curl!(v::VectorData,s::Nodes{Dual},sys::ILMSystem)

The operation ``\\mathbf{v} = C_s s = R_f^T C s``, which maps grid data `s` (like
streamfunction) to vector surface data `v` (like velocity).
This is the adjoint to ``C_s^T``, also given by `surface_curl!`, but with
arguments switched.  Note that
the differential operations are divided either by 1 or by the grid cell size,
depending on whether `sys` has been designated with `IndexScaling` or `GridScaling`,
respectively.
"""
function surface_curl!(v::VectorData,s::Nodes{Dual},cache::BasicILMCache)
  _unscaled_surface_curl!(v,s,cache.E,cache.gdata_cache)
  _scale_derivative!(v,cache)
end

function _unscaled_surface_curl!(v::VectorData{N},s::Nodes{Dual,NX,NY},Ef::InterpolationMatrix,q_cache::Edges{Primal,NX,NY}) where {NX,NY,N}
    fill!(q_cache,0.0)
    curl!(q_cache,s)
    interpolate!(v,q_cache,Ef)
end

"""
    surface_curl_cross!(w::Nodes{Dual},f::ScalarData,cache::BasicILMCache)
    surface_curl_cross!(w::Nodes{Dual},f::ScalarData,sys::ILMSystem)

The operation ``w = \\hat{C}_s^T f = C^T R_f \\mathbf{n}\\times f \\mathbf{e}_z``, which maps scalar surface data `f` (like
a jump in streamfunction, multiplied by out-of-plane unit vector ``\\mathbf{e}_z``) to grid data `w` (like vorticity).
This is the adjoint to ``\\hat{C}_s``, also given by `surface_curl_cross!` (but with arguments switched). Note that
the differential operations are divided either by 1 or by the grid cell size,
depending on whether `sys` has been designated with `IndexScaling` or `GridScaling`,
respectively.
"""
function surface_curl_cross!(w::Nodes{Dual},f::ScalarData,cache::BasicILMCache)
  _unscaled_surface_curl_cross!(w,f,cache.nrm,cache.Rsn,cache.gsnorm_cache,cache.snorm_cache)
  _scale_derivative!(w,cache)
end

function _unscaled_surface_curl_cross!(w::Nodes{Dual,NX,NY},f::ScalarData{N},nrm::VectorData{N},Rf::RegularizationMatrix,q_cache::Edges{Primal,NX,NY},snorm_cache::VectorData{N}) where {NX,NY,N}
    fill!(w,0.0)
    fill!(q_cache,0.0)
    regularize_normal_cross!(q_cache,f,nrm,Rf,snorm_cache)
    curl!(w,q_cache)
end

"""
    surface_curl_cross!(γ::ScalarData,s::Nodes{Dual},cache::BasicILMCache)
    surface_curl_cross!(γ::ScalarData,s::Nodes{Dual},sys::ILMSystem)

The operation ``\\gamma = \\hat{C}_s s = \\mathbf{e}_z\\cdot (\\mathbf{n} \\times R_f^T C s)``, which maps grid data `s` (like
streamfunction) to scalar surface data `γ` (like vorticity in the surface).
This is the adjoint to ``\\hat{C}_s^T``, also given by `surface_curl_cross!`, but with
arguments switched.  Note that
the differential operations are divided either by 1 or by the grid cell size,
depending on whether `sys` has been designated with `IndexScaling` or `GridScaling`,
respectively.
"""
function surface_curl_cross!(vn::ScalarData,s::Nodes{Dual},cache::BasicILMCache)
  _unscaled_surface_curl_cross!(vn,s,cache.nrm,cache.Esn,cache.gsnorm_cache,cache.snorm_cache)
  _scale_derivative!(vn,cache)
end

function _unscaled_surface_curl_cross!(vn::ScalarData{N},s::Nodes{Dual,NX,NY},nrm::VectorData{N},Ef::InterpolationMatrix,q_cache::Edges{Primal,NX,NY},snorm_cache::VectorData{N}) where {NX,NY,N}
    fill!(vn,0.0)
    fill!(q_cache,0.0)
    curl!(q_cache,s)
    normal_cross_interpolate!(vn,q_cache,nrm,Ef,snorm_cache)
end

"""
    surface_divergence!(Θ::Nodes{Primal},f::ScalarData,cache::BasicILMCache)
    surface_divergence!(Θ::Nodes{Primal},f::ScalarData,sys::ILMSystem)

The operation ``\\theta = D_s f = D R_f \\mathbf{n} \\circ f``, which maps surface scalar data `f` (like
jump in scalar potential) to grid data `Θ` (like dilatation, i.e. divergence of velocity).
This is the adjoint of [`surface_grad!`](@ref).
Note that
the differential operations are divided either by 1 or by the grid cell size,
depending on whether `sys` has been designated with `IndexScaling` or `GridScaling`,
respectively.
"""
function surface_divergence!(θ::Nodes{Primal},f::ScalarData,cache::BasicILMCache)
  _unscaled_surface_divergence!(θ,f,cache.nrm,cache.Rsn,cache.gsnorm_cache,cache.snorm_cache)
  _scale_derivative!(θ,cache)
end

function _unscaled_surface_divergence!(θ::Nodes{Primal,NX,NY},f::ScalarData{N},nrm::VectorData{N},Rf::RegularizationMatrix,q_cache::Edges{Primal,NX,NY},snorm_cache::VectorData{N}) where {NX,NY,N}
    fill!(θ,0.0)
    fill!(q_cache,0.0)
    regularize_normal!(q_cache,f,nrm,Rf,snorm_cache)
    divergence!(θ,q_cache)
end

"""
    surface_divergence!(v::Edges{Primal},dv::VectorData,cache::BasicILMCache)
    surface_divergence!(v::Edges{Primal},dv::VectorData,sys::ILMSystem)

The operation ``\\mathbf{v} = D_s d\\mathbf{v} = D R_f (\\mathbf{n} \\circ d\\mathbf{v} + d\\mathbf{v} \\circ \\mathbf{n})``, which maps surface vector data `dv` (like
jump in velocity) to grid data `v` (like velocity). This is the adjoint of [`surface_grad!`](@ref). Note that
the differential operations are divided either by 1 or by the grid cell size,
depending on whether `sys` has been designated with `IndexScaling` or `GridScaling`,
respectively.
"""
function surface_divergence!(θ::Edges{Primal},f::VectorData,cache::BasicILMCache)
  _unscaled_surface_divergence!(θ,f,cache.nrm,cache.Rsn,cache.gsnorm_cache,cache.snorm_cache,cache.snorm2_cache)
  _scale_derivative!(θ,cache)
end

function _unscaled_surface_divergence!(θ::Edges{Primal,NX,NY},f::VectorData{N},nrm::VectorData{N},Rf::RegularizationMatrix,q_cache::EdgeGradient{Primal,Dual,NX,NY},snorm_cache::TensorData{N},snorm2_cache::TensorData{N}) where {NX,NY,N}
    fill!(θ,0.0)
    fill!(q_cache,0.0)
    regularize_normal!(q_cache,f,nrm,Rf,snorm_cache,snorm2_cache)
    divergence!(θ,q_cache)
end

"""
    surface_grad!(vn::ScalarData,ϕ::Nodes{Primal},cache::BasicILMCache)
    surface_grad!(vn::ScalarData,ϕ::Nodes{Primal},sys::ILMSystem)

The operation ``v_n = G_s\\phi = \\mathbf{n} \\cdot R_f^T G\\phi``, which maps grid data `ϕ` (like
scalar potential) to scalar surface data `vn` (like normal component of velocity). This is the adjoint of [`surface_divergence!`](@ref).
Note that
the differential operations are divided either by 1 or by the grid cell size,
depending on whether `sys` has been designated with `IndexScaling` or `GridScaling`,
respectively.
"""
function surface_grad!(vn::ScalarData,ϕ::Nodes{Primal},cache::BasicILMCache)
  _unscaled_surface_grad!(vn,ϕ,cache.nrm,cache.Esn,cache.gsnorm_cache,cache.snorm_cache)
  _scale_derivative!(vn,cache)
end

function _unscaled_surface_grad!(vn::ScalarData{N},ϕ::Nodes{Primal,NX,NY},nrm::VectorData{N},Ef::InterpolationMatrix,q_cache::Edges{Primal,NX,NY},snorm_cache::VectorData{N}) where {NX,NY,N}
    fill!(vn,0.0)
    fill!(q_cache,0.0)
    grad!(q_cache,ϕ)
    normal_interpolate!(vn,q_cache,nrm,Ef,snorm_cache)
end

"""
    surface_grad!(τ::VectorData,v::Edges{Primal},cache::BasicILMCache)
    surface_grad!(τ::VectorData,v::Edges{Primal},sys::ILMSystem)

The operation ``\\mathbf{\\tau} = G_s v = \\mathbf{n} \\cdot R_t^T (G \\mathbf{v} + (G \\mathbf{v})^T)``, which maps grid vector data `v` (like
velocity) to vector surface data `τ` (like traction). This is the adjoint of [`surface_divergence!`](@ref). Note that
the differential operations are divided either by 1 or by the grid cell size,
depending on whether `sys` has been designated with `IndexScaling` or `GridScaling`,
respectively.
"""
function surface_grad!(vn::VectorData,ϕ::Edges{Primal},cache::BasicILMCache)
  _unscaled_surface_grad!(vn,ϕ,cache.nrm,cache.Esn,cache.gsnorm_cache,cache.gsnorm2_cache,cache.snorm_cache)
  _scale_derivative!(vn,cache)
end

function _unscaled_surface_grad!(vn::VectorData{N},ϕ::Edges{Primal,NX,NY},nrm::VectorData{N},Ef::InterpolationMatrix,gsnorm_cache::EdgeGradient{Primal,Dual,NX,NY},gsnorm2_cache::EdgeGradient{Primal,Dual,NX,NY},snorm_cache::TensorData{N}) where {NX,NY,N}
    fill!(vn,0.0)
    fill!(gsnorm_cache,0.0)
    grad!(gsnorm_cache,ϕ)
    normal_interpolate!(vn,gsnorm_cache,nrm,Ef,gsnorm2_cache,snorm_cache)
end

"""
    surface_divergence_cross!(Θ::Nodes{Primal},f::ScalarData,cache::BasicILMCache)
    surface_divergence_cross!(Θ::Nodes{Primal},f::ScalarData,sys::ILMSystem)

The operation ``\\theta = \\hat{D}_s f = D R_f \\mathbf{n} \\times f \\mathbf{e}_z``, which maps surface scalar data `f` (like
jump in streamfunction) to grid data `Θ` (like dilatation, i.e. divergence of velocity).
This is the adjoint of [`surface_grad_cross!`](@ref).
Note that the differential operations are divided either by 1 or by the grid cell size,
depending on whether `sys` has been designated with `IndexScaling` or `GridScaling`,
respectively.
"""
function surface_divergence_cross!(θ::Nodes{Primal},f::ScalarData,cache::BasicILMCache)
  _unscaled_surface_divergence_cross!(θ,f,cache.nrm,cache.Rsn,cache.gsnorm_cache,cache.snorm_cache)
  _scale_derivative!(θ,cache)
end

function _unscaled_surface_divergence_cross!(θ::Nodes{Primal,NX,NY},f::ScalarData{N},nrm::VectorData{N},Rf::RegularizationMatrix,q_cache::Edges{Primal,NX,NY},snorm_cache::VectorData{N}) where {NX,NY,N}
    fill!(θ,0.0)
    fill!(q_cache,0.0)
    regularize_normal_cross!(q_cache,f,nrm,Rf,snorm_cache)
    divergence!(θ,q_cache)
end

"""
    surface_grad_cross!(γ::ScalarData,ϕ::Nodes{Primal},cache::BasicILMCache)
    surface_grad_cross!(γ::ScalarData,ϕ::Nodes{Primal},sys::ILMSystem)

The operation ``\\gamma = \\hat{G}_s\\phi = \\mathbf{e}_z \\cdot (\\mathbf{n} \\times R_f^T G\\phi)``, which maps grid data `ϕ` (like
scalar potential) to scalar surface data `γ` (like surface vorticity). This is the adjoint of [`surface_divergence_cross!`](@ref).
Note that
the differential operations are divided either by 1 or by the grid cell size,
depending on whether `sys` has been designated with `IndexScaling` or `GridScaling`,
respectively.
"""
function surface_grad_cross!(vn::ScalarData,ϕ::Nodes{Primal},cache::BasicILMCache)
  _unscaled_surface_grad_cross!(vn,ϕ,cache.nrm,cache.Esn,cache.gsnorm_cache,cache.snorm_cache)
  _scale_derivative!(vn,cache)
end

function _unscaled_surface_grad_cross!(vn::ScalarData{N},ϕ::Nodes{Primal,NX,NY},nrm::VectorData{N},Ef::InterpolationMatrix,q_cache::Edges{Primal,NX,NY},snorm_cache::VectorData{N}) where {NX,NY,N}
    fill!(vn,0.0)
    fill!(q_cache,0.0)
    grad!(q_cache,ϕ)
    normal_cross_interpolate!(vn,q_cache,nrm,Ef,snorm_cache)
end

"""
    mask(cache::BasicILMCache) -> GridData
    mask(sys::ILMSystem) -> GridData

Create grid data that consist of 1s inside of a surface (i.e., on a side opposite
  the normal vectors) and 0s outside. The grid data are the same type as the
  output data type of `sys`.  Only allows `sys` to have `GridScaling`.
"""
@inline mask(cache::BasicILMCache{N,GridScaling}) where {N} = _mask!(zeros_grid(cache),cache)

"""
    complementary_mask(cache::BasicILMCache) -> GridData
    complementary_mask(sys::ILMSystem) -> GridData

Create grid data that consist of 0s inside of a surface (i.e., on a side opposite
  the normal vectors) and 1s outside. The grid data are the same type as the
  output data type of `cache`.  Only allows `cache` to have `GridScaling`.
"""
@inline complementary_mask(cache::BasicILMCache{N,GridScaling}) where {N} =
          _complementary_mask!(zeros_grid(cache),cache)

"""
    mask!(w::GridData,cache::BasicILMCache)
    mask!(w::GridData,sys::ILMSystem)

Mask the data `w` in place by multiplying it by 1s inside of a surface (i.e., on a side opposite
  the normal vectors) and 0s outside. The grid data `w` must be of the same type as the
  output data type of `cache`. Only allows `cache` to have `GridScaling`.
"""
function mask!(w::T,cache::BasicILMCache{N,GridScaling}) where {T <: GridData,N}
  @unpack gdata_cache = cache
  _mask!(gdata_cache,cache)
  product!(w,gdata_cache,w)
end

"""
    complementary_mask!(w::GridData,cache::BasicILMCache)
    complementary_mask!(w::GridData,sys::ILMSystem)

Mask the data `w` in place by multiplying it by 0s inside of a surface (i.e., on a side opposite
  the normal vectors) and 1s outside. The grid data `w` must be of the same type as the
  output data type of `cache`. Only allows `cache` to have `GridScaling`.
"""
function complementary_mask!(w::T,cache::BasicILMCache{N,GridScaling}) where {T <: GridData,N}
  @unpack gdata_cache = cache
  _complementary_mask!(gdata_cache,cache)
  product!(w,gdata_cache,w)
end

function _mask!(msk,cache)
  @unpack sdata_cache, gdata_cache, L = cache
  typeof(msk) == typeof(gdata_cache) || error("Wrong data type")
  fill!(sdata_cache,1.0)
  surface_divergence!(msk,sdata_cache,cache)
  inverse_laplacian!(msk,cache)
  msk .*= -1.0
end

function _complementary_mask!(msk,cache)
    _mask!(msk,cache)
    msk .= 1.0 - msk
end


text:
type SkewNormalSampler <: Sampleable{Univariate,Continuous}
    ξ::Float64
    ω::Float64
    aux::MvNormal
    function SkewNormalSampler(ξ::Real, ω::Real, α::Real)
        δ = α/(1 + α^2)
        aux = MvNormal([[1 δ]; [δ 1]])
        new(ξ, ω, aux)
    end
end

function rand(dist::SkewNormalSampler)
    x0, x = rand(dist.aux)
    dist.ω * ((x0 > 0) ? x : -x) + dist.ξ
end

type SkewNormal <: ContinuousUnivariateDistribution
    ξ::Float64  # Location
    ω::Float64  # Scale
    α::Float64  # Skewness
end

sampler(dist::SkewNormal) = SkewNormalSampler(dist.ξ, dist.ω, dist.α)

function pdf(dist::SkewNormal, x::Real)
    z = (x - dist.ξ)/dist.ω
    2.0 * (normpdf(z)/dist.ω) * normcdf(dist.α * z)
end

function cdf(dist::SkewNormal, x::Real)
    quadgk(t->pdf(dist,t), -Inf, x)[1]
end

function quantile(dist::SkewNormal, β::Float64)
    newton(x->cdf(dist,x) - β, dist.ξ)
end

minimum(dist::SkewNormal) = -Inf
maximum(dist::SkewNormal) = Inf

mean(dist::SkewNormal) = dist.ω*sqrt(2/π)*(dist.α/sqrt(1.0 + dist.α^2)) + dist.ξ
var(dist::SkewNormal) = (dist.ω^2) * (1.0 - (sqrt(2/π)*(dist.α/sqrt(1.0 + dist.α^2)))^2)

# Cumulant generating function

cgf(dist::SkewNormal, t::Real) = t*dist.ξ + 0.5*t^2*dist.ω^2 + log(2normcdf((dist.ω*t*dist.α)/sqrt(1 + dist.α^2)))
mgf(dist::SkewNormal, t::Real) = exp(cgf(dist, t))


text:
using ControlSystemIdentification, ControlSystems
n = 5
w = exp10.(LinRange(-2, 3, 1000))
av = exp10.(LinRange(-2, 3, n))
signal = randn(1000)

for basis in [kautz(av), laguerre(1,n), laguerre_oo(1,n), adhocbasis(av)]
    p = rand(basislength(basis))
    F = sum_basis(basis, p)
    fr1 = freqresp(F, w) |> vec
    fr2 = freqresp(tf(1), basis, w, p) |> vec
    @test norm(fr1-fr2) < 1e-10
    @test basislength(basis) == n
end
basis = laguerre_oo(1,n)
b2 = add_poles(basis, [-2])
@test basislength(b2) == n+1
@test -2 ∈ pole(b2)

basis = kautz(av, 0.001)
y = filter_bank(basis, signal)
@test size(y) == (1000, n)

## Test that it's possible to identify a basis model in time domain

N = 2000
h = 0.01
t = range(0, step=h, length=N)
u = randn(N)
G = tf(0.8, [1, -0.9], 0.01)
signal = lsim(G, u, t)[1][:]

n = 10
av = exp10.(LinRange(0, log10(80pi), n))
w = exp10.(LinRange(-2, log10(100pi), 300))
# basis = laguerre_oo(1, n)
basis = kautz(av, h)
Y = filter_bank(basis, u)
p = Y\signal
@test norm(signal - Y*p)/norm(signal) < 1e-3
Gh = sum_basis(basis, p)
@test isdiscrete(Gh) == isdiscrete(basis)

@test mean(abs2, abs2.(freqresp(G, w)) - abs2.(freqresp(Gh, w))) < 0.1

@test mean(abs2, (freqresp(G, w)) - (freqresp(Gh, w))) < 0.1

if isinteractive()
    bodeplot(G, w)
    bodeplot!(Gh, w)
end


text:
function latex_title(args...)
    fontsize(20)
    latex(L"E=mc^2", 0, -200)
end

function circ(p = O, color = "black")
    sethue(color)
    circle(p, 25, :fill)
    return p
end

function circ_ret_trans(p = O, color = "black")
    sethue(color)
    circle(p, 25, :fill)
    return Transformation(p, 0.0)
end

function path!(points, pos, color)
    sethue(color)
    push!(points, pos)
    circle.(points, 2, :fill)
end

function rad(p1, p2, color)
    sethue(color)
    line(p1, p2, :stroke)
end

@testset "Dancing circles (gif)" begin
    p1 = Point(100, 0)
    p2 = Point(200, 80)
    from_rot = 0.0
    to_rot = 2π
    path_of_blue = Point[]
    path_of_red = Point[]

    video = Video(500, 500)
    back = Background(1:25, ground, Point(25, 25))

    Object(latex_title)
    red_ball = Object(RFrames(-24:0), (args...) -> circ(O, "red"), p1)
    act!(red_ball, Action(anim_rotate_around(from_rot, to_rot, O)))

    blue_ball = Object(1:25, (args...) -> circ(O, "blue"), p2)
    act!(blue_ball, Action(anim_rotate_around(to_rot, from_rot, red_ball)))
    path_red =
        Object(1:25, (video, args...) -> path!(path_of_red, get_position(red_ball), "red"))
    path_blue =
        Object(:same, (video, args...) -> path!(path_of_blue, pos(blue_ball), "blue"))
    string = Object(1:25, (args...) -> rad(pos(red_ball), pos(blue_ball), "black"))

    render(video; tempdirectory = "images", pathname = "dancing.gif")

    @test_reference "refs/dancing_circles_16.png" load("images/0000000016.png")
    @test isfile("dancing.gif")
    for i in 1:25
        rm("images/$(lpad(i, 10, "0")).png")
    end
    rm("images/palette.png")
    rm("dancing.gif")
end

@testset "Dancing circles (mp4)" begin
    p1 = Point(100, 0)
    p2 = Point(200, 80)
    from_rot = 0.0
    to_rot = 2π
    path_of_blue = Point[]
    path_of_red = Point[]

    video = Video(500, 500)
    back = Background(1:30, ground, Point(25, 25))

    Object(latex_title)
    red_ball = Object((args...) -> circ(O, "red"), p1)
    act!(red_ball, Action(anim_rotate_around(from_rot, to_rot, O)))

    blue_ball = Object(1:30, (args...) -> circ(O, "blue"), p2)
    act!(blue_ball, Action(RFrames(1:30), anim_rotate_around(to_rot, from_rot, red_ball)))
    path_red = Object((video, args...) -> path!(path_of_red, get_position(red_ball), "red"))
    path_blue =
        Object(:same, (video, args...) -> path!(path_of_blue, pos(blue_ball), "blue"))
    string = Object((args...) -> rad(pos(red_ball), pos(blue_ball), "black"))

    render(video; tempdirectory = "images", pathname = "dancing.mp4", framerate = 1)

    # 30 frames with a framerate of 1 should take about 30 seconds ;)
    @test isapprox(VideoIO.get_duration("dancing.mp4"), 30.0, atol = 0.1)
    rm("dancing.mp4")
end

@testset "Dancing circles layered" begin
    p1 = Point(100, 0)
    p2 = Point(200, 80)
    from_rot = 0.0
    to_rot = 2π
    path_of_blue = Point[]
    path_of_red = Point[]

    video = Video(500, 500)
    back = Object(1:25, ground, in_global_layer = true)
    act!(back, Action(anim_rotate(π / 2, π / 2)))
    act!(back, Action(1:1, anim_translate(Point(25, 25))))

    Object(latex_title)
    red_ball = Object(RFrames(-24:0), (args...) -> circ(O, "red"), p1)
    act!(red_ball, Action(anim_rotate_around(from_rot, to_rot, O)))

    blue_ball = Object(@Frames(prev_start(), 25), (args...) -> circ(O, "blue"), p2)
    act!(blue_ball, Action(anim_rotate_around(to_rot, from_rot, red_ball)))
    path_red = Object(
        @Frames(prev_start(), stop = prev_start() + 25 - 1),
        (video, args...) -> path!(path_of_red, get_position(red_ball), "red"),
    )
    path_blue =
        Object(:same, (video, args...) -> path!(path_of_blue, pos(blue_ball), "blue"))
    string = Object(1:25, (args...) -> rad(pos(red_ball), pos(blue_ball), "black"))

    render(video; tempdirectory = "images", pathname = "")

    @test_reference "refs/dancing_circles_16_rot_trans.png" load("images/0000000016.png")
    for i in 1:25
        rm("images/$(lpad(i, 10, "0")).png")
    end
end

@testset "Dancing circles layered return Transformation" begin
    p1 = Point(100, 0)
    p2 = Point(200, 80)
    from_rot = 0.0
    to_rot = 2π
    path_of_blue = Point[]
    path_of_red = Point[]

    video = Video(500, 500)
    back = Object(1:25, ground, in_global_layer = true)
    act!(back, Action(anim_rotate(π / 2, π / 2)))
    act!(back, Action(@Frames(prev_start(), stop = 1), anim_translate(Point(25, 25))))

    Object(latex_title)
    red_ball = Object(RFrames(-24:0), (args...) -> circ_ret_trans(O, "red"), p1)
    act!(red_ball, Action(anim_rotate_around(from_rot, to_rot, O)))

    blue_ball = Object(:all, (args...) -> circ_ret_trans(O, "blue"), p2)
    act!(blue_ball, Action(anim_rotate_around(to_rot, from_rot, red_ball)))
    path_red =
        Object(1:25, (video, args...) -> path!(path_of_red, get_position(red_ball), "red"))
    path_blue =
        Object(:same, (video, args...) -> path!(path_of_blue, pos(blue_ball), "blue"))
    string = Object(1:25, (args...) -> rad(pos(red_ball), pos(blue_ball), "black"))

    render(video; tempdirectory = "images", pathname = "")

    @test_reference "refs/dancing_circles_16_rot_trans.png" load("images/0000000016.png")
    for i in 1:25
        rm("images/$(lpad(i, 10, "0")).png")
    end
end

@testset "Dancing circles layered return Transformation rescale_factor" begin
    p1 = Point(100, 0)
    p2 = Point(200, 80)
    from_rot = 0.0
    to_rot = 2π
    path_of_blue = Point[]
    path_of_red = Point[]

    video = Video(500, 500)
    back = Object(1:25, ground, in_global_layer = true)
    act!(back, Action(anim_rotate(π / 2, π / 2)))
    act!(back, Action(1:1, anim_translate(Point(25, 25))))

    Object(latex_title)
    red_ball = Object(RFrames(-24:0), (args...) -> circ_ret_trans(O, "red"), p1)
    act!(red_ball, Action(anim_rotate_around(from_rot, to_rot, O)))

    blue_ball = Object(:all, (args...) -> circ_ret_trans(O, "blue"), p2)
    act!(blue_ball, Action(anim_rotate_around(to_rot, from_rot, red_ball)))
    path_red =
        Object(1:25, (video, args...) -> path!(path_of_red, get_position(red_ball), "red"))
    path_blue =
        Object(:same, (video, args...) -> path!(path_of_blue, pos(blue_ball), "blue"))
    string = Object(1:25, (args...) -> rad(pos(red_ball), pos(blue_ball), "black"))

    render(video; tempdirectory = "images", pathname = "", rescale_factor = 0.5)

    ref_image = load("refs/dancing_circles_16_rot_trans.png")
    new_size = trunc.(Int, size(ref_image) .* 0.5)
    ref_image = imresize(ref_image, new_size)

    @test load("images/0000000016.png") == ref_image
    for i in 1:25
        rm("images/$(lpad(i, 10, "0")).png")
    end
end

@testset "Drawing grid" begin
    video = Video(500, 500)
    Object(1:40, ground_black_on_white; in_global_layer = true)
    Object(1:10, draw_grid(direction = "BL", line_gap = 25))
    Object(zero_lines(direction = "BL", line_thickness = 10))
    Object(GFrames(11:20), draw_grid(direction = "BR", line_gap = 25))
    Object(zero_lines(direction = "BR", line_thickness = 10))
    Object(RFrames(10), draw_grid(direction = "TL", line_gap = 25))
    Object(zero_lines(direction = "TL", line_thickness = 10))
    Object(RFrames(10), draw_grid(direction = "TR", line_gap = 25))
    Object(zero_lines(direction = "TR", line_thickness = 10))

    render(video; tempdirectory = "images", pathname = "")

    @test_reference "refs/grid_drawing_bl.png" load("images/0000000008.png")
    @test_reference "refs/grid_drawing_br.png" load("images/0000000018.png")
    @test_reference "refs/grid_drawing_tl.png" load("images/0000000028.png")
    @test_reference "refs/grid_drawing_tr.png" load("images/0000000038.png")
    for i in 1:40
        rm("images/$(lpad(i, 10, "0")).png")
    end
end

@testset "Use returned angle" begin
    video = Video(500, 500)
    p = Point(100, 0)

    Background(1:10, ground_black_on_white)
    circ = Object(@Frames(1, stop = prev_end()), (args...) -> circ_ret_trans(), p)
    act!(circ, Action(anim_rotate_around(0.0, 2π, O)))
    Object((args...) -> line(Point(-200, 0), Point(-200, -10 * ang(circ)), :stroke))

    render(video; tempdirectory = "images", pathname = "")

    @test_reference "refs/circle_angle.png" load("images/0000000008.png")
    for i in 1:10
        rm("images/$(lpad(i, 10, "0")).png")
    end
end

function ground_nicholas(args...)
    background("white")
    sethue("blue")
    setline(3)
end

function house_of_nicholas(; p1 = O, width = 100, color = "black")
    # sethue(color)
    #     .p1
    # .p2   .p3
    #
    # .p4   .p5
    width2 = div(width, 2)
    p2 = p1 + Point(-width2, width2)
    p3 = p1 + Point(width2, width2)
    p4 = p2 + Point(0, width)
    p5 = p3 + Point(0, width)
    line(p4, p5, :stroke)
    line(p5, p2, :stroke)
    line(p2, p4, :stroke)
    line(p4, p3, :stroke)
    line(p3, p1, :stroke)
    line(p1, p2, :stroke)
    line(p2, p3, :stroke)
    setline(8)
    line(p3, p5, :stroke)
end

@testset "House of Nicholas line_width" begin
    demo = Video(500, 500)
    Background(1:50, ground_nicholas)
    house = Object((args...) -> house_of_nicholas())
    act!(house, Action(1:25, appear(:fade_line_width)))
    act!(house, Action(RFrames(25), disappear(:fade_line_width)))

    render(demo; tempdirectory = "images", pathname = "")

    @test_reference "refs/nicholas15.png" load("images/0000000015.png")
    @test_reference "refs/nicholas25.png" load("images/0000000025.png")
    @test_reference "refs/nicholas35.png" load("images/0000000035.png")
    for i in 1:50
        rm("images/$(lpad(i, 10, "0")).png")
    end
end

@testset "Squeezing a circle using scale" begin
    demo = Video(500, 500)

    Background(1:150, ground)
    start_scale = Object((args...) -> Javis.Scale(1.0, 1.0))
    circle_obj = Object((args...) -> circ())
    act!(circle_obj, Action(1:25, anim_scale((1.0, 1.5))))
    act!(circle_obj, Action(RFrames(25), anim_scale(Javis.Scale(1.0, 1.5), (2.0, 1.0))))
    act!(circle_obj, Action(RFrames(25), anim_scale(start_scale)))
    act!(circle_obj, Action(RFrames(25), anim_scale(2.0); keep = false))
    act!(circle_obj, Action(126:150, anim_scale(0.0)))
    Object((args...) -> circ(Point(-100, 0)))

    render(demo; tempdirectory = "images", pathname = "")

    @test_reference "refs/squeeze15.png" load("images/0000000015.png")
    @test_reference "refs/squeeze25.png" load("images/0000000025.png")
    @test_reference "refs/squeeze35.png" load("images/0000000035.png")
    @test_reference "refs/squeeze65.png" load("images/0000000065.png")
    @test_reference "refs/squeeze85.png" load("images/0000000085.png")
    @test_reference "refs/squeeze110.png" load("images/0000000110.png")
    @test_reference "refs/squeeze150.png" load("images/0000000150.png")

    for i in 1:150
        rm("images/$(lpad(i, 10, "0")).png")
    end
end

function ground_opacity(args...)
    background("white")
    sethue("black")
    setopacity(0.5)
end

function square_opacity(p1, w)
    setopacity(1.0)
    sethue("red")
    rect(p1, w, w, :fill)
end

@testset "Circle/square appear opacity" begin
    demo = Video(500, 500)

    Background(1:50, ground_opacity)
    circle_obj = Object(1:42, (args...) -> circ())
    act!(circle_obj, Action(1:25, appear(:fade)))
    act!(circle_obj, Action(26:42, disappear(:fade)))

    square_obj = Object(5:50, (args...) -> square_opacity(Point(-100, 0), 60))
    act!(square_obj, Action(1:15, linear(), appear(:fade)))
    act!(square_obj, Action(RFrames(20), linear(), anim_translate(100, 50)))
    act!(square_obj, Action(GFrames(40:44), disappear(:fade)))
    # for global frames 45-50 it should still be disappeared

    render(demo; tempdirectory = "images", pathname = "")
    @test Javis.get_frames(square_obj.actions[1]) == 1:15 # globally 5:19
    @test Javis.get_frames(square_obj.actions[2]) == 16:35 # == 20:39
    @test Javis.get_frames(square_obj.actions[3]) == 36:40 # == 40:44

    @test_reference "refs/circleSquare07opacity.png" load("images/0000000007.png")
    @test_reference "refs/circleSquare25opacity.png" load("images/0000000025.png")
    @test_reference "refs/circleSquare42opacity.png" load("images/0000000042.png")
    # test that the last frame is completely white
    @test sum(load("images/0000000050.png")) ==
          RGB{Float64}(500 * 500, 500 * 500, 500 * 500)
    for i in 1:50
        rm("images/$(lpad(i, 10, "0")).png")
    end
end

@testset "Animations.jl translate()" begin
    video = Video(500, 500)
    circle_anim = Animation(
        [0, 0.3, 0.6, 1], # must go from 0 to 1
        [O, Point(150, 0), Point(150, 150), O],
        [sineio(), polyin(5), expin(8)],
    )

    Background(1:150, ground)
    circle_obj = Object((args...) -> circle(O, 25, :fill))
    act!(circle_obj, Action(1:150, circle_anim, translate()))

    render(video; tempdirectory = "images", pathname = "")
    @test_reference "refs/anim_circle020.png" load("images/0000000020.png")
    @test_reference "refs/anim_circle075.png" load("images/0000000075.png")
    @test_reference "refs/anim_circle142.png" load("images/0000000142.png")
    for i in 1:150
        rm("images/$(lpad(i, 10, "0")).png")
    end
end

@testset "Change with/without interpolation" begin
    video = Video(500, 500)

    Background(1:150, ground)
    circle_obj = Object(
        (args...; radius = 25, color = "blue") ->
            circle_with_color(O, radius, :fill, color),
    )
    act!(circle_obj, Action(1:150, change(:radius, 25 => 100)))
    act!(circle_obj, Action(50:100, change(:color, "red"); keep = false))

    render(video; tempdirectory = "images", pathname = "")
    @test_reference "refs/circle_change020.png" load("images/0000000020.png")
    @test_reference "refs/circle_change075.png" load("images/0000000075.png")
    @test_reference "refs/circle_change101.png" load("images/0000000101.png")
    for i in 1:150
        rm("images/$(lpad(i, 10, "0")).png")
    end
end

@testset "Animations.jl @warn" begin
    video = Video(500, 500)
    circle_anim = Animation(
        [0, 0.3, 0.6, 1.2], # must go from 0 to 1
        [O, Point(150, 0), Point(150, 150), O],
        [sineio(), polyin(5), expin(8)],
    )

    Background(1:2, ground)
    circle_obj = Object((args...) -> circle(O, 25, :fill))
    act!(circle_obj, Action(circle_anim, translate()))

    # warning as animation goes to 1.2 but should go to 1.0
    @test_logs (:warn,) (:warn,) render(video; pathname = "")
end

@testset "Animations.jl rotate, scale, translate" begin
    video = Video(500, 500)
    translate_anim = Animation(
        [0, 1], # must go from 0 to 1
        [O, Point(150, 0)],
        [sineio()],
    )

    translate_back_anim = Animation(
        [0, 1], # must go from 0 to 1
        [O, Point(-150, 0)],
        [sineio()],
    )

    rotate_anim = Animation(
        [0, 1], # must go from 0 to 1
        [0, 2π],
        [linear()],
    )

    color_anim = Animation(
        [0, 0.5, 1], # must go from 0 to 1
        [Lab(colorant"red"), Lab(colorant"cyan"), Lab(colorant"black")],
        [sineio(), sineio()],
    )

    Background(1:150, ground)
    circle_obj = Object((args...) -> circle(O, 25, :fill))
    act!(circle_obj, Action(1:10, sineio(), scale()))
    act!(circle_obj, Action(11:50, translate_anim, translate()))
    act!(circle_obj, Action(51:100, rotate_anim, rotate_around(Point(-150, 0))))
    act!(circle_obj, Action(101:140, translate_back_anim, translate()))
    act!(circle_obj, Action(141:150, rev(sineio()), scale()))
    act!(circle_obj, Action(1:150, color_anim, sethue()))

    bot_circle = Object((args...) -> circle(O, 25, :fill), Point(0, 200))
    act!(bot_circle, Action(sineio(), setopacity()))

    render(video; tempdirectory = "images", pathname = "")

    @test_reference "refs/animations_all_05.png" load("images/0000000005.png")
    @test_reference "refs/animations_all_25.png" load("images/0000000025.png")
    @test_reference "refs/animations_all_65.png" load("images/0000000065.png")
    @test_reference "refs/animations_all_125.png" load("images/0000000125.png")
    @test_reference "refs/animations_all_145.png" load("images/0000000145.png")

    for i in 1:150
        rm("images/$(lpad(i, 10, "0")).png")
    end
end

@testset "Rotate around center animation" begin
    rotate_anim = Animation([0.0, 1.0], [0.0, 2π], [sineio()])

    video = Video(500, 500)
    Background(1:50, ground)
    Background(1:50, (args...) -> scaleto(2))
    circle_obj = Object((args...) -> circ(Point(75, 0)))
    act!(circle_obj, Action(1:50, rotate_anim, rotate()))

    render(video; tempdirectory = "images", pathname = "")

    @test_reference "refs/rotate_center25.png" load("images/0000000025.png")
    @test_reference "refs/rotate_center45.png" load("images/0000000045.png")

    for i in 1:50
        rm("images/$(lpad(i, 10, "0")).png")
    end
end

@testset "Scaling circle" begin
    video = Video(500, 500)

    Background(1:50, ground)
    scale_obj = Object((args...) -> 2)
    circle_obj = Object((args...) -> circ())
    act!(circle_obj, Action(1:15, anim_scale(0.0, scale_obj)))
    act!(circle_obj, Action(36:50, anim_scale(0.0)))

    render(video; tempdirectory = "images", pathname = "")

    @test_reference "refs/scalingCircle07.png" load("images/0000000007.png")
    @test_reference "refs/scalingCircle25.png" load("images/0000000025.png")
    @test_reference "refs/scalingCircle42.png" load("images/0000000042.png")

    # test using appear and disappear
    video = Video(500, 500)

    Background(1:50, ground)
    Background(1:50, (args...) -> scale(2))
    circle_obj = Object((args...) -> circ())
    act!(circle_obj, Action(1:15, appear(:scale)))
    act!(circle_obj, Action(36:50, disappear(:scale)))

    render(video; tempdirectory = "images", pathname = "")

    @test_reference "refs/scalingCircle07.png" load("images/0000000007.png")
    @test_reference "refs/scalingCircle25.png" load("images/0000000025.png")
    @test_reference "refs/scalingCircle42.png" load("images/0000000042.png")
    # test that the last frame is completely white
    @test sum(load("images/0000000050.png")) ==
          RGB{Float64}(500 * 500, 500 * 500, 500 * 500)
    for i in 1:50
        rm("images/$(lpad(i, 10, "0")).png")
    end
end

@testset "animating text" begin
    # does only test that it doesn't fail but I wasn't able to test this on all platforms
    # with using reference tests
    video = Video(400, 300)

    Background(1:100, ground)
    Background(1:100, (args...) -> fontsize(30))
    title = Object(1:100, (args...) -> text("Hello Stream!", -50, 50; halign = :centre))
    act!(title, Action(1:15, sineio(), appear(:draw_text)))
    act!(title, Action(76:100, sineio(), disappear(:draw_text)))
    title2 = Object(
        1:100,
        (args...) -> text("Hello World!", -50, -100; halign = :wrong, valign = :wrong),
    )
    act!(title2, Action(1:15, sineio(), appear(:draw_text)))
    act!(title2, Action(76:100, sineio(), disappear(:draw_text)))

    render(video; tempdirectory = "images", pathname = "")

    img07 = load("images/$(lpad(7, 10, "0")).png")
    img30 = load("images/$(lpad(30, 10, "0")).png")
    img82 = load("images/$(lpad(30, 10, "0")).png")

    # does only test that it doesn't fail but I wasn't able to test this on all platforms
    # with using reference tests
    video = Video(400, 300)

    Background(1:100, ground)
    Background(1:100, (args...) -> fontsize(30))
    title = Object(1:100, (args...) -> text("Hello Stream!", -50, 50; halign = :center))
    act!(title, Action(1:15, sineio(), appear(:draw_text)))
    act!(title, Action(76:100, sineio(), disappear(:draw_text)))
    title2 = Object(1:100, (args...) -> text("Hello World!", -50, -100))
    act!(title2, Action(1:15, sineio(), appear(:draw_text)))
    act!(title2, Action(76:100, sineio(), disappear(:draw_text)))

    render(video; tempdirectory = "images", pathname = "")

    img_other07 = load("images/$(lpad(7, 10, "0")).png")
    img_other30 = load("images/$(lpad(30, 10, "0")).png")
    img_other82 = load("images/$(lpad(30, 10, "0")).png")

    @test img07 == img_other07
    @test img30 == img_other30
    @test img82 == img_other82

    for i in 1:100
        rm("images/$(lpad(i, 10, "0")).png")
    end
end

@testset "Following a path" begin
    video = Video(800, 600)
    Background(1:180, ground)

    anim = Animation([0, 1], [0.0, 2.0], [sineio()])

    color_anim = Animation(
        [0, 0.5, 1], # must go from 0 to 1
        [Lab(colorant"red"), Lab(colorant"cyan"), Lab(colorant"red")],
        [sineio(), sineio()],
    )

    actions = [
        Action(1:150, anim, follow_path(star(O, 100))),
        Action(1:150, color_anim, sethue()),
    ]

    objects = [
        Object(frame_start:(frame_start + 149), (args...) -> star(O, 20, 5, 0.5, 0, :fill)) for frame_start in 1:7:22
    ]

    act!(objects, actions)

    render(video; tempdirectory = "images", pathname = "")

    @test_reference "refs/followPath10.png" load("images/0000000010.png")
    @test_reference "refs/followPath30.png" load("images/0000000030.png")
    @test_reference "refs/followPath100.png" load("images/0000000100.png")
    @test_reference "refs/followPath160.png" load("images/0000000160.png")

    # Following along a bezier path
    function simple_bezier()
        P1 = Point(-300, 0)
        CP1 = Point(-200, -200)
        CP2 = Point(200, -200)
        P2 = Point(300, 0)

        beziersegment = BezierPathSegment(P1, CP1, CP2, P2)
        beziertopoly(beziersegment)
    end

    video = Video(800, 600)
    Background(1:180, ground)

    anim = Animation([0, 1], [0.0, 1.0], [sineio()])

    color_anim = Animation(
        [0, 0.5, 1], # must go from 0 to 1
        [Lab(colorant"red"), Lab(colorant"cyan"), Lab(colorant"red")],
        [sineio(), sineio()],
    )

    actions =
        (frame_start) -> [
            Action(1:10, appear(:fade)),
            Action(
                11:150,
                anim,
                follow_path(
                    simple_bezier() .- (simple_bezier()[1] + Point(0, 3 * frame_start));
                    closed = false,
                ),
            ),
            Action(1:150, color_anim, sethue()),
        ]

    objects = [
        Object(
            frame_start:(frame_start + 149),
            (args...) -> star(O, 20, 5, 0.5, 0, :fill),
            simple_bezier()[1] + Point(0, 3 * frame_start),
        ) for frame_start in 1:7:22
    ]

    for (object, frame_start) in zip(objects, 1:7:22)
        act!(object, actions(frame_start))
    end

    render(video; tempdirectory = "images", pathname = "")

    @test_reference "refs/followPathBezier1.png" load("images/0000000001.png")
    @test_reference "refs/followPathBezier5.png" load("images/0000000005.png")
    @test_reference "refs/followPathBezier10.png" load("images/0000000010.png")
    @test_reference "refs/followPathBezier30.png" load("images/0000000030.png")
    @test_reference "refs/followPathBezier100.png" load("images/0000000100.png")
    @test_reference "refs/followPathBezier160.png" load("images/0000000160.png")

    for i in 1:180
        rm("images/$(lpad(i, 10, "0")).png")
    end
end

@testset "Change" begin
    function object(p, radius, color = "black")
        sethue(color)
        circle(p, radius, :fill)
        return p
    end

    anim =
        Animation([0.0, 0.5, 1.0], [O, Point(100, 0), Point(50, 50)], [sineio(), sineio()])

    rev_anim =
        Animation([0.0, 0.5, 1.0], [Point(50, 50), Point(100, 0), O], [sineio(), sineio()])

    myvideo = Video(500, 500)
    double_rotation = Animation([0.0, 1.0], [0.0, 4π], [linear()])

    Background(1:100, ground)
    circ1 = Object((args...; radius = 25) -> object(Point(100, 0), radius, "red"))
    act!(circ1, Action(double_rotation, rotate()))
    act!(circ1, Action(1:50, change(:radius, 25 => 0)))
    act!(circ1, Action(51:100, change(:radius, 0 => 25)))

    circ2 = Object((args...; point = O) -> object(point, 25, "green"))
    act!(circ2, Action(1:50, anim, change(:point)))
    act!(circ2, Action(51:100, rev_anim, change(:point)))

    render(myvideo; tempdirectory = "images")

    @test_reference "refs/changeKeyword1.png" load("images/0000000001.png")
    @test_reference "refs/changeKeyword25.png" load("images/0000000025.png")
    @test_reference "refs/changeKeyword40.png" load("images/0000000040.png")
    @test_reference "refs/changeKeyword70.png" load("images/0000000070.png")
    @test_reference "refs/changeKeyword90.png" load("images/0000000090.png")

    for i in 1:100
        rm("images/$(lpad(i, 10, "0")).png")
    end
    rm("images/palette.png")
end

@testset "test default kwargs" begin
    video = Video(500, 500)
    Background(1:10, ground)
    Object(1:10, (args...) -> circle(O, 50, :fill))
    pathname = render(video)
    path, ext = splitext(pathname)
    @test ext == ".gif"
    @test isfile(pathname)
    rm(pathname)
end

@testset "test @error .mp3" begin
    video = Video(500, 500)
    Background(1:10, ground)
    Object(1:10, (args...) -> circle(O, 50, :fill))
    @test_logs (:error,) render(video; pathname = "test.mp3")
end

@testset "test layer vs nonlayer actions" begin
    function act_all!(a)
        act!(a, Action(1:70, anim_translate(Point(50, 50))))
        act!(a, Action(71:140, anim_translate(Point(50, 50), Point(-50, -50))))
        act!(a, Action(71:140, anim_rotate(2π)))
        act!(a, Action(140:210, anim_translate(Point(50, 50))))
    end

    n_frames = 257
    nolayer_video = Video(500, 500)
    Background(1:n_frames, ground)
    circ = Object(JBox(O, 20, 20, action = :fill, color = "white"))
    act_all!(circ)
    render(nolayer_video, tempdirectory = "images/without_layer", pathname = "")

    layer_video = Video(500, 500)
    Background(1:n_frames, ground)
    l1 = @JLayer 1:n_frames begin
        Object(JBox(O, 20, 20, action = :fill, color = "white"))
    end
    act_all!(l1)
    render(layer_video, tempdirectory = "images/with_layer", pathname = "")

    for frame in [8, 16, 32, 64, 128, 256]
        png_name = lpad(string(frame), 10, "0")
        @test_reference "images/with_layer/$png_name.png" load(
            "images/without_layer/$png_name.png",
        )
    end

    for i in 1:n_frames
        rm("images/with_layer/$(lpad(i, 10, "0")).png")
        rm("images/without_layer/$(lpad(i, 10, "0")).png")
    end
    rm("images/with_layer/", recursive = true)
    rm("images/without_layer/", recursive = true)

    n_frames = 140
    nolayer_video = Video(500, 500)
    Background(1:n_frames, ground)
    circ = Object(JBox(O, 20, 20, action = :fill, color = "white"))

    act!(circ, Action(1:140, anim_translate(Point(0, 50))))
    act!(circ, Action(1:140, anim_translate(Point(50, 0))))

    render(nolayer_video, tempdirectory = "images/without_layer", pathname = "")

    layer_video = Video(500, 500)
    Background(1:n_frames, ground)
    l1 = @JLayer 1:n_frames begin
        Object(JBox(O, 20, 20, action = :fill, color = "white"))
    end
    act!(l1, Action(1:140, anim_translate(Point(50, 50))))
    render(layer_video, tempdirectory = "images/with_layer", pathname = "")

    for frame in [8, 16, 24, 32, 40, 48]
        png_name = lpad(string(frame), 10, "0")
        @test_reference "images/with_layer/$png_name.png" load(
            "images/without_layer/$png_name.png",
        )
    end

    for i in 1:140
        rm("images/with_layer/$(lpad(i, 10, "0")).png")
        rm("images/without_layer/$(lpad(i, 10, "0")).png")
    end
    rm("images/with_layer/", recursive = true)
    rm("images/without_layer/", recursive = true)
end


text:
using EngEconomics, Roots

# Given
# a = parasailing
# b = kayaking
aInit = -100000
aAnnual = 15000
aN = 15

bInit = -10000
bAnnual = 2000
bN = 15
MARR = 0.10

N = aN

# Determine the better alternative
ΔInit = aInit - bInit
ΔAnnual = aAnnual - bAnnual
ΔIRR(x) = ΔInit + ΔAnnual * seriesPresentAmountFactor(x, N)
ΔIRR_ans = find_zero(ΔIRR, 0.05)

if ΔIRR_ans > MARR
	println("Choose Parasailing")
else
	println("Choose Kayaking")
end


text:
# Write JavaSeis with Julia package TeaSeis.jl

# Run from Julia REPL
# julia> include("C:/Users/xinfa/Documents/code/pieseis/examples/write_file.jl")
# Run from command line (Linux terminal or Windows command prompt)
# C:\Users\xinfa\AppData\Local\Julia-1.0.1\bin\julia.exe C:\Users\xinfa\Documents\code\pieseis\examples\write_file.jl

println("start")

push!(LOAD_PATH, "C:/Users/xinfa/Documents/code/TeaSeis.jl/src")
using TeaSeis
println("done load")

filename = "C:/Users/xinfa/Documents/181116_julia.js"
io = jsopen(filename, "w", axis_lengths=[101, 101, 101], dataformat=Int16)
println("done open")

trcs, hdrs = allocframe(io)

map(i->set!(prop(io, stockprop[:TRC_TYPE]), hdrs, i, tracetype[:live]), 1:size(io,2))
map(i->set!(prop(io, stockprop[:TRACE]   ), hdrs, i, i               ), 1:size(io,2))
map(i->set!(prop(io, stockprop[:FRAME]   ), hdrs, i, 1               ), 1:size(io,2))
println("done write headers")

using Random
rand!(trcs)
writeframe(io, trcs, hdrs)
println("done write traces")

close(io)

max_stars_repo_name: JoeXinfa/pieseis

max_issues_repo_name: JoeXinfa/pieseis

max_forks_repo_name: JoeXinfa/pieseis


text:
module Telegram
include("c_api.jl")
include("telegram_common.jl")
end


text:
push!(LOAD_PATH, joinpath(@__DIR__, "..", ".."))

import ClimaCore.Geometry, LinearAlgebra, UnPack
import ClimaCore:
    Fields,
    Domains,
    Topologies,
    Meshes,
    DataLayouts,
    Operators,
    Geometry,
    Spaces

using OrdinaryDiffEq: OrdinaryDiffEq, ODEProblem, solve, SSPRK33

using Logging: global_logger
using TerminalLoggers: TerminalLogger
global_logger(TerminalLogger())

const FT = Float64

# https://github.com/CliMA/CLIMAParameters.jl/blob/master/src/Planet/planet_parameters.jl#L5
const MSLP = FT(1e5) # mean sea level pressure
const grav = FT(9.8) # gravitational constant
const R_d = FT(287.058) # R dry (gas constant / mol mass dry air)
const γ = FT(1.4) # heat capacity ratio
const C_p = FT(R_d * γ / (γ - 1)) # heat capacity at constant pressure
const C_v = FT(R_d / (γ - 1)) # heat capacit at constant volume
const R_m = FT(R_d) # moist R, assumed to be dry

domain = Domains.IntervalDomain(
    Geometry.ZPoint{FT}(0.0),
    Geometry.ZPoint{FT}(30e3),
    boundary_tags = (:bottom, :top),
)
#mesh = Meshes.IntervalMesh(domain, Meshes.ExponentialStretching(7.5e3); nelems = 30)
mesh = Meshes.IntervalMesh(domain; nelems = 30)

cspace = Spaces.CenterFiniteDifferenceSpace(mesh)
fspace = Spaces.FaceFiniteDifferenceSpace(cspace)

# https://github.com/CliMA/Thermodynamics.jl/blob/main/src/TemperatureProfiles.jl#L115-L155
# https://clima.github.io/Thermodynamics.jl/dev/TemperatureProfiles/#DecayingTemperatureProfile
function decaying_temperature_profile(
    z;
    T_virt_surf = FT(280.0),
    T_min_ref = FT(230.0),
)
    # Scale height for surface temperature
    H_sfc = R_d * T_virt_surf / grav
    H_t = H_sfc

    z′ = z / H_t
    tanh_z′ = tanh(z′)

    ΔTv = T_virt_surf - T_min_ref
    Tv = T_virt_surf - ΔTv * tanh_z′

    ΔTv′ = ΔTv / T_virt_surf
    p =
        MSLP * exp(
            (
                -H_t *
                (z′ + ΔTv′ * (log(1 - ΔTv′ * tanh_z′) - log(1 + tanh_z′) + z′))
            ) / (H_sfc * (1 - ΔTv′^2)),
        )
    ρ = p / (R_d * Tv)
    ρθ = ρ * Tv * (MSLP / p)^(R_m / C_p)
    return (ρ = ρ, ρθ = ρθ)
end

Π(ρθ) = C_p * (R_d * ρθ / MSLP)^(R_m / C_v)
Φ(z) = grav * z

function discrete_hydrostatic_balance!(ρ, w, ρθ, Δz::FT, _grav::FT, Π::Function)
    # compute θ such that
    #   I(θ)[i+1/2] = -g / ∂f(Π(ρθ))
    # discretely, then set
    #   ρ = ρθ/θ
    for i in 1:(length(ρ) - 1)
        #  ρ[i+1] = ρθ[i+1]/(-2Δz*_grav/(Π(ρθ[i+1]) - Π(ρθ[i])) - ρθ[i]/ρ[i])
        ρ[i + 1] =
            ρθ[i + 1] /
            (-2 * _grav / ((Π(ρθ[i + 1]) - Π(ρθ[i])) / Δz) - ρθ[i] / ρ[i])

        ρ[i + 1] =
            ρθ[i + 1] /
            (1 / ((-2 * _grav) * (Π(ρθ[i + 1]) - Π(ρθ[i]))Δz) - ρθ[i] / ρ[i])

        ∂Π∂z = (Π(ρθ[i + 1]) - Π(ρθ[i])) / Δz
    end
end

zc = Fields.coordinate_field(cspace)
Yc = decaying_temperature_profile.(zc.z)
w = Geometry.Cartesian3Vector.(zeros(FT, fspace))

Y_init = copy(Yc)
w_init = copy(w)

Y = (Yc, w)

function tendency!(dY, Y, _, t)
    (Yc, w) = Y.x
    (dYc, dw) = dY.x

    If = Operators.InterpolateC2F()
    ∂ = Operators.DivergenceF2C(
        bottom = Operators.SetValue(Geometry.Cartesian3Vector(zero(FT))),
        top = Operators.SetValue(Geometry.Cartesian3Vector(zero(FT))),
    )
    ∂f = Operators.GradientC2F()
    B = Operators.SetBoundaryOperator(
        bottom = Operators.SetValue(Geometry.Cartesian3Vector(zero(FT))),
        top = Operators.SetValue(Geometry.Cartesian3Vector(zero(FT))),
    )

    @. dYc.ρ = -(∂(w * If(Yc.ρ)))
    @. dYc.ρθ = -(∂(w * If(Yc.ρθ)))
    @. dw = B(
        Geometry.CartesianVector(
            -(If(Yc.ρθ / Yc.ρ) * ∂f(Π(Yc.ρθ))) - ∂f(Φ(zc.z)),
        ),
    )
    return dY
end

using RecursiveArrayTools

Y = ArrayPartition(Yc, w)
dY = tendency!(similar(Y), Y, nothing, 0.0)

Δt = 1.0
ndays = 10

# Solve the ODE operator
prob = ODEProblem(tendency!, Y, (0.0, 60 * 60 * 24 * ndays))
sol = solve(
    prob,
    SSPRK33(),
    dt = Δt,
    saveat = 60 * 60, # save every hour
    progress = true,
    progress_message = (dt, u, p, t) -> t,
);

ENV["GKSwstype"] = "nul"
import Plots
Plots.GRBackend()

dirname = "hydrostatic"
path = joinpath(@__DIR__, "output", dirname)
mkpath(path)

z_centers = parent(Fields.coordinate_field(cspace))
z_faces = parent(Fields.coordinate_field(fspace))

function hydrostatic_plot(u; title = "", size = (1024, 600))
    sub_plt1 = Plots.plot(
        parent(Y_init.ρ),
        z_centers,
        marker = :circle,
        xlabel = "ρ",
        label = "T=0",
    )
    sub_plt1 = Plots.plot!(sub_plt1, parent(u.x[1].ρ), z_centers, label = "T")

    sub_plt2 = Plots.plot(
        parent(w_init),
        z_faces,
        marker = :circle,
        xlim = (-0.2, 0.2),
        xlabel = "ω",
        label = "T=0",
    )
    sub_plt2 = Plots.plot!(sub_plt2, parent(u.x[2]), z_faces, label = "T")

    sub_plt3 = Plots.plot(
        parent(Y_init.ρθ),
        z_centers,
        marker = :circle,
        xlabel = "ρθ",
        label = "T=0",
    )
    sub_plt3 = Plots.plot!(sub_plt3, parent(u.x[1].ρθ), z_centers, label = "T")

    return Plots.plot(
        sub_plt1,
        sub_plt2,
        sub_plt3,
        title = title,
        layout = (1, 3),
        size = size,
    )
end

anim = Plots.@animate for (i, u) in enumerate(sol.u)
    hydrostatic_plot(u, title = "Hour $(i)")
end
Plots.mp4(anim, joinpath(path, "hydrostatic.mp4"), fps = 10)

Plots.png(hydrostatic_plot(sol[end]), joinpath(path, "hydrostatic_end.png"))

function linkfig(figpath, alt = "")
    # buildkite-agent upload figpath
    # link figure in logs if we are running on CI
    if get(ENV, "BUILDKITE", "") == "true"
        artifact_url = "artifact://$figpath"
        print("\033]1338;url='$(artifact_url)';alt='$(alt)'\a\n")
    end
end

linkfig("output/$(dirname)/hydrostatic_end.png", "Hydrostatic End")


text:
# routines related to the ParallelStencil backend

export initialize_backend!

"""
    start_backend(Backend{BackendParallelStencil}; Scalar=Float64, dim=1)

Starts the `ParallelStencil` backend for the `Scalar`, which loads ParallelStencil (and, if requested, MPI) for the correct number of dimensions
"""
function check_backend(b::Backend{BackendParallelStencil}; dim=1, Scalar=Float64)
    # Check whether MPI is loaded
    if (b.mpi==true && !isdefined(Main, :ImplicitGlobalGrid))
        error("ImplicitGlobalGrid is not loaded; ensure it is loaded first with: using ImplicitGlobalGrid")
    end
    if (b.mpi==true && !isdefined(Main, :MPI))
        error("MPI is not loaded; ensure it is loaded first with: using MPI")
    end

    # Check whether ParallelStencil is loaded
    if !isdefined(Main, :ParallelStencil)
        error("ParallelStencil is not loaded; ensure it is loaded first with: using ParallelStencil")
    end

    # Load GPU or CPU backend
    if ((b.arch==:GPU) || (b.arch==:CUDA))
        # Currently, ParallelStencil, only supports CUDA
        if !isdefined(Main, :CUDA)
            error("CUDA is not loaded; ensure it is loaded first with: using CUDA")
        end
    end

    return nothing
end

function initialize_backend!(grid, b::Backend{BackendParallelStencil}, stencilwidth, opts, fields::NamedTuple)

    N_vec = ones(Int,3);
    D = grid.dim

    Nl,dims = localfromglobalsize(backend.mpi, grid.N, opts)      # local size using MPI (w/out halo)
    N_vec[1:D] = collect(Nl)

    initmpi=false
    if !MPI.Initialized()
        initmpi=true
    end

    # Halo
    overlap  = zeros(Int64,3)
    overlap[1:D] .= 2*stencilwidth
    Nl_halo  = N_vec .+ overlap

    # Periodic?
    period  = zeros(Int64,3)
    period[findall(grid.topology.==Periodic)] .= 1

    IGG_data = init_global_grid(Nl_halo[1], Nl_halo[2], Nl_halo[3],
                                                            init_MPI=initmpi,
                                                            quiet=false,
                                                            overlapx=overlap[1], overlapy=overlap[2], overlapz=overlap[3],
                                                            periodx=period[1], periody=period[2], periodz=period[3];
                                                            opts...);

    IGG_data = ImplicitGlobalGrid.get_global_grid()
    grid.IGG = IGG_data;                    # store

    # Store local grid dimensions
    grid.Nl = (Nl_halo[1:D]...,)            # local with halo
    grid.Ng = (IGG_data.nxyz_g[1:D]...,)    # global + halo

    # Compute local indices of grid in global grid (not including the halo)
    ind_s = IGG_data.coords.*(IGG_data.nxyz .- IGG_data.overlaps)
    ind_e = ind_s[1:D] .+ Nl;
    ind_local = ();
    for i=1:D
        ind_local = (ind_local..., ind_s[i]+1:ind_e[i])
    end
    grid.ind_local = ind_local

    # Update corners accordingly, which provides an easy way to iterate over the local portions of the grid, with/without ghost/halo points
    size = ind_e .- ind_s[1:D];
    grid.corners = ( lower=(CartesianIndex( (ones(Int64,D) .+ stencilwidth)...,)),
                    upper=(CartesianIndex( (zeros(Int64,D) .+ stencilwidth + size)...,)),
                    size=(size...,))
    grid.ghostcorners = ( lower=(CartesianIndex( (ones(Int64,D))...,)),
                        upper=(CartesianIndex( (zeros(Int64,D) .+ (2 .* stencilwidth) + size)...,)),
                        size =(Nl_halo[1:D]...,))

    # Update the 1D coordinate vectors
    local_coordinates!(grid, grid.backend)

    return nothing

end

"""
    initialize_fields!(grid, b::Backend{BackendParallelStencil}, fields::NamedTuple)

This initializes fields to a grid structure, in case the `ParallelStencil` is employed.
The names of the fields are a NamedTuple, which specifies the initial value (constant) value of the arrays, such as `fields=(T=0, P=1, Ci=2)`.

Note: typically, this routine is not called directly but as part of setting up the grid

# Example

```
julia> @init_backend(ParallelStencil, Threads, false, Float64);
julia> initialize_fields!(grid, b::Backend{BackendParallelStencil}, fields::NTuple)
```

"""
function initialize_fields!(grid, b::Backend{BackendParallelStencil, FT}, fields::NamedTuple) where FT

    names = keys(fields)
    for ifield=1:length(fields)

        # Created named tuple for new field. We use @eval here, to avoid
        name = names[ifield]
        value = fields[ifield]

        # Create a string with the expression to be evaluated (as calling @ones directly induces a compilation error)
        if b.arch == :Threads
            new_field = ParallelStencil.ParallelKernel.@ones_threads(grid.Nl)*value
        elseif b.arch == :CUDA
            new_field = ParallelStencil.ParallelKernel.@ones_cuda(grid.Nl)*value
        else
            error("Not yet implemented")
        end

        # Add to Tuple
        grid.fields = add_field(grid.fields,name,new_field)
    end
end

"""
    globalfromlocalsize(mpi, localsize::NTuple, opts::NTuple, stencilwidth::Int, b::Backend{BackendParallelStencil, FT})

Compute global from local size in case we use the ParallelStencil backend
"""
function globalfromlocalsize(size, localsize, opts, stencilwidth, topology, b::Backend{BackendParallelStencil, FT}) where FT

    if !b.mpi & !isempty(localsize)
        size = localsize   # not using MPI/IGG
    elseif b.mpi & !isempty(localsize)
        dim = length(localsize);
        dims = zeros(Int64,3)

        dims[dim+1:end] .= 1;       # fix
        if typeof(opts) <: Dict # process command-line options
            dims[1] = get(opts, :dimx, dims[1]);
            if dim>1;   dims[2] = get(opts, :dimy, dims[2]);    end
            if dim>2;   dims[3] = get(opts, :dimz, dims[3]);    end
        end
        size_vec = ones(Int,3);
        size_vec[1:dim] .= localsize;

        # Periodic boundary
        period  = zeros(Int64,3)
        period[findall(topology.==Periodic)] .= 1

        # Use MPI to compute processor partitioning (dims contain the # of processors in each direction)
        overlap = stencilwidth*2

        # Initialize the global grid
        init_global_grid(size_vec[1],size_vec[2],size_vec[3],
                        dimx=dims[1], dimy=dims[2], dimz=dims[3],
                        overlapx=overlap, overlapy=overlap, overlapz=overlap,
                        periodx=period[1],  periody=period[2], periodz=period[3],
                        init_MPI=false, quiet=true)

        size_global = [nx_g(), ny_g(), nz_g()] .- overlap

        size = (size_global[1:dim]...,)
    else
        size = size
    end

    return size
end


text:
export SpeedMappingOpt

struct SpeedMappingOpt end

function __map_optimizer_args(prob::OptimizationProblem, opt::SpeedMappingOpt;
    cb=nothing,
    maxiters::Union{Number, Nothing}=nothing,
    maxtime::Union{Number, Nothing}=nothing,
    abstol::Union{Number, Nothing}=nothing,
    reltol::Union{Number, Nothing}=nothing,
    kwargs...)

    # add optimiser options from kwargs
    mapped_args = (; kwargs...)

    if !(isnothing(maxiters))
        @info "maxiters defines maximum gradient calls for $(opt)"
        mapped_args = (; mapped_args..., maps_limit = maxiters)
    end

    if !(isnothing(maxtime))
        mapped_args = (; mapped_args..., time_limit = maxtime)
    end

    if !isnothing(abstol)
        @warn "common abstol is currently not used by $(opt)"
    end

    if !isnothing(reltol)
        @warn "common reltol is currently not used by $(opt)"
    end

    return mapped_args
end

function __solve(prob::OptimizationProblem, opt::SpeedMappingOpt;
                 maxiters::Union{Number, Nothing} = nothing,
                 maxtime::Union{Number, Nothing} = nothing,
                 abstol::Union{Number, Nothing}=nothing,
                 reltol::Union{Number, Nothing}=nothing,
                 progress = false,
                 kwargs...)
    local x

    maxiters = _check_and_convert_maxiters(maxiters)
    maxtime = _check_and_convert_maxtime(maxtime)

    f = instantiate_function(prob.f,prob.u0,prob.f.adtype,prob.p)

    _loss = function(θ)
        x = f.f(θ, prob.p)
        return first(x)
    end

    if isnothing(f.grad)
        @info "SpeedMapping's ForwardDiff AD backend is used to calculate the gradient information."
    end

    opt_args = _map_optimizer_args(prob,opt, maxiters=maxiters, maxtime=maxtime, abstol=abstol, reltol=reltol; kwargs...)

    t0 = time()
    opt_res = SpeedMapping.speedmapping(prob.u0; f = _loss, g! =f.grad, lower=prob.lb, upper=prob.ub, opt_args...)
    t1 = time()
    opt_ret = Symbol(opt_res.converged)

    SciMLBase.build_solution(prob, opt, opt_res.minimizer, _loss(opt_res.minimizer); original=opt_res, retcode=opt_ret)
end


text:
import Base:start,next,done,iteratorsize,iteratoreltype,eltype,length,size

struct Exterate{I}
    itr::I
    len::Int
    function Exterate(itr::I,len::Int) where I
        new{I}(itr,len)
    end
end

exterate(itr::I,len::Int) where I = Exterate(itr,len)

function eltype(::Type{Exterate{I}}) where I
    eltype(I)
end

done(z::Exterate,state) = state[1] > z.len

_start(itr) = (1,start(itr),start(itr))
start(v::Exterate) = _start(v.itr)

length(z::Exterate) = z.len

@inline function next(z::Exterate{I},state::S)::Tuple{eltype(I),S} where {I,S}
    if done(z.itr,state[2])
        ret = next(z.itr,state[3])
        return ret[1],(state[1]+1,ret[2],state[3])
    end
    ret = next(z.itr,state[2])
    return ret[1],(state[1]+1,ret[2],state[2])
end

max_stars_repo_name: sebastianpech/SpecialIterators

max_issues_repo_name: sebastianpech/SpecialIterators

max_forks_repo_name: sebastianpech/SpecialIterators


text:
function solve(bigMpb::SolverBigM{T}) where T
    pb = bigMpb.model
    e = ones(pb.c)

    ## Compute initial primal dual starting point, and big M initial values
    x0 = get_symmetricprimalpt(pb.A, pb.b)
    z = PointPrimalDual(x0.dims, length(x0.vec), bigMpb.model.m, T, Dense{T})

    printstyled("Begining solve...\n", color=:yellow)

    ## start solve
    it = 0
    keepon = true
    while it < 5 && keepon
        printstyled("--- iteration $it\n", color=:yellow)

        ## Setting problem with new values of M1, M2
        @assert norm(evaluate(pb.A, x0) - pb.b) < 1e-10

        z̄, M1, M2 = derivestartingpoint(bigMpb, x0)
        printstyled("M1, M2 = $M1, $M2\n", color=:yellow)

        update_parameters!(bigMpb, z̄, M1, M2)

        @show min(bigMpb.model_bigM, z̄.x)
        @show min(bigMpb.model_bigM, z̄.s)

        @assert get_primfeaserr(bigMpb.model_bigM, z̄) < 1e-10
        @assert get_dualfeaserr(bigMpb.model_bigM, z̄) < 1e-10

        ## Solving
        z̄opt = solve(bigMpb.model_bigM, z̄)

        ξ1 = z̄opt.x.vec[end-1]
        ξ2 = z̄opt.x.vec[end]
        η = z̄opt.y[end]
        σ1 = z̄opt.s.vec[end-1]
        σ2 = z̄opt.s.vec[end]

        println("ξ1 = ", ξ1)
        println("ξ2 = ", ξ2)
        println("η = ", η)
        println("σ1 = ", σ1)
        println("σ2 = ", σ2)

        extractinitpbsol!(z, z̄opt)
        printstyled("prim feaserr: ", get_primfeaserr(bigMpb.model, z), "\n", color=:yellow)
        printstyled("dual feaserr: ", get_primfeaserr(bigMpb.model, z), "\n", color=:yellow)
        printstyled("prim obj: ", get_primobj(bigMpb.model, z), "\n", color=:yellow)
        printstyled("dual obj: ", get_dualobj(bigMpb.model, z), "\n", color=:yellow)
        printstyled("min_K(x): ", min(bigMpb.model_bigM, z̄.x), "\n", color=:yellow)
        printstyled("min_K(s): ", min(bigMpb.model_bigM, z̄.s), "\n", color=:yellow)

        if ξ2 < 1e-10 && η < 1e-10
            keepon = false
            printstyled("Optimal solution found.\n", color=:yellow)
            return z

        elseif bigMpb.M1 > 1e5
            keepon = false
            printstyled("Dual problem is likely to be infeasible.\n", color=:yellow)
            return z

        elseif bigMpb.M2 > 1e5
            keepon = false
            printstyled("Primal problem is likely to be infeasible.\n", color=:yellow)
            return z
        end

        x0 = z.x
        # if η > 0
        #     M1 *= 2
        # end
        #
        # if ξ2 > 0
        #     M2 *= 2
        # end

        it += 1
    end

    return false
end


text:
#!/usr/bin/env bash
    #=
    exec julia --project="$(realpath $(dirname $0))/../" "${BASH_SOURCE[0]}" "$@" -e 'include(popfirst!(ARGS))' \
    "${BASH_SOURCE[0]}" "$@"
    =#

#=
	IntegralArray{T, N, A} <: AbstractArray{T, N}

Rectangle features can be computed very rapidly using an intermediate representation for the image, which we call the integral image.
The integral image at location $x,y$ contains the sum of the pixels above and to the left of $x,y$ inclusive.
Original    Integral
+--------   +------------
| 1 2 3 .   | 1  3  6 .
| 4 5 6 .   | 5 12 21 .
| . . . .   | . . . . .
=#
struct IntegralArray{T, N, A} <: AbstractArray{T, N}
	data::A
end

"""
	to_integral_image(img_arr::AbstractArray) -> AbstractArray

Calculates the integral image based on this instance's original image data.

# Arguments

- `img_arr::AbstractArray`: Image source data

# Returns

 - `integral_image_arr::AbstractArray`: Integral image for given image
 """
function to_integral_image(img_arr::AbstractArray{T, N}) where {T, N}
	array_size = size(img_arr)
    integral_image_arr = Array{Images.accum(eltype(img_arr))}(undef, array_size)
    sd = coords_spatial(img_arr)
    cumsum!(integral_image_arr, img_arr; dims=first(sd))
    for i = 2:length(sd)
        cumsum!(integral_image_arr, integral_image_arr; dims=sd[i])
    end

    return Array{T, N}(integral_image_arr)
end

LinearIndices(A::IntegralArray) = Base.LinearFast()
@inline size(A::IntegralArray) = size(A.data)
@inline getindex(A::IntegralArray, i::Int...) = A.data[i...]
@inline getindex(A::IntegralArray, ids::Tuple...) = getindex(A, first(ids)...)

"""
	sum_region(
		integral_image_arr::AbstractArray,
		top_left::Tuple{Int,Int},
		bottom_right::Tuple{Int,Int}
	) -> Number

# Arguments

- `integral_image_arr::AbstractArray`: The intermediate Integral Image
- `top_left::Tuple{Integer, Integer}`: (x,y) of the rectangle's top left corner
- `bottom_right::Tuple{Integer, Integer}`: (x,y) of the rectangle's bottom right corner

# Returns

- `sum::Number` The sum of all pixels in the given rectangle defined by the parameters top_left and bottom_right
"""
function sum_region(
	integral_image_arr::AbstractArray,
	top_left::Tuple{T, T},
	bottom_right::Tuple{T, T}
) where T <: Integer
    _1 = one(T)
	_0 = zero(0)
	sum = integral_image_arr[last(bottom_right), first(bottom_right)]
    sum -= first(top_left) > _1 ? integral_image_arr[last(bottom_right), first(top_left) - _1] : _0
    sum -= last(top_left) > _1 ? integral_image_arr[last(top_left) - _1, first(bottom_right)] : _0
    sum += last(top_left) > _1 && first(top_left) > _1 ? integral_image_arr[last(top_left) - _1, first(top_left) - _1] : _0
    return sum
end


text:
mutable struct AddConsoleTab <: GtkScrolledWindow
    handle::Ptr{Gtk.GObject}
    textview::GtkTextView
    worker_idx::Int

    function AddConsoleTab(main_window)
        vbox = GtkBox(:v)
        sc = GtkScrolledWindow(vbox)

        b = GtkButton("Add Console")
        b.expand[Bool] = false
        signal_connect(add_console_button_press_cb, b, "button-press-event",
        Cint, (Ptr{Gtk.GdkEvent},), false, main_window)
        push!(vbox, b)

        buffer = GtkTextBuffer()
        textview = GtkTextView(buffer)
        textview.editable[Bool] = false
        push!(vbox, textview)

        n = new(sc.handle, textview, -1)#invalid index so we don't think it's a Console
        Gtk.gobject_move_ref(n, sc)
    end
end

function add_plus_button(console_mng::ConsoleManager)
    tab = AddConsoleTab(console_mng.main_window)
    push!(console_mng, tab)
    Gtk.GAccessor.tab_label_text(console_mng, tab, "+")
end

@guarded (INTERRUPT) function add_console_button_press_cb(widgetptr::Ptr, eventptr::Ptr, user_data)
    button = convert(GtkButton, widgetptr)
    main_window = user_data
    add_remote_console(main_window, main_window.console_manager.top_module)
    return INTERRUPT
end

function on_console_mng_switch_page(cm::ConsoleManager, tab::AddConsoleTab)
    #
    b = tab.textview.buffer[GtkTextBuffer]
    t = String[]
    for i = 1:length(cm)
        c = get_tab(cm, i)
        if typeof(c) <: Console
            push!(t," ID: $(c.worker_idx) \t Port: $(c.worker_port) \t Status: $(c.busy)\n")
        end
    end
    b.text[String] = string("Consoles:\n",t...)

end


text:
module TestUtils

using Base.Test
using RLEVectors

@testset begin

# rep
@test RLEVectors.rep([4,5,6], each=2) == [4,4,5,5,6,6]
@test RLEVectors.rep([4,5,6], times=2) == [4,5,6,4,5,6]
@test RLEVectors.rep([4,5,6], each=3, times=2) == [4,4,4,5,5,5,6,6,6,4,4,4,5,5,5,6,6,6]
@test RLEVectors.rep(["p","d","q"], each=3, times=2) == ["p","p","p","d","d","d","q","q","q","p","p","p","d","d","d","q","q","q"]
@test_throws ArgumentError RLEVectors.rep([1,2,3], each=[4,5])

# searchsortedfirst
@test searchsortedfirst( [3, 6, 9, 12, 15 ], 6, 1, 5 ) == 2
@test searchsortedfirst( [3, 6, 9, 12, 15 ], 7, 1, 5 ) == 3
@test searchsortedfirst( [3, 6, 9, 12, 15 ], 7, 0, 22 ) == 3
@test searchsortedfirst( [3, 6, 9, 12, 15 ], 40, 1, 5 ) == 6
@test searchsortedfirst( [3, 6, 9, 12, 15 ], -2, 1, 5 ) == 1
@test searchsortedfirst( [3, 6, 9, 12, 15 ], 14, 4, 5 ) == 5
@test searchsortedfirst( [3, 6, 9, 12, 15 ], 3, 4, 5 ) == 4

@test searchsortedfirst( [0,5,10,15], [-3,2,3,7,22] ) == [1, 2, 2, 3, 5]
@test searchsortedfirst( [0,5,10,15], [-3,2,-3,7,22] ) == [1, 2, 1, 3, 5]

end # testset

end # module


text:
function send_stream(rd::IO, sock::TCPSocket, id::Int)
    nb = nb_available(rd)
    if nb > 0
        d = read(rd, nb)
        s = String(copy(d))

        if !isempty(s)
            # remotecall_fetch(print_to_console_remote,sock,s,id)

            info("sending $s to $(string(sock)) with id: $id")

            remotecall_fetch(include_string, sock,"
                eval(GtkIDE,:(
                    print_to_console_remote(\"$(s)\", $(id))
                ))
            ")

        end
    end
end

function watch_stream(rd::IO, sock::TCPSocket, id::Int)
    while !eof(rd) # blocks until something is available
        send_stream(rd,sock,id)
        sleep(0.01) # a little delay to accumulate output
    end
end

function print_to_console_remote(s,idx::Integer)
    info("received print data with index $idx" )
    #print the output to the right console
    for i = 1:length(main_window.console_manager)
        c = get_tab(main_window.console_manager,i)
        if c.worker_idx == idx
            write(c.stdout_buffer,s)
        end
    end
end


text:
# ---
# title: 394. Decode String
# id: problem394
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Stack, Depth-first Search
# link: <https://leetcode.com/problems/decode-string/description/>
# hidden: true
# ---
#
# Given an encoded string, return its decoded string.
#
# The encoding rule is: `k[encoded_string]`, where the _encoded_string_ inside
# the square brackets is being repeated exactly _k_ times. Note that _k_ is
# guaranteed to be a positive integer.
#
# You may assume that the input string is always valid; No extra white spaces,
# square brackets are well-formed, etc.
#
# Furthermore, you may assume that the original data does not contain any digits
# and that digits are only for those repeat numbers, _k_. For example, there
# won't be input like `3a` or `2[4]`.
#
#
#
# **Example 1:**
#
#
#
#     Input: s = "3[a]2[bc]"
#     Output: "aaabcbc"
#
#
# **Example 2:**
#
#
#
#     Input: s = "3[a2[c]]"
#     Output: "accaccacc"
#
#
# **Example 3:**
#
#
#
#     Input: s = "2[abc]3[cd]ef"
#     Output: "abcabccdcdcdef"
#
#
# **Example 4:**
#
#
#
#     Input: s = "abc3[cd]xyz"
#     Output: "abccdcdcdxyz"
#
#
#
#
# **Constraints:**
#
#   * `1 <= s.length <= 30`
#   * `s` consists of lowercase English letters, digits, and square brackets `'[]'`.
#   * `s` is guaranteed to be **a valid** input.
#   * All the integers in `s` are in the range `[1, 300]`.
#
#
## @lc code=start
using LeetCode

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


text:
"""
Convert angle from cartesian coordinate system to meterological coordinate system
"""

function tip_cart2met(wind_direction_cart)
    # convert from cart. polar system (CCW, 0 rad.=E) to meteorological polar system (CW, 0 rad.=N)
    wind_direction_met = 3.0*pi/2.0 - wind_direction_cart

    # adjust to be between 0 and 2pi
    if wind_direction_met < 0.0
        wind_direction_met += 2.0*pi
    end

    return wind_direction_met
end

function tip_wake_combination_model(deficit, wind_speed, old_deficit_sum)
    # Katic et al. 1986
    new_deficit_sum = sqrt(old_deficit_sum^2 + (wind_speed*deficit)^2)
    if new_deficit_sum > wind_speed
        new_deficit_sum = wind_speed
    end

    return new_deficit_sum

end

# get angle for two given turbines
function tip_get_angle(dx, dy)
    angle = atan(dy,dx)
end

# get probability for a given angle
function tip_get_probability_spline(wind_rose, step)
    pre_spline = Akima(wind_rose[:,1], wind_rose[:,3])
    p_spline(x) = pre_spline(x)*length(wind_rose[:,1])/(2*pi/step)
    return p_spline
end

# get speed for a given angle
function tip_get_speed_spline(wind_rose)
    speed_spline = Akima(wind_rose[:,1], wind_rose[:,2])
    return speed_spline
end

# get effective velocity for a given turbine
function tip_get_effvelocity(turbid, x, y, uave, probability_spline, speed_spline, ct_spline, pstep; r0=0.5, alpha=0.1)
    nturbines = length(x)
    deficit_sum = 0.0
    for j = 1:nturbines
        dx = x[j] - x[turbid]
        dy = y[j] - y[turbid]
        dr = sqrt(dx^2+dy^2)
        angle = tip_cart2met(tip_get_angle(dx, dy)) # radians

        wake_range = angle-atan(alpha):pstep:angle+atan(alpha)
        p = sum(probability_spline(wake_range))
        wind_speed = sum(speed_spline(wake_range))/length(wake_range)
        # counter = 0
        # for ang in (angle-beta*0.0):1:(angle+beta*0.0)
        #     p += probability_spline(ang)
        #     wind_speed += speed_spline(ang)
        #     counter += 1
        # end
        # wind_speed /= counter
        a = tip_get_ai(wind_speed, ct_spline)
        deficit = p*2.0*a*(r0/(r0+alpha*dr))^2

        deficit_sum = tip_wake_combination_model(deficit, wind_speed, deficit_sum)
    end

    ueff = uave - deficit_sum
    return ueff
end

# get effective cp for a turbine in a given direction
function tip_get_cp_spline(cpdata)
    cp_spline = Akima(cpdata[:,1], cpdata[:,2])
    return cp_spline
end

# create spline for thrust coefficient
function tip_get_ct_spline(ctdata)
    ct_spline = Akima(ctdata[:,1], ctdata[:,2])
    return ct_spline
end

# calculate ai given wind speed and a ct curve
function tip_get_ai(speed, ct_spline)
    ct = ct_spline(speed)
    if (ct > 0.96) # Glauert condition
        ai = 0.143 + sqrt(0.0203-0.6427*(0.889 - ct))
    else
        ai = 0.5(1 - sqrt(1-ct))
    end
    return ai
end

# get average wind speed
function tip_get_ave_wind_speed(wind_rose)
    wind_speeds = wind_rose[:,2]
    probabilities = wind_rose[:,3]
    ave_speed = sum(wind_speeds.*probabilities)
    return ave_speed
end

# get average power for a given turbine
function tip_get_effpower(turbid, x, y, uave, probability_spline, speed_spline, cp_spline, ct_spline, pstep; prated=5E6, u_cut_in=3.0, u_cut_out=25.0, density=1.225, r0=0.5, alpha=0.1)
    area = pi*r0^2
    ueff = tip_get_effvelocity(turbid, x, y, uave, probability_spline, speed_spline, ct_spline, pstep, r0=r0, alpha=alpha)
    peff = 0.0
    if ueff < u_cut_in
        return peff
    elseif ueff > u_cut_out
        return peff
    else
        cp = cp_spline(ueff)
        peff = 0.5*density*area*cp*ueff^3
    end

    if peff > prated
        peff = prated
    end
    return peff
end

# get annual energy production
function tip_get_aep(x, y, wind_rose, cpdata, ctdata; pstep=0.01, prated=5E6, u_cut_in=3.0, u_cut_out=25.0, density=1.225, r0=0.5, alpha=0.1)
    nturbines = length(x)
    uave = tip_get_ave_wind_speed(wind_rose)
    power = zeros(nturbines)
    probability_spline = tip_get_probability_spline(wind_rose,pstep)
    speed_spline = tip_get_speed_spline(wind_rose)
    cp_spline = tip_get_cp_spline(cpdata)
    ct_spline = tip_get_ct_spline(ctdata)
    for i = 1:nturbines
        power[i] = tip_get_effpower(i, x, y, uave, probability_spline, speed_spline, cp_spline, ct_spline, pstep, prated=prated, u_cut_in=u_cut_in, u_cut_out=u_cut_out, density=density, r0=r0, alpha=alpha)
    end
    aep = (365.25)*(24.0)*sum(power)
    return aep
end


text:
#------------------------Grad CAM--------------------------

struct GradCAM
  model::Chain
end

function save_gradcam(gradient, original_image_path, grad_file_name, heatmap_file_name, combined_file_name)
  gradient = max.(gradient, zero(gradient))
  gradient = permutedims((gradient - minimum(gradient))/(maximum(gradient) - minimum(gradient)), (2, 1))
  img = Gray.(gradient)
  img = imresize(img, 224, 224)
  display(img)
  h1 = plot(heatmap(float.(augment(img, FlipY())), color = :rainbow))
  try
    display(h1)
  catch
    info("The heatmap could not be displayed. The file will be saved at $heatmap_file_name")
  end
  mapped = float.(channelview(imresize(load(original_image_path), (224, 224)))) * 0.9 .+ reshape(float.(img), 1, 224, 224) * 3
  mapped -= minimum(mapped)
  mapped /= maximum(mapped)
  h2 = plot(heatmap(float.(Gray.(augment(colorview(RGB{eltype(mapped)}, mapped), FlipY()))), color = :rainbow))
  try
    display(h2)
  catch
    info("The heatmap could not be displayed. The file will be saved at $combined_file_name")
  end
  save(grad_file_name, img)
  try
    savefig(h1, heatmap_file_name)
    savefig(h2, combined_file_name)
  catch
    info("Encountered Errors while trying to save heatmaps. So the heatmaps are being returned")
    info("They need to be saved manually")
    (h1, h2)
  end
end

function save_gradient(x)
  global target_output = x |> cpu
  x
end

save_gradient(x::TrackedArray) = Tracker.track(save_gradient, x)

normalize_grad(grad) = grad / (sqrt(mean(abs2.(grad))) + eps())

function Tracker.back(::typeof(save_gradient), Δ, x)
  global guided_gradients = Δ |> cpu
  Tracker.@back(x, Δ)
end

function GradCAM(model::Chain, target::Int)
  length(model.layers) - 1 < target && error("target should be less than the model length")
  updated_model = []
  for (i, l) in enumerate(model.layers[1:end-1])
    push!(updated_model, l)
    i == target && push!(updated_model, save_gradient)
  end
  GradCAM(Chain(updated_model...) |> gpu)
end

function (m::GradCAM)(img, top = 1)
  if !Tracker.istracked(img)
    img = param(img) |> gpu
  end
  preds = m.model(img)
  probs = softmax(preds)
  prob, inds = get_topk(probs.data, top)
  grads = []
  for (i, idx) in enumerate(inds)
    Flux.back!(preds, one_hot_encode(preds, idx))
    weights = reshape(maximum(normalize_grad(guided_gradients), [1, 2, 4]), 1, 1, size(target_output)[3], 1)
    cam = ones(size(target_output)[1:2])
    cam += squeeze(sum(weights .* target_output, [3, 4]), dims = (3, 4))
    push!(grads, (cam, prob[i], idx))
    img.grad .= 0.0
  end
  grads
end

#--------------------Guided Grad CAM-----------------------

struct GuidedGradCAM
  gcam::GradCAM
  backprop::Backprop
end

GuidedGradCAM(model::Chain, target::Int) = GuidedGradCAM(GradCAM(model, target), GuidedBackprop(model))

function (m::GuidedGradCAM)(img, top = 1)
  # Ideally we should not be doing the forward pass twice
  gcam_grads = m.gcam(img, top)
  backprop_grads = m.backprop(img, top)
  grads = []
  for i in 1:top
    normalized_camgrads = gcam_grads[i][1] - minimum(gcam_grads[i][1])
    normalized_camgrads /= maximum(normalized_camgrads)
    normalized_camgrads = reshape(float.(imresize(Gray.(normalized_camgrads), 224, 224)), 224, 224, 1, 1)
    normalized_bpropgrads = backprop_grads[i][1] - minimum(backprop_grads[i][1])
    normalized_bpropgrads /= maximum(normalized_bpropgrads)
    push!(grads, (normalized_camgrads .* normalized_bpropgrads, gcam_grads[i][2], gcam_grads[i][3]))
  end
  grads
end


text:
"""
# A flux of a scalar quantity is itself a quantity
# """
type Flux{S<:Real,T<:Quantity{S}} <: Quantity{S}
    ts::Vector{DateTime}
    q::Vector{S}
end

function Base.:*{S<:Real,T<:Quantity{S}}(a::T,b::Discharge)
    times(a) == times(b) || error("Time series are not sampled at the same times. Consider interpolating")
    Flux{S,T}(times(a),quantity(a).*quantity(b))
end

function Base.:*{S<:Real,T<:Quantity{S}}(a::Discharge,b::T)
    times(a) == times(b) || error("Time series are not sampled at the same times. Consider interpolating")
    Flux{S,T}(times(a),quantity(a).*quantity(b))
end


text:
# =================================================
# Platooning model
# See https://easychair.org/publications/paper/3QLs
# =================================================
using Reachability, HybridSystems, MathematicalSystems, MathematicalPredicates,
      SparseArrays
using LazySets: HalfSpace  # resolve name-space conflicts with Polyhedra

"""
    platooning(;
               [deterministic_switching]::Bool=true,
               [time_horizon]::Float64=20.,
               [allowed_distance]::Float64=50.)

Construct the platooning model consisting of three vehicles.

### Input

- `deterministic_switching` -- (optional, default: `true`) flag to use
                               deterministic transitions
- `time_horizon`            -- (optional, default: `20.`) time horizon
- `allowed_distance`        -- (optional, default: `50.`) specification of
                               allowed distance between the vehicles
"""
function platooning(;
                    deterministic_switching::Bool=true,
                    time_horizon::Float64=20.,
                    allowed_distance::Float64=50.)
    # three variables for each vehicle, (ei, d(et)/dt, ai) for
    # (spacing error, relative velocity, speed), and the last dimension is time
    n = 9 + 1

    # constants
    c1 = c2 = 5.   # clock constraints
    tb = 10.       # lower bound for loss of communication
    tc = tr = 20.  # upper bound for loss of communication (tc) and reset time (tr)

    # additional bloating to get a non-flat guard intersection
    guard_bloating = sqrt(eps(Float64))

    # transition graph
    automaton = LightAutomaton(2)
    add_transition!(automaton, 1, 2, 1)
    add_transition!(automaton, 2, 1, 2)

    # common inputs
    B = sparse([2, n], [1, 2], [1., 1.], n, 2)
    U = Hyperrectangle(low=[-9., 1.], high=[1., 1.])  # acceleration of the lead vehicle + time

    # mode 1 ("connected")
    A = zeros(n, n)
    for (i, j) in [(1, 2), (4, 5), (5, 3), (7, 8), (8, 6)]
        A[i, j] = 1.
    end
    for (i, j) in [(2, 3), (5, 6), (8, 9)]
        A[i, j] = -1.
    end
    A[3, :] = [1.6050, 4.8680, -3.5754, -0.8198, 0.4270, -0.0450, -0.1942,  0.3626, -0.0946, 0.]
    A[6, :] = [0.8718, 3.8140, -0.0754,  1.1936, 3.6258, -3.2396, -0.5950,  0.1294, -0.0796, 0.]
    A[9, :] = [0.7132, 3.5730, -0.0964,  0.8472, 3.2568, -0.0876,  1.2726,  3.0720, -3.1356, 0.]
    invariant = deterministic_switching ?
        HalfSpace(sparsevec([n], [1.], n), c1) :
        Universe(n)
    m_1 = ConstrainedLinearControlContinuousSystem(A, B, invariant, U)

    # mode 2 ("not connected/connection broken")
    A = copy(A)
    A[3, 4:9] = zeros(6)
    A[6, 1:3] = A[6, 7:9] = zeros(3)
    invariant = deterministic_switching ?
        HalfSpace(sparsevec([n], [1.], n), c2) :
        Universe(n)
    m_2 = ConstrainedLinearControlContinuousSystem(A, B, invariant, U)

    # modes
    modes = [m_1, m_2]

    # common reset
    reset = Dict(n => 0.)

    # transition l1 -> l2
    # (using a hyperplane in the deterministic case causes floating-point issues)
    guard = deterministic_switching ?
        HalfSpace(sparsevec([n], [-1.], n), -c1 + guard_bloating) :
        HPolyhedron([HalfSpace(sparsevec([n], [-1.], n), -tb),
                     HalfSpace(sparsevec([n], [1.], n), tc)])
    t1 = ConstrainedResetMap(n, guard, reset)

    # transition l2 -> l1
    guard = deterministic_switching ?
        HalfSpace(sparsevec([n], [-1.], n), -c2 + guard_bloating) :
        HalfSpace(sparsevec([n], [1.], n), tr)
    t2 = ConstrainedResetMap(n, guard, reset)

    # transition annotations
    resetmaps = [t1, t2]

    # switching
    switchings = [AutonomousSwitching()]

    ℋ = HybridSystem(automaton, modes, resetmaps, switchings)

    # initial condition in mode 1
    X0 = Singleton(zeros(n))
    initial_condition = [(1, X0)]

    system = InitialValueProblem(ℋ, initial_condition)

    # safety property
    d1 = zeros(n); d1[1] = -1.  # x1 >= -dmin
    d4 = zeros(n); d4[4] = -1.  # x4 >= -dmin
    d7 = zeros(n); d7[7] = -1.  # x7 >= -dmin
    property = Conjunction(
        [is_contained_in(HalfSpace(d, allowed_distance)) for d in [d1, d4, d7]])

    # default options
    options = Options(:T=>time_horizon, :property=>property)

    return (system, options)
end

function run_platooning(system, options)
    opC = BFFPSV18(:δ => 0.01, :assume_sparse => true)
    opD = LazyDiscretePost(:lazy_R⋂I => true, :lazy_R⋂G => false)
    options[:mode] = "check"
    options[:plot_vars] = [0, 1]

    solve(system, options, opC, opD)
end

max_stars_repo_name: JuliaReach/ReachabilityBenchmarks

max_issues_repo_name: JuliaReach/ReachabilityBenchmarks

max_forks_repo_name: JuliaReach/ReachabilityBenchmarks


text:
println("starting!!!")
using DelimitedFiles
using TicketToRide
using Flux
using Flux.Optimise
using TicketToRide: train!

Ls = [4, 6, 8]
Nlayers = collect(1:50)
repeats = [i for i in 1:5]

arg = 1
 #parse(Int, ARGS[1])

list = Iterators.product(Ls, Nlayers, repeats) |> collect
println("length of array jobs is:\t", length(list))

n, nlayers, r = list[arg]

for L in Ls
    for nlayers in Nlayers
        for repeats in 1:10
            @info "\nL is $n"
            @info "nlayers is $nlayers"
            @info "r is $r"
            # make directory if none exists
            dir = "/scratch/mbeach/new_new_tickettoride/L-$n/"
            mkpath(dir)
            file = dir * "layers-$nlayers-r-$r.txt"
            @info file

            # pre-compile once
            opt = Optimise.ADAMW()
            circuit = variational_circuit(n, nlayers);
            ham = heisenberg1D(n)
            dispatch!(circuit, :random);

            ind = findall(x -> x == n, Ls)
            println(ind)

            opt = Optimise.ADAM()
            # @time history = train!(opt, 1, ham, circuit; verbose=true)
            history = train!(opt, 2000, ham, circuit; verbose=false)
            writedlm(file, history)

            opt = Optimise.Descent(1e-3)
            history = train!(opt, 2000, ham, circuit; verbose=false)
            writedlm(file, history)

            opt = Optimise.Descent(1e-4)
            history = train!(opt, 20000, ham, circuit; verbose=false)
            writedlm(file, history)

        end
    end
end