text:
"""
    AnisotropicMinimumDissipation{FT} <: AbstractTurbulenceClosure

Parameters for the "anisotropic minimum dissipation" turbulence closure for large eddy simulation
proposed originally by [Rozema15](@cite) and [Abkar16](@cite), and then modified
by [Verstappen18](@cite), and finally described and validated for by [Vreugdenhil18](@cite).
"""
struct AnisotropicMinimumDissipation{FT, PK, PN, K} <: AbstractEddyViscosityClosure
    Cν :: PN
    Cκ :: PK
    Cb :: FT
     ν :: FT
     κ :: K

    function AnisotropicMinimumDissipation{FT}(Cν, Cκ, Cb, ν, κ) where FT
        return new{FT, typeof(Cκ), typeof(Cν), typeof(κ)}(Cν, Cκ, Cb, ν, convert_diffusivity(FT, κ))
    end
end

const AMD = AnisotropicMinimumDissipation

Base.show(io::IO, closure::AMD{FT}) where FT =
    print(io, "AnisotropicMinimumDissipation{$FT} turbulence closure with:\n",
              "           Poincaré constant for momentum eddy viscosity Cν: ", closure.Cν, '\n',
              "    Poincaré constant for tracer(s) eddy diffusivit(ies) Cκ: ", closure.Cκ, '\n',
              "                        Buoyancy modification multiplier Cb: ", closure.Cb, '\n',
              "                Background diffusivit(ies) for tracer(s), κ: ", closure.κ, '\n',
              "             Background kinematic viscosity for momentum, ν: ", closure.ν)

"""
    AnisotropicMinimumDissipation(FT=Float64; C=1/12, Cν=nothing, Cκ=nothing,
                                            Cb=0.0, ν=ν₀, κ=κ₀)

Returns parameters of type `FT` for the `AnisotropicMinimumDissipation`
turbulence closure.

Keyword arguments
=================
    - `C`  : Poincaré constant for both eddy viscosity and eddy diffusivities. `C` is overridden
             for eddy viscosity or eddy diffusivity if `Cν` or `Cκ` are set, respecitvely.
    - `Cν` : Poincaré constant for momentum eddy viscosity.
    - `Cκ` : Poincaré constant for tracer eddy diffusivities. If one number or function, the same
             number or function is applied to all tracers. If a `NamedTuple`, it must possess
             a field specifying the Poncaré constant for every tracer.
    - `Cb` : Buoyancy modification multiplier (`Cb = 0` turns it off, `Cb = 1` turns it on)
    - `ν`  : Constant background viscosity for momentum.
    - `κ`  : Constant background diffusivity for tracer. If a single number, the same background
             diffusivity is applied to all tracers. If a `NamedTuple`, it must possess a field
             specifying a background diffusivity for every tracer.

By default: `C = Cν = Cκ` = 1/12, which is appropriate for a finite-volume method employing a
second-order advection scheme, `Cb` = 0, which terms off the buoyancy modification term,
the molecular viscosity of seawater at 20 deg C and 35 psu is used for `ν`, and
the molecular diffusivity of heat in seawater at 20 deg C and 35 psu is used for `κ`.

`Cν` or `Cκ` may be constant numbers, or functions of `x, y, z`.

Example
=======

julia> pretty_diffusive_closure = AnisotropicMinimumDissipation(C=1/2)
AnisotropicMinimumDissipation{Float64} turbulence closure with:
           Poincaré constant for momentum eddy viscosity Cν: 0.5
    Poincaré constant for tracer(s) eddy diffusivit(ies) Cκ: 0.5
                        Buoyancy modification multiplier Cb: 0.0
                Background diffusivit(ies) for tracer(s), κ: 1.46e-7
             Background kinematic viscosity for momentum, ν: 1.05e-6

julia> const Δz = 0.5; # grid resolution at surface

julia> surface_enhanced_tracer_C(x, y, z) = 1/12 * (1 + exp((z + Δz/2) / 8Δz))
surface_enhanced_tracer_C (generic function with 1 method)

julia> fancy_closure = AnisotropicMinimumDissipation(Cκ=surface_enhanced_tracer_C)
AnisotropicMinimumDissipation{Float64} turbulence closure with:
           Poincaré constant for momentum eddy viscosity Cν: 0.08333333333333333
    Poincaré constant for tracer(s) eddy diffusivit(ies) Cκ: surface_enhanced_tracer_C
                        Buoyancy modification multiplier Cb: 0.0
                Background diffusivit(ies) for tracer(s), κ: 1.46e-7
             Background kinematic viscosity for momentum, ν: 1.05e-6

julia> tracer_specific_closure = AnisotropicMinimumDissipation(Cκ=(c₁=1/12, c₂=1/6))
AnisotropicMinimumDissipation{Float64} turbulence closure with:
           Poincaré constant for momentum eddy viscosity Cν: 0.08333333333333333
    Poincaré constant for tracer(s) eddy diffusivit(ies) Cκ: (c₁ = 0.08333333333333333, c₂ = 0.16666666666666666)
                        Buoyancy modification multiplier Cb: 0.0
                Background diffusivit(ies) for tracer(s), κ: 1.46e-7
             Background kinematic viscosity for momentum, ν: 1.05e-6

References
==========
Vreugdenhil C., and Taylor J. (2018), "Large-eddy simulations of stratified plane Couette
    flow using the anisotropic minimum-dissipation model", Physics of Fluids 30, 085104.

Verstappen, R. (2018), "How much eddy dissipation is needed to counterbalance the nonlinear
    production of small, unresolved scales in a large-eddy simulation of turbulence?",
    Computers & Fluids 176, pp. 276-284.
"""
function AnisotropicMinimumDissipation(FT=Float64; C=1/12, Cν=nothing, Cκ=nothing,
                                                 Cb=0.0, ν=ν₀, κ=κ₀)
    Cν = Cν === nothing ? C : Cν
    Cκ = Cκ === nothing ? C : Cκ

    return AnisotropicMinimumDissipation{FT}(Cν, Cκ, Cb, ν, κ)
end

function with_tracers(tracers, closure::AnisotropicMinimumDissipation{FT}) where FT
    κ = tracer_diffusivities(tracers, closure.κ)
    Cκ = tracer_diffusivities(tracers, closure.Cκ)
    return AnisotropicMinimumDissipation{FT}(closure.Cν, Cκ, closure.Cb, closure.ν, κ)
end

#####
##### Diffusive fluxes
#####

function diffusive_flux_x(i, j, k, grid, clock, closure::AnisotropicMinimumDissipation, c,
                          ::Val{tracer_index}, diffusivities, args...) where tracer_index

    κₑ = diffusivities.κₑ[tracer_index]
    return diffusive_flux_x(i, j, k, grid, clock, κₑ, c)
end

function diffusive_flux_y(i, j, k, grid, clock, closure::AnisotropicMinimumDissipation, c,
                          ::Val{tracer_index}, diffusivities, args...) where tracer_index

    κₑ = diffusivities.κₑ[tracer_index]
    return diffusive_flux_y(i, j, k, grid, clock, κₑ, c)
end

function diffusive_flux_z(i, j, k, grid, clock, closure::AnisotropicMinimumDissipation, c,
                          ::Val{tracer_index}, diffusivities, args...) where tracer_index

    κₑ = diffusivities.κₑ[tracer_index]
    return diffusive_flux_z(i, j, k, grid, clock, κₑ, c)
end

#####
##### Kernel functions
#####

# Dispatch on the type of the user-provided AMD model constant.
# Only numbers, arrays, and functions supported now.
@inline Cᴾᵒⁱⁿ(i, j, k, grid, C::Number) = C
@inline Cᴾᵒⁱⁿ(i, j, k, grid, C::AbstractArray) = @inbounds C[i, j, k]
@inline Cᴾᵒⁱⁿ(i, j, k, grid, C::Function) = C(xnode(Center, i, grid), ynode(Center, j, grid), znode(Center, k, grid))

@inline function νᶜᶜᶜ(i, j, k, grid::AbstractGrid{FT}, closure::AMD, buoyancy, U, C) where FT
    ijk = (i, j, k, grid)
    q = norm_tr_∇uᶜᶜᶜ(ijk..., U.u, U.v, U.w)

    if q == 0 # SGS viscosity is zero when strain is 0
        νˢᵍˢ = zero(FT)
    else
        r = norm_uᵢₐ_uⱼₐ_Σᵢⱼᶜᶜᶜ(ijk..., closure, U.u, U.v, U.w)
        ζ = norm_wᵢ_bᵢᶜᶜᶜ(ijk..., closure, buoyancy, U.w, C) / Δᶠzᶜᶜᶜ(ijk...)
        δ² = 3 / (1 / Δᶠxᶜᶜᶜ(ijk...)^2 + 1 / Δᶠyᶜᶜᶜ(ijk...)^2 + 1 / Δᶠzᶜᶜᶜ(ijk...)^2)
        νˢᵍˢ = - Cᴾᵒⁱⁿ(i, j, k, grid, closure.Cν) * δ² * (r - closure.Cb * ζ) / q
    end

    return max(zero(FT), νˢᵍˢ) + closure.ν
end

@inline function κᶜᶜᶜ(i, j, k, grid::AbstractGrid{FT}, closure::AMD, c, ::Val{tracer_index},
                       U) where {FT, tracer_index}

    ijk = (i, j, k, grid)

    @inbounds κ = closure.κ[tracer_index]
    @inbounds Cκ = closure.Cκ[tracer_index]

    σ =  norm_θᵢ²ᶜᶜᶜ(i, j, k, grid, c)

    if σ == 0 # denominator is zero: short-circuit computations and set subfilter diffusivity to zero.
        κˢᵍˢ = zero(FT)
    else
        ϑ =  norm_uᵢⱼ_cⱼ_cᵢᶜᶜᶜ(ijk..., closure, U.u, U.v, U.w, c)
        δ² = 3 / (1 / Δᶠxᶜᶜᶜ(ijk...)^2 + 1 / Δᶠyᶜᶜᶜ(ijk...)^2 + 1 / Δᶠzᶜᶜᶜ(ijk...)^2)
        κˢᵍˢ = - Cᴾᵒⁱⁿ(i, j, k, grid, Cκ) * δ² * ϑ / σ
    end

    return max(zero(FT), κˢᵍˢ) + κ
end

function calculate_diffusivities!(K, arch, grid, closure::AnisotropicMinimumDissipation, buoyancy, U, C)
    workgroup, worksize = work_layout(grid, :xyz)

    barrier = Event(device(arch))

    viscosity_kernel! = calculate_viscosity!(device(arch), workgroup, worksize)
    diffusivity_kernel! = calculate_tracer_diffusivity!(device(arch), workgroup, worksize)

    viscosity_event = viscosity_kernel!(K.νₑ, grid, closure, buoyancy, U, C, dependencies=barrier)

    events = [viscosity_event]

    for (tracer_index, κₑ) in enumerate(K.κₑ)
        @inbounds c = C[tracer_index]
        event = diffusivity_kernel!(κₑ, grid, closure, c, Val(tracer_index), U, dependencies=barrier)
        push!(events, event)
    end

    wait(device(arch), MultiEvent(Tuple(events)))

    return nothing
end

@kernel function calculate_viscosity!(νₑ, grid, closure::AnisotropicMinimumDissipation, buoyancy, U, C)
    i, j, k = @index(Global, NTuple)
    @inbounds νₑ[i, j, k] = νᶜᶜᶜ(i, j, k, grid, closure, buoyancy, U, C)
end

@kernel function calculate_tracer_diffusivity!(κₑ, grid, closure, c, tracer_index, U)
    i, j, k = @index(Global, NTuple)
    @inbounds κₑ[i, j, k] = κᶜᶜᶜ(i, j, k, grid, closure, c, tracer_index, U)
end

#####
##### Filter width at various locations
#####

# Recall that filter widths are 2x the grid spacing in AMD
@inline Δᶠxᶜᶜᶜ(i, j, k, grid) = 2 * Oceananigans.Operators.Δx(i, j, k, grid)
@inline Δᶠyᶜᶜᶜ(i, j, k, grid) = 2 * Oceananigans.Operators.Δy(i, j, k, grid)
@inline Δᶠzᶜᶜᶜ(i, j, k, grid) = 2 * Oceananigans.Operators.ΔzC(i, j, k, grid)

for loc in (:ccf, :fcc, :cfc, :ffc, :cff, :fcf), ξ in (:x, :y, :z)
    Δ_loc = Symbol(:Δᶠ, ξ, :_, loc)
    Δᶜᶜᶜ = Symbol(:Δᶠ, ξ, :ᶜᶜᶜ)
    @eval begin
        const $Δ_loc = $Δᶜᶜᶜ
    end
end

#####
##### The *** 30 terms *** of AMD
#####

@inline function norm_uᵢₐ_uⱼₐ_Σᵢⱼᶜᶜᶜ(i, j, k, grid, closure, u, v, w)
    ijk = (i, j, k, grid)
    uvw = (u, v, w)
    ijkuvw = (i, j, k, grid, u, v, w)

    uᵢ₁_uⱼ₁_Σ₁ⱼ = (
         norm_Σ₁₁(ijkuvw...) * norm_∂x_u(ijk..., u)^2
      +  norm_Σ₂₂(ijkuvw...) * ℑxyᶜᶜᵃ(ijk..., norm_∂x_v², uvw...)
      +  norm_Σ₃₃(ijkuvw...) * ℑxzᶜᵃᶜ(ijk..., norm_∂x_w², uvw...)

      +  2 * norm_∂x_u(ijkuvw...) * ℑxyᶜᶜᵃ(ijk..., norm_∂x_v_Σ₁₂, uvw...)
      +  2 * norm_∂x_u(ijkuvw...) * ℑxzᶜᵃᶜ(ijk..., norm_∂x_w_Σ₁₃, uvw...)
      +  2 * ℑxyᶜᶜᵃ(ijk..., norm_∂x_v, uvw...) * ℑxzᶜᵃᶜ(ijk..., norm_∂x_w, uvw...)
           * ℑyzᵃᶜᶜ(ijk..., norm_Σ₂₃, uvw...)
    )

    uᵢ₂_uⱼ₂_Σ₂ⱼ = (
      + norm_Σ₁₁(ijkuvw...) * ℑxyᶜᶜᵃ(ijk..., norm_∂y_u², uvw...)
      + norm_Σ₂₂(ijkuvw...) * norm_∂y_v(ijk..., v)^2
      + norm_Σ₃₃(ijkuvw...) * ℑyzᵃᶜᶜ(ijk..., norm_∂y_w², uvw...)

      +  2 * norm_∂y_v(ijkuvw...) * ℑxyᶜᶜᵃ(ijk..., norm_∂y_u_Σ₁₂, uvw...)
      +  2 * ℑxyᶜᶜᵃ(ijk..., norm_∂y_u, uvw...) * ℑyzᵃᶜᶜ(ijk..., norm_∂y_w, uvw...)
           * ℑxzᶜᵃᶜ(ijk..., norm_Σ₁₃, uvw...)
      +  2 * norm_∂y_v(ijkuvw...) * ℑyzᵃᶜᶜ(ijk..., norm_∂y_w_Σ₂₃, uvw...)
    )

    uᵢ₃_uⱼ₃_Σ₃ⱼ = (
      + norm_Σ₁₁(ijkuvw...) * ℑxzᶜᵃᶜ(ijk..., norm_∂z_u², uvw...)
      + norm_Σ₂₂(ijkuvw...) * ℑyzᵃᶜᶜ(ijk..., norm_∂z_v², uvw...)
      + norm_Σ₃₃(ijkuvw...) * norm_∂z_w(ijk..., w)^2

      +  2 * ℑxzᶜᵃᶜ(ijk..., norm_∂z_u, uvw...) * ℑyzᵃᶜᶜ(ijk..., norm_∂z_v, uvw...)
           * ℑxyᶜᶜᵃ(ijk..., norm_Σ₁₂, uvw...)
      +  2 * norm_∂z_w(ijkuvw...) * ℑxzᶜᵃᶜ(ijk..., norm_∂z_u_Σ₁₃, uvw...)
      +  2 * norm_∂z_w(ijkuvw...) * ℑyzᵃᶜᶜ(ijk..., norm_∂z_v_Σ₂₃, uvw...)
    )

    return uᵢ₁_uⱼ₁_Σ₁ⱼ + uᵢ₂_uⱼ₂_Σ₂ⱼ + uᵢ₃_uⱼ₃_Σ₃ⱼ
end

#####
##### trace(∇u) = uᵢⱼ uᵢⱼ
#####

@inline function norm_tr_∇uᶜᶜᶜ(i, j, k, grid, uvw...)
    ijk = (i, j, k, grid)

    return (
        # ccc
        norm_∂x_u²(ijk..., uvw...)
      + norm_∂y_v²(ijk..., uvw...)
      + norm_∂z_w²(ijk..., uvw...)

        # ffc
      + ℑxyᶜᶜᵃ(ijk..., norm_∂x_v², uvw...)
      + ℑxyᶜᶜᵃ(ijk..., norm_∂y_u², uvw...)

        # fcf
      + ℑxzᶜᵃᶜ(ijk..., norm_∂x_w², uvw...)
      + ℑxzᶜᵃᶜ(ijk..., norm_∂z_u², uvw...)

        # cff
      + ℑyzᵃᶜᶜ(ijk..., norm_∂y_w², uvw...)
      + ℑyzᵃᶜᶜ(ijk..., norm_∂z_v², uvw...)
    )
end

@inline function norm_wᵢ_bᵢᶜᶜᶜ(i, j, k, grid, closure, buoyancy, w, C)
    ijk = (i, j, k, grid)

    wx_bx = (ℑxzᶜᵃᶜ(ijk..., norm_∂x_w, w)
             * Δᶠxᶜᶜᶜ(ijk...) * ℑxᶜᵃᵃ(ijk..., ∂xᶠᵃᵃ, buoyancy_perturbation, buoyancy, C))

    wy_by = (ℑyzᵃᶜᶜ(ijk..., norm_∂y_w, w)
             * Δᶠyᶜᶜᶜ(ijk...) * ℑyᵃᶜᵃ(ijk..., ∂yᵃᶠᵃ, buoyancy_perturbation, buoyancy, C))

    wz_bz = (norm_∂z_w(ijk..., w)
             * Δᶠzᶜᶜᶜ(ijk...) * ℑzᵃᵃᶜ(ijk..., ∂zᵃᵃᶠ, buoyancy_perturbation, buoyancy, C))

    return wx_bx + wy_by + wz_bz
end

@inline function norm_uᵢⱼ_cⱼ_cᵢᶜᶜᶜ(i, j, k, grid, closure, u, v, w, c)
    ijk = (i, j, k, grid)

    cx_ux = (
                  norm_∂x_u(ijk..., u) * ℑxᶜᵃᵃ(ijk..., norm_∂x_c², c)
        + ℑxyᶜᶜᵃ(ijk..., norm_∂x_v, v) * ℑxᶜᵃᵃ(ijk..., norm_∂x_c, c) * ℑyᵃᶜᵃ(ijk..., norm_∂y_c, c)
        + ℑxzᶜᵃᶜ(ijk..., norm_∂x_w, w) * ℑxᶜᵃᵃ(ijk..., norm_∂x_c, c) * ℑzᵃᵃᶜ(ijk..., norm_∂z_c, c)
    )

    cy_uy = (
          ℑxyᶜᶜᵃ(ijk..., norm_∂y_u, u) * ℑyᵃᶜᵃ(ijk..., norm_∂y_c, c) * ℑxᶜᵃᵃ(ijk..., norm_∂x_c, c)
        +         norm_∂y_v(ijk..., v) * ℑyᵃᶜᵃ(ijk..., norm_∂y_c², c)
        + ℑxzᶜᵃᶜ(ijk..., norm_∂y_w, w) * ℑyᵃᶜᵃ(ijk..., norm_∂y_c, c) * ℑzᵃᵃᶜ(ijk..., norm_∂z_c, c)
    )

    cz_uz = (
          ℑxzᶜᵃᶜ(ijk..., norm_∂z_u, u) * ℑzᵃᵃᶜ(ijk..., norm_∂z_c, c) * ℑxᶜᵃᵃ(ijk..., norm_∂x_c, c)
        + ℑyzᵃᶜᶜ(ijk..., norm_∂z_v, v) * ℑzᵃᵃᶜ(ijk..., norm_∂z_c, c) * ℑyᵃᶜᵃ(ijk..., norm_∂y_c, c)
        +         norm_∂z_w(ijk..., w) * ℑzᵃᵃᶜ(ijk..., norm_∂z_c², c)
    )

    return cx_ux + cy_uy + cz_uz
end

@inline norm_θᵢ²ᶜᶜᶜ(i, j, k, grid, c) = (
      ℑxᶜᵃᵃ(i, j, k, grid, norm_∂x_c², c)
    + ℑyᵃᶜᵃ(i, j, k, grid, norm_∂y_c², c)
    + ℑzᵃᵃᶜ(i, j, k, grid, norm_∂z_c², c)
)


text:
using Revise
using AdFem
using PyPlot
using LinearAlgebra

mmesh = Mesh(10, 10, 0.1)
u = rand(2mmesh.ndof)
ε = eval_strain_on_gauss_pts(u, mmesh)
K = [1.0 0.0 0.0;0.0 1.0 0.0; 0.0 0.0 0.5]
σ = zeros(size(ε)...)
for i = 1:get_ngauss(mmesh)
    σ[i,:] = K * ε[i,:]
end
e1 = compute_strain_energy_term(σ, mmesh)

Q = compute_fem_stiffness_matrix(K, mmesh)
e2 = Q * u

@info norm(e1 - e2 )

mmesh = Mesh(10, 10, 0.1, degree=2)
u = rand(2mmesh.ndof)
ε = eval_strain_on_gauss_pts(u, mmesh)
K = [1.0 0.0 0.0;0.0 1.0 0.0; 0.0 0.0 0.5]
σ = zeros(size(ε)...)
for i = 1:get_ngauss(mmesh)
    σ[i,:] = K * ε[i,:]
end
e1 = compute_strain_energy_term(σ, mmesh)

Q = compute_fem_stiffness_matrix(K, mmesh)
e2 = Q * u

@info norm(e1 - e2 )

ext: jl


text:
using Documenter, DiffEqBase, SciMLBase, DiffEqProblemLibrary

ODEProblemLibrary = DiffEqProblemLibrary.ODEProblemLibrary
ODEProblemLibrary.importodeproblems()

SDEProblemLibrary = DiffEqProblemLibrary.SDEProblemLibrary
SDEProblemLibrary.importsdeproblems()

DDEProblemLibrary = DiffEqProblemLibrary.DDEProblemLibrary
DDEProblemLibrary.importddeproblems()

DAEProblemLibrary = DiffEqProblemLibrary.DAEProblemLibrary
DAEProblemLibrary.importdaeproblems()

makedocs(modules=[DiffEqBase,SciMLBase,DiffEqProblemLibrary,ODEProblemLibrary,SDEProblemLibrary,DDEProblemLibrary,DAEProblemLibrary],
         doctest=false, clean=true,
         format = Documenter.HTML(analytics = "UA-90474609-3",
                                  assets = ["assets/favicon.ico"],
                                  canonical="https://diffeq.sciml.ai/stable/"),
         sitename="DifferentialEquations.jl",
         authors="Chris Rackauckas",
         pages = Any[
         "DifferentialEquations.jl: Scientific Machine Learning (SciML) Enabled Simulation and Estimation" => "index.md",
         "Tutorials" => Any[
           "tutorials/ode_example.md",
           "tutorials/advanced_ode_example.md",
           "tutorials/sde_example.md",
           "tutorials/rode_example.md",
           "tutorials/dde_example.md",
           "tutorials/dae_example.md",
           "tutorials/discrete_stochastic_example.md",
           "tutorials/jump_diffusion.md",
           "tutorials/bvp_example.md",
           "tutorials/additional.md"
         ],
         "Basics" => Any[
           "basics/overview.md",
           "basics/common_solver_opts.md",
           "basics/solution.md",
           "basics/plot.md",
           "basics/integrator.md",
           "basics/problem.md",
           "basics/faq.md",
           "basics/compatibility_chart.md"
         ],
         "Problem Types" => Any[
           "types/discrete_types.md",
           "types/ode_types.md",
           "types/nonautonomous_linear_ode.md",
           "types/dynamical_types.md",
           "types/split_ode_types.md",
           "types/steady_state_types.md",
           "types/bvp_types.md",
           "types/sde_types.md",
           "types/sdae_types.md",
           "types/rode_types.md",
           "types/dde_types.md",
           "types/sdde_types.md",
           "types/dae_types.md",
           "types/jump_types.md",
         ],
         "Solver Algorithms" => Any[
           "solvers/discrete_solve.md",
           "solvers/ode_solve.md",
           "solvers/nonautonomous_linear_ode.md",
           "solvers/dynamical_solve.md",
           "solvers/split_ode_solve.md",
           "solvers/steady_state_solve.md",
           "solvers/bvp_solve.md",
           "solvers/jump_solve.md",
           "solvers/sde_solve.md",
           "solvers/sdae_solve.md",
           "solvers/rode_solve.md",
           "solvers/dde_solve.md",
           "solvers/sdde_solve.md",
           "solvers/dae_solve.md",
           "solvers/benchmarks.md"
         ],
         "Additional Features" => Any[
           "features/performance_overloads.md",
           "features/diffeq_arrays.md",
           "features/diffeq_operator.md",
           "features/noise_process.md",
           "features/linear_nonlinear.md",
           "features/callback_functions.md",
           "features/callback_library.md",
           "features/ensemble.md",
           "features/io.md",
           "features/low_dep.md",
           "features/progress_bar.md"
         ],
         "Analysis Tools" => Any[
           "analysis/parameterized_functions.md",
           "analysis/parameter_estimation.md",
           "analysis/bifurcation.md",
           "analysis/sensitivity.md",
           "analysis/global_sensitivity.md",
           "analysis/uncertainty_quantification.md",
           "analysis/neural_networks.md",
           "analysis/dev_and_test.md"
         ],
         "Domain Modeling Tools" => Any[
             "models/multiscale.md",
             "models/physical.md",
             "models/financial.md",
             "models/chemical_reactions.md",
             "models/external_modeling.md"
         ],
         "Extra Details" => Any[
             "extras/timestepping.md",
             "extras/sensitivity_math.md",
         ]
         ])

#Redirect old links
# cd(joinpath(@__DIR__, "build")) do
#     for (root, dirs, files) in walkdir(".")
#         for file in files
#             path = relpath(joinpath(root, file), ".")
#             m = match(r"(.+)/index\.html$", path)
#             m === nothing && continue
#             redirect = "$(m[1]).html"
#             @info "Adding redirect for $(m[1]) from $(redirect)"
#             isfile(redirect) && (@warn "$redirect exists, skip"; continue)
#             open(redirect, "w") do io
#                 write(io, """
#                 <meta http-equiv="refresh" content="0; url=$(basename(m[1]))/"/>
#                 """)
#             end
#         end
#     end
# end

deploydocs(
   repo = "github.com/SciML/DiffEqDocs.jl.git"
)


text:
module SampledVectors

export SampledVector, sampled, sampledindexes, capacity, setcapacity!, factor, setfactor!, downsample!

mutable struct SampledVector{T, ST<:Integer, LT<:Integer, CT<:Integer, FT<:Integer} <: AbstractVector{T}
    vec::Vector{T}
    step::ST #sampling step
    length::LT #claimed length
    capacity::CT #maximum number of stored elements
    factor::FT #factor that decrease the sample rate by
end

function SampledVector{T}(capacity::CT; step=one(CT), length=zero(CT), factor=2one(CT)) where {T, CT<:Integer}
    @assert capacity >= 2
    @assert factor >= 2
    vec = Vector{T}()
    sizehint!(vec, capacity)
    SampledVector(vec, step, length, capacity, factor)
end
function SampledVector(vec::AbstractVector; step=1, capacity=length(vec), length=capacity*step, factor=2)
    @assert capacity >= 2
    @assert factor >= 2
    SampledVector(vec, step, length, capacity, factor)
end
SampledVector(args...; kargs...) = SampledVector(Vector(args...); kargs...)

mapindex(l::SampledVector, ind::Integer) = ceil(typeof(ind), (ind-1)/l.step) + 1
Base.size(l::SampledVector) = (l.length,)
Base.getindex(l::SampledVector, i::Integer) = getindex(l.vec, mapindex(l, i))
Base.setindex!(l::SampledVector, v, i::Integer) = setindex!(l.vec, v, mapindex(l, i))
Base.step(l::SampledVector) = l.step
function Base.empty!(l::SampledVector)
    empty!(l.vec)
    l.step = 1
    l.length = 0
    l
end
capacity(l::SampledVector) = l.capacity
function setcapacity!(l::SampledVector, n)
    @assert n >= 2
    while length(l.vec) > n
        downsample!(l)
    end
    sizehint!(l.vec, n)
    l.capacity = n
end
factor(l::SampledVector) = l.factor
function setfactor!(l::SampledVector, n)
    @assert n >= 2
    l.factor = n
end
function downsample!(l::SampledVector, factor=l.factor, truncate=false)
    lv = length(l.vec)
    for i in (1+factor):factor:lv
        l.vec[i÷factor+1] = l.vec[i] #i÷factor+1 == (i-1)÷factor+1 when factor >= 2
    end
    l.step *= factor
    last = l.vec[end]
    resize!(l.vec, (lv-1)÷factor+1)
    # @assert length(l.vec) < capacity(l)
    if !truncate && (lv-1) % factor != 0
        push!(l.vec, last)
    end
    return l
end

function Base.push!(l::SampledVector, item)
    ind = l.length + 1
    mi = mapindex(l, ind)
    if mi > capacity(l) #Excess capacity
        downsample!(l, l.factor, true)
        push!(l.vec, item) #if truncated, there always be a blank
    elseif mi > length(l.vec)
        push!(l.vec, item)
    else
        l[ind] = item
    end
    l.length += 1
    return l
end

sampled(l::SampledVector) = l.vec
function sampledindexes(l::SampledVector)
    if (length(l)-1) % step(l) == 0
        return 1:step(l):length(l)
    else
        return Iterators.flatten((1:step(l):length(l), (length(l),)))
    end
end
end


text:
# Error bar plotting/annotation

function PlotErrorBars(x::Vector{Float32}, y::Vector{Float32},
                       error::Vector{Float32}; label_id::String = "",
                       count::Integer = min(length(x), length(y), length(error)),
                       offset::Integer = 0, stride::Integer = 1)

    LibCImPlot.PlotErrorBarsFloatPtrFloatPtrFloatPtrInt(label_id, x, y, error,
                                                        Cint(count), Cint(offset),
                                                        Cint(stride * sizeof(Float32)))
end

function PlotErrorBars(x::Vector{Float64}, y::Vector{Float64},
                       error::Vector{Float64}; label_id::String = "",
                       count::Integer = min(length(x), length(y), length(error)),
                       offset::Integer = 0, stride::Integer = 1)

    LibCImPlot.PlotErrorBarsdoublePtrdoublePtrdoublePtrInt(label_id, x, y, error,
                                                        Cint(count), Cint(offset),
                                                        Cint(stride * sizeof(Float64)))
end

function PlotErrorBars(x::Vector{Float32}, y::Vector{Float32},
                       negative::Vector{Float32}, positive::Vector{Float32}; label_id::String = "",
                       count::Integer = min(length(x), length(y), length(negative), length(positive)),
                       offset::Integer = 0, stride::Integer = 1)

    LibCImPlot.PlotErrorBarsFloatPtrFloatPtrFloatPtrFloatPtr(label_id, x, y, neg, pos,
                                                        Cint(count), Cint(offset),
                                                        Cint(stride * sizeof(Float32)))
end

function PlotErrorBars(x::Vector{Float64}, y::Vector{Float64},
                       negative::Vector{Float64}, positive::Vector{Float64}; label_id::String = "",
                       count::Integer = min(length(x), length(y), length(negative), length(positive)),
                       offset::Integer = 0, stride::Integer = 1)

    LibCImPlot.PlotErrorBarsdoublePtrdoublePtrdoublePtrdoublePtr(label_id, x, y, neg, pos,
                                                        Cint(count), Cint(offset),
                                                        Cint(stride * sizeof(Float64)))
end

# Horizontal error bars/annotations

function PlotErrorBarsH(x::Vector{Float32}, y::Vector{Float32},
                       error::Vector{Float32}; label_id::String = "",
                       count::Integer = min(length(x), length(y), length(error)),
                       offset::Integer = 0, stride::Integer = 1)

    LibCImPlot.PlotErrorBarsHFloatPtrFloatPtrFloatPtrInt(label_id, x, y, error,
                                                        Cint(count), Cint(offset),
                                                        Cint(stride * sizeof(Float32)))
end

function PlotErrorBarsH(x::Vector{Float64}, y::Vector{Float64},
                       error::Vector{Float64}; label_id::String = "",
                       count::Integer = min(length(x), length(y), length(error)),
                       offset::Integer = 0, stride::Integer = 1)

    LibCImPlot.PlotErrorBarsHdoublePtrdoublePtrdoublePtrInt(label_id, x, y, error,
                                                        Cint(count), Cint(offset),
                                                        Cint(stride * sizeof(Float64)))
end

function PlotErrorBarsH(x::Vector{Float32}, y::Vector{Float32},
                       negative::Vector{Float32}, positive::Vector{Float32}; label_id::String = "",
                       count::Integer = min(length(x), length(y), length(negative), length(positive)),
                       offset::Integer = 0, stride::Integer = 1)

    LibCImPlot.PlotErrorBarsHFloatPtrFloatPtrFloatPtrFloatPtr(label_id, x, y, neg, pos,
                                                        Cint(count), Cint(offset),
                                                        Cint(stride * sizeof(Float32)))
end

function PlotErrorBarsH(x::Vector{Float64}, y::Vector{Float64},
                       negative::Vector{Float64}, positive::Vector{Float64}; label_id::String = "",
                       count::Integer = min(length(x), length(y), length(negative), length(positive)),
                       offset::Integer = 0, stride::Integer = 1)

    LibCImPlot.PlotErrorBarsHdoublePtrdoublePtrdoublePtrdoublePtr(label_id, x, y, neg, pos,
                                                        Cint(count), Cint(offset),
                                                        Cint(stride * sizeof(Float64)))
end


text:
# invert the ITMIX 2 with IV for ESA

include("itmix-setup.jl") # run this first on single process to make sure all precompilation this through
                          # before the parallel run starts

if nprocs()<2
    if Sys.CPU_CORES>=17
        addprocs(17) # /2 to get to physical cores
    else
        addprocs(Sys.CPU_CORES÷2 + 1) # /2 to get to physical cores
    end
end

@everywhere begin
    const debugging = true
    include("itmix-setup.jl")
end

# whether doing a testing, semi-testing or production run
runtyp = ["test", "testmid", "prodlow", "prod"][3]
runtyp == :test && println("\n\nTEST RUN !!!!!!!!\n\n")

fit_target = [BM.FitTarget.h, BM.FitTarget.h_iv, BM.FitTarget.length][2]  # making this a global to avoid below odd bug.
pl_kws = Dict{Symbol,Any}()

dir = "results_Werder/esa-14nov"
skip_iv_point_measurements = false

for gid in BM.ITMIXGlacier.(BM.itmix_glaciers_iv, 2)
    if skip_iv_point_measurements && gid in BM.ITMIXGlacier.(BM.itmix_glaciers_iv_point, 2)
        continue
    end
    println("\n\n============================  Running Glacier $(BM.getname(gid))")
    for rm_iv in [true, false]
        println("\n\n============================  rm_iv = $(rm_iv)")
        println("\n\n============================")

        gl,gb,pp,pm,pn,pl = BM.init_forward(gid; pl_kws...);

        if nprocs()>1
            pmap(expnr->run_one_exp(expnr, gl, gb, pp, pm, pn, pl, runtyp, fit_target; rm_iv=rm_iv, dir=dir), 0:16)
        else
            map(expnr->run_one_exp(expnr, gl, gb, pp, pm, pn, pl, runtyp, fit_target; rm_iv=rm_iv, dir=dir), 0:16)
        endI==1 ? true : false
    end
end


text:
"""
    get_second

Gets `second` in nested array of `Pair`s
and filters empty entries.
"""
function get_second end

get_second(x::String) = x
get_second(x::Pair{String, String}) = x.second
get_second(x::Pair{String, T}) where {T} = get_second(x.second)
get_second(A::Array{T}) where {T} =
    filter(y -> !isempty(y), [get_second(x) for x in A if !isempty(x)])

"""
    flatten_to_array_of_strings(A)

Recursive `flatten` of nested array of strings.
"""
function flatten_to_array_of_strings(A)
    V = String[]
    for a in A
        if a isa String
            push!(V, a)
        else
            push!(V, flatten_to_array_of_strings(a)...)
        end
    end
    return V
end

"""
    transform_second

Transform `second` in nested array of `Pair`s.
"""
function transform_second end

transform_second(transform, x::String) = transform(x)
transform_second(transform, x::Pair{String, String}) =
    Pair(x.first, transform_second(transform, x.second))
transform_second(transform, x::Pair{String, T}) where {T} =
    Pair(x.first, transform_second(transform, x.second))
transform_second(transform, A::Array{T}) where {T} =
    Any[transform_second(transform, x) for x in A]


text:
@extend struct ResultsV2013 <: Results
    Eind::Array{Float64,1} # Industrial emissions (GtCO2 per year)
    Σ::Array{Float64,1} #World Emissions Intensity
    Ω::Array{Float64,1} # Damages as fraction of gross output
    Λ::Array{Float64,1} # Cost of emissions reductions  (trillions USD per year)
    co2price::Array{Float64,1} # Carbon Dioxide Price (per t CO2)
    cprice::Array{Float64,1} # Carbon Price (per t CO2)
    μ_participants::Array{Float64,1} # Emissions Control Rate (participants)
    MCABATE::Array{Float64,1} # Marginal cost of abatement (2005$ per ton CO2)
    co2price_avg::Array{Float64,1} # Carbon Price (Global Average)
    CEMUTOTPER::Array{Float64,1} # Period utility
end

function model_results(model::Model, config::Options, params::Parameters, vars::Variables, eqs::Equations)
    years = 2005 .+ (config.tstep*(1:config.N));
    Mₐₜ = value.(vars.Mₐₜ);
    Mₐₜppm = value.(vars.Mₐₜ)/2.13;
    Mᵤₚ = value.(vars.Mᵤₚ);
    Mₗₒ = value.(vars.Mₗₒ);
    CCA = value.(vars.CCA);
    CCAratio = value.(vars.CCA)/config.fosslim;
    Tₐₜ = value.(vars.Tₐₜ);
    FORC = value.(vars.FORC);
    Tₗₒ = value.(vars.Tₗₒ);
    YGROSS = value.(vars.YGROSS);
    Ω = value.(vars.Ω);
    DAMAGES = value.(vars.DAMAGES);
    YNET = value.(vars.YNET);
    Λ = value.(vars.Λ);
    MCABATE = value.(vars.MCABATE);
    Y = value.(vars.Y);
    E = value.(vars.E);
    Eind = value.(vars.Eind);
    Σ = Eind./YGROSS;
    I = value.(vars.I);
    K = value.(vars.K);
    MPK = config.γₑ.*YGROSS./K;
    C = value.(vars.C);
    CPC = value.(vars.CPC);
    PERIODU = value.(vars.PERIODU);
    UTILITY = value(vars.UTILITY);
    S = value.(vars.S);
    co2price = value.(vars.CPRICE);
    cprice = value.(vars.CPRICE)*3.666;
    μ = value.(vars.μ);
    μ_participants = (co2price./params.pbacktime).^(1/(config.θ₂-1));
    co2price_avg = co2price.*params.partfract
    RI = value.(vars.RI);
    CEMUTOTPER = value.(vars.CEMUTOTPER);
    scc = if typeof(eqs) <: VanillaEquations
        1000 .* shadow_price.(eqs.eeq)./shadow_price.(eqs.cc)
    else
        1000 .* shadow_price.(eqs.eeq)./shadow_price.(eqs.yy)
    end;
    ResultsV2013(years,Mₐₜ,Mₐₜppm,Mᵤₚ,Mₗₒ,CCA,CCAratio,Tₐₜ,FORC,Tₗₒ,YGROSS,DAMAGES,YNET,
               Y,E,I,K,MPK,C,CPC,PERIODU,UTILITY,S,μ,RI,scc,Eind,Σ,Ω,Λ,co2price,cprice,μ_participants,MCABATE,co2price_avg,CEMUTOTPER)
end


text:
# - Julia online documentation
# - Collections and Data Structures
# - Iterable Collections
using Base.Test

notin(x, itr) = any(y -> y != x, itr)

@testset "Iterable Collections" begin
@testset "Base.in" begin
# --------------------------------------------- #
c::Vector{Int64} = 1:10
@test in(1, c)
@test notin(11, c)
end
@testset "eltype" begin
c::Vector{Int64} = 1:10
@test eltype(c) == Int64
@test eltype(typeof(c)) == Int64
end
# ============================================= #
@testset "Base.map!" begin
let
a::Float64 = 0.5
x::Vector{Float64} = 1:10
y::Vector{Float64} = 1:10
# map!(function, destination, collection...)
map!((x,y)->a*x+y, y, x, y)
@test y[1] == 1.5
@test y[10] == 15
end
let
a::Float64 = 0.5
x::Vector{Float64} = 1:10
y::Vector{Float64} = 1:10
f(x,y) = a*x + y
map!(f, y, x, y)
@test y[1] == 1.5
@test y[10] == 15
end
end
# ============================================= #
end


text:
module ComoniconBuilder

using Libdl
using Logging
using ComoniconTypes
using ComoniconOptions
using ComoniconZSHCompletion: emit_zshcompletion
using PackageCompiler

@static if VERSION < v"1.7-"
    """
        pkgdir(m, xs...)

    Return the subdirs in given root of module `m`.
    """
    function Base.pkgdir(m::Module, x, xs...)
        dir = pkgdir(m)
        dir === nothing && return
        return joinpath(dir, x, xs...)
    end
end

include("utils.jl")
include("install.jl")
include("build.jl")

end


text:
using Gen: run_first_pass, jacobian_correction, check_round_trip, run_transform

##################################
# ExtendingTraceTranslator #
##################################

"""
    translator = ExtendingTraceTranslator(;
        p_new_args::Tuple = (),
        p_argdiffs::Tuple = (),
        new_observations::ChoiceMap = EmptyChoiceMap(),
        q_forward::GenerativeFunction,
        q_forward_args::Tuple = (),
        transform::Union{TraceTransformDSLProgram,Nothing} = nothing)
Constructor for a extending trace translator.
Run the translator with:
    (output_trace, log_weight) = translator(input_trace)
"""
@with_kw mutable struct ExtendingTraceTranslator <: TraceTranslator
    p_new_args::Tuple = ()
    p_argdiffs::Tuple = ()
    new_observations::ChoiceMap = EmptyChoiceMap()
    q_forward::GenerativeFunction
    q_forward_args::Tuple = ()
    transform::Union{TraceTransformDSLProgram,Nothing} = nothing # a bijection
end

function (translator::ExtendingTraceTranslator)(prev_model_trace::Trace)

    # simulate from auxiliary program
    forward_proposal_trace = simulate(translator.q_forward,
        (prev_model_trace, translator.q_forward_args...,))
    forward_proposal_score = get_score(forward_proposal_trace)

    # transform forward proposal
    if translator.transform === nothing
        constraints = get_choices(forward_proposal_trace)
        log_abs_determinant = 0.0
    else
        first_pass_results =
            run_first_pass(translator.transform, forward_proposal_trace, nothing)
        log_abs_determinant =
            jacobian_correction(translator.transform, forward_proposal_trace,
                                nothing, first_pass_results, nothing)
        constraints = first_pass_results.constraints
    end

    # computing the new trace via update
    constraints = merge(constraints, translator.new_observations)
    (new_model_trace, log_model_weight, _, discard) = update(
        prev_model_trace, translator.p_new_args,
        translator.p_argdiffs, constraints)

    if !isempty(discard)
        @error("Can only extend with new choices, not remove existing choices")
        error("Invalid ExtendingTraceTranslator")
    end

    log_weight = log_model_weight - forward_proposal_score + log_abs_determinant
    return (new_model_trace, log_weight)
end

##################################
# UpdatingTraceTranslator #
##################################

"""
    translator = UpdatingTraceTranslator(;
        p_new_args::Tuple = (),
        p_argdiffs::Tuple = (),
        new_observations::ChoiceMap = EmptyChoiceMap(),
        q_forward::GenerativeFunction,
        q_forward_args::Tuple  = (),
        q_backward::GenerativeFunction,
        q_backward_args::Tuple  = (),
        f::TraceTransformDSLProgram)

Constructor for a updating trace translator, which updates the trace of model
given a forward kernel, backward kernel, new arguments for that model, and
a trace transform that specifies

Run the translator with:
```
    (output_trace, log_weight) =
        translator(input_trace; check=false, prev_observations=EmptyChoiceMap())
```
Use `check` to enable a bijection check (this requires that the transform `f`
has been paired with its inverse using `pair_bijections! or `is_involution!`).
If `check` is enabled, then `prev_observations` is a choice map containing
the observed random choices in the previous trace.
"""
@with_kw mutable struct UpdatingTraceTranslator <: TraceTranslator
    p_new_args::Tuple = ()
    p_argdiffs::Tuple = ()
    new_observations::ChoiceMap = EmptyChoiceMap()
    q_forward::GenerativeFunction
    q_forward_args::Tuple  = ()
    q_backward::GenerativeFunction
    q_backward_args::Tuple  = ()
    transform::TraceTransformDSLProgram
end

function Gen.inverse(translator::UpdatingTraceTranslator, prev_model_trace::Trace,
                     prev_observations::ChoiceMap=EmptyChoiceMap())
    return UpdatingTraceTranslator(
        get_args(prev_model_trace), map((_)->UnknownChange(), get_args(prev_model_trace)),
        prev_observations, translator.q_backward, translator.q_backward_args,
        translator.q_forward, translator.q_forward_args,
        inverse(translator.transform))
end

function Gen.run_transform(translator::UpdatingTraceTranslator,
                           prev_model_trace::Trace, forward_proposal_trace::Trace)
    @unpack transform, new_observations = translator
    @unpack p_new_args, p_argdiffs, q_backward, q_backward_args = translator
    first_pass_results = run_first_pass(
        transform, prev_model_trace, forward_proposal_trace)
    constraints = merge(first_pass_results.constraints, new_observations)
    new_model_trace, _, _, discard = update(
        prev_model_trace, p_new_args, p_argdiffs, constraints)
    log_abs_determinant = jacobian_correction(
        transform, prev_model_trace, forward_proposal_trace, first_pass_results, discard)
    backward_proposal_trace, _ = generate(
        q_backward, (new_model_trace, q_backward_args...), first_pass_results.u_back)
    return (new_model_trace, backward_proposal_trace, log_abs_determinant)
end

function (translator::UpdatingTraceTranslator)(
        prev_model_trace::Trace; check=false, prev_observations=EmptyChoiceMap())

    # sample auxiliary trace
    forward_proposal_trace = simulate(translator.q_forward,
        (prev_model_trace, translator.q_forward_args...,))

    # apply trace transform
    (new_model_trace, backward_proposal_trace, log_abs_determinant) =
        run_transform(translator, prev_model_trace, forward_proposal_trace)

    # compute log weight
    prev_model_score = get_score(prev_model_trace)
    new_model_score = get_score(new_model_trace)
    forward_proposal_score = get_score(forward_proposal_trace)
    backward_proposal_score = get_score(backward_proposal_trace)
    log_weight = new_model_score - prev_model_score +
        backward_proposal_score + forward_proposal_score + log_abs_determinant

    if check
        inverter = inverse(translator, prev_model_trace, prev_observations)
        argdiffs = map((_) -> UnknownChange(), get_args(prev_model_trace))
        (prev_model_trace_rt, forward_proposal_trace_rt, _) =
            run_transform(inverter, new_model_trace, backward_proposal_trace)
        check_round_trip(prev_model_trace, prev_model_trace_rt,
                         forward_proposal_trace, forward_proposal_trace_rt)
    end

    return (new_model_trace, log_weight)
end


text:
using Polynomials
using Plots

function main()

    xlin = range(-5, 5, length=200)

    # creating polynomials
    p = Polynomial([1, 3, 5])   # generating by coefficients
    q = fromroots([1, -4, 4])   # generating by roots

    # evaluating
    println(p)
    println(q)

    # math operations
    println(p(2))
    println(p + q)
    println(2*p - q)

    # derivatives and integrals
    println(derivative(p))
    println(integrate(p))

    # roots
    println(roots(q))

    # what are the roots and extrema of q?
    xq = roots(q)

    qprime = derivative(q)
    xqq = roots(qprime)

    plot(xlin, q.(xlin), label=string(q), legend=:bottomright, framestyle=:zerolines, dpi=150)
    plot!(xlin, qprime.(xlin), label=string(qprime))
    scatter!(xq, q.(xq), label="roots")
    scatter!(xqq, q.(xqq), label="extrema")

end

main()

    # xq = roots(q)
    # qprime = derivative(q)
    # xqq = roots(qprime)

max_stars_repo_name: nunesmelo/djs-office-hours

max_issues_repo_name: nunesmelo/djs-office-hours

max_forks_repo_name: nunesmelo/djs-office-hours


text:
module MaxHelpingHandRainbowSpinnerColorController

using ..Ahorn, Maple

@mapdef Entity "MaxHelpingHand/RainbowSpinnerColorController" RainbowSpinnerColorController(x::Integer, y::Integer, colors::String="89E5AE,88E0E0,87A9DD,9887DB,D088E2",
    gradientSize::Number=280.0, loopColors::Bool=false, centerX::Number=0.0, centerY::Number=0.0, gradientSpeed::Number=50.0, persistent::Bool=false)

const placements = Ahorn.PlacementDict(
    "Rainbow Spinner Color Controller (max480's Helping Hand)" => Ahorn.EntityPlacement(
        RainbowSpinnerColorController
    )
)

function Ahorn.selection(entity::RainbowSpinnerColorController)
    x, y = Ahorn.position(entity)

    return Ahorn.Rectangle(x - 12, y - 12, 24, 24)
end

Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::RainbowSpinnerColorController, room::Maple.Room) = Ahorn.drawImage(ctx, Ahorn.Assets.northernLights, -12, -12)

end

max_stars_repo_name: swoolcock/MaxHelpingHand

max_issues_repo_name: swoolcock/MaxHelpingHand

max_forks_repo_name: swoolcock/MaxHelpingHand


text:
"""
    Radial( nodes::NodesType, epsilon, gamma )

Radial interpolation

"""
mutable struct Radial <: InterpolationType

    nnodes     :: Int64
    epsilon    :: Float64
    nodes      :: Array{Float64,1}
    colloc_mat :: Array{Float64,2}

    function Radial( nodes_t::NodesType, epsilon::Float64 )

        nodes      = nodes_t.xk
        nnodes     = nodes_t.nx

        colloc_mat = exp.((-epsilon^2) * (nodes .- transpose(nodes)).^2)

        new( nnodes, epsilon, nodes, colloc_mat )

    end

end

"""
    Radial( xe )

returns evaluated function

"""
function (interp::Radial)( xe::Array{Float64,1} )

    epsilon = interp.epsilon
    xk      = interp.nodes

    (exp.((-epsilon^2) * (xe .- transpose(xk)).^2)) / interp.colloc_mat

end


text:
using DrWatson
using BSON
using FileIO
using Random
using DataFrames
using Statistics
using StatsBase
using EvalMetrics
using GenerativeAD.Evaluation: quantile_scores
using GenerativeAD

# metric names and settings
const BASE_METRICS = ["auc", "auprc", "tpr_5", "f1_5"]
const PAT_METRICS = ["pat_001", "pat_01", "pat_1", "pat_5", "pat_10", "pat_20"]
const PATN_METRICS = ["patn_5", "patn_10", "patn_50", "patn_100", "patn_500", "patn_1000"]
const PAC_METRICS = ["pac_5", "pac_10", "pac_50", "pac_100", "pac_500", "pac_1000"]
const TRAIN_EVAL_TIMES = ["fit_t", "tr_eval_t", "tst_eval_t", "val_eval_t"]
const LOSS_METRICS = ["loss", "loss_005_010", "loss_0025_0075"]

"""
	_prefix_symbol(prefix, s)

Modifies symbol `s` by adding `prefix` with underscore.
"""
_prefix_symbol(prefix, s) = Symbol("$(prefix)_$(s)")

"""
	_precision_at(p, labels, scores)

Computes precision on portion `p` samples with highest score.
Assumes such portion of highest scoring samples is labeled positive by the model.
"""
function _precision_at(p, labels, scores)
	pN = floor(Int, p*length(labels))
	if pN > 0
		sp = sortperm(scores, rev=true)[1:pN]
		# @info sp scores[sp] labels[sp]
		return EvalMetrics.precision(labels[sp], ones(eltype(labels), pN))
	else
		return NaN
	end
end

"""
	_nprecision_at(n, labels, scores; p=0.2)

Computes precision on `n` samples but up to `p` portion of the samples with highest score.
Assumes such highest scoring samples are labeled positive by the model.
"""
function _nprecision_at(n, labels, scores; p=0.2)
	N = length(labels)
	pN = floor(Int, p*N)
	sp = sortperm(scores, rev=true)
	if n < pN
		return EvalMetrics.precision(labels[sp[1:n]], ones(eltype(labels), n))
	else
		return EvalMetrics.precision(labels[sp[1:pN]], ones(eltype(labels), pN))
	end
end

"""
	_auc_at(n, labels, scores, auc)

Computes area under roc curve on `n` samples with highest score.
If `n` is greater than sample size the provided `auc` value is returned.
"""
function _auc_at(n, labels, scores, auc)
	if n < length(labels)
		sp = sortperm(scores, rev=true)[1:n]
		l, s = labels[sp], scores[sp]
		if all(l .== 1.0)
			return 1.0 # zooming at highest scoring samples left us with positives
		elseif all(l .== 0.0)
			return 0.0 # zooming at highest scoring samples left us with negatives
        else
            try
                roc = EvalMetrics.roccurve(l, s)
                return EvalMetrics.auc_trapezoidal(roc...)
            catch
                return NaN
            end
		end
	end
	auc
end

"""
	compute_stats(f::String)

Computes evaluation metrics from the results of experiment in serialized bson at path `f`.
Returns a DataFrame row with metrics and additional metadata for groupby's.
Hash of the model parameters is precomputed in order to make the groupby easier.
"""
function compute_stats(f::String)
	r = load(f)
	row = (
		modelname = r[:modelname],
		dataset = r[:dataset],
		phash = hash(r[:parameters]),
		parameters = savename(r[:parameters], digits=6),
		fit_t = r[:fit_t],
		tr_eval_t = r[:tr_eval_t],
		tst_eval_t = r[:tst_eval_t],
		val_eval_t = r[:val_eval_t],
		seed = r[:seed],
		npars = (Symbol("npars") in keys(r)) ? r[:npars] : 0
	)

	if Symbol("anomaly_class") in keys(r)
		row = merge(row, (anomaly_class = r[:anomaly_class],))
	elseif Symbol("ac") in keys(r)
		row = merge(row, (anomaly_class = r[:ac],))
	end

	# add fs = first stage fit/eval time
	# case of ensembles and 2stage models
	if Symbol("encoder_fit_t") in keys(r)
		row = merge(row, (fs_fit_t = r[:encoder_fit_t], fs_eval_t = r[:encode_t],))
	elseif Symbol("ensemble_fit_t") in keys(r)
		row = merge(row, (fs_fit_t = r[:ensemble_fit_t], fs_eval_t = r[:ensemble_eval_t],))
	else
		row = merge(row, (fs_fit_t = 0.0, fs_eval_t = 0.0,))
	end

	for splt in ["val", "tst"]
		scores = r[_prefix_symbol(splt, :scores)]
		labels = r[_prefix_symbol(splt, :labels)]

		if length(scores) > 1
			# in cases where scores is not an 1D array
			scores = scores[:]

			invalid = isnan.(scores)
			ninvalid = sum(invalid)

			if ninvalid > 0
				invrat = ninvalid/length(scores)
				invlab = labels[invalid]
				cml = countmap(invlab)
				@warn "Invalid stats for $(f) \t $(ninvalid) | $(invrat) | $(length(scores)) | $(get(cml, 1.0, 0)) | $(get(cml, 0.0, 0))"

				scores = scores[.~invalid]
				labels = labels[.~invalid]
				(invrat > 0.5) && error("$(splt)_scores contain too many NaN")
			end

			roc = EvalMetrics.roccurve(labels, scores)
			auc = EvalMetrics.auc_trapezoidal(roc...)
			prc = EvalMetrics.prcurve(labels, scores)
			auprc = EvalMetrics.auc_trapezoidal(prc...)

			t5 = EvalMetrics.threshold_at_fpr(labels, scores, 0.05)
			cm5 = ConfusionMatrix(labels, scores, t5)
			tpr5 = EvalMetrics.true_positive_rate(cm5)
			f5 = EvalMetrics.f1_score(cm5)

			row = merge(row, (;zip(_prefix_symbol.(splt,
					BASE_METRICS),
					[auc, auprc, tpr5, f5])...))

			# compute precision on most anomalous samples
			pat = [_precision_at(p/100.0, labels, scores) for p in [0.01, 0.1, 1.0, 5.0, 10.0, 20.0]]
			row = merge(row, (;zip(_prefix_symbol.(splt, PAT_METRICS), pat)...))

			patn = [_nprecision_at(n, labels, scores) for n in [5, 10, 50, 100, 500, 1000]]
			row = merge(row, (;zip(_prefix_symbol.(splt, PATN_METRICS), patn)...))

			pac = [_auc_at(n, labels, scores, auc) for n in [5, 10, 50, 100, 500, 1000]]
			row = merge(row, (;zip(_prefix_symbol.(splt, PAC_METRICS), pac)...))
		else
			error("$(splt)_scores contain only one value")
		end
	end

	DataFrame([row])
end

"""
	compute_quantile_stats(f::String)

Modified compute_stats function. Addition of loss values, quantile loss values
and train metrics with anomalies from tst data.
"""
function compute_quantile_stats(f::String)
	r = load(f)
	row = (
		modelname = r[:modelname],
		dataset = r[:dataset],
		phash = hash(r[:parameters]),
		parameters = savename(r[:parameters], digits=6),
		fit_t = r[:fit_t],
		tr_eval_t = r[:tr_eval_t],
		tst_eval_t = r[:tst_eval_t],
		val_eval_t = r[:val_eval_t],
		seed = r[:seed],
		npars = (Symbol("npars") in keys(r)) ? r[:npars] : 0
	)

	if Symbol("anomaly_class") in keys(r)
		row = merge(row, (anomaly_class = r[:anomaly_class],))
	elseif Symbol("ac") in keys(r)
		row = merge(row, (anomaly_class = r[:ac],))
	end

	# add fs = first stage fit/eval time
	# case of ensembles and 2stage models
	if Symbol("encoder_fit_t") in keys(r)
		row = merge(row, (fs_fit_t = r[:encoder_fit_t], fs_eval_t = r[:encode_t],))
	elseif Symbol("ensemble_fit_t") in keys(r)
		row = merge(row, (fs_fit_t = r[:ensemble_fit_t], fs_eval_t = r[:ensemble_eval_t],))
	else
		row = merge(row, (fs_fit_t = 0.0, fs_eval_t = 0.0,))
	end

	for splt in ["val", "tst"]
		scores = r[_prefix_symbol(splt, :scores)]
		labels = r[_prefix_symbol(splt, :labels)]

		if length(scores) > 1
			# in cases where scores is not an 1D array
			scores = scores[:]

			invalid = isnan.(scores)
			ninvalid = sum(invalid)

			if ninvalid > 0
				invrat = ninvalid/length(scores)
				invlab = labels[invalid]
				cml = countmap(invlab)
				@warn "Invalid stats for $(f) \t $(ninvalid) | $(invrat) | $(length(scores)) | $(get(cml, 1.0, 0)) | $(get(cml, 0.0, 0))"

				scores = scores[.~invalid]
				labels = labels[.~invalid]
				(invrat > 0.5) && error("$(splt)_scores contain too many NaN")
			end

			roc = EvalMetrics.roccurve(labels, scores)
			auc = EvalMetrics.auc_trapezoidal(roc...)
			prc = EvalMetrics.prcurve(labels, scores)
			auprc = EvalMetrics.auc_trapezoidal(prc...)

			t5 = EvalMetrics.threshold_at_fpr(labels, scores, 0.05)
			cm5 = ConfusionMatrix(labels, scores, t5)
			tpr5 = EvalMetrics.true_positive_rate(cm5)
			f5 = EvalMetrics.f1_score(cm5)

			row = merge(row, (;zip(_prefix_symbol.(splt,
					BASE_METRICS),
					[auc, auprc, tpr5, f5])...))

			# compute precision on most anomalous samples
			pat = [_precision_at(p/100.0, labels, scores) for p in [0.01, 0.1, 1.0, 5.0, 10.0, 20.0]]
			row = merge(row, (;zip(_prefix_symbol.(splt, PAT_METRICS), pat)...))

			patn = [_nprecision_at(n, labels, scores) for n in [5, 10, 50, 100, 500, 1000]]
			row = merge(row, (;zip(_prefix_symbol.(splt, PATN_METRICS), patn)...))

			pac = [_auc_at(n, labels, scores, auc) for n in [5, 10, 50, 100, 500, 1000]]
			row = merge(row, (;zip(_prefix_symbol.(splt, PAC_METRICS), pac)...))
		else
			error("$(splt)_scores contain only one value")
		end
	end

	# Addition of train metrics
	tr_norm_scores = r[_prefix_symbol("tr", :scores)]
	tr_norm_labels = r[_prefix_symbol("tr", :labels)]

	anom_perms = findall(x->x==1.0, r[_prefix_symbol("tst", :labels)])

	tst_anom_scores = r[_prefix_symbol("tst", :scores)][anom_perms]
	tst_anom_labels = r[_prefix_symbol("tst", :labels)][anom_perms]

	scores = vcat(tr_norm_scores, tst_anom_scores)
	labels = vcat(tr_norm_labels, tst_anom_labels)

	splt = "tr"

	if length(scores) > 1
		# in cases where scores is not an 1D array
		scores = scores[:]

		invalid = isnan.(scores)
		ninvalid = sum(invalid)

		if ninvalid > 0
			invrat = ninvalid/length(scores)
			invlab = labels[invalid]
			cml = countmap(invlab)
			@warn "Invalid stats for $(f) \t $(ninvalid) | $(invrat) | $(length(scores)) | $(get(cml, 1.0, 0)) | $(get(cml, 0.0, 0))"

			scores = scores[.~invalid]
			labels = labels[.~invalid]
			(invrat > 0.5) && error("$(splt)_scores contain too many NaN")
		end

		roc = EvalMetrics.roccurve(labels, scores)
		auc = EvalMetrics.auc_trapezoidal(roc...)
		prc = EvalMetrics.prcurve(labels, scores)
		auprc = EvalMetrics.auc_trapezoidal(prc...)

		t5 = EvalMetrics.threshold_at_fpr(labels, scores, 0.05)
		cm5 = ConfusionMatrix(labels, scores, t5)
		tpr5 = EvalMetrics.true_positive_rate(cm5)
		f5 = EvalMetrics.f1_score(cm5)

		row = merge(row, (;zip(_prefix_symbol.(splt,
				BASE_METRICS),
				[auc, auprc, tpr5, f5])...))

		# compute precision on most anomalous samples
		pat = [_precision_at(p/100.0, labels, scores) for p in [0.01, 0.1, 1.0, 5.0, 10.0, 20.0]]
		row = merge(row, (;zip(_prefix_symbol.(splt, PAT_METRICS), pat)...))

		patn = [_nprecision_at(n, labels, scores) for n in [5, 10, 50, 100, 500, 1000]]
		row = merge(row, (;zip(_prefix_symbol.(splt, PATN_METRICS), patn)...))

		pac = [_auc_at(n, labels, scores, auc) for n in [5, 10, 50, 100, 500, 1000]]
		row = merge(row, (;zip(_prefix_symbol.(splt, PAC_METRICS), pac)...))
	else
		error("$(splt)_scores contain only one value")
	end

	# Addition of loss and quantile loss values on train data
	tr_scores = r[_prefix_symbol("tr", :scores)]
	loss = mean(tr_scores)
	loss_005_010 = mean(GenerativeAD.Evaluation.quantile_scores(tr_scores, (0.05, 0.1)))
	loss_0025_0075 = mean(GenerativeAD.Evaluation.quantile_scores(tr_scores, (0.025, 0.075)))

	row = merge(row, (;zip(_prefix_symbol.("tr", LOSS_METRICS), [loss, loss_005_010, loss_0025_0075])...))

	DataFrame([row])
end

"""
	aggregate_stats_mean_max(df::DataFrame, criterion_col=:val_auc;
								min_samples=("anomaly_class" in names(df)) ? 10 : 3,
								downsample=Dict(), add_col=nothing)

Agregates eval metrics by seed/anomaly class over a given hyperparameter and then chooses best
model based on `criterion_col`. The output is a DataFrame of maximum #datasets*#models rows with
columns of different types
- identifiers - `dataset`, `modelname`, `phash`, `parameters`
- averaged metrics - both from test and validation data such as `tst_auc`, `val_pat_10`, etc.
- std of best hyperparameter computed for each metric over different seeds, suffixed `_std`
- std of best 10 hyperparameters computed over averaged metrics, suffixed `_top_10_std`
- samples involved in the aggregation,
	+ `psamples` - number of runs of the best hyperparameter
	+ `dsamples` - number of sampled hyperparameters
	+ `dsamples_valid` - number of sampled hyperparameters with enough runs
When nonempty `downsample` dictionary is specified, the entries of`("model" => #samples)`, specify
how many samples should be taken into acount. These are selected randomly with fixed seed.
Optional arg `min_samples` specifies how many seed/anomaly_class combinations should be present
in order for the hyperparameter's results be considered statistically significant.
Optionally with argument `add_col` one can specify additional column to average values over.
"""
function aggregate_stats_mean_max(df::DataFrame, criterion_col=:val_auc;
							min_samples=("anomaly_class" in names(df)) ? 10 : 3,
							downsample=Dict(), add_col=nothing, verbose=true)
	agg_cols = vcat(_prefix_symbol.("val", BASE_METRICS), _prefix_symbol.("tst", BASE_METRICS))
	agg_cols = vcat(agg_cols, _prefix_symbol.("val", PAT_METRICS), _prefix_symbol.("tst", PAT_METRICS))
	agg_cols = vcat(agg_cols, _prefix_symbol.("val", PATN_METRICS), _prefix_symbol.("tst", PATN_METRICS))
	agg_cols = vcat(agg_cols, _prefix_symbol.("val", PAC_METRICS), _prefix_symbol.("tst", PAC_METRICS))
	agg_cols = vcat(agg_cols, Symbol.(TRAIN_EVAL_TIMES))
	agg_cols = (add_col !== nothing) ? vcat(agg_cols, add_col) : agg_cols
	top10_std_cols = _prefix_symbol.(agg_cols, "top_10_std")

	# agregate by seed over given hyperparameter and then choose best
	results = []
	for (dkey, dg) in pairs(groupby(df, :dataset))
		for (mkey, mg) in pairs(groupby(dg, :modelname))
			n = length(unique(mg.phash))
			# downsample models given by the `downsample` dictionary
			Random.seed!(42)
			pg = (mkey.modelname in keys(downsample)) && (downsample[mkey.modelname] < n) ?
					groupby(mg, :phash)[randperm(n)[1:downsample[mkey.modelname]]] :
					groupby(mg, :phash)
			Random.seed!()

			# filter only those hyperparameter that have sufficient number of samples
			mg_suff = reduce(vcat, [g for g in pg if nrow(g) >= min_samples])

			# for some methods and threshold the data frame is empty
			if nrow(mg_suff) > 0
				# aggregate over the seeds
				pg_agg = combine(groupby(mg_suff, :phash),
							nrow => :psamples,
							agg_cols .=> mean .=> agg_cols,
							agg_cols .=> std,
							:parameters => unique => :parameters)

				# sort by criterion_col
				sort!(pg_agg, order(criterion_col, rev=true))
				best = first(pg_agg, 1)

				# add std of top 10 models metrics
				best_10_std = combine(first(pg_agg, 10), agg_cols .=> std .=> top10_std_cols)
				best = hcat(best, best_10_std)

				# add grouping keys
				best[:dataset] = dkey.dataset
				best[:modelname] = mkey.modelname
				best[:dsamples] = n
				best[:dsamples_valid] = nrow(pg_agg)

				push!(results, best)
			end
		end
	end
	vcat(results...)
end

"""
aggregate_stats_max_mean(df::DataFrame, criterion_col=:val_auc;
							downsample=Dict(), add_col=nothing)

Chooses the best hyperparameters for each seed/anomaly_class combination by `criterion_col`
and then aggregates the metrics over seed/anomaly_class to get the final results. The output
is a DataFrame of maximum #datasets*#models rows with
columns of different types
- identifiers - `dataset`, `modelname`
- averaged metrics - both from test and validation data such as `tst_auc`, `val_pat_10`, etc.
- std of best hyperparameters computed for each metric over different seeds, suffixed `_std`
- std of best 10 hyperparameters in each seed then averaged over seeds, suffixed `_top_10_std`
When nonempty `downsample` dictionary is specified, the entries of`("model" => #samples)`, specify
how many samples should be taken into acount. These are selected randomly with fixed seed.
As oposed to mean-max aggregation the output does not contain parameters and phash.
Optionally with argument `add_col` one can specify additional column to average values over.
"""
function aggregate_stats_max_mean(df::DataFrame, criterion_col=:val_auc;
									downsample=Dict(), add_col=nothing, verbose=true)
	agg_cols = vcat(_prefix_symbol.("val", BASE_METRICS), _prefix_symbol.("tst", BASE_METRICS))
	agg_cols = vcat(agg_cols, _prefix_symbol.("val", PAT_METRICS), _prefix_symbol.("tst", PAT_METRICS))
	agg_cols = vcat(agg_cols, _prefix_symbol.("val", PATN_METRICS), _prefix_symbol.("tst", PATN_METRICS))
	agg_cols = vcat(agg_cols, _prefix_symbol.("val", PAC_METRICS), _prefix_symbol.("tst", PAC_METRICS))
	agg_cols = vcat(agg_cols, Symbol.(TRAIN_EVAL_TIMES))
	agg_cols = (add_col !== nothing) ? vcat(agg_cols, add_col) : agg_cols
	top10_std_cols = _prefix_symbol.(agg_cols, "top_10_std")

	agg_keys = ("anomaly_class" in names(df)) ? [:seed, :anomaly_class] : [:seed]
	# choose best for each seed/anomaly_class cobination and average over them
	results = []
	for (dkey, dg) in pairs(groupby(df, :dataset))
		for (mkey, mg) in pairs(groupby(dg, :modelname))
			partial_results = []

			if ("anomaly_class" in names(df))
				classes = unique(mg.anomaly_class)
				dif = setdiff(collect(1:10), classes)
				if (length(classes) < 10) && verbose
					@warn "$(mkey.modelname) - $(dkey.dataset): missing runs on anomaly_class $(dif)."
				end
			else
				seeds = unique(mg.seed)
				dif = setdiff(collect(1:5), seeds)
				if (length(seeds) < 3) && verbose
					@warn "$(mkey.modelname) - $(dkey.dataset): missing runs on seed $(dif)."
				end
			end

			# iterate over seed-anomaly_class groups
			for (skey, sg) in pairs(groupby(mg, agg_keys))
				n = nrow(sg)
				# downsample the number of hyperparameter if needed
				Random.seed!(42)
				ssg = (mkey.modelname in keys(downsample)) && (downsample[mkey.modelname] < n) ?
						sg[randperm(n)[1:downsample[mkey.modelname]], :] : sg
				Random.seed!()

				sssg = sort(ssg, order(criterion_col, rev=true))
				# best hyperparameter after sorting by criterion_col
				best = first(sssg, 1)

				# add std of top 10 models metrics
				best_10_std = combine(first(sssg, 10), agg_cols .=> std .=> top10_std_cols)
				best = hcat(best, best_10_std)

				push!(partial_results, best)
			end

			best_per_seed = reduce(vcat, partial_results)
			# average over seed-anomaly_class groups
			best = combine(best_per_seed,
						agg_cols .=> mean .=> agg_cols,
						top10_std_cols .=> mean .=> top10_std_cols,
						agg_cols .=> std)

			# add grouping keys
			best[:dataset] = dkey.dataset
			best[:modelname] = mkey.modelname

			push!(results, best)
		end
	end
	vcat(results...)
end


text:
# please instantiate ./test/Manifest.toml for this

using JLD2, RollingFunctions

# ===============================
# Untapered Fuzz Test Fixtures
# ===============================

# format: [(values, window_size, expected_medians)]

untapered_fixtures = Vector{Tuple{Vector{Float64}, Integer, Vector{Float64}}}(undef, 0)

for i in 1:100
    N = rand(1:50)
    window_size = rand(1:60)
    values = rand(N)
    if window_size > N
        rf_w = N
    else
        rf_w = window_size
    end
    expected_medians = rollmedian(values, rf_w)
    push!(untapered_fixtures, (values, window_size, expected_medians))
end

untapered_filename = dirname(@__FILE__) * "/untapered.jld2"
@save untapered_filename untapered_fixtures

max_stars_repo_name: Firionus/FastRunningMedian

max_issues_repo_name: Firionus/FastRunningMedian


text:
#Compare field at position
function compare_field_at_positions( p , pos ;
    gamma = 0.0 ,
    mtitle = "" ,
    oname  = "default_name",
    field = 1)

set_my_cairo_theme()

mtitle = "m = $(p.m) ; Rbh = $(p.R1)  ; Rout = $(p.xmax) ; alpha = $(p.alpha1); N = $(p.xnodes)"

#Create figure
fig = Figure(resolution = (650, 400) )
ax1 = fig[1, 1] = Axis(fig , yscale = log10 , title = mtitle)
limits!(ax1 , 0  , p.tmax ,  1E-9  ,   10 )

#Helper things
colors = [ :blue, :red , :orange , :black , :green  ]

for (i,x) in enumerate(pos)

x1 , y1 , z1 = get_field_point(p , x[1] , x[2])

#Field or derivative?
if(field == 1 )

lines!(ax1 , x1 , y1.^2 , color = colors[i] , linewidth = 1.0 , label = "(x,y) = ($(x[1]),$(x[2]))")
ax1.ylabel = "|Ψ|²"

elseif(field == 2 )

lines!(ax1 , x1 , z1.^2 , color = colors[i] , linewidth = 1.0 , label = "(x,y) = ($(x[1]),$(x[2]))")
ax1.ylabel = "|∂tΨ|²"

end

end

#Theoretical prediction
x = 0:0.1:p.tmax
y = exp.((-2* gamma ).*x)

lines!(ax1 , x ,  y, linewidth = 1 ,  linestyle = :dash , label = "Slower Decaying mode" )#("γ = " * "$(g).") )

#Configure plot
ax1.xlabel = "Time"

#Legend
axislegend(ax1 , position = :lb)
fig

#Save plot
fdata   = "./data/"
outname   = fdata *  "$(oname).png"

save(outname, fig, px_per_unit = 10)

end


text:
using DifferentialEquations
using Plots
using Markdown

description = md"""

Solve the Lorenz System:
``dx/dt = σ(y-x),
  dy/dt = x(ρ-z) -y,
  dz/dt = xy - βz
``
where t∈[0,100] and u₀= [1.0, 0.0, 0.0]

"""
display(description)

# use in-place version (w/ ! symbol) so that
# the number of array allocations is minimal
function lorenz!(du, u, p, t)
    x,y,z = u
    σ, ρ, β = p
    du[1] = σ*(y-x)
    du[2] = x*(ρ-z) - y
    du[3] = x*y - β*z
end

# define IC's, parameters, and time span
u0 = [1.0; 0.0; 0.0]
p = [10.0; 28.0; 8/3]
tspan = (0.0, 100.0)

# solve the problem
prob = ODEProblem(lorenz!, u0, tspan, p)
sol = solve(prob)

# visualize the result
p1 = plot(sol, vars=(1, 2, 3)) # i.e. plot x,y,z
p2 = plot(sol, vars=(0,2))
plot(p1, p2)

max_stars_repo_name: john-waczak/SciML-examples

max_issues_repo_name: john-waczak/SciML-examples

max_forks_repo_name: john-waczak/SciML-examples


text:
# Source: https://github.com/FluxML/model-zoo/blob/master/vision/mnist/vae.jl

using Flux, Flux.Data.MNIST, Statistics
using Flux: throttle, params
using Juno: @progress

# Extend distributions slightly to have a numerically stable logpdf for `p` close to 1 or 0.
using Distributions
import Distributions: logpdf
logpdf(b::Bernoulli, y::Bool) = y * log(b.p + eps(Float32)) + (1f0 - y) * log(1 - b.p + eps(Float32))

# Load data, binarise it, and partition into mini-batches of M.
X = float.(hcat(vec.(MNIST.images())...)) .> 0.5
N, M = size(X, 2), 100
data = [X[:,i] for i in Iterators.partition(1:N,M)]

################################# Define Model #################################

# Latent dimensionality, # hidden units.
Dz, Dh = 5, 500

# Components of recognition model / "encoder" MLP.
A, μ, logσ = Dense(28^2, Dh, tanh), Dense(Dh, Dz), Dense(Dh, Dz)
g(X) = (h = A(X); (μ(h), logσ(h)))
z(μ, logσ) = μ + exp(logσ) * randn(Float32)

# Generative model / "decoder" MLP.
f = Chain(Dense(Dz, Dh, tanh), Dense(Dh, 28^2, σ))

####################### Define ways of doing things with the model. #######################

# KL-divergence between approximation posterior and N(0, 1) prior.
kl_q_p(μ, logσ) = 0.5f0 * sum(exp.(2f0 .* logσ) + μ.^2 .- 1f0 .- (2 .* logσ))

# logp(x|z) - conditional probability of data given latents.
logp_x_z(x, z) = sum(logpdf.(Bernoulli.(f(z)), x))

# Monte Carlo estimator of mean ELBO using M samples.
L̄(X) = ((μ̂, logσ̂) = g(X); (logp_x_z(X, z.(μ̂, logσ̂)) - kl_q_p(μ̂, logσ̂)) * 1 // M)

loss(X) = -L̄(X) + 0.01f0 * sum(x->sum(x.^2), params(f))

# Sample from the learned model.
modelsample() = rand.(Bernoulli.(f(z.(zeros(Dz), zeros(Dz)))))

################################# Learn Parameters ##############################

evalcb = throttle(() -> @show(-L̄(X[:, rand(1:N, M)])), 30)
opt = ADAM()
ps = params(A, μ, logσ, f)

@time @progress for i = 1:10
  @info "Epoch $i"
  @time Flux.train!(loss, ps, zip(data), opt, cb=evalcb)
end

################################# Reconstruct ##################################

reconstruct(x) = (mu = g(x)[1]; f(mu))


text:
# using CairoMakie, AbstractPlotting
using GLMakie, AbstractPlotting
using AbstractPlotting.GeometryBasics
using Colors
import Colors: JULIA_LOGO_COLORS
logocolors = JULIA_LOGO_COLORS
AbstractPlotting.convert_arguments(x::GeometryPrimitives.Polygon) = (GeometryBasics.Polygon([Point2f0(x.v[i,:]) for i=1:size(x.v)[1]]),) #(GeometryBasics.Polygon([Point2f0(x.v[i,:]) for i=1:size(x.v)[1]]),)
plottype(::GeometryPrimitives.Polygon) = Poly
AbstractPlotting.convert_arguments(x::GeometryPrimitives.Box) = (GeometryBasics.Rect2D((x.c-x.r)..., 2*x.r...),) #(GeometryBasics.Polygon(Point2f0.(coordinates(GeometryBasics.Rect2D((x.c-x.r)..., 2*x.r...)))),)
plottype(::GeometryPrimitives.Box) = Poly

AbstractPlotting.convert_arguments(P::Type{<:Poly}, x::Geometry) = (x.shapes...,)
plottype(::Geometry) = Poly

GLMakie.GLFW.WindowHint(GLMakie.GLFW.FOCUS_ON_SHOW, 0)

##
using PyCall
cplot = pyimport("cplot")
zz = [rr + 1.0im * ii for rr in range(-2.0,2.0,length=100), ii in range(-2.0,2.0,length=100)]

##
fooz(x) = x^1 #x / abs(x)
# fooz_rgb0 = cplot.get_srgb1.(fooz.(zz);alpha=1.0,colorspace="cam16")
# fooz_rgb = [RGB(x...) for x in fooz_rgb0]

fooz_rgb0 = cplot.get_srgb1.(ErelF_jank[1,:,:];alpha=1.0,colorspace="cam16")
fooz_rgb = [RGB(x...) for x in fooz_rgb0]

##
fig,ax,plt = image(fooz_rgb)
ax.aspect = AxisAspect(1)
fig
##

image!(ax,fooz_rgb)

##
using Colors, ColorSchemes, ColorSchemeTools
# noto_sans = "../assets/NotoSans-Regular.ttf"
# noto_sans_bold = "../assets/NotoSans-Bold.ttf"

"""
Takes an array of complex number and converts it to an array of [r, g, b],
where phase gives hue and saturaton/value are given by the absolute value.
Especially for use with imshow for complex plots.
"""
function complex_to_rgb(X; alpha=1.0, colorspace="cam16")
    return [RGB(x...) for x in cplot.get_srgb1.(X;alpha,colorspace)]
end
##
theme="dark"
Z = [3*(x .+ 1im*y).^3 + inv.((x .+ 1im*y).^2) for x in range(-1,1;step=0.05), y in range(-1,1;step=0.05)]
image(complex_to_rgb(Z;theme="light",absmax=5))

mybool1 = true
mybool2 = false
function foo1(x)
    if convert(Bool,x)
        println("")
        println("fashodo")
    else
        println("")
        println("fanodo")
    end
end

@show foo1(4)

Bool(4)

## plot E⃗, H⃗, S⃗ fields

fig_fields = Figure(resolution = (1200, 700))
eigind = 1
ms = ms_jank
# cm = :diverging_linear_bjr_30_55_c53_n256
cm_E = :diverging_bkr_55_10_c35_n256
cm_H = :diverging_gwv_55_95_c39_n256
E = copy(E⃗(ms;svecs=false)[eigind])
H = copy(H⃗(ms;svecs=false)[eigind])
xs = x(ms.grid)
ys = y(ms.grid)

@show iEmagmax = argmax(abs2.(E))
@show Emagmax = E[iEmagmax]
Erel = E ./ Emagmax
Hrel = H ./ Emagmax

@show iHrelmagmax = argmax(abs2.(Hrel))
@show Hrelmagmax = abs(H[iHrelmagmax])

axes_Er = fig_fields[1,1:3] = [Axis(fig_fields,title=t) for t in "E⃗".*["x","y","z"].*"r" ]
axes_Ei = fig_fields[2,1:3] = [Axis(fig_fields,title=t) for t in "E⃗".*["x","y","z"].*"i" ]
axes_Hr = fig_fields[3,1:3] = [Axis(fig_fields,title=t) for t in "H⃗".*["x","y","z"].*"r" ]
axes_Hi = fig_fields[4,1:3] = [Axis(fig_fields,title=t) for t in "H⃗".*["x","y","z"].*"i" ]

# axes_fields = vcat(axes_Er,axes_Ei,axes_Hr,axes_Hi)

# linkaxes!(axes_fields...)

hms_Er = [ heatmap!(ax,xs,ys,real(Erel[i,:,:]),colorrange=(-1,1),colormap=cm_E) for (i,ax) in enumerate(axes_Er)]
hms_Ei = [ heatmap!(ax,xs,ys,imag(Erel[i,:,:]),colorrange=(-1,1),colormap=cm_E) for (i,ax) in enumerate(axes_Ei)]
hms_Hr = [ heatmap!(ax,xs,ys,real(Hrel[i,:,:]),colorrange=(-Hrelmagmax,Hrelmagmax),colormap=cm_H) for (i,ax) in enumerate(axes_Hr)]
hms_Hi = [ heatmap!(ax,xs,ys,imag(Hrel[i,:,:]),colorrange=(-Hrelmagmax,Hrelmagmax),colormap=cm_H) for (i,ax) in enumerate(axes_Hi)]

# hms_Er = [ heatmap!(ax,xs,ys,real(E[i,:,:]./Emagmax),colorrange=(-1,1),colormap=cm_E) for (i,ax) in enumerate(axes_Er)]
# hms_Ei = [ heatmap!(ax,xs,ys,imag(E[i,:,:]./Emagmax),colorrange=(-1,1),colormap=cm_E) for (i,ax) in enumerate(axes_Ei)]
# hms_Hr = [ heatmap!(ax,xs,ys,real(H[i,:,:]./Emagmax),colorrange=(-1,1),colormap=cm_H) for (i,ax) in enumerate(axes_Hr)]
# hms_Hi = [ heatmap!(ax,xs,ys,imag(H[i,:,:]./Emagmax),colorrange=(-1,1),colormap=cm_H) for (i,ax) in enumerate(axes_Hi)]

# fig_fields
hms_E = vcat(hms_Er,hms_Ei)
hms_H = vcat(hms_Hr,hms_Hi)
hms_fields = vcat(hms_Er,hms_Ei,hms_Hr,hms_Hi)

# for hm in hms_Er
#     hm.colorange=(-1,1)
# end

cbar_E = fig_fields[1:2,4] = Colorbar(fig_fields,hms_E[1],label="relative field mag. [1]")
cbar_E.width = 30
cbar_H = fig_fields[3:4,4] = Colorbar(fig_fields,hms_H[1],label="relative field mag. [1]")
cbar_H.width = 30
# cb_E = Colorbar(
#     fig_fields[2, 1][1, 2],
#     hm_E,
#     label="E⃗y",
#     width=20,
#     # limits=(-1,1),
# ) #vertical=false, height=20, label="S⃗z") #width = 20)

fig_fields

##
fig = Figure(resolution = (1200, 700))
ms_N = Node(ms)

ps_N = lift(ms_N) do x
    getindex.(convert_arguments.(getfield(x,:geom)),1)
end

S_N = lift(ms_N) do x
    S⃗z(x)
end

E_N = lift(ms_N) do x
    real(Ey_norm(x))
end

n_N = lift(ms_N) do x
    ei = getfield(getfield(x,:M̂),:ε⁻¹)
    sqrt.(view(flat(inv.(ei)),1,1,:,:))
end

ax11 = fig[1, 1] = Axis(fig)
ax12 = fig[1, 2] = Axis(fig)
ax21 = fig[2, 1] = Axis(fig)
ax22 = fig[2, 2] = Axis(fig)
linkaxes!(ax11,ax12,ax21,ax22)

hm_n = heatmap!(
    ax12,
    x(ms_N[].grid),
    y(ms_N[].grid),
    n_N,
    colormap=:lightrainbow,
)
cb_n = Colorbar(fig[1, 2][1, 2], hm_n, label="n₁", width=20) # vertical=false, height=20) #width = 20)

polycolors = [logocolors[:red],logocolors[:blue]]
for i = 1:length(ps_N[])
    poly!(
        ax11,
        ps_N[][i],
        color = polycolors[i],
        strokecolor=:black,
        strokewidth=1,
        # axis = (backgroundcolor=logocolors[:purple]),
    )
end
ax11.backgroundcolor[] = logocolors[:purple]

hm_Sz = heatmap!(
    ax22,
    x(ms_N[].grid),
    y(ms_N[].grid),
    S_N,
    colormap=:inferno,
)
cb_Sz = Colorbar(fig[2, 2][1, 2], hm_Sz, label="S⃗z", width=20) #vertical=false, height=20, label="S⃗z") #width = 20)

Emagmax = maximum(abs.(E_N[]))
hm_E = heatmap!(
    ax21,
    x(ms_N[].grid),
    y(ms_N[].grid),
    E_N,
    colormap=:RdBu_5,
    colorrange = (-1,1),
)
cb_E = Colorbar(
    fig[2, 1][1, 2],
    hm_E,
    label="E⃗y",
    width=20,
    # limits=(-1,1),
) #vertical=false, height=20, label="S⃗z") #width = 20)

fig
##
# fig[1,1] = laxis = LAxis(scene, title="workspace")
# mouseposition(laxis.scene)

mp11 = on(ax11.scene.events.mouseposition) do mp
    area = pixelarea(ax11.scene)[]
    mp = Point2f0(mp) .- minimum(area)
    r = [Point2f0(AbstractPlotting.to_world(ax11.scene, mp))]
    scatter!(ax11, mp11; markersize=3)
    # return r
end

# scatter!(ax11, mp11; markersize=3)

fig
##
geom_N = lift(x->getfield(x,:geom),ms_N)
ℋ_N = lift(x->getfield(x,:H⃗),ms_N)
ε⁻¹_N = lift(x->getfield(getfield(x,:M̂),:ε⁻¹),ms_N)
grid_N = lift(x->getfield(x,:grid),ms_N)
xs_N = lift(x,grid_N)
ys_N = lift(y,grid_N)
Ex_N = lift(E⃗x,ms_N)
Sx_N = lift(S⃗x,ms_N)

# zs_N[]
zs_n = lift(ε⁻¹_N) do ei
     [sqrt(eemm[1,1]) for eemm in inv.(ei) ]
 end

zs_E = lift()

ax_n, hm_n = heatmap(fig[1, 1],
    xs_N,
    ys_N,
    zs_n,
    colormap=:rainbow,
)

ax_S, hm_S = heatmap(fig[2, 2],
    xs_N,
    ys_N,
    zs_n,
    colormap=:rainbow,
)

##
esm1N = Node(esm1)
grN = Node(gr)
##
fig = Figure(resolution = (1200, 700), backgroundcolor = RGBf0(0.98, 0.98, 0.98) ) # , font = noto_sans)
ax_shapes = fig[2, 1] = Axis(fig, title = "ε")
ax_eps = fig[2, 2] = Axis(fig, title = "ε")
# ax_Esq = fig[3,1]
linkaxes!(ax_shapes,ax_eps)
geom_sublayout = GridLayout()
fig[2,1:2] = geom_sublayout
geom_sublayout[:h] = [ax_shapes,ax_eps]

hm_eps = epsilonheatmap!(ax_eps,grN,esm1N)

poly!(ax_shapes,geom2.shapes[1],color=:red,strokecolor=:black,strokewidth=1)
poly!(ax_shapes,geom2.shapes[2],color=:blue,strokecolor=:black,strokewidth=1)

xlims!(ax_eps,extrema(x(grN[])))
ylims!(ax_eps,extrema(y(grN[])))
ax_eps.aspect = DataAspect()
ax_shapes.aspect[] = ax_eps.aspect[] #( grN[].Δx / grN[].Δy )

cbar = fig[1,2] = Colorbar(fig, hm_eps.plots[1], label = "ε₁₁", vertical=false, height=30)

fig
##

poly!(pgn1,color=:blue,strokecolor=:black,strokewidth=1)
poly!(rect1,color=:red,strokecolor=:black,strokewidth=1)

#(geom.shapes[1])

rect1 = geom.shapes[2])
R1 = Rect2D((rect1.c-rect1.r)..., 2*rect1.r...)
poly!(ax2,R1,color=:orange,strokecolor=:black,strokewidth=1)

poly(geom2.shapes[1])

# ax_eps.autolimitaspect = ( grN[].Δx / grN[].Δy )

# esm2 = copy(esm1)

esm1N[] = copy(esm2)
esm1N[] = [SMatrix{3,3}(rand(3,3)) for i=1:128,j=1:128]

fig
##
@recipe(MSplot,ms) do scene
    Attributes(
        colormap_n = :lightrainbow,
        colormap_E = :RdBu_5,
        colormap_I = :inferno,
    )
end

function AbstractPlotting.plot!(msp::MSplot{<:Tuple{<:ModeSolver}})
    # Nodes
    ms = msp[:ms]
    # geom = lift(x->getfield(x,:geom),ms)
    # ℋ = lift(x->getfield(x,:H⃗),ms)
    # ε⁻¹ = getfield(getfield(ms,:M̂),:ε⁻¹)
    # grid = lift(x->getfield(x,:grid),ms)
    # xs = lift(x,gr)
    # ys = lift(y,gr)

    # update methods
    # function update_grid(gr,ε)
    #
    #     xs[] = lift(x,gr)
    #     ys[] = lift(y,gr)[]
    #     # zs = lift(flat,hm[:ε])
    #     zs_n[] = [sqrt(eemm[1,1]) for eemm in inv.(ε⁻¹[]) ]
    # end

    # function update_plots(gr,ε⁻¹)
    #     zs_n[]
    #     # xs[] = lift(x,gr)[]
    #     # ys[] = lift(y,gr)[]
    #     # zs = lift(flat,hm[:ε])
    #     zs_n[] = [sqrt(eemm[1,1]) for eemm in inv.(ε⁻¹[]) ]
    # end
    #
    # AbstractPlotting.Observables.onany(update_plot, geom, ε⁻¹, ℋ)

    ps = lift(ms) do x
        getindex.(convert_arguments.(getfield(x,:geom)),1)
    end

    S = lift(ms) do x
        S⃗z(x)
    end

    E = lift(ms) do x
        real(Ey_norm(x))
    end

    n = lift(ms) do x
        ei = getfield(getfield(x,:M̂),:ε⁻¹)
        sqrt.(view(flat(inv.(ei)),1,1,:,:))
    end

    # fig = Figure(resolution = (1200, 700))
    ax11 = msp[1, 1] = Axis(msp)
    ax12 = msp[1, 2] = Axis(msp)
    ax21 = msp[2, 1] = Axis(msp)
    ax22 = msp[2, 2] = Axis(msp)
    linkaxes!(ax11,ax12,ax21,ax22)

    hm_n = heatmap!(
        ax12,
        x(ms[].grid),
        y(ms[].grid),
        n,
        colormap=msp.colormap_n,
    )
    cb_n = Colorbar(msp[1, 2][1, 2], hm_n, label="n₁", width=20) # vertical=false, height=20) #width = 20)

    polycolors = [logocolors[:red],logocolors[:blue]]
    for i = 1:length(ps[])
        poly!(
            ax11,
            ps[][i],
            color = polycolors[i],
            strokecolor=:black,
            strokewidth=1,
            # axis = (backgroundcolor=logocolors[:purple]),
        )
    end
    ax11.backgroundcolor[] = logocolors[:purple]

    hm_Sz = heatmap!(
        ax22,
        x(ms[].grid),
        y(ms[].grid),
        S,
        colormap=msp.colormap_I,
    )
    cb_Sz = Colorbar(msp[2, 2][1, 2], hm_Sz, label="S⃗z", width=20) #vertical=false, height=20, label="S⃗z") #width = 20)

    Emagmax = maximum(abs.(E[]))
    hm_E = heatmap!(
        ax21,
        x(ms[].grid),
        y(ms[].grid),
        E,
        colormap=msp.colormap_E,
        colorrange = (-1,1),
    )
    cb_E = Colorbar(
        msp[2, 1][1, 2],
        hm_E,
        label="E⃗y",
        width=20,
        # limits=(-1,1),
    ) #vertical=false, height=20, label="S⃗z") #width = 20)
    msp
end
##

@recipe(EpsilonHeatmap, gr, ε) do scene
    Attributes(
        colormap = :viridis,
    )
end

function AbstractPlotting.plot!(hm::EpsilonHeatmap{<:Tuple{<:Grid{2,<:Real},<:Matrix{<:SMatrix{3,3,<:Real,9}}}})
    gr = hm[:gr]
    ε = hm[:ε]
    xs = lift(x,gr)
    ys = lift(y,gr)
    xcs = lift(xc,gr)
    ycs = lift(yc,gr)
    xlm = extrema(xs[])
    ylm = extrema(ys[])
    zs = lift(ε) do ems
        [eemm[1,1] for eemm in ems ]
    end
    #Node{AbstractMatrix}(rand(Nx,Ny))

    function update_plot(gr,ε)
        zs[]
        # xs[] = lift(x,gr)[]
        # ys[] = lift(y,gr)[]
        # zs = lift(flat,hm[:ε])
        new_zs = [eemm[1,1] for eemm in ε ] # lift(ε) do ems
        #     [eemm[1,1] for eemm in ems ]
        # end
        zs[] = new_zs
    end
    AbstractPlotting.Observables.onany(update_plot, gr, ε)
    update_plot(gr[], ε[])
    heatmap!(hm,
        xs,
        ys,
        zs,
        # colormap=sc.colormap,
    )
    # for xx in xcs[]
    #     lines!(hm,[xx,xx],[ylm...],color= :white)
    # end
    # for yy in ycs[]
    #     lines!(hm,[xlm...],[yy,yy],color= :white)
    # end
    hm
end

@recipe(GeomPoly,geom) do scene
    Attributes(
        colormap = :viridis,
    )
end

function AbstractPlotting.plot!(p::GeomPoly{<:Tuple{<:Geometry}})
    geom = p[:geom]
    shapes = geom[].shapes
    bnds = bounds.shapes
    xmin = minimum([bb[1][1] for bb in bnds])
    ymin = minimum([bb[1][2] for bb in bnds])
    xmax = maximum([bb[2][1] for bb in bnds])
    xmax = maximum([bb[2][2] for bb in bnds])

    #Node{AbstractMatrix}(rand(Nx,Ny))

    function update_plot(gr,ε)
        zs[]
        # xs[] = lift(x,gr)[]
        # ys[] = lift(y,gr)[]
        # zs = lift(flat,hm[:ε])
        new_zs = [eemm[1,1] for eemm in ε ] # lift(ε) do ems
        #     [eemm[1,1] for eemm in ems ]
        # end
        zs[] = new_zs
    end
    AbstractPlotting.Observables.onany(update_plot, gr, ε)
    update_plot(gr[], ε[])
    heatmap!(hm,
        xs,
        ys,
        zs,
        # colormap=sc.colormap,
    )
    # for xx in xcs[]
    #     lines!(hm,[xx,xx],[ylm...],color= :white)
    # end
    # for yy in ycs[]
    #     lines!(hm,[xlm...],[yy,yy],color= :white)
    # end
    hm
end

pgn1 = geom.shapes[1]
pgn_pts = [Point2f0(pgn1.v[i,:]) for i=1:size(pgn1.v)[1]]
poly!(ax2,pgn_pts,color=:red,strokecolor=:black,strokewidth=1)

ax2.limits[] = ax1.limits[]

fig

rect1 = geom.shapes[2])
R1 = Rect2D((rect1.c-rect1.r)..., 2*rect1.r...)
poly!(ax2,R1,color=:orange,strokecolor=:black,strokewidth=1)

hm1 = myheatmap(gr,esm1)

hm2 = heatmap(x(gr),
    y(gr),
    [eemm[1,1] for eemm in esm1],
    overdraw=true,
    axis = ( xlims = extrema(x(gr)),
        ylims = extrema(y(gr)),
        aspect = DataAspect(),
        ygridcolor=:white,
        xgridcolor=:white,
        yminorgridcolor=:white,
        xminorgridcolor=:white,
        backgroundcolor=:transparent,
     )
    )

hm2.axis.backgroundcolor[] = :transparent

hm2.axis.ygridcolor=:red

xlims!(extrema(x(grN[])))
ylims!(extrema(y(grN[])))
hm2.axis.aspect = DataAspect()
hm2.axis.attributes.ygridlinestyle = :dash
hm2.axis.attributes.yminorticksvisible = true
# CairoMakie.heatmap!(ax1,esm1_11); fig
#
#
# CairoMakie.heatmap!(ax1,
# # CairoMakie.heatmap(
#     x(gr),
#     y(gr),
#     flat(esm1)[1,1,:,:],
# )
#; fig

#
# data1 = randn(50, 2) * [1 2.5; 2.5 1] .+ [10 10]
#
# line1 = lines!(ax1, 5..15, x -> x, color = :red, linewidth = 2)
# scat1 = scatter!(ax1, data1,
#     color = (:red, 0.3), markersize = 15px, marker = '■')
#
# # fig
#
# ax2, line2 = lines(fig[1, 2], 7..17, x -> -x + 26,
#     color = :blue, linewidth = 2,
#     axis = (title = "C0v!d",))
#
# # fig
#
# data2 = randn(50, 2) * [1 -2.5; -2.5 1] .+ [13 13]
#
# scat2 = scatter!(data2,
#     color = (:blue, 0.3), markersize = 15px, marker = '▲')
#
# # fig
#
# linkaxes!(ax1,ax2)
# hideydecorations!(ax2, grid = false)
#
# ax1.xlabel = "time (minutes)"
# ax2.xlabel = "time (years)"
# ax1.ylabel = "\$\$\$\$"
#
# leg = fig[1, end+1] = Legend(fig,
#     [line1, scat1, line2, scat2],
#     ["f(x) = x", "Data", "f(x) = -x + 26", "Data"])
#
# fig

##
using Plots: plot, plot!, heatmap, @layout, cgrad, grid, heatmap!
using LaTeXStrings
cmap_n=cgrad(:viridis)
cmap_e=cgrad(:plasma)

# p0 = [
#     1.7,                #   top ridge width         `w_top`         [μm]
#     0.7,                #   ridge thickness         `t_core`        [μm]
#     π / 14.0,           #   ridge sidewall angle    `θ`             [radian]
# ]
# # p_lower = [0.4, 0.3, 0.]
# # p_upper = [2., 1.8, π/4.]
# p_opt0 = [1.665691811699148, 0.36154202879652847, 0.2010932097703251]
# #
# # optimized using ωs = collect(0.625:0.025:0.7) => λs0 = [1.6, 1.538, 1.481, 1.428]
# # plot with higher density of frequency points,
# ωs = collect(0.53:0.02:0.8); λs = 1. ./ ωs
# x = ms.M̂.x; y = ms.M̂.y;

p = [
                   1.7,                #   top ridge width         `w_top`         [μm]
                   0.7,                #   ridge thickness         `t_core`        [μm]
                   π / 14.0,           #   ridge sidewall angle    `θ`             [radian]
                   2.4,                #   core index              `n_core`        [1]
                   1.4,                #   substrate index         `n_subs`        [1]
                   0.5,                #   vacuum gap at boundaries `edge_gap`     [μm]
               ];
Δx,Δy,Δz,Nx,Ny,Nz = 6.0, 4.0, 1.0, 128, 128, 1
rwg(p) = ridge_wg(p[1],p[2],p[3],p[6],p[4],p[5],Δx,Δy)
ms = ModeSolver(1.45, rwg(p), Δx, Δy, Δz, Nx, Ny, Nz);
mg = MaxwellGrid(Δx,Δy,Nx,Ny)
ω = 1/1.55

##

ng = solve_n(ms,ω,rwg(p))
k,H = solve_k(ms,ω,rwg(p))
mn = vcat(reshape(ms.M̂.m,(1,3,Nx,Ny,Nz)),reshape(ms.M̂.n,(1,3,ms.M̂.Nx,ms.M̂.Ny,ms.M̂.Nz)))
e = ε⁻¹_dot( fft( kx_tc(reshape(H,(2,ms.M̂.Nx,ms.M̂.Ny,ms.M̂.Nz)),mn,ms.M̂.mag), (2:4) ), ms.M̂.ε⁻¹ )
enorm = e ./ e[argmax(abs2.(e))]
x = ms.M̂.x; y = ms.M̂.y;

ε⁻¹ = ms.M̂.ε⁻¹
ε = [ inv(ε⁻¹[:,:,ix,iy,1]) for ix=1:Nx, iy=1:Ny ]
n₁ = [ √ε[ix,iy,1][1,1] for ix=1:Nx, iy=1:Ny ]
hm_e = heatmap(
        x,
        y,
        real(enorm[1,:,:,1])',
        c=cmap_e, #cgrad(:cherry),
        aspect_ratio=:equal,
        legend=false,
        colorbar = true,
        clim = (0,1),
        xlabel = "x (μm)",
        ylabel = "y (μm)",
        colorbar_title = "|Eₓ|²",
    )
hm_n = heatmap(
    x,
    y,
    transpose(n₁),
    c=cmap_n, #cgrad(:cherry),
    aspect_ratio=:equal,
    legend=false,
    colorbar = true,
    clim = (1,2.5),
    xlabel = "x (μm)",
    ylabel = "y (μm)",
    colorbar_title = "nₓ",
    title = "optimized params",
)

##
nng = solve_n(ms,ωs,rwg(p))
k,H = solve_k(ms,ωs[14],rwg(p_opt0))
mn = vcat(reshape(ms.M̂.m,(1,3,Nx,Ny,Nz)),reshape(ms.M̂.n,(1,3,ms.M̂.Nx,ms.M̂.Ny,ms.M̂.Nz)))
e = ε⁻¹_dot( fft( kx_tc(reshape(H,(2,ms.M̂.Nx,ms.M̂.Ny,ms.M̂.Nz)),mn,ms.M̂.mag), (2:4) ), ms.M̂.ε⁻¹ )
enorm = e ./ e[argmax(abs2.(e))]
x = ms.M̂.x; y = ms.M̂.y;

ε⁻¹ = ms.M̂.ε⁻¹
ε = [ inv(ε⁻¹[:,:,ix,iy,1]) for ix=1:Nx, iy=1:Ny ]
n₁ = [ √ε[ix,iy,1][1,1] for ix=1:Nx, iy=1:Ny ]

nng0 = solve_n(ms,ωs,rwg(p0))
k0,H0 = solve_k(ms,ωs[14],rwg(p0))
mn0 = vcat(reshape(ms.M̂.m,(1,3,Nx,Ny,Nz)),reshape(ms.M̂.n,(1,3,ms.M̂.Nx,ms.M̂.Ny,ms.M̂.Nz)))
e0 = ε⁻¹_dot( fft( kx_tc(reshape(H0,(2,ms.M̂.Nx,ms.M̂.Ny,ms.M̂.Nz)),mn0,ms.M̂.mag), (2:4) ), ms.M̂.ε⁻¹ )
enorm0 = e0 ./ e0[argmax(abs2.(e0))]

ε⁻¹0 = ms.M̂.ε⁻¹
ε0 = [ inv(ε⁻¹0[:,:,ix,iy,1]) for ix=1:Nx, iy=1:Ny ]
n₁0 = [ √ε0[ix,iy,1][1,1] for ix=1:Nx, iy=1:Ny ]

##
pyplot()
ylim_ng = (2.485,2.545)

hm_e = heatmap(
        x,
        y,
        real(enorm[1,:,:,1])',
        c=cmap_e, #cgrad(:cherry),
        aspect_ratio=:equal,
        legend=false,
        colorbar = true,
        clim = (0,1),
        xlabel = "x (μm)",
        ylabel = "y (μm)",
        colorbar_title = "|Eₓ|²",
    )
hm_n = heatmap(
    x,
    y,
    transpose(n₁),
    c=cmap_n, #cgrad(:cherry),
    aspect_ratio=:equal,
    legend=false,
    colorbar = true,
    clim = (1,2.5),
    xlabel = "x (μm)",
    ylabel = "y (μm)",
    colorbar_title = "nₓ",
    title = "optimized params",
)
hm_n0 = heatmap(
        x,
        y,
        transpose(n₁0),
        c=cmap_n, #cgrad(:cherry),
        aspect_ratio=:equal,
        legend=false,
        colorbar = true,
        clim = (1,2.5),
        xlabel = "x (μm)",
        ylabel = "y (μm)",
        colorbar_title = "nₓ",
        title = "intial params",
    )
hm_e0 = heatmap(
    x,
    y,
    real(enorm0[1,:,:,1])',
    c=cmap_e, #cgrad(:cherry),
    aspect_ratio=:equal,
    legend=false,
    colorbar = true,
    clim = (0,1),
    xlabel = "x (μm)",
    ylabel = "y (μm)",
    colorbar_title = "|Eₓ|²",

)

plt_n = plot(
    λs,nng0[1],
    xlabel="λ (μm)",
    ylabel="effective index n",
    label="init. params",
    legend=:bottomleft,
    m=:dot,
    msize=2,
    )
plot!(plt_n,λs,
    nng[1],
    label="opt. params",
    m=:dot,
    msize=2,
    )
plt_ng = plot(λs,
    nng0[2],
    xlabel="λ (μm)",
    ylabel="group index ng",
    legend=false,
    m=:dot,
    msize=2,
    ylim=ylim_ng,
    )
plot!(plt_ng,
    λs,
    nng[2],
    m=:dot,
    msize=2,
    fillalpha=0.2,
    )
annot_str1 = latexstring("\$ \\mathcal{L}(\\vec{p}) = \$")
annot_str2 = latexstring("\$ \\langle n_g^2 \\rangle - \\langle n_g \\rangle^2 \$")
annot_str = "minimize\n" * annot_str1 * "\n" * annot_str2 * "\n here"
plot!(plt_ng,
    [1.428,1.6],
    [2.52,2.52],
    linecolor=nothing,
    fill_between=(ylim_ng[1],ylim_ng[2]),
    fillcolor=:blue,
    fillalpha=0.2,
    annotations = (1.51,2.53, annot_str)
    )

l = @layout [   a   b
                c   d
                e   f   ]
plot(plt_n,
        plt_ng,
        hm_n0,
        hm_n,
        hm_e0,
        hm_e,
        layout=l,
        size= (900,700),
        )


text:
function foo()
    file = open("data/1.txt", "r")
    freq = 0
    x = Set{Int64}(freq)
    adjs = parse.(Int64, readlines(file))
    while true
        for a in adjs
            freq += a
            if freq in x
                return freq
            else
                push!(x, freq)
            end
        end
    end
end

println(foo())


text:
"""
    d_func()

Lorem ipsum dolor sit amet.
"""
function d_func()
  1
end

"""
    d_func(cur_obj::MType)

Lorem ipsum dolor sit amet.
"""
function d_func(cur_obj::MType)
  2
end

"""
    d_func(cur_obj::SType)

Lorem ipsum dolor sit amet.
"""
function d_func(cur_obj::SType)
  3
end