text:
mutable struct FEFSJacobians{dim,T,MM}
    detJ::T
    Jinv::Tensors.Tensor{2,dim,T,MM}
end

struct FEMFunctionSpace{dim,T<:Real,FE<:FiniteElement,mtype <: AbstractPolytopalMesh,MM} <: AbstractFEMFunctionSpace{dim,T,FE}
    N::Matrix{T}
    dNdx::Matrix{Vec{dim,T}}
    dNdξ::Matrix{Vec{dim,T}}
    M::Matrix{T}
    dMdξ::Matrix{Vec{dim,T}}
    qr_weights::Vector{T}
    fe::FE
    mesh::mtype
    jacobians::FEFSJacobians{dim,T,MM}
    components::Int
end

getelement(fs::FEMFunctionSpace) = fs.fe

function Base.show(io::IO, fs::FEMFunctionSpace)
    if fs.components == 1
        println(io, "Discrete scalar function space")
    else
        println(io, "Discrete vectorial function space (components = $(fs.components))")
    end
    println(io, "Finite Element: \n", fs.fe)
    println(io, "Domain mesh: \n", fs.mesh)
    println(io, "Number of local dofs: ", getnlocaldofs(fs))
  end
"""
function FEMFunctionSpace(mesh::AbstractPolytopalMesh, felem::FiniteElement{dim,shape,order,gorder}, components::Int;
    quad_degree::Int = order+1) where {dim, shape, order,gorder}
    Build FEM finite element of type `felem` FunctionSpace with `components` and domain `mesh`
"""
function FEMFunctionSpace(mesh::AbstractPolytopalMesh, felem::FiniteElement{dim,shape,order,gorder}, components::Int;
    quad_degree::Int = order+1) where {dim, shape, order,gorder}
    quad_rule = QuadratureRule{shape}(DefaultQuad(), quad_degree)
    _scalar_fs(Float64, mesh, quad_rule, felem, components)
end

function _scalar_fs(::Type{T}, mesh::AbstractPolytopalMesh{dim,T}, quad_rule::QuadratureRule{shape},
    felem::FiniteElement{dim,shape,order,1}, components::Int) where {dim,T,shape<:Shape,order}
    n_qpoints = length(getpoints(quad_rule))
    n_cells = getncells(mesh)

    # Function interpolation
    n_func_basefuncs = getnbasefunctions(felem)
    N    = fill(zero(T)          * T(NaN), n_func_basefuncs, n_qpoints)
    dNdx = fill(zero(Vec{dim,T}) * T(NaN), n_func_basefuncs, n_qpoints)
    dNdξ = fill(zero(Vec{dim,T}) * T(NaN), n_func_basefuncs, n_qpoints)

    # Geometry interpolation
    n_geom_basefuncs = getngeombasefunctions(felem)
    M    = fill(zero(T)          * T(NaN), n_geom_basefuncs, n_qpoints)
    dMdξ = fill(zero(Vec{dim,T}) * T(NaN), n_geom_basefuncs, n_qpoints)

    for (qp, ξ) in enumerate(quad_rule.points)
        for i in 1:n_func_basefuncs
            dNdξ[i, qp], N[i, qp]  = gradient(ξ -> value(felem, i, ξ), ξ, :all)
        end
        for i in 1:n_geom_basefuncs
            dMdξ[i, qp], M[i, qp] = gradient(ξ -> geom_value(felem, i, ξ), ξ, :all)
        end
    end

    detJ = T(NaN)
    Jinv = zero(Tensors.Tensor{2,dim,T,2*dim})

    MM = Tensors.n_components(Tensors.get_base(typeof(Jinv)))

    jacobians = FEFSJacobians{dim,T,MM}(detJ,Jinv)
    FEMFunctionSpace{dim,T,typeof(felem),typeof(mesh),MM}(N, dNdx, dNdξ,
    M, dMdξ, getweights(quad_rule), felem, mesh,jacobians,components)
end

function reinit!(fs::FEMFunctionSpace{dim}, x::AbstractVector{Vec{dim,T}}) where {dim,T}
    n_geom_basefuncs = getngeombasefunctions(fs.fe)
    n_func_basefuncs = getnbasefunctions(fs.fe)
    @assert length(x) == n_geom_basefuncs
    fecv_J = zero(Tensors.Tensor{2,dim})
    for j in 1:n_geom_basefuncs
        fecv_J += x[j] ⊗ fs.dMdξ[j, 1]
    end
    detJ = det(fecv_J)
    detJ > 0.0 || throw(ArgumentError("det(J) is not positive: det(J) = $(detJ)"))
    fs.jacobians.detJ = detJ
    fs.jacobians.Jinv = inv(fecv_J)

    for j in 1:n_func_basefuncs, i in 1:length(fs.qr_weights)
        fs.dNdx[j, i] = fs.dNdξ[j, i] ⋅ fs.jacobians.Jinv
    end
end

########### Data Functions
getngeobasefunctions(fs::AbstractFEMFunctionSpace) = size(fs.M, 1)
getnquadpoints(fs::AbstractFEMFunctionSpace) = length(fs.qr_weights)
getnbasefunctions(fs::AbstractFEMFunctionSpace) = size(fs.N,1)*fs.components
getdetJdV(fs::FEMFunctionSpace, q_point::Int) = fs.jacobians.detJ*fs.qr_weights[q_point]
geometric_value(fs::AbstractFEMFunctionSpace, q_point::Int, base_func::Int) = fs.M[base_func, q_point]
getdim(::AbstractFEMFunctionSpace{dim}) where {dim} = dim
reference_coordinate(fs::AbstractFEMFunctionSpace{dim,T},cell::Int, mesh, x::Vec{dim,T}) where {dim,T} = fs.jacpbians.Jinv⋅(x-mesh.nodes[mesh.cells[cell].nodes[1]].x)
getfiniteelement(fs::AbstractFEMFunctionSpace) = fs.fe
getnlocaldofs(fs::AbstractFEMFunctionSpace) = getnbasefunctions(fs)
getmesh(fs::AbstractFEMFunctionSpace) = fs.mesh
getncomponents(fs::AbstractFEMFunctionSpace) = fs.components

function shape_value(fs::AbstractFEMFunctionSpace{dim,T}, q_point::Int, base_func::Int) where {dim,T}
    if fs.components == 1
        fs.N[base_func, q_point]
    else
        N_comp = zeros(T, fs.components)
        n = size(fs.N,1)
        N_comp[div(base_func,n+1)+1] = fs.N[mod1(base_func,n),q_point]
        return N_comp
    end
end

function shape_gradient(fs::AbstractFEMFunctionSpace{dim,T}, q_point::Int, base_func::Int)  where {dim,T}
    if fs.components == 1
        fs.dNdx[base_func, q_point]
    else
        dN_comp = zeros(T, fs.components, fs.components)
        n = size(fs.N,1)
        dN_comp[div(base_func,n+1)+1, :] = fs.dNdξ[mod1(base_func,n), q_point]
        return Tensors.Tensor{2,dim,T}((dN_comp...,)) ⋅ fs.jacobians.Jinv
    end
end

function shape_divergence(fs::AbstractFEMFunctionSpace, q_point::Int, base_func::Int)
    if fs.components == 1
        sum(fs.dNdx[base_func, q_point])
    else
        tr(shape_gradient(fs, q_point, base_func))
    end
end

"""
function function_value(f::Function, fs::FEMFunctionSpace{dim,T}, cell::Int,q_point::Int) where {dim,T}
Evaluate function f(x): x ∈ ℝⁿ ↦ ℝ with x given by cell number `cell` and quadrature point with index `qpoint`
    using fs FunctionSpace data
"""
function function_value(f::Function, fs::FEMFunctionSpace{dim,T}, cell::Int,q_point::Int) where {dim,T}
    coords = getverticescoords(getmesh(fs),cell)
    return f(_spatial_coordinate(fs, q_point, coords))
end

"""
function _spatial_coordinate(fs::ScalarFunctionSpace{dim}, q_point::Int, x::AbstractVector{Vec{dim,T}})
Map coordinates of quadrature point `q_point` of Function Space `fs`
into domain simplex with vertices `x`
"""
function _spatial_coordinate(fs::FEMFunctionSpace{dim}, q_point::Int, x::AbstractVector{Vec{dim,T}}) where {dim,T}
    n_base_funcs = getngeobasefunctions(fs)
    @assert length(x) == n_base_funcs
    vec = zero(Vec{dim,T})
    @inbounds for i in 1:n_base_funcs
        vec += geometric_value(fs, q_point, i) * x[i]
    end
    return vec
end

max_stars_repo_name: Paulms/jFEMToools

max_issues_repo_name: Paulms/jFEMToools

max_forks_repo_name: Paulms/jFEMToools


text:
@deprecate read(t::AITask, precision::DataType, num_samples_per_chan::Integer = -1) read(t::AITask, num_samples_per_chan, precision)
@deprecate read(t::DITask, precision::DataType, num_samples_per_chan::Integer = -1) read(t::DITask, num_samples_per_chan, precision)
@deprecate analog_input_ranges(device::String)  analog_voltage_input_ranges(device::String)
@deprecate analog_output_ranges(device::String) analog_voltage_output_ranges(device::String)

# see also `analog_input` and `analog_input_configs`

ext: jl


text:
using Statistics
using Logging
using Printf
using SpecialFunctions
using DataDeps
using JLD2

using Oceananigans
using Oceananigans.Units

using Dates: now, Second, format
using Oceananigans.Models.HydrostaticFreeSurfaceModels: VerticalVorticityField
using Oceananigans.Diagnostics: accurate_cell_advection_timescale
using Oceananigans.CubedSpheres: CubedSphereFaces, inject_cubed_sphere_exchange_boundary_conditions

Logging.global_logger(OceananigansLogger())

#####
##### Progress monitor
#####

mutable struct Progress
    interval_start_time :: Float64
end

function (p::Progress)(sim)
    wall_time = (time_ns() - p.interval_start_time) * 1e-9
    progress = sim.model.clock.time / sim.stop_time
    ETA = (1 - progress) / progress * sim.run_time
    ETA_datetime = now() + Second(round(Int, ETA))

    @info @sprintf("[%06.2f%%] Time: %s, iteration: %d, max(|u⃗|): (%.2e, %.2e) m/s, extrema(η): (min=%.2e, max=%.2e), CFL: %.2e",
                   100 * progress,
                   prettytime(sim.model.clock.time),
                   sim.model.clock.iteration,
                   maximum(abs, sim.model.velocities.u),
                   maximum(abs, sim.model.velocities.v),
                   minimum(sim.model.free_surface.η),
                   maximum(sim.model.free_surface.η),
                   sim.parameters.cfl(sim.model))

    @info @sprintf("           ETA: %s (%s), Δ(wall time): %s / iteration",
                   format(ETA_datetime, "yyyy-mm-dd HH:MM:SS"),
                   prettytime(ETA),
                   prettytime(wall_time / sim.iteration_interval))

    p.interval_start_time = time_ns()

    return nothing
end

#####
##### Loading grids
#####

ENV["DATADEPS_ALWAYS_ACCEPT"] = "true"

Logging.global_logger(OceananigansLogger())

dd32 = DataDep("cubed_sphere_32_grid",
    "Conformal cubed sphere grid with 32×32 grid points on each face",
    "https://github.com/CliMA/OceananigansArtifacts.jl/raw/main/cubed_sphere_grids/cubed_sphere_32_grid.jld2",
    "b1dafe4f9142c59a2166458a2def743cd45b20a4ed3a1ae84ad3a530e1eff538" # sha256sum
)

dd96 = DataDep("cubed_sphere_96_grid",
    "Conformal cubed sphere grid with 96×96 grid points on each face",
    "https://github.com/CliMA/OceananigansArtifacts.jl/raw/main/cubed_sphere_grids/cs96/cubed_sphere_96_grid.jld2"
)

DataDeps.register(dd32)
DataDeps.register(dd96)

#####
##### Utils
#####

function diagnose_velocities_from_streamfunction(ψ, grid)
    ψᶠᶠᶜ = Field(Face, Face,   Center, CPU(), grid)
    uᶠᶜᶜ = Field(Face, Center, Center, CPU(), grid)
    vᶜᶠᶜ = Field(Center, Face, Center, CPU(), grid)

    for (f, grid_face) in enumerate(grid.faces)
        Nx, Ny, Nz = size(grid_face)

        ψᶠᶠᶜ_face = ψᶠᶠᶜ.data.faces[f]
        uᶠᶜᶜ_face = uᶠᶜᶜ.data.faces[f]
        vᶜᶠᶜ_face = vᶜᶠᶜ.data.faces[f]

        for i in 1:Nx+1, j in 1:Ny+1
            ψᶠᶠᶜ_face[i, j, 1] = ψ(grid_face.λᶠᶠᵃ[i, j], grid_face.φᶠᶠᵃ[i, j])
        end

        for i in 1:Nx+1, j in 1:Ny
            uᶠᶜᶜ_face[i, j, 1] = (ψᶠᶠᶜ_face[i, j, 1] - ψᶠᶠᶜ_face[i, j+1, 1]) / grid.faces[f].Δyᶠᶜᵃ[i, j]
        end

        for i in 1:Nx, j in 1:Ny+1
            vᶜᶠᶜ_face[i, j, 1] = (ψᶠᶠᶜ_face[i+1, j, 1] - ψᶠᶠᶜ_face[i, j, 1]) / grid.faces[f].Δxᶜᶠᵃ[i, j]
        end
    end

    return uᶠᶜᶜ, vᶜᶠᶜ, ψᶠᶠᶜ
end

function cubed_sphere_surface_momentum_flux_bcs(τx, τy, grid; ε=0.01)

    for f in 1:length(grid.faces)
        @. τx.data.faces[f] += ε * randn()
        @. τy.data.faces[f] += ε * randn()
    end

    Nx, Ny, Nz = size(grid.faces[1])

    u_bcs_faces = Tuple(
        inject_cubed_sphere_exchange_boundary_conditions(
            FieldBoundaryConditions(face_grid, (Face, Center, Center); top=FluxBoundaryCondition(τx.data.faces[face_number][1:Nx+1, 1:Ny, 1])),
            face_number,
            grid.face_connectivity
        )
        for (face_number, face_grid) in enumerate(grid.faces)
    )

    v_bcs_faces = Tuple(
        inject_cubed_sphere_exchange_boundary_conditions(
            FieldBoundaryConditions(face_grid, (Center, Face, Center); top=FluxBoundaryCondition(τy.data.faces[face_number][1:Nx, 1:Ny+1, 1])),
            face_number,
            grid.face_connectivity
        )
        for (face_number, face_grid) in enumerate(grid.faces)
    )

    u_bcs = CubedSphereFaces{typeof(u_bcs_faces[1]), typeof(u_bcs_faces)}(u_bcs_faces)
    v_bcs = CubedSphereFaces{typeof(v_bcs_faces[1]), typeof(v_bcs_faces)}(v_bcs_faces)

    return u_bcs, v_bcs
end

#####
##### Script starts here
#####

function cubed_sphere_eddying_aquaplanet(grid_filepath)

    ## Grid setup

    H = 100meters
    grid = ConformalCubedSphereGrid(grid_filepath, Nz=1, z=(-H, 0))

    ## "Tradewind-like" zonal wind stress pattern where -π/2 ≤ φ ≤ π/2
    τx(φ) = 0.1 * exp(-10 * (φ - π/3)^2) + 0.1 * exp(-10 * (φ + π/3)^2) - 0.2 * exp(-8φ^2) + 0.19 * exp(-20φ^2)

    coriolis = HydrostaticSphericalCoriolis(scheme=VectorInvariantEnstrophyConserving())
    R = grid.faces[1].radius
    Ω = coriolis.rotation_rate

    # Streamfunction integrated using Mathematica
    # Scale by Ω^2 * R^2 * sin(φ) to get units of momentum flux (m²/s²) on the sphere.
    # Minus sign because a negative flux drives currents in the positive x-direction.
    # Hmmm, I think I messed up the magnitude but the resulting wind stress patterns should do the trick for now. It's a little too strong (stronger than the 0.1 N/m² we expect from τx).
    # I'm further scaling it down by a factor so the wind stress is even weaker (for stability/debugging).
    ψ̃(λ, φ) = - 1e-4 * Ω^2 * R^2 * sin(φ) * (-0.02802495608 * erf(1.054092553 * (π - 3φ)) - 0.06266570687 * erf(2.828427125φ) + 0.03765160933 * erf(4.472135955φ) + 0.02802495608 * erf(1.054092553 * (π + 3φ)))

    λ̃(λ) = (λ + 180)/ 360 * 2π
    φ̃(φ) = φ / 180 * π
    ψ(λ, φ) = ψ̃(λ̃(λ), φ̃(φ))

    u_top_flux, v_top_flux, ψ₀ = diagnose_velocities_from_streamfunction(ψ, grid)

    u_bcs, v_bcs = cubed_sphere_surface_momentum_flux_bcs(u_top_flux, v_top_flux, grid, ε=0)

    # Linear damping so the wind stress doesn't keep accelerating the fluid.
    # @inline linear_damping(λ, φ, z, t, u, μ) = - μ * u

    # μ = 1e-3
    # u_forcing = Forcing(linear_damping, parameters=μ, field_dependencies=:u)
    # v_forcing = Forcing(linear_damping, parameters=μ, field_dependencies=:v)

    # Since the continuous forcing seems to slow down the model by a factor of ~2x.
    @inline linear_damping_u(i, j, k, grid, clock, model_fields, μ) = @inbounds - μ * model_fields.u[i, j, k]
    @inline linear_damping_v(i, j, k, grid, clock, model_fields, μ) = @inbounds - μ * model_fields.v[i, j, k]

    μ = 1/year
    u_forcing = Forcing(linear_damping_u, parameters=μ, discrete_form=true)
    v_forcing = Forcing(linear_damping_v, parameters=μ, discrete_form=true)

    ## Model setup

    model = HydrostaticFreeSurfaceModel(
               architecture = CPU(),
                       grid = grid,
         momentum_advection = VectorInvariant(),
               free_surface = ExplicitFreeSurface(gravitational_acceleration=0.5),
                   coriolis = nothing,
                    closure = HorizontallyCurvilinearAnisotropicDiffusivity(νh=200),
        boundary_conditions = (u=u_bcs, v=v_bcs),
                  # forcing = (u=u_forcing, v=v_forcing),
                    tracers = nothing,
                   buoyancy = nothing
    )

    # Some random noise to get things going.
    ε(λ, φ, z) = 0.1 * randn()
    Oceananigans.set!(model, u=ε, v=ε)

    ## Simulation setup

    Δt = 10minutes

    g = model.free_surface.gravitational_acceleration
    gravity_wave_speed = √(g * H)
    min_spacing = filter(!iszero, grid.faces[1].Δyᶠᶠᵃ) |> minimum
    wave_propagation_time_scale = min_spacing / gravity_wave_speed
    gravity_wave_cfl = Δt / wave_propagation_time_scale
    @info @sprintf("Gravity wave CFL = %.4f", gravity_wave_cfl)

    deformation_radius_45°N = √(g * H) / (2Ω*sind(45))
    @info @sprintf("Deformation radius @ 45°N: %.2f km", deformation_radius_45°N / 1000)

    cfl = CFL(Δt, accurate_cell_advection_timescale)

    simulation = Simulation(model,
                        Δt = Δt,
                 stop_time = 5years,
        iteration_interval = 20,
                  progress = Progress(time_ns()),
                parameters = (; cfl)
    )

    output_fields = merge(model.velocities, (η=model.free_surface.η, ζ=VerticalVorticityField(model)))

    simulation.output_writers[:fields] =
        JLD2OutputWriter(model, output_fields,
            schedule = TimeInterval(1day),
              prefix = "cubed_sphere_eddying_aquaplanet",
               force = true)

    run!(simulation)

    return simulation
end

include("animate_on_map_projection.jl")

function run_cubed_sphere_eddying_aquaplanet()

    simulation = cubed_sphere_eddying_aquaplanet(datadep"cubed_sphere_96_grid/cubed_sphere_96_grid.jld2")

    projections = [
        ccrs.NearsidePerspective(central_longitude=0, central_latitude=30),
        ccrs.NearsidePerspective(central_longitude=180, central_latitude=-30)
    ]

    animate_eddying_aquaplanet(projections=projections)

    return simulation
end


text:
immutable ContinuousSpace <: Space{RealBasis,PiecewiseInterval,1}
    domain::PiecewiseInterval
end

Space(d::PiecewiseInterval)=ContinuousSpace(d)

isperiodic(C::ContinuousSpace)=isperiodic(domain(C))

spacescompatible(a::ContinuousSpace,b::ContinuousSpace)=domainscompatible(a,b)
conversion_rule{CD<:Tuple{Vararg{ChebyshevDirichlet{1,1}}}}(a::ContinuousSpace,b::PiecewiseSpace{CD,RealBasis})=a

function transform(S::ContinuousSpace,vals::Vector)
    n=length(vals)
    d=domain(S)
    K=numpieces(d)
    k=div(n,K)

    PT=promote_type(eltype(eltype(d)),eltype(vals))
    if k==0
        vals
    elseif isperiodic(d)
        ret=Array(PT,max(K,n-K))
        r=n-K*k

        for j=1:r
            cfs=transform(ChebyshevDirichlet{1,1}(d[j]),vals[(j-1)*(k+1)+1:j*(k+1)])
            if j==1
                ret[1]=cfs[1]-cfs[2]
                ret[2]=cfs[1]+cfs[2]
            elseif j < K
                ret[j+1]=cfs[1]+cfs[2]
            end
            ret[K+j:K:end]=cfs[3:end]
        end

        for j=r+1:K
            cfs=transform(ChebyshevDirichlet{1,1}(d[j]),vals[r*(k+1)+(j-r-1)*k+1:r*(k+1)+(j-r)*k])
            if length(cfs)==1 && j <K
                ret[j+1]=cfs[1]
            elseif j==1
                ret[1]=cfs[1]-cfs[2]
                ret[2]=cfs[1]+cfs[2]
            elseif j < K
                ret[j+1]=cfs[1]+cfs[2]
            end
            ret[K+j:K:end]=cfs[3:end]
        end

        ret
    else
        ret=Array(PT,n-K+1)
        r=n-K*k

        for j=1:r
            cfs=transform(ChebyshevDirichlet{1,1}(d[j]),vals[(j-1)*(k+1)+1:j*(k+1)])
            if j==1
                ret[1]=cfs[1]-cfs[2]
            end

            ret[j+1]=cfs[1]+cfs[2]
            ret[K+j+1:K:end]=cfs[3:end]
        end

        for j=r+1:K
            cfs=transform(ChebyshevDirichlet{1,1}(d[j]),vals[r*(k+1)+(j-r-1)*k+1:r*(k+1)+(j-r)*k])
            if j==1
                ret[1]=cfs[1]-cfs[2]
            end
            ret[j+1]=cfs[1]+cfs[2]
            ret[K+j+1:K:end]=cfs[3:end]
        end

        ret
    end
end

canonicalspace(S::ContinuousSpace)=PiecewiseSpace(map(ChebyshevDirichlet{1,1},pieces(domain(S))))

## pieces

Base.vec{T}(f::Fun{ContinuousSpace,T},j::Integer)=vec(Fun(f,canonicalspace(f)),j)
Base.vec{T}(f::Fun{ContinuousSpace,T})=vec(Fun(f,canonicalspace(space(f))))
pieces{T}(f::Fun{ContinuousSpace,T})=vec(f)

function points(f::Fun{ContinuousSpace})
    n=ncoefficients(f)
    d=domain(f)
    K=numpieces(d)

    m=isperiodic(d)?max(K,n+2K-1):n+K
    points(f.space,m)
end

## Conversion

coefficients(cfsin::Vector,A::ContinuousSpace,B::PiecewiseSpace)=defaultcoefficients(cfsin,A,B)

# We implemnt conversion between continuous space and PiecewiseSpace with Chebyshev dirichlet
Conversion{CD<:Tuple{Vararg{ChebyshevDirichlet{1,1}}},DD}(ps::PiecewiseSpace{CD,RealBasis,DD,1},cs::ContinuousSpace)=
                ConcreteConversion(ps,cs)

Conversion{CD<:Tuple{Vararg{ChebyshevDirichlet{1,1}}},DD}(cs::ContinuousSpace,ps::PiecewiseSpace{CD,RealBasis,DD,1})=
                ConcreteConversion(cs,ps)

bandinds{CD<:Tuple{Vararg{ChebyshevDirichlet{1,1}}},DD}(C::ConcreteConversion{PiecewiseSpace{CD,RealBasis,DD,1},ContinuousSpace})=-1,numpieces(domain(rangespace(C)))

function getindex{T,DD,CD<:Tuple{Vararg{ChebyshevDirichlet{1,1}}}}(C::ConcreteConversion{PiecewiseSpace{CD,RealBasis,DD,1},ContinuousSpace,T},k::Integer,j::Integer)
    d=domain(rangespace(C))
    K=numpieces(d)
    if isperiodic(d)
        if k==j==1
            one(T)
        elseif k==1 && j==K+1
            -one(T)
        elseif 2≤k≤K && (j==k-1 || j==K+k-1)
            one(T)
        elseif K<k && j==k+K
            one(T)
        else
            zero(T)
        end
    else
        if k==j==1
            one(T)
        elseif k==1 && j==K+1
            -one(T)
        elseif 2≤k≤K+1 && (j==k-1 || j==K+k-1)
            one(T)
        elseif K+1<k && j==k+K-1
            one(T)
        else
            zero(T)
        end
    end
end

bandinds{CD<:Tuple{Vararg{ChebyshevDirichlet{1,1}}},DD}(C::ConcreteConversion{ContinuousSpace,PiecewiseSpace{CD,RealBasis,DD,1}})=isperiodic(domainspace(C))?(1-2numpieces(domain(rangespace(C))),1):(-numpieces(domain(rangespace(C))),1)

function getindex{T,CD<:Tuple{Vararg{ChebyshevDirichlet{1,1}}},DD}(C::ConcreteConversion{ContinuousSpace,PiecewiseSpace{CD,RealBasis,DD,1},T},k::Integer,j::Integer)
    d=domain(domainspace(C))
    K=numpieces(d)
    if isperiodic(d)
        if k < K && (j==k || j==k+1)
            one(T)/2
        elseif k==K && (j==K || j==1)
            one(T)/2
        elseif K+1≤k≤2K && j==k-K
            -one(T)/2
        elseif K+1≤k<2K && j==k-K+1
            one(T)/2
        elseif k==2K && j==1
            one(T)/2
        elseif k>2K && j==k-K
            one(T)
        else
            zero(T)
        end
    else
        if k≤K && (j==k || j==k+1)
            one(T)/2
        elseif K+1≤k≤2K && j==k-K
            -one(T)/2
        elseif K+1≤k≤2K && j==k-K+1
            one(T)/2
        elseif k>2K && j==k-K+1
            one(T)
        else
            zero(T)
        end
    end
end


text:
function solve()
    x = parse(Int, readline())
    100 - (x % 100)
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
# nb: nothing in this file should really be in this package, but it's necessary for now for things to work

# TODO: figure out how to use rrule instead
@adjoint function Base.adjoint(A::LinearMap)
      pullback(Ā) = (Ā.lmap, )
      return A', pullback
end

# TODO: change to clunkier but more Julian subtyping
function Identity(size)
    LinearMap(x -> x, x -> x, size, size)
end

## rrules for LinearMaps

# CompositeMap's

rrule(::typeof(Base.adjoint), A::LinearMap) = A', (∂A -> (NoTangent(), ∂A.lmap))

function rrule(config::RuleConfig{>:HasReverseMode}, ::typeof(*), A::CompositeMap, x::AbstractVector)
    function safemul(A::CompositeMap, x::AbstractVector)
        n = length(A.maps)
        for i in 1:n
            x = A.maps[i] * x
        end
        x
    end
    rrule_via_ad(config, safemul, A, x)
end

# TODO: generalize to adjoint using metaprogramming, as in https://github.com/Jutho/LinearMaps.jl/blob/master/src/blockmap.jl
function rrule(config::RuleConfig{>:HasReverseMode}, ::typeof(*), At::TransposeMap{<:Any, <:CompositeMap}, x::AbstractVector)
    function safemul(At::TransposeMap{<:Any, <:CompositeMap}, x::AbstractVector)
        A = At.lmap
        n = length(A.maps)
        for i in n:-1:1
            x = A.maps[i]' * x
        end
        x
    end
    rrule_via_ad(config, safemul, A, x)
end

# TODO: investigate list comprehensions efficiency in backwards pass
# TODO: Ideally, I'd like to implement A' * ∂y rule for all LinearMap's;
# but I don't see how to avoid the undesirable default of
# having not implemented w.r.t. the map itself
function rrule(::typeof(*), A::Union{BlockMap, TransposeMap{<:Any, <:BlockMap}}, x::AbstractVector)
    pullback(∂y) = NoTangent(), NoTangent(), A' * ∂y
    # @not_implemented("Derivative of BlockMap multiplication w.r.t. the map itself is not implemented") # TODO: use this, while avoiding error when Zygote backpropagates not implemented for no reason
    A * x, pullback
end

# fixes bug in LinearMaps.jl: _unsafe_mul! (used, for example, in CompositeMap's) does not fall back to mul! for transpose maps
_unsafe_mul!(y::AbstractVecOrMat, A::TransposeMap, x::AbstractVector) =
                issymmetric(A.lmap) ?
               _unsafe_mul!(y, A.lmap, x) : LinearMaps._generic_mapvec_mul!(y, A, x, 1, 0)

_unsafe_mul!(y::AbstractVecOrMat, A::AdjointMap, x::AbstractVector) =
                ishermitian(A.lmap) ?
               _unsafe_mul!(y, A.lmap, x) : LinearMaps._generic_mapvec_mul!(y, A, x, 1, 0)


text:
module VulkanExamples

using MLStyle
using Meshes
using ColorTypes
using Vulkan

function __init__()
    debug_callback_c[] = @cfunction(
        default_debug_callback,
        UInt32,
        (
            DebugUtilsMessageSeverityFlagEXT,
            DebugUtilsMessageTypeFlagEXT,
            Ptr{vk.VkDebugUtilsMessengerCallbackDataEXT},
            Ptr{Cvoid},
        )
    )
end

include("init.jl")
include("memory.jl")
include("vertex.jl")

export init,

    # memory
    find_memory_type,
    buffer_size,
    upload_data,
    download_data,

    # vertex
    PosColor,
    PosUV,
    Point4f,
    invert_y_axis

end


text:
#  -- updated drivers -- #
using Distributions, JSON, DistributedTensorConstruction, Random, Distributed
@everywhere using LambdaTAME

driver_path = "./"
all_results_path = "results/"
sparse_matrices_path = "../data/sparse_matrices/LVGNA/"
sparse_tensors_path = "../data/sparse_tensors/LVGNA/"

function save_results(file_path,data)
    open(file_path,"w") do f
        JSON.print(f,data)
    end
end

include(driver_path*"max_rank_experiments.jl")        # Figure  5.1 (a)
include(driver_path*"Increasing_Clique_experiments.jl") # Figures 5.1 (b), 5.2 (b), SM3
include(driver_path*"LVGNA_experiments.jl")           # Figures 5.1 (a) 5.3 (a) & (b), SM5
include(driver_path*"RG_Dup_noise_experiments.jl")    # Figures 5.2 (a), SM 4 (a) & (b)
include(driver_path*"RG_ER_noise_experiments.jl")     # Figures 5.2 (a), SM 4 (a) & (b)
include(driver_path*"singular_values_experiments.jl") # Figure  A.1

function run_all()
    # All functions will run in testing mode unless
    # 'testing=false' is passed in.

    #  -- max_rank_experiments.jl functions --  #
    LVGNA_max_rank_experiments()
    LVGNA_shift_experiment()
    random_graph_experiments()

    #  -- Increasing_Clique_experiments.jl functions --  #
    TAME_LRT_increasing_cliques()
    increased_clique_size_duplication_noise_experiments()

    #  -- LVGNA_experiments.jl functions --  #
    LVGNA_alignments()

    #  -- RG_Dup_noise_experiments.jl functions --  #
    RandomGraph_duplication_noise_experiments()

    #  -- RG_ER_noise_experiments.jl functions --  #
    RandomGraph_ER_noise_experiments()

    #  -- singular_values_experiments.jl functions --  #
    singular_values_experiments()
end

max_stars_repo_name: chuckcol/LambdaTAME

max_issues_repo_name: chuckcol/LambdaTAME

max_forks_repo_name: chuckcol/LambdaTAME


text:
@testset "Hessian of Lagrangian" begin

    # horizon
    T = 3

    # acrobot
    nx = 4
    nu = 1
    nw = 0
    w_dim = [nw for t = 1:T]

    function acrobot(x, u, w)
        # dimensions
        n = 4
        m = 1
        d = 0

        # link 1
        mass1 = 1.0
        inertia1 = 0.33
        length1 = 1.0
        lengthcom1 = 0.5

        # link 2
        mass2 = 1.0
        inertia2 = 0.33
        length2 = 1.0
        lengthcom2 = 0.5

        gravity = 9.81
        friction1 = 0.1
        friction2 = 0.1

        # mass matrix
        function M(x, w)
            a = (inertia1 + inertia2 + mass2 * length1 * length1
                + 2.0 * mass2 * length1 * lengthcom2 * cos(x[2]))

            b = inertia2 + mass2 * length1 * lengthcom2 * cos(x[2])

            c = inertia2

        return [a b; b c]
        end

        # dynamics bias
        function τ(x, w)
            a = (-1.0 * mass1 * gravity * lengthcom1 * sin(x[1])
                - mass2 * gravity * (length1 * sin(x[1])
                + lengthcom2 * sin(x[1] + x[2])))

            b = -1.0 * mass2 * gravity * lengthcom2 * sin(x[1] + x[2])

            return [a; b]
        end

        function C(x, w)
            a = -2.0 * mass2 * length1 * lengthcom2 * sin(x[2]) * x[4]
            b = -1.0 * mass2 * length1 * lengthcom2 * sin(x[2]) * x[4]
            c = mass2 * length1 * lengthcom2 * sin(x[2]) * x[3]
            d = 0.0

            return [a b; c d]
        end

        # input Jacobian
        function B(x, w)
            [0.0; 1.0]
        end

        # dynamics
        q = view(x, 1:2)
        v = view(x, 3:4)

        qdd = M(q, w) \ (-1.0 * C(x, w) * v
                + τ(q, w) + B(q, w) * u[1] - [friction1; friction2] .* v)

        return [x[3]; x[4]; qdd[1]; qdd[2]]
    end

    function midpoint_implicit(y, x, u, w)
        h = 0.05 # timestep
        y - (x + h * acrobot(0.5 * (x + y), u, w))
    end

    dt = Dynamics(midpoint_implicit, nx, nx, nu, nw=nw, eval_hess=true)
    dyn = [dt for t = 1:T-1]
    model = DynamicsModel(dyn, w_dim=w_dim)

    # initial state
    x1 = [0.0; 0.0; 0.0; 0.0]

    # goal state
    xT = [0.0; π; 0.0; 0.0]

    # objective
    ot = (x, u, w) -> 0.1 * dot(x[3:4], x[3:4]) + 0.1 * dot(u, u)
    oT = (x, u, w) -> 0.1 * dot(x[3:4], x[3:4])
    objt = Cost(ot, nx, nu, nw, [t for t = 1:T-1], eval_hess=true)
    objT = Cost(oT, nx, 0, nw, [T], eval_hess=true)
    obj = [objt, objT]

    # constraints
    x_init = Bound(nx, nu, [1], xl=x1, xu=x1)
    x_goal = Bound(nx, 0, [T], xl=xT, xu=xT)
    ct = (x, u, w) -> [-5.0 * ones(nu) - cos.(u) .* sum(x.^2); cos.(x) .* tan.(u) - 5.0 * ones(nx)]
    cT = (x, u, w) -> sin.(x.^3.0)
    cont = StageConstraint(ct, nx, nu, nw, [t for t = 1:T-1], :inequality, eval_hess=true)
    conT = StageConstraint(cT, nx, 0, nw, [T], :equality, eval_hess=true)
    cons = ConstraintSet([x_init, x_goal], [cont, conT])

    # problem
    trajopt = TrajectoryOptimizationProblem(obj, model, cons)
    s = Solver(trajopt, eval_hess=true)

    # Lagrangian
    function lagrangian(z)
        x1 = z[1:nx]
        u1 = z[nx .+ (1:nu)]
        x2 = z[nx + nu .+ (1:nx)]
        u2 = z[nx + nu + nx .+ (1:nu)]
        x3 = z[nx + nu + nx + nu .+ (1:nx)]
        λ1_dyn = z[nx + nu + nx + nu + nx .+ (1:nx)]
        λ2_dyn = z[nx + nu + nx + nu + nx + nx .+ (1:nx)]

        λ1_stage = z[nx + nu + nx + nu + nx + nx + nx .+ (1:(nu + nx))]
        λ2_stage = z[nx + nu + nx + nu + nx + nx + nx + nu + nx .+ (1:(nu + nx))]
        λ3_stage = z[nx + nu + nx + nu + nx + nx + nx + nu + nx + nu + nx .+ (1:nx)]

        L = 0.0
        L += ot(x1, u1, zeros(nw))
        L += ot(x2, u2, zeros(nw))
        L += oT(x3, zeros(0), zeros(nw))
        L += dot(λ1_dyn, midpoint_implicit(x2, x1, u1, zeros(nw)))
        L += dot(λ2_dyn, midpoint_implicit(x3, x2, u2, zeros(nw)))
        L += dot(λ1_stage, ct(x1, u1, zeros(nw)))
        L += dot(λ2_stage, ct(x2, u2, zeros(nw)))
        L += dot(λ3_stage, cT(x3, zeros(0), zeros(nw)))
        return L
    end

    nz = nx + nu + nx + nu + nx + nx + nx + nu + nx + nu + nx + nx
    np = nx + nu + nx + nu + nx
    nd = nx + nx + nu + nx + nu + nx + nx
    @variables z[1:nz]
    L = lagrangian(z)
    Lxx = Symbolics.hessian(L, z[1:np])
    Lxx_sp = Symbolics.sparsehessian(L, z[1:np])
    spar = [findnz(Lxx_sp)[1:2]...]
    Lxx_func = eval(Symbolics.build_function(Lxx, z)[1])
    Lxx_sp_func = eval(Symbolics.build_function(Lxx_sp.nzval, z)[1])

    z0 = rand(nz)
    nh = length(s.p.sp_hess_lag)
    h0 = zeros(nh)

    σ = 1.0
    fill!(h0, 0.0)
    DirectTrajectoryOptimization.trajectory!(s.p.trajopt.x, s.p.trajopt.u, z0[1:np],
        s.p.trajopt.model.idx.x, s.p.trajopt.model.idx.u)
    DirectTrajectoryOptimization.duals!(s.p.trajopt.λ_dyn, s.p.trajopt.λ_stage, z0[np .+ (1:nd)], s.p.idx.dyn_con, s.p.idx.stage_con)
    DirectTrajectoryOptimization.eval_obj_hess!(h0, s.p.idx.obj_hess, s.p.trajopt.obj, s.p.trajopt.x, s.p.trajopt.u, s.p.trajopt.w, σ)
    DirectTrajectoryOptimization.eval_hess_lag!(h0, s.p.idx.dyn_hess, s.p.trajopt.model.dyn, s.p.trajopt.x, s.p.trajopt.u, s.p.trajopt.w, s.p.trajopt.λ_dyn)
    DirectTrajectoryOptimization.eval_hess_lag!(h0, s.p.idx.stage_hess, s.p.trajopt.con.stage, s.p.trajopt.x, s.p.trajopt.u, s.p.trajopt.w, s.p.trajopt.λ_stage)

    sp_obj_hess = DirectTrajectoryOptimization.sparsity_hessian(obj, model.dim.x, model.dim.u)
    sp_dyn_hess = DirectTrajectoryOptimization.sparsity_hessian(dyn, model.dim.x, model.dim.u)
    sp_con_hess = DirectTrajectoryOptimization.sparsity_hessian(cons.stage, model.dim.x, model.dim.u)
    sp_hess = collect([sp_obj_hess..., sp_dyn_hess..., sp_con_hess...])
    sp_key = sort(unique(sp_hess))

    idx_obj_hess = DirectTrajectoryOptimization.hessian_indices(obj, sp_key, model.dim.x, model.dim.u)
    idx_dyn_hess = DirectTrajectoryOptimization.hessian_indices(dyn, sp_key, model.dim.x, model.dim.u)
    idx_con_hess = DirectTrajectoryOptimization.hessian_indices(cons.stage, sp_key, model.dim.x, model.dim.u)

    # indices
    @test sp_key[vcat(idx_obj_hess...)] == sp_obj_hess
    @test sp_key[vcat(idx_dyn_hess...)] == sp_dyn_hess
    @test sp_key[vcat(idx_con_hess...)] == sp_con_hess

    # Hessian
    h0_full = zeros(np, np)
    for (i, h) in enumerate(h0)
        h0_full[sp_key[i][1], sp_key[i][2]] = h
    end
    @test norm(h0_full - Lxx_func(z0)) < 1.0e-8
    @test norm(norm(h0 - Lxx_sp_func(z0))) < 1.0e-8

    # h0 = zeros(nh)
    # MOI.eval_hessian_lagrangian(s.p, h0, z0[1:np], σ, z0[np .+ (1:nd)])
    # @test norm(norm(h0 - Lxx_sp_func(z0))) < 1.0e-8

    # a = z0[1:np]
    # b = z0[np .+ (1:nd)]
    # info = @benchmark MOI.eval_hessian_lagrangian($s.p, $h0, $a, $σ, $b)
end


text:
"""
    exponents = get_exponents(variables, dmax::Int)

Compute the set of all exponents in `variables` variables, of degree up to
`dmax`.
"""
function compute_exponents(variables::Set{Variable}, dmax::Int; compute_conj=false)
    cur_order = Set{Exponent}([Exponent()])
    result = deepcopy(cur_order)
    prev_order = Set{Exponent}()
    for i=1:dmax
        prev_order = deepcopy(cur_order)
        cur_order = Set{Exponent}()
        for var in variables
            if compute_conj
                union!(cur_order, Set([product(conj(var), elt) for elt in prev_order]))
            else
                union!(cur_order, Set([product(Exponent(var), elt) for elt in prev_order]))
            end
        end
        union!(result, cur_order)
    end
    return result
end

"""
    ishomo = is_homogeneous(p, kind)

    Check wether `p`, of kind `:Real` or `:Complex` is homogeneous or not.
"""
function is_homogeneous(p::Polynomial, kind::Symbol)
    ishomo = true
    for (expo, λ) in p
        ishomo = ishomo && is_homogeneous(expo, kind)
    end
    return ishomo
end

"""
    ishomo = is_homogeneous(expo, kind)

    Check wether `expo`, of kind `:Real` or `:Complex` is homogeneous or not.
"""
function is_homogeneous(expo::Exponent, kind::Symbol)
    explsum, conjsum = get_sumdegs(expo)
    if kind == :Real
        (conjsum != 0) && error(LOGGER, "is_homogeneous(): Exponent $expo has conjugated variables for a real hierarchy.")
        return explsum % 2 == 0
    elseif kind == :Complex
        return explsum == conjsum
    else
        error(LOGGER, "is_homogeneous(expo, kind): kind should be either :Real or :Complex ($kind here).")
    end
end

"""
    explsum, conjsum = get_sumdegs(expo)

    Compute `|α|`, `|β|` the sum of the real variables exponents and conjugated variables exponents.
"""
function get_sumdegs(expo::Exponent)
    explsum = conjsum = 0

    state1 = start(expo)
    while !done(expo, state1)
        (i1, state1) = next(expo, state1)
        explsum += last(i1).explvar
        conjsum += last(i1).conjvar
    end

    return explsum, conjsum
end

"""
    cstrtype = get_cstrtype(cstr::Constraint)

    Return a cstraint type among `:eq`, `:ineqhi`, `:ineqlo`, `:ineqdouble`.
"""
function get_cstrtype(cstr::Constraint)
    if cstr.lb == cstr.ub && isfinite(cstr.ub)
        return :eq
    elseif (cstr.lb == -Inf-im*Inf) && (cstr.ub != Inf+im*Inf)
        return :ineqhi
    elseif (cstr.lb != -Inf-im*Inf) && (cstr.ub == Inf+im*Inf)
        return :ineqlo
    elseif (cstr.lb != -Inf-im*Inf) && (cstr.ub != Inf+im*Inf)
        return :ineqdouble
    else
        error(LOGGER, "get_cstrtype(): unknown constraint type.\nConstraint is $cstr")
    end
end

function compute_degree(expo::Exponent)
    expldeg = conjdeg = 0
    for (var, deg) in expo.expo
        expldeg = max(expldeg, deg.explvar)
        conjdeg = max(conjdeg, deg.conjvar)
    end
    return Degree(expldeg, conjdeg)
end

function compute_degree(p::Polynomial)
    expldeg = conjdeg = 0
    for (expo, λ) in p
        for (var, deg) in expo
            expldeg = max(expldeg, deg.explvar)
            conjdeg = max(conjdeg, deg.conjvar)
        end
    end
    return Degree(expldeg, conjdeg)
end

function update_degree!(expo::Exponent)
    updeg = compute_degree(expo)
    expo.degree.explvar = updeg.explvar
    expo.degree.conjvar = updeg.conjvar
end

function update_degree!(p::Polynomial)
    updeg = compute_degree(p)
    p.degree.explvar = updeg.explvar
    p.degree.conjvar = updeg.conjvar
end


text:
import NaiveNASlib
import NaiveNASlib:reset_in!, reset_out!, reset!
import InteractiveUtils:subtypes

@testset "Vertex mutation operations" begin

    expectedtype(t::Type{<:MutationOp}) = Integer
    expectedtype(t::Type{IoIndices}) = AbstractArray{<:Integer,1}

    @testset "Method contracts $subtype" for subtype in implementations(MutationOp)
        @test hasmethod(Δnin, (subtype,expectedtype(subtype)))
        @test hasmethod(Δnout, (subtype, expectedtype(subtype)))
        @test hasmethod(clone, (subtype,))
    end

    @testset "Method contracts $subtype" for subtype in implementations(MutationState)
        @test hasmethod(nin, (subtype,))
        @test hasmethod(nout, (subtype,))
        @test hasmethod(in_inds, (subtype,))
        @test hasmethod(out_inds, (subtype,))
    end

    @testset "IoSize" begin
        size = IoSize(3, 4)
        @test nin(size) == [3]
        @test nout(size) == 4

        Δnin(size, -1)
        @test nin(size) == [2]
        @test nout(size) == 4

        Δnout(size, 2)
        @test nin(size) == [2]
        @test nout(size) == 6

        @test in_inds(size) == [1:2]
        @test out_inds(size) == 1:6

        size = IoSize([2, 3, 4], 5)
        @test nin(size) == [2, 3, 4]

        Δnin(size, 1,-2, 3)
        @test nin(size) == [3, 1 ,7]

        Δnin(size, missing,4, -2)
        @test nin(size) == [3, 5 ,5]

        @test in_inds(size) == [1:3, 1:5, 1:5]

        @test issame(size, clone(size))
    end

    @testset "IoIndices" begin
        inds = IoIndices([3,4], 5)

        @test nin(inds) == [3, 4]
        @test nout(inds) == 5

        Δnin(inds, [1,2], [2,3,4])
        @test nin(inds) == [2, 3]
        @test nout(inds) == 5
        @test in_inds(inds) == [[1, 2], [2,3,4]]

        reset_in!(inds)
        @test in_inds(inds) == [[1, 2], [1,2,3]]

        Δnout(inds, [2,4])
        @test nin(inds) == [2, 3]
        @test nout(inds) == 2
        @test out_inds(inds) == [2, 4]

        reset_out!(inds)
        @test out_inds(inds) == [1, 2]

        Δnin(inds, missing, [-1, 2, -1])
        @test in_inds(inds) == [[1, 2], [-1,2,-1]]

        @test issame(inds, clone(inds))
    end

    @testset "IoChange" begin
        s = IoChange(3, 4)

        @test nin(s) == [3]
        @test nout(s) == 4
        @test in_inds(s) == [1:3]
        @test out_inds(s) == 1:4

        Δnin(s, -1)
        @test nin(s) == [2]
        @test nout(s) == 4
        @test in_inds(s) == [1:2]

        s = IoChange([3,4], 5)

        @test nin(s) == [3, 4]
        @test nout(s) == 5
        @test in_inds(s) == [1:3, 1:4]
        @test out_inds(s) == 1:5

        Δnin(s, 1, -2)
        @test nin(s) == [4, 2]
        @test nout(s) == 5
        @test in_inds(s) == [vcat(1:3, -1), 1:2]
        @test out_inds(s) == 1:5

        Δnout(s, 2)
        @test nin(s) == [4, 2]
        @test nout(s) == 7
        @test in_inds(s) == [vcat(1:3, -1), 1:2]
        @test out_inds(s) == vcat(1:5, [-1, -1])

        Δnin(s, [1, -1, 2, -1], [2, -1])
        @test nin(s) == [4, 2]
        @test nout(s) == 7
        @test in_inds(s) == [[1, -1, 2, -1], [2, -1]]
        @test out_inds(s) == vcat(1:5, [-1, -1])

        Δnout(s, [-1, 2, 3, -1, -1, 5, 7])
        @test nin(s) == [4, 2]
        @test nout(s) == 7
        @test in_inds(s) == [[1, -1, 2, -1], [2, -1]]
        @test out_inds(s) == [-1, 2, 3, -1, -1, 5, 7]

        Δnin(s, [1, 3], [-1, 2, 1])
        @test nin(s) == [2, 3]
        @test in_inds(s) == [[1, 3], [-1, 2, 1]]

        Δnout(s, [2, -1, 3, -1])
        @test nout(s) == 4
        @test out_inds(s) == [2, -1, 3, -1]

        reset_in!(s)
        @test nin_org(s) == nin(s) == [2, 3]
        @test in_inds(s) == [1:2, 1:3]
        @test nout_org(s) != nout(s)
        @test out_inds(s) == [2, -1, 3, -1]

        reset_out!(s)
        @test nout_org(s) == nout(s) == 4
        @test out_inds(s) == 1:4

        Δnin(s, 2, -1)
        @test nin(s) == [4, 2]

        Δnin(s, missing, [1, 2])
        @test nin(s) == [4, 2]

        s = IoChange([4 ,2], 4)
        Δnin(s, -1, 3)
        reset_in!(s)

        @test in_inds(s) == [[1,2,3], [1,2,3,4,5]]

        Δnout(s, 3)
        reset_out!(s)
        @test out_inds(s) == 1:7

    end
end


text:
abstract type Chronology end

struct Synchronize <: Chronology
    from_steps::Int64
    to_steps::Int64
end

struct RecedingHorizon <: Chronology end

struct Sequential <: Chronology end


text:
@testset "transition rate matrix checks" begin
    @test_throws ArgumentError TransitionRateMatrix(ones(3, 4))  # non-square
    @test_throws DomainError TransitionRateMatrix(-ones(3, 3))   # negative
end

@testset "transition matrix creation (dense)" begin
    for i in 1:100
        n = rand(5:10)
        A = abs.(randn(n, n))
        # copying constructor
        Q = TransitionRateMatrix(A)
        @test size(Q) ≡ (n, n)
        @test Base.IndexStyle(typeof(Q)) ≡ Base.IndexStyle(typeof(A))
        for i in axes(Q, 1)
            @test sum(Q[i, :]) ≈ 0 atol = n*10*eps()
        end

        # non-copying constructor
        Q2 = TransitionRateMatrix!(A)
        @test size(Q2) ≡ (n, n)
        @test Base.IndexStyle(typeof(Q2)) ≡ Base.IndexStyle(typeof(A))
        for i in axes(Q2, 1)
            @test sum(Q2[i, :]) ≈ 0 atol = n*10*eps()
        end

        π = stationary_distribution(Q)
        π2 = stationary_distribution(Q2)
        @test π == π2           # exact equivalence
        @test sum(π) ≈ 1
        @test all(π .≥ 0)
        @test Q'*π ≈ zero(π) atol = 10*eps()
    end
end

@testset "transition matrix creation (sparse)" begin
    for i in 1:100
        n = rand(5:10)
        A = abs.(sprandn(n, n, 1/n))
        # copying constructor
        Q = TransitionRateMatrix(A)
        @test size(Q) ≡ (n, n)
        @test Base.IndexStyle(typeof(Q)) ≡ Base.IndexStyle(typeof(A))
        for i in axes(Q, 1)
            @test sum(Q[i, :]) ≈ 0 atol = n*10*eps()
        end

        # non-copying constructor
        Q2 = TransitionRateMatrix!(A)
        @test size(Q2) ≡ (n, n)
        @test Base.IndexStyle(typeof(Q2)) ≡ Base.IndexStyle(typeof(A))
        for i in axes(Q2, 1)
            @test sum(Q2[i, :]) ≈ 0 atol = n*10*eps()
        end

        # FIXME tests commented out, cf
        # https://github.com/tpapp/ContinuousTimeMarkov.jl/issues/2

        # π = stationary_distribution(Q)
        # π2 = stationary_distribution(Q2)
        # @test π == π2           # exact equivalence
        # @test sum(π) ≈ 1
        # @test all(π .≥ 0)
        # @test Q'*π ≈ zero(π) atol = 10*eps()
    end
end


text:
include("features.jl")
include("optimization.jl")

function initialize()
    coord, locinfo = get_mapinfo()
    config = Config()
    feat = Features(coord, locinfo)
    bookkeep = BookKeep(coord.dim, feat.dim)

    fill_diag!(coord.pairwise, floatmax(Float64))

    return coord, locinfo, config, bookkeep, feat
end

function main(latitude=51.355468, longitude=11.100790; state=nothing, output=true)
    coord, locinfo, config, bookkeep, feat = isnothing(state) ? initialize() : state

    # precomputing pointwise distance to provided coordinates
    dist_pointwise!(bookkeep.pointwise, degtorad(longitude), degtorad(latitude), coord.longitudes, coord.latitudes)
    update!(feat, bookkeep.pointwise)
    Random.seed!(config.rng, 1234) # reset seed to make results reproductable

    run_optimization(coord, locinfo, config, bookkeep, feat, output)
end


text:
module TestUtils

using ..Memento
using Test

import Test: @test_warn, @test_throws

export @test_log, @test_nolog

const EMPTY_MATCH = nothing

# These function are just a copy of `ismatch_warn` on 0.6 or `contains_warn` on 0.7`.
# contains_msg(output, s) = contains_warn(output, s)
occursin_msg(s::AbstractString, output) = occursin(s, output)
occursin_msg(s::Regex, output) = occursin(s, output)
occursin_msg(s::Function, output) = s(output)
occursin_msg(S::Union{AbstractArray,Tuple}, output) = all(s -> occursin_msg(s, output), S)

"""
    @test_log(logger, level, msg, expr)

Test that the expression `expr` emits a record in the `logger` with the specified `level`
string and contains the `msg` string or matches the `msg` regular expression.
If `msg` is a boolean function, tests whether `msg(output)` returns true.
If `msg` is a tuple or array, checks that the output contains/matches each item in `msg`.
Returns the result of evaluating `expr`.

This will temporarily add a test handler to the `logger` which will always be removed after
executing the expression.

See also [`@test_nolog`](@ref) to check for the absence of a record.

## Example
```jldoctest
julia> using Memento, Memento.TestUtils

julia> logger = getlogger("test_log");

julia> m = "Hello World!";

julia> @test_log logger "info" "Hello" info(logger, m)

julia> @test_log logger "info" r"^Hello" info(logger, m)

julia> @test_log logger "info" ("Hello", r"World!\$") info(logger, m)  # All elements occursin
```
"""
macro test_log(logger, level, msg, expr)
    quote
        handler = TestHandler($(esc(level)), $(esc(msg)))
        handlers = copy(gethandlers($(esc(logger))))

        try
            $(esc(logger)).handlers = Dict{Any, Handler}()
            push!($(esc(logger)), handler)

            setpropagating!($(esc(logger)), false) do
                ret = $(esc(expr))
                @test handler.found[1] == $(esc(level))
                @test occursin_msg($(esc(msg)), handler.found[2])
                ret
            end
        finally
            $(esc(logger)).handlers = handlers
        end
    end
end

"""
    @test_nolog(logger, level, msg, expr)

Test that the expression `expr` does not emit a record in the `logger` with the specified
`level` string containing the `msg` string or matching the `msg` regular expression.

See also [`@test_log`](@ref) for further details.
"""
macro test_nolog(logger, level, msg, expr)
    quote
        handler = TestHandler($(esc(level)), $(esc(msg)))
        handlers = copy(gethandlers($(esc(logger))))

        try
            $(esc(logger)).handlers = Dict{Any, Handler}()
            push!($(esc(logger)), handler)

            setpropagating!($(esc(logger)), false) do
                ret = $(esc(expr))
                @test handler.found[1] == EMPTY_MATCH && handler.found[2] == EMPTY_MATCH
                ret
            end
        finally
            $(esc(logger)).handlers = handlers
        end
    end
end

"""
    @test_warn(logger, msg, expr)

Convenience macro that calls [`@test_log(logger, "warn", msg, expr)`](@ref @test_log).
"""
macro test_warn(logger, msg, expr)
    quote
        @test_log($(esc(logger)), "warn", $(esc(msg)), $(esc(expr)))
    end
end

"""
    @test_throws(logger, extype, expr)

Disables the `logger` and calls [`@test_throws extype expr`](https://docs.julialang.org/en/v1/stdlib/Test/#Test.@test_throws).
"""
macro test_throws(logger, extype, expr)
    quote
        setlevel!($(esc(logger)), "not_set") do
            setpropagating!($(esc(logger)), false) do
                @test_throws $(esc(extype)) $(esc(expr))
            end
        end
    end
end

mutable struct TestHandler{F, O} <: Handler{F, O}
    level::String
    msg
    found::Tuple
end

function TestHandler(level, msg)
    TestHandler{DefaultFormatter, IOBuffer}(
        String(level), msg, (EMPTY_MATCH, EMPTY_MATCH)
    )
end

function Base.log(handler::TestHandler, rec::Record)
    level = getlevel(rec)

    # Uncomment the lines below to debug issues with the `TestHandler`
    # println(string("Record: ", level, ", ", rec.msg))
    # println(string("Search: ", handler.level, ", ", handler.msg))
    # println(string("Found: ", level == handler.level && occursin_msg(handler.msg, String(rec.msg))))

    if level == handler.level && occursin_msg(handler.msg, String(rec.msg))
        handler.found = (level, rec.msg)
    end
end

end


text:
using CellularInfo, JSON, Plots, Plots.PlotMeasures

function plotit(basename)
    μ = open(JSON.parse, "$basename.json")
    class1 = μ[ECA_CLASS_I .+ 1]
    class2 = μ[ECA_CLASS_II .+ 1]
    class3 = μ[ECA_CLASS_III .+ 1]
    class4 = μ[ECA_CLASS_IV .+ 1]

    c12 = max(maximum(class1), maximum(class2))

    @show maximum(class1)
    @show maximum(class2)
    @show maximum(class3)
    @show maximum(class4)

    p = scatter(ECA_CLASS_I, class1, markershape=:cross, markercolor=:black, label="Class I", left_margin=1cm, bottom_margin=8mm)
    scatter!(p, ECA_CLASS_II, class2, markershape=:xcross, markercolor=:black, label="Class II")
    scatter!(p, ECA_CLASS_III, class3, markershape=:rect, markercolor=:green, label="Class III")
    scatter!(p, ECA_CLASS_IV, class4, markershape=:rect, markercolor=:orange, label="Class IV")
    hline!(p, [c12], linestyle=:dash, label="Class I & II Limit")
    xlabel!(p, "Rule")
    ylabel!(p, "Mean TE")
    title!(p, "Mean TE vs Rule by Wolfram Class")
    savefig(p, "$basename.png")
end

plotit("onepoint")
plotit("rand")

max_stars_repo_name: ELIFE-ASU/dglmoore

max_issues_repo_name: ELIFE-ASU/dglmoore

max_forks_repo_name: ELIFE-ASU/dglmoore


text:
@testset "Gevrey functions" begin

    t₀ = 0.0
    Tf = 10.0 # Final time
    Δt = 0.1  # Time discretization
    tspan = t₀ : Δt : Tf # Time span

    w = 1.5 # Steepness

    try
        y = pulsebell(t₀, Tf; ω = w)
        @test true
    catch
        @test false
    end

    try
        y = smoothjump(t₀, Tf; ω = w)
        @test true
    catch
        @test false
    end

    try
        y = smoothjump.(tspan, Tf; ω = w)
        @test true
    catch
        @test false
    end

end


text:
multidims = zeros(Int8, 2, 3)
tupdims = zeros(Int8, (2, 3))
defaulttype = zeros((2, 3))

println(multidims)
println(tupdims)
println(defaulttype)

max_stars_repo_name: JulialangOrgCN/howtojulia

max_issues_repo_name: JulialangOrgCN/howtojulia

max_forks_repo_name: JulialangOrgCN/howtojulia


text:
using Test
using QuantumLattices.Mathematics.Combinatorics

@testset "Combinations" begin
    @test Combinations{0}(1:3)|>collect==[()]
    @test Combinations{1}(1:3)|>collect==[(1,),(2,),(3,)]
    @test Combinations{2}(1:3)|>collect==[(1,2),(1,3),(2,3)]
    @test Combinations{3}(1:3)|>collect==[(1,2,3)]
end

@testset "DulCombinations" begin
    @test DulCombinations{0}(1:3)|>collect==[()]
    @test DulCombinations{1}(1:3)|>collect==[(1,),(2,),(3,)]
    @test DulCombinations{2}(1:3)|>collect==[(1,1),(1,2),(1,3),(2,2),(2,3),(3,3)]
    @test DulCombinations{3}(1:3)|>collect==[(1,1,1),(1,1,2),(1,1,3),(1,2,2),(1,2,3),(1,3,3),(2,2,2),(2,2,3),(2,3,3),(3,3,3)]
end

@testset "Permutations" begin
    @test Permutations{0}(1:3)|>collect==[()]
    @test Permutations{1}(1:3)|>collect==[(1,),(2,),(3,)]
    @test Permutations{2}(1:3)|>collect==[(1,2),(1,3),(2,1),(2,3),(3,1),(3,2)]
    @test Permutations{3}(1:3)|>collect==[(1,2,3),(1,3,2),(2,1,3),(2,3,1),(3,1,2),(3,2,1)]
end

@testset "DulPermutations" begin
    @test DulPermutations{0}(1:3)|>collect==[()]
    @test DulPermutations{1}(1:3)|>collect==[(1,),(2,),(3,)]
    @test DulPermutations{2}(1:3)|>collect==[(1,1),(1,2),(1,3),(2,1),(2,2),(2,3),(3,1),(3,2),(3,3)]
    @test DulPermutations{3}(1:3)|>collect==[(i,j,k) for i=1:3 for j=1:3 for k=1:3]
end


text:
# L2,1 norm/Sum of norms of columns or rows (times a constant)

"""
  NormL21(λ::Real=1.0, dim=1)

Returns the function `g(X) = λsum(||x_i||)` for a nonnegative parameter `λ ⩾ 0`,
where x_i are the columns of X if dim == 1, or the rows of X if dim == 2.
In words, it is the (weighted) sum of the Euclidean norm of the columns (rows) of X.
"""

immutable NormL21{R <: Real, I <: Integer} <: ProximableConvex
  lambda::R
  dim::I
  function NormL21(lambda::R, dim::I)
    if lambda < 0
      error("parameter λ must be nonnegative")
    else
      new(lambda, dim)
    end
  end
end

NormL21{R <: Real, I <: Integer}(lambda::R=1.0, dim::I=1) = NormL21{R, I}(lambda, dim)

function (f::NormL21){T <: RealOrComplex}(X::AbstractArray{T,2})
  nslice = 0.0
  n21X = 0.0
  if f.dim == 1
    for j = 1:size(X,2)
      nslice = 0.0
      for i = 1:size(X,1)
      	nslice += abs(X[i,j])^2
      end
      n21X += sqrt(nslice)
    end
  elseif f.dim == 2
    for i = 1:size(X,1)
      nslice = 0.0
      for j = 1:size(X,2)
      	nslice += abs(X[i,j])^2
      end
      n21X += sqrt(nslice)
    end
  end
  return f.lambda*n21X
end

function prox!{T <: RealOrComplex}(f::NormL21, X::AbstractArray{T,2}, Y::AbstractArray{T,2}, gamma::Real=1.0)
  gl = gamma*f.lambda
  nslice = zero(Float64)
  n21X = zero(Float64)
  if f.dim == 1
    for j = 1:size(X,2)
      nslice = 0.0
      for i = 1:size(X,1)
      	nslice += abs(X[i,j])^2
      end
      nslice = sqrt(nslice)
      scal = 1-gl/nslice
      scal = scal <= 0.0 ? 0.0 : scal
      for i = 1:size(X,1)
        Y[i,j] = scal*X[i,j]
      end
      n21X += scal*nslice
    end
  elseif f.dim == 2
    for i = 1:size(X,1)
      nslice = 0.0
      for j = 1:size(X,2)
      	nslice += abs(X[i,j])^2
      end
      nslice = sqrt(nslice)
      scal = 1-gl/nslice
      scal = scal <= 0.0 ? 0.0 : scal
      for j = 1:size(X,2)
        Y[i,j] = scal*X[i,j]
      end
      n21X += scal*nslice
    end
  end
  return f.lambda*n21X
end

fun_name(f::NormL21) = "sum of Euclidean norms"
fun_dom(f::NormL21) = "AbstractArray{Real,2}, AbstractArray{Complex,2}"
fun_expr(f::NormL21) = "x ↦ λsum(||x_i||)"
fun_params(f::NormL21) = "λ = $(f.lambda), dim = $(f.dim)"

function prox_naive{T <: RealOrComplex}(f::NormL21, X::AbstractArray{T,2}, gamma::Real=1.0)
  Y = max.(0, 1-f.lambda*gamma./sqrt.(sum(abs.(X).^2, f.dim))).*X
  return Y, f.lambda*sum(sqrt.(sum(abs.(Y).^2, f.dim)))
end


text:
using LEDStrip

const RGB_COLORS = [
    0xff0000,
    0xff0000,
    0x00ff00,
    0x00ff00,
    0x0000ff,
    0x0000ff,
    0x000000,
    0x000000,
]

const SNAKE_COLORS = [
    0xffffff, 0xffffff, 0xffffff, 0xffffff, 0xffffff, 0xffffff, 0xffffff,
    0xff0000, 0x00ff00, 0x0000ff,
    0x000000, 0x000000,
]

const RAINBOW_COLORS = [
    0xFF0000,   # red
    0xFF7F00,   # orange
    0xFFFF00,   # yellow
    0x00FF00,   # green
    0x0000FF,   # lightblue
    0x4B0082,   # blue
    0x9400D3,   # violet
]

const PIXEL_COUNT = 100
s = StripSPI(1; pixel_count = PIXEL_COUNT)

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

@info "R-R-G-G-B-B..."

set_pixels!(s, RGB_COLORS, replicate = false)
show_pixels(s)
sleep(5.)
hide_pixels(s)
clean_pixels!(s)
sleep(1.)

@info "R-R-G-G-B-B repeated..."
set_pixels!(s, RGB_COLORS, replicate = true)
show_pixels(s)
sleep(5.)
hide_pixels(s)
clean_pixels!(s)

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

@info "All RED... All GREEN... All BLUE... All WHITE..."

set_pixels!(s, [0xff0000]; replicate = true)
show_pixels(s)
sleep(5)
set_pixels!(s, [0x00ff00]; replicate = true)
show_pixels(s)
sleep(5)
set_pixels!(s, [0x0000ff]; replicate = true)
show_pixels(s)
sleep(5)
set_pixels!(s, [0xffffff]; replicate = true)
show_pixels(s)
sleep(5)

hide_pixels(s)
clean_pixels!(s)
sleep(1.)

@info "Rainbow 7 colors..."

set_pixels!(s, RAINBOW_COLORS)
show_pixels(s)
sleep(5)
hide_pixels(s)
clean_pixels!(s)
sleep(1.)

@info "W-W-W-W-W-W-W-W-R-G-B snake moves forward and back..."

set_pixels!(s, SNAKE_COLORS)

for i in 1:90
    show_pixels(s)
    shift_forward!(s)
    sleep(0.1)
end

for i in 1:90
    show_pixels(s)
    shift_backward!(s)
    sleep(0.1)
end

hide_pixels(s)
clean_pixels!(s)

@info "Stop."


text:
using Test, AdvancedHMC
include("common.jl")

ϵ = 0.01
lf = Leapfrog(ϵ)

θ_init = randn(D)
h = Hamiltonian(UnitEuclideanMetric(D), logπ, ∂logπ∂θ)
r_init = AdvancedHMC.rand_momentum(h)

n_steps = 10

@testset "step(::Leapfrog) against steps(::Leapfrog)" begin
    θ_step, r_step = copy(θ_init), copy(r_init)

    t_step = @elapsed for i = 1:n_steps
        θ_step, r_step, _ = AdvancedHMC.step(lf, h, θ_step, r_step)
    end

    t_steps = @elapsed θ_steps, r_steps, _ = AdvancedHMC.steps(lf, h, θ_init, r_init, n_steps)

    @info "Performance of step() v.s. steps()" n_steps t_step t_steps t_step / t_steps

    @test θ_step ≈ θ_steps atol=DETATOL
    @test r_step ≈ r_steps atol=DETATOL
end

using Turing: Inference
@testset "steps(::Leapfrog) against Turing.Inference._leapfrog()" begin
    t_Turing = @elapsed θ_Turing, r_Turing, _ = Inference._leapfrog(θ_init, r_init, n_steps, ϵ, x -> (nothing, ∂logπ∂θ(x)))
    t_AHMC = @elapsed θ_AHMC, r_AHMC, _ = AdvancedHMC.steps(lf, h, θ_init, r_init, n_steps)
    @info "Performance of leapfrog of AdvancedHMC v.s. Turing" n_steps t_Turing t_AHMC t_Turing / t_AHMC

    @test θ_Turing ≈ θ_AHMC atol=DETATOL
    @test r_Turing ≈ r_AHMC atol=DETATOL
end

using LinearAlgebra: dot
using Statistics: mean
@testset "Eq (2.11) from Neal (2011)" begin
    D = 1
    negU(q::AbstractVector{T}) where {T<:Real} = -dot(q, q) / 2
    ∂negU∂q = q -> gradient(negU, q)

    ϵ = 0.01
    lf = Leapfrog(ϵ)

    q_init = randn(D)
    h = Hamiltonian(UnitEuclideanMetric(D), negU, ∂negU∂q)
    p_init = AdvancedHMC.rand_momentum(h)

    q, p = copy(q_init), copy(p_init)

    n_steps = 10_000
    qs = zeros(n_steps)
    ps = zeros(n_steps)
    Hs = zeros(n_steps)
    for i = 1:n_steps
        q, p, _ = AdvancedHMC.step(lf, h, q, p)
        qs[i] = q[1]
        ps[i] = p[1]
        Hs[i] = AdvancedHMC.hamiltonian_energy(h, q, p)
    end

    # Throw first 1_000 steps
    qs = qs[1_000:end]
    ps = ps[1_000:end]
    Hs = Hs[1_000:end]

    # Check if all points located at a cirle centered at the origin
    rs = sqrt.(qs.^2 + ps.^2)
    @test all(x-> abs(x - mean(rs)) < 2e-3, rs)

    # Check if the Hamiltonian energy is stable
    @test all(x-> abs(x - mean(Hs)) < 2e-3, Hs)
end


text:
using Gtk

@testset "list" begin

@testset "v3 populate with push!" begin
    cbstr = GtkComboBoxText()
    push!(cbstr, "1", "abc")
    push!(cbstr, "2", "xyz")

    set_gtk_property!(cbstr, :active, 0)
    @test get_gtk_property(cbstr, "active-id", String) == "1"

    set_gtk_property!(cbstr, :active, 1)
    @test get_gtk_property(cbstr, "active-id", String) == "2"

    cbsym = GtkComboBoxText()
    push!(cbsym, :a, "abc")
    push!(cbsym, :b, "xyz")

    set_gtk_property!(cbsym, :active, 0)
    @test get_gtk_property(cbsym, "active-id", String) == "a"

    set_gtk_property!(cbsym, :active, 1)
    @test get_gtk_property(cbsym, "active-id", String) == "b"
end

@testset "v3 populate with pushfirst!" begin
    cbstr = GtkComboBoxText()
    pushfirst!(cbstr, "1", "abc")
    pushfirst!(cbstr, "2", "xyz")

    set_gtk_property!(cbstr, :active, 0)
    @test get_gtk_property(cbstr, "active-id", String) == "2"

    set_gtk_property!(cbstr, :active, 1)
    @test get_gtk_property(cbstr, "active-id", String) == "1"

    cbsym = GtkComboBoxText()
    pushfirst!(cbsym, :a, "abc")
    pushfirst!(cbsym, :b, "xyz")

    set_gtk_property!(cbsym, :active, 0)
    @test get_gtk_property(cbsym, "active-id", String) == "b"

    set_gtk_property!(cbsym, :active, 1)
    @test get_gtk_property(cbsym, "active-id", String) == "a"
end

@testset "v3 populate with insert!" begin
    cbstr = GtkComboBoxText()
    push!(cbstr, "1", "abc")
    push!(cbstr, "3", "mno")
    insert!(cbstr, 2, "2", "xyz")

    set_gtk_property!(cbstr, :active, 0)
    @test get_gtk_property(cbstr, "active-id", String) == "1"

    set_gtk_property!(cbstr, :active, 1)
    @test get_gtk_property(cbstr, "active-id", String) == "2"

    set_gtk_property!(cbstr, :active, 2)
    @test get_gtk_property(cbstr, "active-id", String) == "3"

    cbsym = GtkComboBoxText()
    push!(cbsym, :a, "abc")
    push!(cbsym, :c, "mno")
    insert!(cbsym, 2, :b, "xyz")

    set_gtk_property!(cbsym, :active, 0)
    @test get_gtk_property(cbsym, "active-id", String) == "a"

    set_gtk_property!(cbsym, :active, 1)
    @test get_gtk_property(cbsym, "active-id", String) == "b"

    set_gtk_property!(cbsym, :active, 2)
    @test get_gtk_property(cbsym, "active-id", String) == "c"
end

end

max_stars_repo_name: adamwillisMastery/DeepLearner

max_issues_repo_name: adamwillisMastery/DeepLearner

max_forks_repo_name: adamwillisXanax/DeepLearner


text:
module StaticTrafficAssignment

using CSV: read
using DataFrames: DataFrameRow, nrow
using LightGraphs: AbstractGraph, AbstractEdge
using Distributed: @distributed
using DataStructures: PriorityQueue, dequeue!
using LinearAlgebra: Diagonal

import Base: show

# i/o
export readtntpdata

# interface
export
        AbstractNework, AbstractLink, AbstractZone, AbstractBush,
        has_link, add_link!, rem_link!, add_zone!, rem_zone!,
        idx, numnodes, numzones, numlinks, numsources, numsinks,
        sources, sinks, links, zones,
        outneighbors, inneighbors, outdegree, indegree, fstar, bstar,
        upn, dwn, props,
        id, issource, issink, throughflowallowed, src, net, topo_order

# cost functions
export CostFunction, CostFunctionUE, CostFunctionSO, TimeFunctionContainer, BPR

# algorithms
export
        dijkstra, topologicalorder, orderednodes, acyclic,
        allornothing, msa, frankwolfe, conjugatefrankwolfe, algorithmB

# simple types
export SimpleNetwork, SimpleLink, SimpleZone, SimpleBush

include("utils.jl")

# i/o
include("data_parsers/utils.jl")
include("data_parsers/tntp_parser.jl")

# interface
include("interface.jl")

# costs
include("costfunctions/utils.jl")
include("costfunctions/time_function_container.jl")
include("costfunctions/generator.jl")
include("costfunctions/bpr.jl")

# algorithms
include("algorithms/utils.jl")
include("algorithms/shortestpaths/dijkstra.jl")
include("algorithms/topologicalorder.jl")
include("algorithms/shortestpaths/acyclic.jl")
include("algorithms/allornothing.jl")

include("algorithms/linkbased/msa.jl")
include("algorithms/linkbased/frankwolfe.jl")
include("algorithms/linkbased/conjugatefrankwolfe.jl")

include("algorithms/bushbased/algorithmB.jl")

# simple types
include("SimpleNetworks/SimpleNetworks.jl")
using .SimpleNetworks

end # module StaticTrafficAssignment


text:
"""
    NeuralNetworkKernel()

Neural network kernel function.

```math
    κ(x, y) =  asin(x' * y / sqrt[(1 + x' * x) * (1 + y' * y)])
```
# Significance
Neal (1996) pursued the limits of large models, and showed that a Bayesian neural network
becomes a Gaussian process with a **neural network kernel** as the number of units
approaches infinity. Here, we give the neural network kernel for single hidden layer
Bayesian neural network with erf (Error Function) as activation function.

# References:
- [GPML Pg 105](http://www.gaussianprocess.org/gpml/chapters/RW4.pdf)
- [Neal(1996)](https://www.cs.toronto.edu/~radford/bnn.book.html)
- [Andrew Gordon's Thesis Pg 45](http://www.cs.cmu.edu/~andrewgw/andrewgwthesis.pdf)
"""
struct NeuralNetworkKernel <: BaseKernel end

function (κ::NeuralNetworkKernel)(x, y)
    return asin(dot(x, y) / sqrt((1 + sum(abs2, x)) * (1 + sum(abs2, y))))
end

Base.show(io::IO, κ::NeuralNetworkKernel) = print(io, "Neural Network Kernel")


text:
# deps/jl needs these to be available
libclang_path() = libclang.Clang_jll.libclang_path
libclang_version() = libclang_version(libclang_path())
function libclang_version(path::String)
	dir = joinpath(dirname(dirname(path)), "lib", "clang")
	for entry in readdir(dir)
		m = match(r"^(\d+)\.\d+\.\d+$", entry)
		!isnothing(m) && return m.captures[1]
	end
	error("Failed to determine the version of libclang")
end

module libclang
	import Clang_jll
	using ..C

	include(joinpath(dirname(@__DIR__), "deps/libclang-$(C.libclang_version(Clang_jll.libclang_path)).jl"))
end

using .libclang

function _string(func, args...)
	cxstr = nothing
	try
		cxstr = func(args...)
		ptr = clang_getCString(cxstr)
		return ptr == C_NULL ? nothing : unsafe_string(ptr)
	finally
		isnothing(cxstr) || clang_disposeString(cxstr)
	end
end

Base.string(cursor::CXType) = _string(clang_getTypeSpelling, cursor)
Base.string(cursor::CXCursor) = _string(clang_getCursorSpelling, cursor)
Base.string(diag::CXDiagnostic) = _string(clang_formatDiagnostic, diag, clang_defaultDiagnosticDisplayOptions())
Base.string(tu::CXTranslationUnit, token::CXToken) = _string(clang_getTokenSpelling, tu, token)

function children(cursor::CXCursor)
	function visitor(c::CXCursor, p::CXCursor, cs::CXClientData)::CXChildVisitResult
		push!(unsafe_pointer_to_objref(cs)::Vector{CXCursor}, c)
		return CXChildVisit_Continue
	end

	cursors = CXCursor[]
	Bool(clang_visitChildren(cursor, visitor, pointer_from_objref(cursors))) && error("Failed to visit cursor children")
	return cursors
end

function tokenize(tu::CXTranslationUnit, cursor::CXCursor)
	tokens = Ref{Cptr{CXToken}}(C_NULL)
	count = Ref{Cuint}(0)
	try
		rng = clang_getCursorExtent(cursor)
		cursrange = getlocation(rng)

		clang_tokenize(tu, rng, tokens, count)

		# tokens are wrong when cursor involves a macro expansion, so filter down to the cursor's lines
		return filter(map(i -> unsafe_load(tokens[], i), 1:count[])) do tok
			rng = clang_getTokenExtent(tu, tok)
			loc = getlocation(rng)

			isnothing(loc) && return false
			isnothing(cursrange) && return false
			first(loc).file == first(cursrange).file || return false

			return first(loc).line >= first(cursrange).line && last(loc).line <= last(cursrange).line
		end
	finally
		tokens[] == C_NULL || clang_disposeTokens(tu, tokens[], count[])
	end
end

getlocation(cursor::CXCursor) = getlocation(clang_getCursorExtent(cursor))

function getlocation(rng::CXSourceRange)
	Bool(clang_Range_isNull(rng)) && return nothing

	start = clang_getRangeStart(rng)
	stop  = clang_getRangeEnd(rng)
	Bool(clang_Location_isInSystemHeader(start)) && return nothing
	Bool(clang_Location_isInSystemHeader(stop))  && return nothing

	locs = map(getlocation, (start, stop))
	any(isnothing, locs) && return nothing

	return locs
end

function getlocation(loc::CXSourceLocation)
	cxfile = Ref{CXFile}(C_NULL)
	line = Ref{Cuint}(0)
	col = Ref{Cuint}(0)
	offset = Ref{Cuint}(0)
	clang_getSpellingLocation(loc, cxfile, line, col, offset)

	return cxfile[] == C_NULL ? nothing : CodeLocation(_string(clang_getFileName, cxfile[]), Int(line[]), Int(col[]))
end


text:
function paste(sep, s...)
    output = ""
    for i in s
        if(isa(i,String))
            output *= i * sep
        elseif isa(i,Array{ASCIIString,1})
            output *= paste(sep,i) * sep
        end
    end
    return output
end
function paste(sep, s::Array{ASCIIString,1})
    output = ""
    for i in s
        output *= i * sep
    end
    return output
end

heaviside(x) = 0.5 * (sign(x) + 1)
dirac(x) = abs(x) <= eps(1.0)  ? 0 : Inf
min(n::Number) = n
max(n::Number) = n

# These functions are reserved in order to have a basic mathematical lingo to write the models
# It also includes operators, reserved words such as time and control flow functions
const list_of_functions = ["abs", "abs2", "acos", "acosd", "acosh", "acot", "acotd",
                     "acoth", "acsc", "acscd", "acsch", "asec", "asecd", "asech",
                     "asin", "asind", "asinh", "atan", "atan2", "atand", "atanh",
                     "cbrt", "cos", "cosd", "cosh", "cot", "cotd", "coth", "csc",
                     "cscd", "csch", "exp", "exp2", "expm1", "log", "log", "log10",
                     "log1p", "log2", "sec", "secd", "sech", "sin", "sind", "sinh",
                     "sqrt", "tan", "tand", "tanh", "hypot","*","+","-","/","^","%","⋅",
                     "√","ifelse","min","max","<",">","=","<=",">=","time"];


text:
abstract type VegasDamping end
struct NoDamping <: VegasDamping end

"""
Damping for quantile estimation as proposed in original Vegas Paper by Lepage.
"""
@qstruct LepageDamping{T}(;
    alpha::T=1) <: VegasDamping

# TODO renames, this is not really Vegas algorithm
@settable @qstruct Vegas{RNG,REG}(
        neval::Int=10^4;
        rng::RNG=GLOBAL_PAR_RNG,
        regularization::REG=LepageDamping(),
    ) do
        @argcheck neval > 2
end <: MCAlgorithm

@qstruct VegasGrid{N,T}(
    boundaries::NTuple{N,Vector{T}}
   ) do
    for i in 1:N
        @argcheck length(boundaries[i]) >= 2
        @argcheck issorted(boundaries[i])
    end
end

function domain(iq::VegasGrid)
    lo = map(first, iq.boundaries)
    hi = map(last , iq.boundaries)
    Domain(SVector(lo), SVector(hi))
end

function volume(iq::VegasGrid)
    volume(domain(iq))
end

function uniform(rng, iq::VegasGrid)
    uniform(rng, domain(iq))
end

@qstruct CDF{P,V}(
    positions::Vector{P},
    values::Vector{V},
    ) do
    @argcheck first(values) == 0.
    @argcheck issorted(values)
    @argcheck issorted(positions)
    @argcheck length(positions) == length(values)
    @argcheck length(positions) >= 2
end

function linterpol(x, (x_min, x_max), (y_min, y_max))
    @argcheck x_min <= x <= x_max
    if x == x_min
        y_min
    elseif x == x_max
        y_max
    else
        w1 = x_max - x
        w2 = x - x_min
        w = w1 + w2
        @assert w1 + w2 ≈ w
        y_min * w1/w + y_max * w2/w
    end
end

function quantiles(cdf::CDF, nwalls)
    ret = [first(cdf.positions)]
    qvals = range(first(cdf.values), stop=last(cdf.values), length=nwalls)
    i = 1
    for qval in qvals[2:end]
        while cdf.values[i] < qval
            i += 1
        end
        xs = cdf.values[i-1], cdf.values[i]
        ys = cdf.positions[i-1], cdf.positions[i]
        q = linterpol(qval, xs, ys)
        push!(ret, q)
    end
    @assert length(ret) == nwalls
    ret
end

function Base.size(s::VegasGrid)
    map(s.boundaries) do b
        length(b) - 1
    end
end
Base.length(s::VegasGrid) = prod(size(s))
Base.eachindex(s::VegasGrid) = CartesianIndices(s)
Base.axes(s::VegasGrid) = map(Base.OneTo, size(s))
function Base.CartesianIndices(s::VegasGrid)
    CartesianIndices(axes(s))
end
function Base.LinearIndices(s::VegasGrid)
    LinearIndices(axes(s))
end

function rand_index(rng, s::VegasGrid)
    li = rand(rng, LinearIndices(s))
    CartesianIndices(s)[li]
end

function Base.getindex(s::VegasGrid, index::CartesianIndex)
    lower = map(s.boundaries, Tuple(index)) do walls, i
        walls[i]
    end
    upper = map(s.boundaries, Tuple(index)) do walls, i
        walls[i+1]
    end
    Domain(SVector(lower), SVector(upper))
end

function draw_x_index_cell(rng, s::VegasGrid)
    i = rand_index(rng,s)
    cell = s[i]
    x = uniform(rng, cell)
    (position=x, index=i, cell=cell)
end

function draw(rng, iq::VegasGrid)
    s = draw_x_index_cell(rng, iq)
    wt = volume(s.cell) * length(iq)
    (position=s.position, weight = wt)
end

function estimate_pdf(histogram, ::NoDamping)
    pdf = map(/, histogram.sums, histogram.counts)
    pdf ./= sum(pdf)
    pdf
end

function estimate_pdf(histogram, r::LepageDamping)
    pdf = estimate_pdf(histogram, NoDamping())
    updf = map(pdf) do p
        ((p - 1) / log(p))^r.alpha
    end
    normalize!(updf, 1)
end

function estimate_cdf(histogram, alg::Vegas)
    pdf = estimate_pdf(histogram, alg.regularization)
    cdf_vals = cumsum(pdf)
    z = zero(eltype(cdf_vals))
    pushfirst!(cdf_vals, z)
    CDF(histogram.walls, cdf_vals)
end

function tuneonce(f, iq::VegasGrid,
              alg::Vegas=Vegas();
              neval=1000,
              outsize=size(iq),
            )
    s = draw_x_index_cell(alg.rng, iq)
    y = float(norm(f(s.position)))

    hists = map(iq.boundaries) do xs
        nbins = length(xs) - 1
        (
            walls = xs,
            counts = fill(0, nbins),
            sums = fill(zero(typeof(y)), nbins),
        )
    end

    for _ in 1:neval
        s = draw_x_index_cell(alg.rng, iq)
        y = norm(f(s.position))
        cell = s.cell
        for i in eachindex(iq.boundaries)
            h = hists[i]
            wt = cell.upper[i] - cell.lower[i]
            h.counts[s.index[i]] += 1
            h.sums[s.index[i]] += y*wt
        end
    end

    bdries_new = map(hists, outsize) do h, ncells
        nwalls = ncells + 1
        cdf = estimate_cdf(h, alg)
        quantiles(cdf, nwalls)
    end

    VegasGrid(bdries_new)
end

function integral_kernel(f, dom::VegasGrid, alg::Vegas)
    mc_kernel(f, alg.rng, dom; neval=alg.neval)
end

function default_grid_size(dom::Domain{N}) where {N}
    ntuple(_ -> 100, Val(N))
end

function equidistant_grid(
        dom::Domain,
        size=default_grid_size(dom)
    )

    bdries = map(Tuple(dom.lower), Tuple(dom.upper), size) do lo, hi, nbins
        nwalls = nbins + 1
        r = range(lo, stop=hi, length=nwalls)
        collect(r)
    end
    VegasGrid(bdries)
end

function initvr(f, dom::Domain, alg::Vegas)
    equidistant_grid(dom)
end

function tune(f, iq::VegasGrid, alg::Vegas)
    iq = tuneonce(f, iq, alg, neval=2000)
    iq = tuneonce(f, iq, alg, neval=2000)
    iq = tuneonce(f, iq, alg, neval=2000)
    iq = tuneonce(f, iq, alg, neval=2000)
    iq = tuneonce(f, iq, alg, neval=2000)
end

function tune(f, dom::Domain, alg)
    iq = initvr(f, dom, alg)
    tune(f, iq, alg)
end

function integral_kernel(f, dom::Domain, alg::Vegas)
    iq = tune(f, dom, alg)
    integral_kernel(f, iq, alg)
end

function canonicalize(f, dom::VegasGrid, alg::Vegas)
    f, dom, alg
end


text:
#
# This file is a part of KronLinInv. License is MIT
# Copyright (c) 2019 Andrea Zunino
#

##==========================================================
module KronLinInv

export CovMats,FwdOps,KLIFactors
export calcfactors,posteriormean,blockpostcov

using Distributed
using LinearAlgebra

include("kronlininv_parallel.jl")

include("kronlininv_serial.jl")

##==========================================================
end # module
#########################################


text:
@safetestset "isbefore" begin
    using Validator

    # date
    @test Validator.isbefore("2022-02-28", "2022-02-27") === false
    @test Validator.isbefore("2022-02-28", "2022-02-28") === false
    @test Validator.isbefore("2022-02-28", "2022-03-01") === true

    # datetime
    @test Validator.isbefore("2022-02-28T12:00:00", "2022-02-27T12:00:00") === false
    @test Validator.isbefore("2022-02-28T12:00:00", "2022-02-28T09:00:00") === false
    @test Validator.isbefore("2022-02-28T12:00:00", "2022-02-28T12:00:00") === false
    @test Validator.isbefore("2022-02-28T12:00:00", "2022-03-01T12:00:00") === true

    # default
    @test Validator.isbefore("2022-02-27") === true
    @test Validator.isbefore("2050-02-28") === false
    @test Validator.isbefore("2022-02-28T12:00:00") === true
    @test Validator.isbefore("2050-02-28T09:00:00") === false
end


text:
module Merlin

using Base.Threads
@info "# CPU threads: $(nthreads())"

const GDICT = []

include("cuda/CUDA.jl")
using .CUDA

using Markdown
import LinearAlgebra.BLAS: scal!, axpy!, gemv, gemv!, gemm, gemm!
export gemv, gemm
#export Parametric
#abstract type Parametric end

const UniArray{T,N} = Union{Array{T,N},CuArray{T,N}}
const UniVector{T} = UniArray{T,1}
const UniMatrix{T} = UniArray{T,2}

include("add.jl")
include("var.jl")
include("graph.jl")
include("dataset.jl")
include("check.jl")
include("config.jl")
include("device.jl")

for name in [
    "activation/crelu",
    "activation/elu",
    "activation/leaky_relu",
    "activation/relu",
    "activation/selu",
    "activation/sigmoid",
    "activation/swish",
    "activation/tanh",

    "cnn/conv1d",
    # "cnn/conv2d",

    "loss/crossentropy",
    "loss/l2",
    "loss/mse",
    "loss/softmax_crossentropy",
    "math/arithmetic",
    "math/broadcast",

    "reduction/max",
    #"reduction/mean",
    #"reduction/sum",

    "rnn/lstm",

    "blas",
    "concat",
    "dropout",
    "getindex",
    "linear",
    "lookup",
    "pack",
    "reshape",
    "softmax",
    "split",
    "window1d"
    ]
    include("functions/$name.jl")
end

include("initializers/fill.jl")
include("initializers/normal.jl")
include("initializers/orthogonal.jl")
include("initializers/orthonormal.jl")
include("initializers/uniform.jl")
include("initializers/xavier.jl")

include("optimizers/adagrad.jl")
include("optimizers/adam.jl")
include("optimizers/sgd.jl")

#include("datasets/Datasets.jl")
#include("caffe/Caffe.jl")

end


text:
Base.@kwdef mutable struct NFFTParams{T}
  m::Int = 4
  σ::T = 2.0
  reltol::T = 1e-9
  window::Symbol = :kaiser_bessel
  LUTSize::Int64 = 0
  precompute::PrecomputeFlags = POLYNOMIAL
  sortNodes::Bool = false
  storeApodizationIdx::Bool = false
  blocking::Bool = true
end

mutable struct NFFTPlan{T,D,R} <: AbstractNFFTPlan{T,D,R}
    N::NTuple{D,Int64}
    NOut::NTuple{R,Int64}
    M::Int64
    x::Matrix{T}
    n::NTuple{D,Int64}
    dims::UnitRange{Int64}
    params::NFFTParams{T}
    forwardFFT::FFTW.cFFTWPlan{Complex{T},-1,true,D,UnitRange{Int64}}
    backwardFFT::FFTW.cFFTWPlan{Complex{T},1,true,D,UnitRange{Int64}}
    tmpVec::Array{Complex{T},D}
    tmpVecHat::Array{Complex{T},D}
    # Caches for deconvolve
    deconvolveIdx::Array{Int64,1}
    windowHatInvLUT::Vector{Vector{T}}
    # Cache for precompute = LUT
    windowLinInterp::Vector{T}
    # Cache for precompute = POLYNOMIAL
    windowPolyInterp::Matrix{T}
    # Caches for blocking
    blocks::Array{Array{Complex{T},D},D}
    nodesInBlock::Array{Vector{Int64},D}
    blockOffsets::Array{NTuple{D,Int64},D}
    idxInBlock::Array{Matrix{Tuple{Int,T}},D}
    windowTensor::Array{Array{T,3},D}
    # Cache for precompute = FULL
    B::SparseMatrixCSC{T,Int64}
end

function Base.copy(p::NFFTPlan{T,D,R}) where {T,D,R}
    tmpVec = similar(p.tmpVec)
    tmpVecHat = similar(p.tmpVecHat)
    deconvolveIdx = copy(p.deconvolveIdx)
    windowLinInterp = copy(p.windowLinInterp)
    windowPolyInterp = copy(p.windowPolyInterp)
    windowHatInvLUT = copy(p.windowHatInvLUT)
    B = copy(p.B)
    blocks = deepcopy(p.blocks)
    nodesInBlock = deepcopy(p.nodesInBlock)
    blockOffsets = copy(p.blockOffsets)
    idxInBlock = copy(p.idxInBlock)
    windowTensor = copy(p.windowTensor)
    x = p.x

    FP = plan_fft!(tmpVec, p.dims; flags = p.forwardFFT.flags)
    BP = plan_bfft!(tmpVec, p.dims; flags = p.backwardFFT.flags)

    return NFFTPlan{T,D,R}(p.N, p.NOut, p.M, x, p.n, p.dims, p.params, FP, BP, tmpVec,
                           tmpVecHat, deconvolveIdx, windowHatInvLUT, windowLinInterp, windowPolyInterp,
                           blocks, nodesInBlock, blockOffsets, idxInBlock, windowTensor, B)
end

################
# constructor
################

function NFFTPlan(x::Matrix{T}, N::NTuple{D,Int}; dims::Union{Integer,UnitRange{Int64}}=1:D,
                 fftflags=nothing, kwargs...) where {T,D}

    checkNodes(x)

    params, N, NOut, M, n, dims_ = initParams(x, N, dims; kwargs...)

    if length(NOut) > 1
      params.precompute = LINEAR
    end

    tmpVec = Array{Complex{T},D}(undef, n)

    fftflags_ = (fftflags != nothing) ? (flags=fftflags,) : NamedTuple()
    FP = plan_fft!(tmpVec, dims_; fftflags_...)
    BP = plan_bfft!(tmpVec, dims_; fftflags_...)

    calcBlocks = (params.precompute == LINEAR ||
                  params.precompute == TENSOR ||
                  params.precompute == POLYNOMIAL ) &&
                     params.blocking && length(dims_) == D

    blocks, nodesInBlocks, blockOffsets, idxInBlock, windowTensor = precomputeBlocks(x, n, params, calcBlocks)

    windowLinInterp, windowPolyInterp, windowHatInvLUT, deconvolveIdx, B =
            precomputation(x, N[dims_], n[dims_], params)

    U = params.storeApodizationIdx ? N : ntuple(d->0,D)
    tmpVecHat = Array{Complex{T},D}(undef, U)

    NFFTPlan(N, NOut, M, x, n, dims_, params, FP, BP, tmpVec, tmpVecHat,
                       deconvolveIdx, windowHatInvLUT, windowLinInterp, windowPolyInterp,
                       blocks, nodesInBlocks, blockOffsets, idxInBlock, windowTensor, B)
end

function AbstractNFFTs.nodes!(p::NFFTPlan{T}, x::Matrix{T}) where {T}
    checkNodes(x)

    # Sort nodes in lexicographic way
    if p.params.sortNodes
        x .= sortslices(x, dims=2)
    end

    calcBlocks = (p.params.precompute == LINEAR ||
                  p.params.precompute == TENSOR ||
                  p.params.precompute == POLYNOMIAL ) &&
                     p.params.blocking && length(p.dims) == length(p.N)

    blocks, nodesInBlocks, blockOffsets, idxInBlock, windowTensor = precomputeBlocks(x, p.n, p.params, calcBlocks)

    windowLinInterp, windowPolyInterp, windowHatInvLUT, deconvolveIdx, B =
       precomputation(x, p.N, p.n, p.params)

    p.blocks = blocks
    p.nodesInBlock = nodesInBlocks
    p.blockOffsets = blockOffsets
    p.idxInBlock = idxInBlock
    p.windowTensor = windowTensor

    p.M = size(x, 2)
    p.windowLinInterp = windowLinInterp
    p.windowPolyInterp = windowPolyInterp
    p.windowHatInvLUT = windowHatInvLUT
    p.deconvolveIdx = deconvolveIdx
    p.B = B
    p.x = x

    return p
end

function Base.show(io::IO, p::NFFTPlan{T,D,R}) where {T,D,R}
    print(io, "NFFTPlan with ", p.M, " sampling points for an input array of size",
           p.N, " and an output array of size", p.NOut, " with dims ", p.dims)
end

AbstractNFFTs.size_in(p::NFFTPlan) = p.N
AbstractNFFTs.size_out(p::NFFTPlan) = p.NOut

################
# nfft functions
################

function LinearAlgebra.mul!(fHat::StridedArray, p::NFFTPlan{T,D,R}, f::AbstractArray;
               verbose=false, timing::Union{Nothing,TimingStats} = nothing) where {T,D,R}
    consistencyCheck(p, f, fHat)

    fill!(p.tmpVec, zero(Complex{T}))
    t1 = @elapsed @inbounds deconvolve!(p, f, p.tmpVec)
    t2 = @elapsed p.forwardFFT * p.tmpVec
    t3 = @elapsed @inbounds convolve!(p, p.tmpVec, fHat)
    if verbose
        @info "Timing: apod=$t1 fft=$t2 conv=$t3"
    end
    if timing != nothing
      timing.conv = t3
      timing.fft = t2
      timing.apod = t1
    end
    return fHat
end

function LinearAlgebra.mul!(f::StridedArray, pl::Adjoint{Complex{T},<:NFFTPlan{T}}, fHat::AbstractArray;
                       verbose=false, timing::Union{Nothing,TimingStats} = nothing) where {T}
    p = pl.parent
    consistencyCheck(p, f, fHat)

    t1 = @elapsed @inbounds convolve_transpose!(p, fHat, p.tmpVec)
    t2 = @elapsed p.backwardFFT * p.tmpVec
    t3 = @elapsed @inbounds deconvolve_transpose!(p, p.tmpVec, f)
    if verbose
        @info "Timing: conv=$t1 fft=$t2 apod=$t3"
    end
    if timing != nothing
      timing.conv_adjoint = t1
      timing.fft_adjoint = t2
      timing.apod_adjoint = t3
    end
    return f
end


text:
@testset "DeviceModel Tests" begin
    @test_throws ArgumentError DeviceModel(ThermalGen, ThermalStandardUnitCommitment)
    @test_throws ArgumentError DeviceModel(ThermalStandard, PSI.AbstractThermalFormulation)
    @test_throws ArgumentError NetworkModel(PM.AbstractPowerModel)
end

@testset "NetworkModel Tests" begin
    @test_throws ArgumentError NetworkModel(PM.AbstractPowerModel)
end

@testset "ServiceModel Tests" begin
    @test_throws ArgumentError ServiceModel(AGC, PSI.AbstractServiceFormulation, "TestName")
    @test_throws ArgumentError ServiceModel(
        VariableReserve{PSY.ReserveUp},
        PSI.AbstractServiceFormulation,
        "TestName2",
    )
end

@testset "Feedforward Struct Tests" begin
    ffs = [
        UpperBoundFeedforward(
            component_type=RenewableDispatch,
            source=ActivePowerVariable,
            affected_values=[ActivePowerVariable],
        ),
        LowerBoundFeedforward(
            component_type=RenewableDispatch,
            source=ActivePowerVariable,
            affected_values=[ActivePowerVariable],
        ),
        SemiContinuousFeedforward(
            component_type=ThermalMultiStart,
            source=OnVariable,
            affected_values=[ActivePowerVariable, ReactivePowerVariable],
        ),
        EnergyLimitFeedforward(
            component_type=GenericBattery,
            source=EnergyVariable,
            affected_values=[EnergyVariable],
            number_of_periods=10,
        ),
    ]

    for ff in ffs
        for av in PSI.get_affected_values(ff)
            @test isa(av, PSI.VariableKey)
        end
    end

    ff = FixValueFeedforward(
        component_type=HydroDispatch,
        source=OnVariable,
        affected_values=[OnStatusParameter],
    )

    for av in PSI.get_affected_values(ff)
        @test isa(av, PSI.ParameterKey)
    end

    @test_throws ErrorException EnergyLimitFeedforward(
        component_type=GenericBattery,
        source=EnergyVariable,
        affected_values=[OnStatusParameter],
        number_of_periods=10,
    )

    @test_throws ErrorException UpperBoundFeedforward(
        component_type=RenewableDispatch,
        source=ActivePowerVariable,
        affected_values=[OnStatusParameter],
    )

    @test_throws ErrorException LowerBoundFeedforward(
        component_type=RenewableDispatch,
        source=ActivePowerVariable,
        affected_values=[OnStatusParameter],
    )

    @test_throws ErrorException SemiContinuousFeedforward(
        component_type=ThermalMultiStart,
        source=OnVariable,
        affected_values=[ActivePowerVariable, OnStatusParameter],
    )
end


text:
using ForwardDiff

#thresh = 1e-12

@testset "Derivative checks" begin

    ###################################
    # Test Jacobians
    ###################################

    @testset "Jacobian checks" begin

        # Quaternion to rotation matrix
        @testset "Jacobian (Quat -> RotMatrix)" begin
            for i = 1:10    # do some repeats
                q = rand(Quat)  # a random quaternion

                # test jacobian to a Rotation matrix
                R_jac = Rotations.jacobian(RotMatrix, q)
                FD_jac = ForwardDiff.jacobian(x -> SVector{9}(Quat(x[1], x[2], x[3], x[4])),
                                              SVector(q.w, q.x, q.y, q.z))

                # compare
                @test FD_jac ≈ R_jac
            end
        end

        # SPQuat to Quat
        @testset "Jacobian (SPQuat -> Quat)" begin
            for i = 1:10    # do some repeats
                spq = rand(SPQuat)  # a random SPQUat
                # TODO: should test the full domain of SPQuats, not just the one with ||.|| < 1

                # test jacobian to a Rotation matrix
                R_jac = Rotations.jacobian(Quat, spq)
                FD_jac = ForwardDiff.jacobian(x -> (q = Quat(SPQuat(x[1],x[2],x[3]));
                                                    SVector(q.w, q.x, q.y, q.z)),
                                              SVector(spq.x, spq.y, spq.z))

                # compare
                @test FD_jac ≈ R_jac
            end
        end

        # SPQuat to Quat
        @testset "Jacobian (Quat -> SPQuat)" begin
            for i = 1:10    # do some repeats
                q = rand(Quat)  # a random Quat

                # test jacobian to a Rotation matrix
                R_jac = Rotations.jacobian(SPQuat, q)
                FD_jac = ForwardDiff.jacobian(x -> (spq = SPQuat(Quat(x[1],x[2],x[3],x[4]));
                                                    SVector(spq.x, spq.y, spq.z)),
                                              SVector(q.w, q.x, q.y, q.z))

                # compare
                @test FD_jac ≈ R_jac
            end
        end

        # SPQuat to rotation matrix
        @testset "Jacobian (SPQuat -> RotMatrix)" begin
            for i = 1:10    # do some repeats
                spq = rand(SPQuat)  # a random quaternion

                # test jacobian to a Rotation matrix
                R_jac = Rotations.jacobian(RotMatrix, spq)
                FD_jac = ForwardDiff.jacobian(x -> SVector{9}(SPQuat(x[1], x[2], x[3])),
                                              SVector(spq.x, spq.y, spq.z))

                # compare
                @test FD_jac ≈ R_jac
            end
        end

        ##############################
        # Jacobians for rotating stuff
        ##############################
#=
        # Quaternion multiplication
        @testset "Jacobian (Quaternion muliplication w.r.t. the right quaternion)" begin

            for i = 1:10    # do some repeats

                ql = quatrand()    # a random quaternion (should work for non-unit quaternions)
                qr = quatrand()    # a random quaternion (should work for non-unit quaternions)

                # transformation function
                J(X) = vec(ql * Quaternion(X...))

                # test jacobian to a Rotation matrix
                R_jac = Rotations.jacobian(ql, qr)
                FD_jac = ForwardDiff.jacobian(J, vec(qr))

                # compare
                @test all(abs(Matrix(R_jac) - FD_jac) .< thresh)

            end
        end

        @testset "Jacobian (Quaternion muliplication w.r.t. the left quaternion)" begin

            for i = 1:10    # do some repeats

                ql = quatrand()    # a random quaternion (should work for non-unit quaternions)
                qr = quatrand()    # a random quaternion (should work for non-unit quaternions)

                # transformation function
                J(X) = vec(Quaternion(X...) * qr)

                # test jacobian to a Rotation matrix
                R_jac = Rotations.jacobian(ql, qr, Val{false})
                FD_jac = ForwardDiff.jacobian(J, vec(ql))

                # compare
                @test all(abs(Matrix(R_jac) - FD_jac) .< thresh)

            end
        end =#

        # rotate a point by a RotMatrix
        @testset "Jacobian (RotMatrix rotation)" begin
            for i = 1:10    # do some repeats
                r = rand(RotMatrix{Float64})    # a random quaternion
                v = randn(SVector{3,Float64})

                # test jacobian to a Rotation matrix
                R_jac = Rotations.jacobian(r, v)
                FD_jac = ForwardDiff.jacobian(x -> RotMatrix(x)*v,
                                              SVector(r))

                # compare
                @test FD_jac ≈ R_jac
            end
        end

        # rotate a point by a quaternion
        @testset "Jacobian (Quat rotation)" begin
            for i = 1:10    # do some repeats
                q = rand(Quat)    # a random quaternion
                v = randn(SVector{3,Float64})

                # test jacobian to a Rotation matrix
                R_jac = Rotations.jacobian(q, v)
                FD_jac = ForwardDiff.jacobian(x -> Quat(x[1], x[2], x[3], x[4])*v,
                                              SVector(q.w, q.x, q.y, q.z))

                # compare
                @test FD_jac ≈ R_jac
            end
        end

        # rotate a point by a SPQuat
        @testset "Jacobian (SPQuat rotation)" begin
            for i = 1:10    # do some repeats
                spq = rand(SPQuat)    # a random quaternion
                v = randn(SVector{3,Float64})

                # test jacobian to a Rotation matrix
                R_jac = Rotations.jacobian(spq, v)
                FD_jac = ForwardDiff.jacobian(x -> SPQuat(x[1], x[2], x[3])*v,
                                              SVector(spq.x, spq.y, spq.z))

                # compare
                @test FD_jac ≈ R_jac
            end
        end
#=
        # rotate a point by an SpQuat
        @testset "Jacobian (SpQuat rotation)" begin

            for i = 1:10    # do some repeats

                spq = SpQuat(nquatrand())    # a random quaternion
                X = randn(Vec{3,Float64})

                # transformation function
                J(spq) = Vector(rotate(SpQuat(spq...), X))

                # test jacobian to a Rotation matrix
                R_jac = Rotations.jacobian(spq, X)
                FD_jac = ForwardDiff.jacobian(J, vec(spq))

                # compare
                @test all(abs(Matrix(R_jac) - FD_jac) .< thresh)

            end
        end
=#

    end

    ###################################
    # Test Hessians
    ###################################
    #=
    @testset "Hessian checks" begin

        # SpQuat to Quaternion
        @testset "Hessian (SpQuat -> Quaternion)" begin

            for i = 1:10    # do some repeats

                q = nquatrand()                           # a random quaternion
                spq = Rotations.quat_to_spquat_naive(q)   # I want to test the full domain of SpQuats, not just the one with ||.|| < 1

                # test jacobian to a Rotation matrix
                R_hess = Rotations.hessian(Quaternion, spq)
                for d = 1:4
                    H(X) = vec(Quaternion(SpQuat(X...)))[d]   # transformation function
                    FD_hess = ForwardDiff.hessian(H, vec(spq))

                    # compare
                    @test all(abs(Matrix(R_hess[d]) - FD_hess) .< thresh)
                end
            end
        end

        # Quaternion to SpQuat
        @testset "Hessian (Quaternion -> SpQuat)" begin

            for i = 1:10    # do some repeats

                q = nquatrand()                           # a random quaternion

                # test jacobian to a Rotation matrix
                R_hess = Rotations.hessian(SpQuat, q)
                for d = 1:3
                    H(X) = vec(SpQuat(Quaternion(X...)))[d]   # transformation function
                    FD_hess = ForwardDiff.hessian(H, vec(q))

                    # compare
                    @test all(abs(Matrix(R_hess[d]) - FD_hess) .< thresh)
                end
            end
        end

        @testset "Hessian (Quaternion rotation)" begin

            for i = 1:10   # do some repeats

                q = nquatrand()    # a random quaternion
                X = randn(Vec{3,Float64})

                # transformation function
                R_hess = Rotations.hessian(q, X)

                for d = 1:3
                    H(q) = Vector(rotate(Quaternion(q...), X))[d]

                    # test
                    FD_hess = ForwardDiff.hessian(H, vec(q))

                    # compare
                    @test all(abs(Matrix(R_hess[d]) - FD_hess) .< thresh)
                end
            end
        end

        @testset "Hessian (SpQuat rotation)" begin

            for i = 1:10    # do some repeats

                q = nquatrand()    # a random quaternion
                spq = SpQuat(q)
                X = randn(Vec{3,Float64})

                # transformation function
                R_hess = Rotations.hessian(spq, X)

                for d = 1:3
                    H(spq) = Vector(rotate(SpQuat(spq...), X))[d]

                    # test
                    FD_hess = ForwardDiff.hessian(H, vec(spq))

                    # compare
                    @test all(abs(Matrix(R_hess[d]) - FD_hess) .< thresh)
                end
            end
        end
    end =#
end


text:
include("air.jl")
UPPE.permittivity(w) = permittivity(w)
UPPE.permeability(w) = permeability(w)

grid = GridT(
    CPU();
    tu = 1e-15,   # [s] unit of time
    tmin = -200,   # [tu] left boundary of t grid
    tmax = 200,   # [tu] right boundary of t grid
    Nt = 2^14,   # number of points on t grid
    tguard = 20,   # [tu] the width of the lossy slab at the ends of t grid
    wguard = 1e16,   # [1/s] the cut-off angular frequency
)

lam0 = 0.8e-6   # [m] central wavelength
tau0 = 35e-15 / (2 * sqrt(log(2)))  # [s] initial pulse duration
I0 = 1e12 * 1e4   # [W/cm^2] initial intensity

w0 = 2 * pi * C0 / lam0   # central frequency
n0 = refractive_index(w0)
Iu = I0   # [W/m**2] unit of intensity
Eu = intensity2field(Iu, n0)

E = @. sqrt(I0 / Iu) *
    exp(-0.5 * (grid.t * grid.tu)^2 / tau0^2) * cos(w0 * grid.t * grid.tu)

field = Field(grid, w0, E; Eu, Iu)

ne, Pnl, J = zeros(size(field.E)), zero(field.E), zero(field.E)
response =  Response(ne, Pnl, J, nothing, nothing)

model = Model(
    grid, field, response;
    zu = 1,   # [m] unit of space in z direction
    z = 0,   # [m] initial propagation distance
    zmax = 21,   # [zu] propagation distance
    kparaxial = true,
    nonlinearity = false,
)

simulation = Simulation(
    model, n2, N0;
    prefix = "results/",
    dz0 = 1,   # [zu] initial z step
    dzout = model.zmax,   # [zu] z step for writing field into output file
)

run!(simulation)

# Compare to analytic solution:
(; tu, Nt, t) = grid

zmax = model.zmax * model.zu
k0 = k_func(w0)
vg = group_velocity(w0)
zdisp = dispersion_length(w0, tau0)
taud = tau0 * sqrt(1 - 1im * zmax / zdisp)

Ea = @. tau0 / taud *
     exp(-0.5 * (t * tu)^2 / taud^2) *
     exp(1im * k0 * zmax - 1im * w0 * (t * tu + zmax / vg))

@test isapprox(field.E, Ea; rtol=1e-2)