text:
module CoNLL

function read(path::String, columns::Int...)
    doc = []
    sent = []
    lines = open(readlines, path)
    for line in lines
        line = chomp(line)
        if length(line) == 0
            length(sent) > 0 && push!(doc, sent)
            sent = []
        else
            items = split(line, '\t')
            data = map(c -> String(items[c]), columns)
            push!(sent, data)
        end
    end
    length(sent) > 0 && push!(doc, sent)
    doc
end

end

max_stars_repo_name: hshindo/JukaiNLP

max_issues_repo_name: hshindo/JukaiNLP

max_forks_repo_name: hshindo/JukaiNLP


text:
using RoME: MeanMaxPPE
#=
Load a DEM as a Point3 set in a factor graph
=#

using Images
using FileIO
using Caesar
using RoME
using ProgressMeter

using TensorCast
using Gadfly
Gadfly.set_default_plot_size(35cm, 25cm)

##

prjPath = dirname(dirname(pathof(Caesar)))

# load dem (18x18km span, ~17m/px)
img = load(joinpath(prjPath, "examples","dev","scalar","dem.png")) .|> Gray
img = Float64.(img)
img = img[1:100, 1:100]

x = range(-9e3, 9e3, length = size(img)[1]) # North
y = range(-9e3, 9e3, length = size(img)[2]) # East

# point uncertainty - 2.5m horizontal, 1m vertical
# horizontal uncertainty chosen so that 3sigma is approx half the resolution
sigma = diagm([2.5, 2.5, 1.0])

function loadDEM!(  fg::AbstractDFG,
                    dem::Matrix, # assume grayscale image for now
                    x::AbstractVector,
                    y::AbstractVector;
                    solvable::Int=0,
                    marginalized::Bool=true,
                    meshEdgeSigma=diagm([1;1;1]),
                    refKey::Symbol = :simulated )
    #
    # assume regular grid
    dx, dy = x[2]-x[1], y[2]-y[1]
    @showprogress 1 "loading rows..." for i in 1:length(x)
        for j in 1:length(y)
            s = Symbol("pt$(i)_$(j)") # unique identifier
            v_ij = addVariable!(fg, s, Point3, solvable=solvable)
            setMarginalized!(v_ij, marginalized) # assume solveKey=:default
            refVal = [x[i];y[j];dem[i,j]]
            simPPE = DFG.MeanMaxPPE(refKey, refVal, refVal, refVal)
            setPPE!(v_ij, refKey, typeof(simPPE), simPPE)

            # Regular grid triangulation
            # add factor to (i-1,j)
            # add factor to (i, j-1)
            # add factor to (i-1, j-1)

            # skip first corner
            # i == 1 && j == 1 ? continue : @warn("skip corner")
            # eg i=2, j=1
            # @show i,j
            if i>1
                dVal1 = dem[i,j]-dem[i-1,j]
                f = Point3Point3(MvNormal([dx, 0, dVal1], meshEdgeSigma))
                addFactor!(fg, [Symbol("pt$(i-1)_$(j)"), s], f, solvable=solvable, graphinit=false )
            end

            # skip out on first column
            # j == 1 ? continue : nothing
            if j>1
                dVal2 = dem[i,j]-dem[i,j-1]
                f = Point3Point3(MvNormal([0, dy, dVal2], meshEdgeSigma))
                addFactor!(fg, [Symbol("pt$(i)_$(j-1)"),s], f, solvable=solvable, graphinit=false )
            end

            # skip out on first row
            # i == 1 ? continue : nothing

            if i>1 && j>1
                dVal3 = dem[i,j]-dem[i-1,j-1]
                f = Point3Point3(MvNormal([dx, dy, dVal3], meshEdgeSigma))
                addFactor!(fg,[Symbol("pt$(i-1)_$(j-1)"),s], f, solvable=solvable, graphinit=false )
            end

        end
    end

    nothing
end

##

# function getTM(dem::Matrix)
#     # get z transition samples over x and y
#     # NOTE: this works over a regular grid, so take care when using
#     # with non-uniformly spaced points
#     s = size(dem)
#     dc = img[2:(s[1]-1), 2:(s[2]-1)]

#     dx1 = dc - img[1:(s[1]-2), 2:(s[2]-1)]
#     dx2 = dc - img[3:s[1],     2:(s[2]-1)]
#     dy1 = dc - img[2:(s[1]-1), 1:(s[1]-2)]
#     dy2 = dc - img[2:(s[1]-1), 3:s[2]]
#     [[dx1[:];dx2[:]],[dy1[:];dy2[:]]]
# end

fg = initfg()

##

loadDEM!(fg, img, (x), (y));

##

#callback
cb(g, lbl) = @show lbl

generateCanonicalFG_Helix2DSlew!(100, posesperturn=30, radius=500, dfg=fg, graphinit=false, slew_y=2/3, postpose_cb=cb)

##

getPPE(fg, :x0, :simulated).suggested

getPPE(fg, :pt50_50, :simulated).suggested

## fetch and plot pose variables

traj = ls(fg, r"x\d+") |> sortDFG

traj_xy = traj .|> x->getPPE(fg, x, :simulated).suggested[1:2]

@cast _xy[i,j] := traj_xy[i][j]

plot(
    layer(x=_xy[:,1],y=_xy[:,2], Geom.point),
    layer(x=_xy[:,1],y=_xy[:,2], Geom.path)
)

##

ext: jl


text:
col_altro = "orange"
col_ipopt = "blue"
col_snopt = "green!80!black"
col_alilqr = "cyan"

function trajectory_plot(prob::Problem; kwargs...)
    x = [x[1] for x in prob.X]
    y = [x[2] for x in prob.X]
    PGF.Plots.Linear(x,y; kwargs...)
end

function trajectory_plot_flip(prob::Problem; kwargs...)
    x = [x[2] for x in prob.X]
    y = [x[1] for x in prob.X]
    PGF.Plots.Linear(x,y; kwargs...)
end

max_issues_repo_name: MehrdadZareian/GMARA

max_forks_repo_name: MehrdadZareian/GMARA


text:
module Day9
export solvePart1, solvePart2

numberOfPlayers = 418
lastMarblePoints = 71339

mutable struct Marble
    next
    prev
    number::Int64
end

function indexCounterClockwise(ind, numPlaces, len)
    i = ind - numPlaces
    while i < 0
        i += len
    end
    return i
end

function indexClockwise(ind, numPlaces, len)
    i = ind + numPlaces
    while i > len
        i -= len
    end
    return i
end

function removeMarbleSevenPlacesCounterClockwise(current::Marble)
    toRemove = current
    for i in 1:7
        toRemove = toRemove.prev
    end
    before = toRemove.prev
    after = toRemove.next
    before.next = after
    after.prev = before
    return toRemove
end

function winningElfScore(numberOfPlayers, lastNumberScore)
    # Start with a marble whose previous and next is itself
    firstMarble = Marble(nothing, nothing, 0)
    firstMarble.next = firstMarble
    firstMarble.prev = firstMarble

    whoseTurn = 1
    points = Dict()
    current = firstMarble
    for nextMarbleNumber in 1:lastNumberScore
        if nextMarbleNumber % 23 == 0
            removed = removeMarbleSevenPlacesCounterClockwise(current)
            toGain = nextMarbleNumber + removed.number
            points[whoseTurn] = get(points, whoseTurn, 0) + toGain
            current = removed.next
        else
            before = current.next
            after = before.next
            toInsert = Marble(before, after, nextMarbleNumber)
            before.next = toInsert
            after.prev = toInsert
            toInsert.prev = before
            toInsert.next = after
            current = toInsert
        end

        whoseTurn += 1
        if whoseTurn > numberOfPlayers
            whoseTurn = 1
        end
    end
    return max(values(points)...)
end

function solvePart1()
    winningElfScore(numberOfPlayers, lastMarblePoints)
end

function solvePart2()
    winningElfScore(numberOfPlayers, lastMarblePoints*100)
end
end


text:
md"""
# Axes

We define two singleton types LengthAxis and WidthAxis with the common supertype Axis.

Note that the instances of these Axis types are LengthAxis() and WidthAxis().

Measurements are made along an axis.

Cuts are made across an axis.
"""

abstract type Axis end

struct LengthAxis <: Axis end
struct WidthAxis <: Axis end

"""
    other(axis)
returns the axis that is perpendicular to the specified axis.
"""
function other end

function other(axis::LengthAxis)::Axis
    WidthAxis()
end

function other(axis::WidthAxis)::Axis
    LengthAxis()
end

@assert LengthAxis() isa Axis
@assert WidthAxis() isa Axis
@assert other(WidthAxis()) == LengthAxis()
@assert other(LengthAxis()) == WidthAxis()

pretty(axis::Axis) = typeof(axis)

unit_vector(::LengthAxis) = [1 0]
unit_vector(::WidthAxis) = [0 1]

const zerovector = [0u"inch" 0u"inch"]

area(v::Array{<:Number, 2}) = v[1] * v[2]

@assert area([2 3]) == 6

const zeroarea = area(zerovector)

export Axis, LengthAxis, WidthAxis, other
export unit_vector, zerovector, zeroarea


text:
# Implements conversion from some high-level Flux ops to their equivalent nGraph
# representations when running under Cassette.
#include("lib.jl")

# Some tuple utilities
astuple(x::Tuple) = x
astuple(x) = (x,)

untuple(x::Tuple) = x
untuple(x::Tuple{T}) where {T} = first(x)

# NOTE: Register was required in the past for some nodes like `BatchNorm` that left
# dangling nodes in the computation graph that then led ngraph to segfault
#
# TODO: Does this still happen?
"""
    __register(x::Node) -> Nothing

By default, becomes a no-op. When executed under [`TraceCtx`](@ref), registers `x` as a
hidden output of a ngraph function.
"""
__register(x::Node) = nothing

#####
##### Cassette Magic
#####

Cassette.@context TraceCtx

struct TraceMetadata
    training::Bool

    # Mapping of a Parameter array to nGraph Node
    array_to_node::IdDict{Any, Node}

    # Keep track of the constants we have created.
    constants::IdDict{Any,Node}

    # Implicitly created outputs
    implicit_outputs::Vector{Node}
end

function TraceMetadata(array_to_node, training::Bool)
    return TraceMetadata(
        training,
        array_to_node,
        IdDict{Any,Node}(),
        Node[],
    )
end

# Flag to indicate if we are training or not.
#
# If so, we need to add another pass to compute the gradients of all parameters.
istraining(x::TraceCtx) = x.metadata.training

#####
##### Cassette Overdubs
#####

# Hijack Node constructors from Arrays
function Cassette.overdub(ctx::TraceCtx, f::Type{Node{T,N}}, x::AbstractArray) where {T,N}
    # Check if this array is a parameter.
    # If so, get our cached input for it.
    node = get(ctx.metadata.array_to_node, x, nothing)
    if !isnothing(node)
        @assert isa(node, Node{T,N})
        return node
    end

    # Check if we've already made a constant for this object.
    node = get(ctx.metadata.constants, x, nothing)
    if isnothing(node)
        node = constant(x)
        ctx.metadata.constants[x] = node
    end
    return node
end

# Do not hijack creating nodes from nodes.
Cassette.overdub(ctx::TraceCtx, f::Type{<:Node}, x::Node) = f(x)

# Register implicit outputs
# I'm looking at YOU batch-norm training!
function Cassette.overdub(ctx::TraceCtx, ::typeof(__register), x::Node)
    return push!(ctx.metadata.secondary, x)
end

# Get around Cassette bug with `reverse`
Cassette.overdub(::TraceCtx, ::typeof(reverse), args...) = reverse(args...)

# # Hijack these layers
# Cassette.overdub(ctx::TraceCtx, f::Flux.Dense, args...) =
#     Cassette.overdub(ctx, _dense_impl, f, args...)
#
# Cassette.overdub(ctx::TraceCtx, f::Flux.Conv, args...) =
#     Cassette.overdub(ctx, _conv_impl, f, args...)
#
# Cassette.overdub(ctx::TraceCtx, f::Flux.CrossCor, args...) =
#     Cassette.overdub(ctx, _conv_impl, f, args...)
#
# Cassette.overdub(ctx::TraceCtx, f::Flux.BatchNorm, args...) =
#     Cassette.overdub(ctx, _batchnorm_impl, f, args...)

# Skip recursing initialization calls - recursing turns out to take a very, very long time.
Cassette.overdub(ctx::TraceCtx, f::typeof(rand), args...) = f(args...)
#Cassette.overdub(ctx::TraceCtx, f::typeof(Flux.glorot_normal), args...) = f(args...)
#Cassette.overdub(ctx::TraceCtx, f::typeof(Flux.glorot_uniform), args...) = f(args...)

#####
##### Main `compile` entrypoint
#####

"""
    compile(backend, f, x...; [training]) -> Executable

Trace and compile a Flux model `f` with `args`.
"""
function compile(backend::Backend, f, x...; kw...)
    # Build the nGraph computation graph
    trace, _ = snoop(f, x...; kw...)

    # pass the trace as well as the optimizer arguments to create the executable.
    return make(backend, trace; kw...)
end

"""
Traced parameters from a compiled function: `f(args...)`.

* `inputs`: The nodes that are explicitly used as an input to `f`. Corresponds to
    `args...` in the signature of `f`.
* `outputs`: The nodes that represent the results of `f(args...)`.
* `implicit_outputs::Vector{Node}`: Outputs that were implicitly created during tracing -
    usually the result of BatchNorms being wierd.

"""
struct Trace
    inputs::Vector{Node}
    outputs::Vector{Node}
    implicit_outputs::Vector{Node}

    # Record of what is a parameter
    parameters::Flux.Params
    parameter_nodes::Vector{Node}

    array_to_node::IdDict{Any, Node}
    node_to_array::IdDict{Node, Any}

    args::Vector
end

# Step 1 of compilation - trace the provided function with the provided arguments to
# construct an nGraph graph - apply the requested optimizer to finish graph construction.
#
# Returns an argument named tuple that is given to the next stage of compilation.
snoop(f, x...; kw...) = snoop(f, Flux.Params(), x...; kw...)
function snoop(f, parameters::Flux.Params, x...; training = false, kw...)
    # Extract the inputs and other parameters
    inputs = collect(Node.(x))

    parameter_nodes = isempty(parameters) ? Node[] : Node.(collect(parameters))

    array_to_node = IdDict{Any,Node}()
    node_to_array = IdDict{Node,Any}()

    # Construct parameter-to-node mappings.
    for (node, array) in zip(parameter_nodes, parameters)
        array_to_node[array] = node
        node_to_array[node] = array
    end

    metadata = TraceMetadata(array_to_node, training)
    ctx = TraceCtx(metadata = metadata)

    # Perform traced execution on the function.
    outputs = collect(astuple(Cassette.overdub(ctx, f, inputs...)))
    @assert all(x -> isa(x, Node), outputs)

    trace = Trace(
        inputs,
        outputs,
        ctx.metadata.implicit_outputs,
        parameters,
        parameter_nodes,
        array_to_node,
        node_to_array,
        collect(x),
    )

    return (trace, metadata)
end

function make(backend::Backend, trace; training = false, kw...)
    # TODO: Reimplement training

    # Create an nGraph Executable
    allinputs = [trace.inputs; trace.parameter_nodes]
    alloutputs = [trace.outputs; trace.implicit_outputs]

    ex = compile(
        backend,
        allinputs,
        alloutputs,
        kw...
    )

    # Create TensorViews for each of the inputs and outputs
    input_tensors = Tensor.(Ref(backend), trace.args)
    output_tensors = Tuple(Tensor.(Ref(backend), trace.outputs))
    implicit_output_tensors = Tensor.(Ref(backend), trace.implicit_outputs)

    return CallableFunction(
        ex,
        input_tensors,
        output_tensors,
        implicit_output_tensors,
    )
end

#####
##### Executable
#####

struct CallableFunction{O}
    ex::Executable
    inputs::Vector{Tensor}
    outputs::NTuple{O,Tensor}
    implicit_outputs::Vector{Tensor}
end

allinputs(x::CallableFunction) = x.inputs
alloutputs(x::CallableFunction) = collect(Iterators.flatten((x.outputs, x.implicit_outputs)))

function (ex::CallableFunction)()
    inputs = allinputs(ex)
    outputs = alloutputs(ex)

    # Since we're passing wrapped type to C++, we have to cast them to Any's which is
    # kind of gross - but w/e
    ex.ex(inputs, outputs)
    return untuple(ex.outputs)
end


text:
# haskey
# in
# summary
# show
# _truncate_at_width_or_chars: Not exported
# showdict: Not exported
# showkv: Not exported
# length: interfaces.jl
# isempty: interfaces.jl
# _tt1: Not exported
# _tt2: Not exported
# eltype: interfaces.jl
# start: interfaces.jl
# done: interfaces.jl
# next: TODO
# keys: TODO
# eachindex: interfaces.jl
# values: TODO
# copy: interfaces.jl
# merge!: overwriting
# keytype
# valtype
# merge
# filter!: overwriting
# filter: interfaces.jl
# isequal: interfaces.jl
# ==: interfaces.jl
# hash: works for Node
# _tablesz: Not exported
# getindex: interfaces.jl
# setindex!: overwriting
# push!: overwriting
# ObjectIdDict
# similar: interfaces.jl
# get: TODO
# pop!: overwriting
# delete!: overwriting
# empty!: overwriting
# _oidd_nextind: Not exported
# get!: overwriting
# SerializationState
# Dict
# dict_with_eltype: Not exported
# convert
# hashindex: Not exported
# isslotempty: Not exported
# isslotfilled: Not exported
# isslotmissing: Not exported
# rehash!: Not exported
# sizehint!
# ht_keyindex: Not exported
# ht_keyindex2: Not exported
# _setindex!: Not exported
# getkey
# _pop!: Not exported
# _delete!: Not exported
# skip_deleted: Not exported
# WeakKeyDict
# weak_key_delete!: Not exported


text:
using Dolo

model = Dolo.yaml_import("examples/models/rbc.yaml")


text:
using AbstractGPs
using TemporalGPs

# Doing approximate inference in space-time GPs using pseudo-points is, as always, similar
# to performing exact inference. The primary difference is that `rand` isn't available,
# and the `approx_posterior_marginals` function must be used to query the approximate
# posterior. This is hopefully a temporary solution, so should change at some point in the
# future.

# Load up the separable kernel from TemporalGPs. You need to use this to tell TemporalGPs
# that you're using a separable kernel (it's not enough just to use a kernel which
# happens to be separable).
# RegularInTime is a data structure for inputs which allows for different spatial locations
# at each point in time, and can be used with the approximate inference scheme presented
# here.
using TemporalGPs: Separable, RegularInTime, approx_posterior_marginals

# Specify a separable kernel.
# The first argument is always the kernel over space, the second the kernel over time.
# You can also use sums of separble kernels.
k = Separable(SEKernel(), Matern52Kernel());

# Build a GP, and convert it to an SDE as per usual.
# Use `ArrayStorage`, not `SArrayStorage`, for these kinds of GPs.
f = to_sde(GP(k), ArrayStorage(Float64));

# Construct inputs. Spatial locations change at each point in time.
# Also works with RectilinearGrids of inputs.
# Times must be increasing, points in space can be anywhere.
N = 50;
T = 1_000;
points_in_space = [randn(N) for _ in 1:T];
points_in_time = RegularSpacing(0.0, 0.1, T);
x = RegularInTime(points_in_time, points_in_space);

# Since it's not straightforward to generate samples from this GP at `x`, use a known
# function, under a bit of iid noise.
xs = collect(x);
y = sin.(first.(xs)) .+ cos.(last.(xs)) + sqrt.(0.1) .* randn(length(xs));

# Spatial pseudo-point locations.
z_r = range(-3.0, 3.0; length=5);

# Locations in space at which to make predictions. Assumed to be the same at each point in
# time, but this assumption could easily be relaxed.
N_pr = 150;
x_r_pr = range(-5.0, 5.0; length=N_pr);

# Compute the approximate posterior marginals.
f_post_marginals = approx_posterior_marginals(dtc, f(x, 0.1), y, z_r, x_r_pr);
m_post_marginals = mean.(f_post_marginals);
σ_post_marginals = std.(f_post_marginals);

# Visualise the posterior marginals. We don't do this during in CI because it causes
# problems.
if get(ENV, "TESTING", "FALSE") == "FALSE"
    using Plots
    savefig(
        plot(
            heatmap(reshape(m_post_marginals, N_pr, T)),
            heatmap(reshape(σ_post_marginals, N_pr, T));
            layout=(1, 2),
        ),
        "posterior.png",
    );
end


text:
@testset "GQLLocation" begin
    str = """{"line": 1, "column": 10}"""
    loc = JSON3.read(str, GraphQLClient.GQLLocation)
    @test loc.line == 1
    @test loc.column == 10
    @test sprint(show, loc) == "Line 1 Column 10"
end

@testset "GQLError" begin
    str = """{"message": "Message text"}"""
    err = JSON3.read(str, GraphQLClient.GQLError)
    @test err.message == "Message text"
    @test isnothing(err.locations)
    str = """{"message": "Message text", "locations":[{"line": 1, "column": 10}]}"""
    err = JSON3.read(str, GraphQLClient.GQLError)
    @test err.message == "Message text"
    @test !isempty(err.locations)
end

@testset "GQLResponse" begin
    # Successful query, no errors, default parametric type
    response="""
    {
        "data": {
            "queryName1": {
                "fieldName": "fieldValue"
            },
            "queryName2": {
                "fieldName": 0
            }
        }
    }
    """
    resp = JSON3.read(response, GraphQLClient.GQLResponse{Any})
    @inferred JSON3.read(response, GraphQLClient.GQLResponse{Any})
    @test isnothing(resp.errors)
    @test !isempty(resp.data)
    @test haskey(resp.data, "queryName1")
    @test haskey(resp.data, "queryName2")

    # Unuccessful query, error before execution, default parametric type
    response = """
    {
        "errors": [
            {
                "message": "Failed before execution",
                "locations": [
                    {"line": 1, "column": 10},
                    {"line": 2, "column": 12}
                ]
            },
            {
                "message": "No locations"
            }
        ]
    }
    """
    resp = JSON3.read(response, GraphQLClient.GQLResponse{Any})
    @test length(resp.errors) == 2
    @test isnothing(resp.data)

    # Unuccessful query, error during execution, default parametric type
    response = """
    {
        "errors": [
            {
                "message": "Failed before execution",
                "locations": [
                    {"line": 1, "column": 10},
                    {"line": 2, "column": 12}
                ]
            },
            {
                "message": "No locations"
            }
        ],
        "data": {
            "queryName": null
        }
    }
    """
    resp = JSON3.read(response, GraphQLClient.GQLResponse{Any})
    @test length(resp.errors) == 2
    @test !isempty(resp.data)
    @test isnothing(resp.data["queryName"])

    # Custom type
    response = """
    {
        "data": {
            "queryName": {
                "fieldName1": "fieldValue",
                "fieldName2": 1.0
            }
        }
    }"""
    struct QueryNameResponse
        fieldName1::String
        fieldName2::Float64
    end
    StructTypes.StructType(::Type{QueryNameResponse}) = StructTypes.OrderedStruct()
    resp = JSON3.read(response, GraphQLClient.GQLResponse{QueryNameResponse})
    @inferred JSON3.read(response, GraphQLClient.GQLResponse{QueryNameResponse})
    @test resp.data["queryName"] isa QueryNameResponse

    # failed custom type
    response = """{"data": {"queryName": null}}"""
    resp = JSON3.read(response, GraphQLClient.GQLResponse{QueryNameResponse})
    @test isnothing(resp.data["queryName"])
end


text:
# See interface docs
@inline inbounds(data::Union{GridData,AbstractSimData}, I::Tuple) = inbounds(data, I...)
@inline inbounds(data::Union{GridData,AbstractSimData}, I...) =
    _inbounds(boundary(data), gridsize(data), I...)

@inline function _inbounds(boundary::BoundaryCondition, size::Tuple, i1, i2)
    a, inbounds_a = _inbounds(boundary, size[1], i1)
    b, inbounds_b = _inbounds(boundary, size[2], i2)
    (a, b), inbounds_a & inbounds_b
end
@inline _inbounds(::Remove, size::Number, i::Number) = i, _isinbounds(size, i)
@inline function _inbounds(::Wrap, size::Number, i::Number)
    if i < oneunit(i)
        size + rem(i, size), true
    elseif i > size
        rem(i, size), true
    else
        i, true
    end
end

# See interface docs
@inline isinbounds(data::Union{GridData,AbstractSimData}, I::Tuple) = isinbounds(data, I...)
@inline isinbounds(data::Union{GridData,AbstractSimData}, I...) = _isinbounds(gridsize(data), I...)

@inline _isinbounds(size::Tuple, I...) = all(map(_isinbounds, size, I))
@inline _isinbounds(size, i) = i >= one(i) && i <= size

# _updateboundary
# Reset or wrap boundary where required. This allows us to ignore
# bounds checks on neighborhoods and still use a wraparound grid.
_updateboundary!(grids::Tuple) = map(_updateboundary!, grids)
function _updateboundary!(g::GridData{S,R}) where {S<:Tuple{Y,X},R} where {Y,X}
    R < 1 && return g
    return _updateboundary!(g, boundary(g))
end
function _updateboundary!(g::GridData{S,R}, ::Remove) where {S<:Tuple{Y,X},R} where {Y,X}
    src = parent(source(g))
    # Left
    @inbounds src[1:Y, 1:R] .= Ref(padval(g))
    # Right
    @inbounds src[1:Y, X+R+1:X+2R] .= Ref(padval(g))
    # Top middle
    @inbounds src[1:R, R+1:X+R] .= Ref(padval(g))
    # Bottom middle
    @inbounds src[Y+R+1:Y+2R, R+1:X+R] .= Ref(padval(g))
    status = sourcestatus(g)
    return g
end
function _updateboundary!(g::GridData{S,R}, ::Wrap) where {S<:Tuple{Y,X},R} where {Y,X}
    src = parent(source(g))
    nrows, ncols = gridsize(g)
    startpadrow = startpadcol = 1:R
    endpadrow = nrows+R+1:nrows+2R
    endpadcol = ncols+R+1:ncols+2R
    startrow = startcol = 1+R:2R
    endrow = nrows+1:nrows+R
    endcol = ncols+1:ncols+R
    rows = 1+R:nrows+R
    cols = 1+R:ncols+R

    # Sides ---
    # Left
    @inbounds copyto!(src, CartesianIndices((rows, startpadcol)),
                      src, CartesianIndices((rows, endcol)))
    # Right
    @inbounds copyto!(src, CartesianIndices((rows, endpadcol)),
                      src, CartesianIndices((rows, startcol)))
    # Top
    @inbounds copyto!(src, CartesianIndices((startpadrow, cols)),
                      src, CartesianIndices((endrow, cols)))
    # Bottom
    @inbounds copyto!(src, CartesianIndices((endpadrow, cols)),
                      src, CartesianIndices((startrow, cols)))

    # Corners ---
    # Top Left
    @inbounds copyto!(src, CartesianIndices((startpadrow, startpadcol)),
                      src, CartesianIndices((endrow, endcol)))
    # Top Right
    @inbounds copyto!(src, CartesianIndices((startpadrow, endpadcol)),
                      src, CartesianIndices((endrow, startcol)))
    # Botom Left
    @inbounds copyto!(src, CartesianIndices((endpadrow, startpadcol)),
                      src, CartesianIndices((startrow, endcol)))
    # Botom Right
    @inbounds copyto!(src, CartesianIndices((endpadrow, endpadcol)),
                      src, CartesianIndices((startrow, startcol)))

    _wrapstatus!(sourcestatus(g))
    return g
end

# _wrapstatus!
# Copies status from opposite sides/corners in Wrap boundary mode
_wrapstatus!(status::Nothing) = nothing
function _wrapstatus!(status::AbstractArray)
    # !!! The end row/column is always empty !!!
    # Its padding for block opt. So we work with end-1 and end-2

    # This could be further optimised by not copying the end-2
    # block column/row when the blocks are aligned at both ends.

    # Sides
    status[1, :] .|= status[end-1, :] .| status[end-2, :]
    status[:, 1] .|= status[:, end-1] .| status[:, end-2]
    status[end-1, :] .|= status[1, :]
    status[:, end-1] .|= status[:, 1]
    status[end-2, :] .|= status[1, :]
    status[:, end-2] .|= status[:, 1]

    # Corners
    status[1, 1] |= status[end-1, end-1] | status[end-2, end-1] |
                    status[end-1, end-2] | status[end-2, end-2]
    status[end-1, 1] |= status[1, end-1] | status[1, end-2]
    status[end-2, 1] |= status[1, end-1] | status[1, end-2]
    status[1, end-1] |= status[end-1, 1] | status[end-2, 1]
    status[1, end-2] |= status[end-1, 1] | status[end-2, 1]
    status[end-1, end-1] |= status[1, 1]
    status[end-2, end-2] |= status[1, 1]
end


text:
"""
UA scheme for 4d Vlasov in Fluid-scaling with b(x)
External E
finite difference solver
with new IMEX1 and IMEX2

"""
function efd( ntau, nbpart )

    tau   = zeros(Float64, ntau)
    ltau  = zeros(Float64, ntau)
    et    = zeros(Float64, (2,ntau))

    pl    = zeros(ComplexF64, ntau)
    ql    = zeros(ComplexF64, ntau)
    xf    = zeros(ComplexF64, (ntau,2))
    yf    = zeros(ComplexF64, (ntau,2))
    gx    = zeros(ComplexF64, (ntau,2))
    gy    = zeros(ComplexF64, (ntau,2))
    temp  = zeros(ComplexF64, (ntau,2))
    tilde = zeros(ComplexF64, (ntau,2))
    h     = zeros(ComplexF64, (ntau,2))
    r     = zeros(ComplexF64, (ntau,2))
    fx    = zeros(ComplexF64, (ntau,2))
    fy    = zeros(ComplexF64, (ntau,2))
    xt    = zeros(ComplexF64, (ntau,2))
    yt    = zeros(ComplexF64, (ntau,2))

    ave  = zeros(Float64, 2)
    ave2 = zeros(Float64, 2)
    ave3 = zeros(Float64, 2)

    ftau = plan_fft(tilde, 1)

    nx = 128
    ny = 64

    alpha = 0.05
    kx    = 0.50
    ky    = 1.0
    dimx  = 2pi/kx
    dimy  = 2pi/ky
    poids = dimx * dimy
    dx    = dimx / nx
    dy    = dimy / ny

    dt     = pi/2/(2.0^3)
    tfinal = pi/2
    nstep  = trunc(Int64,tfinal/dt)
    time   = 0.
    ep     = 0.001
    dtau   = 2pi/ntau

    m    = ntau÷2
    ltau = [0:m-1; -m:-1]
    for i=1:ntau
      tau[i] = (i-1)*dtau
    end

    xmin = 0.0
    xmax = dimx
    ymin = 0.0
    ymax = dimy

    mesh = Mesh( xmin, xmax, nx, ymin, ymax, ny )
    f    = MeshFields( mesh )
    p    = read_particles( "particles.dat", mesh )

    @show sum(p.x[1,:]), sum(p.x[2,:]), sum(p.v[1,:]), sum(p.v[2,:])

    for m=1:nbpart

        x1 = p.x[1,m]
        x2 = p.x[2,m]
        v1 = p.v[1,m]
        v2 = p.v[2,m]

        time  = 0.
        bx    = 1+0.5*sin(x1)*sin(x2)
        ds    = dt*bx

        pl[1] = ds
        ql[1] = ds^2/2.0

        for i=2:ntau
            pl[i] = ep*im*(exp(-im*ltau[i]*ds/ep)-1)/ltau[i]
            ql[i] = ep*(ep*(1-exp(-im*ltau[i]*ds/ep))-im*ltau[i]*ds)/ltau[i]^2
        end

        xt[:,1] .= x1
        xt[:,2] .= x2
        yt[:,1] .= v1
        yt[:,2] .= v2

        #--preparation initial data--
        temp[1,1] = v1/bx
        temp[1,2] = v2/bx
        for n=1:ntau
            h[n,1] = ep*(sin(tau[n])*temp[1,1]-cos(tau[n])*temp[1,2])
            h[n,2] = ep*(sin(tau[n])*temp[1,2]+cos(tau[n])*temp[1,1])
        end

        xt[:,1] .= x1 .+ h[:,1] .- h[1,1]
        xt[:,2] .= x2 .+ h[:,2] .- h[1,2]

        p.e[1,m] = (0.5*cos(x1/2)*sin(x2))*(1+0.5*sin(time))
        p.e[2,m] = (sin(x1/2)*cos(x2))*(1+0.5*sin(time))

        @show bx
        for n=1:ntau
            interv = (1+0.5*sin(real(xt[n,1]))*sin(real(xt[n,2]))-bx)/bx
            r[n,1] =  interv*v2
            r[n,2] = -interv*v1
        end
        @show sum(r)

        mul!(tilde, ftau, r)

        @show sum(tilde)
        ave .= real(tilde[1,:])/ep

        @show ave

        for n=2:ntau
            tilde[n,:] .= -im*tilde[n,:]/ltau[n]
        end

        tilde[1,:] .= 0.0
        ldiv!(r, ftau, tilde)

        for n=1:ntau
            r[n,1]  = ep*(sin(tau[n])*p.e[1,m]+cos(tau[n])*p.e[2,m])/bx+r[n,1]
            r[n,2]  = ep*(sin(tau[n])*p.e[2,m]-cos(tau[n])*p.e[1,m])/bx+r[n,2]
            yt[n,1] = v1 + (r[n,1] - r[1,1])
            yt[n,2] = v2 + (r[n,2] - r[1,2])
        end

        @show sum(yt[:,1]), sum(yt[:,2])

        #--more preparation--

        for n=1:ntau
            temp[n,1]=ep*(cos(tau[n])*yt[n,1]+sin(tau[n])*yt[n,2])/bx
            temp[n,2]=ep*(cos(tau[n])*yt[n,2]-sin(tau[n])*yt[n,1])/bx
        end

        mul!(tilde, ftau, temp)

        for n=2:ntau
            tilde[n,:] .= -im*tilde[n,:]/ltau[n]
        end

        tilde[1,:] .= 0.0

        ldiv!(h, ftau, tilde)

        for n=1:ntau
            h[n,1] = h[n,1]-ep^2/bx*(-cos(tau[n])*ave[1]-sin(tau[n])*ave[2])
            h[n,2] = h[n,2]-ep^2/bx*(-cos(tau[n])*ave[2]+sin(tau[n])*ave[1])
        end

        xt[:,1] .= x1 .+ h[:,1] .- h[1,1]
        xt[:,2] .= x2 .+ h[:,2] .- h[1,2]

        @show sum(xt[:,1]), sum(xt[:,2])

        p.e[1,m] = (0.5*cos(x1/2)*sin(x2))*0.5*cos(time)
        p.e[2,m] = (sin(x1/2)*cos(x2))*0.5*cos(time)

        @show sum(p.e[1,1]), sum(p.e[2,1])

        for n=1:ntau

            interv = (1+0.5*sin(real(xt[n,1]))*sin(real(xt[n,2]))-bx)/bx

            fx[n,1] =  interv*ave[2]
            fx[n,2] = -interv*ave[1]

            fy[n,1] = ep/bx*(sin(tau[n])*ave[1]-cos(tau[n])*ave[2])
            fy[n,2] = ep/bx*(cos(tau[n])*ave[1]+sin(tau[n])*ave[2])

            interv = cos(x1)*sin(x2)*fy[n,1]+sin(x1)*cos(x2)*fy[n,2]

            fy[n,1] =  interv/bx/2 * v2 + fx[n,1]
            fy[n,2] = -interv/bx/2 * v1 + fx[n,2]

            fx[n,1] = ep/bx^2*(-sin(tau[n])*p.e[2,m]+cos(tau[n])*p.e[1,m])
            fx[n,2] = ep/bx^2*( sin(tau[n])*p.e[1,m]+cos(tau[n])*p.e[2,m])

        end

        temp .= fy .+ fx

        @show sum(fx[:,1]), sum(fx[:,2])
        @show sum(fy[:,1]), sum(fy[:,2])
        @show sum(temp)
        @show ep, bx, ave
        return

        tilde .= fft(temp,1)

        for n=2:ntau
            fx[n,:]    .= -im*tilde[n,:] ./ ltau[n]
            tilde[n,:] .= -tilde[n,:] ./ ltau[n]^2
        end

        fx[1,:]    .= 0.0
        tilde[1,:] .= 0.0

        temp .= ifft(tilde,1)

        r .= - ep * temp

        fy .= ifft(fx, 1)

        gx .= fy

        for n=1:ntau
            tilde[n,1]=(cos(tau[n])*fy[n,1]+sin(tau[n])*fy[n,2])/bx
            tilde[n,2]=(cos(tau[n])*fy[n,2]-sin(tau[n])*fy[n,1])/bx
        end

        temp .= fft(tilde, 1)

        ave2 .= temp[1,:] ./ ntau

        for n=1:ntau

            p.e[1,m]=(0.5*cos(real(xt[n,1])/2.)*sin(real(xt[n,2])))*(1+0.5*sin(time))
            p.e[2,m]=(sin(real(xt[n,1])/2.)*cos(real(xt[n,2])))*(1+0.5*sin(time))

            interv=(1+0.5*sin(real(xt[n,1]))*sin(real(xt[n,2]))-bx)/bx

            temp[n,1] =  interv*yt[n,2]+ep/bx*(-sin(tau[n])*p.e[2,m]+cos(tau[n])*p.e[1,m])
            temp[n,2] = -interv*yt[n,1]+ep/bx*( sin(tau[n])*p.e[1,m]+cos(tau[n])*p.e[2,m])
        end

        tilde .= fft(temp, 1)

        xf[1,:] .= tilde[1,:]/ep

        for n=2:ntau
            tilde[n,:] .= -im*tilde[n,:]/ltau[n]
        end
        tilde[1,:] .= 0.0
        temp .= ifft(tilde, 1)
        r .= r .+ temp

        yt[:,1] .= v1 .+ r[:,1] .- r[1,1]
        yt[:,2] .= v2 .+ r[:,2] .- r[1,2]

        #----end more preparation
        #--even more preparation--

        for n=1:ntau
            temp[n,1]=(cos(tau[n])*r[n,1]+sin(tau[n])*r[n,2])/bx
            temp[n,2]=(cos(tau[n])*r[n,2]-sin(tau[n])*r[n,1])/bx
        end

        tilde .= fft( temp, 1)

        temp[1,:] .= tilde[1,:]
        ave3      .= temp[1,:]
        interv     = cos(x1)*sin(x2)*temp[1,1]+sin(x1)*cos(x2)*temp[1,2]

        fx[1,1]    =  interv/ep/bx*v2/2
        fx[1,2]    = -interv/ep/bx*v1/2

        for n=1:ntau
            temp[n,1]=(1+0.5*sin(real(xt[n,1]))*sin(real(xt[n,2]))-bx)/bx
        end

        tilde[:,1] .= fft(temp[:,1], 1)

        fx[1,1] = fx[1,1] + tilde[1,1]/ep*ave[2]
        fx[1,2] = fx[1,2] - tilde[1,1]/ep*ave[1]

        for n=1:ntau
            yf[n,:] .= xf[1,:] .+ fy[n,:]
            temp[n,1] = (cos(tau[n])*yf[n,1]+sin(tau[n])*yf[n,2])
            temp[n,2] = (cos(tau[n])*yf[n,2]-sin(tau[n])*yf[n,1])
        end

        tilde .= fft(temp, 1)

        for n=2:ntau
            tilde[n,:] .= -im*tilde[n,:]/ltau[n]
        end

        tilde[1,:] .= 0.0

        temp .= ifft( tilde,1)

        for n=1:ntau
            fy[n,1]=temp[n,1]*ep/bx-ep^2/bx*(-cos(tau[n])*fx[1,1]-sin(tau[n])*fx[1,2])
            fy[n,2]=temp[n,2]*ep/bx-ep^2/bx*(-cos(tau[n])*fx[1,2]+sin(tau[n])*fx[1,1])
        end

        tilde .= fft(fy,1)

        for n=2:ntau
            tilde[n,:] .= -im * tilde[n,:]/ltau[n]
        end

        tilde[1,:] .= 0.0

        temp .= ifft(tilde,1)

        h .= -ep .* temp

        for n=1:ntau
            temp[n,1]=ep*(cos(tau[n])*yt[n,1]+sin(tau[n])*yt[n,2])/bx
            temp[n,2]=ep*(cos(tau[n])*yt[n,2]-sin(tau[n])*yt[n,1])/bx
        end

        tilde .= fft(temp, 1)

        for n=2:ntau
            tilde[n,:] .= -im * tilde[n,:]/ltau[n]
        end

        tilde[1,:] .= 0.0

        temp = ifft(tilde, 1)

        h .= h .+ temp

        xt[:,1] .= x1 .+ h[:,1] .- h[1,1]
        xt[:,2] .= x2 .+ h[:,2] .- h[1,2]

        #--end even more
        #--iteration
        for istep = 1:nstep

            #---imex2 New---
            for n=1:ntau

                et[1,n] = (0.5*cos(real(xt[n,1])/2)*sin(real(xt[n,2])))*(1+0.5*sin(time))
                et[2,n] = (cos(real(xt[n,2]))*sin(real(xt[n,1])/2))*(1+0.5*sin(time))

                interv = (1+0.5*sin(real(xt[n,1]))*sin(real(xt[n,2]))-bx)/bx/ep

                temp[n,1] = (cos(tau[n])*et[1,n]-sin(tau[n])*et[2,n])/bx
                temp[n,2] = (cos(tau[n])*et[2,n]+sin(tau[n])*et[1,n])/bx

                fy[n,1] = temp[n,1]+interv*yt[n,2]
                fy[n,2] = temp[n,2]-interv*yt[n,1]

            end

            gy .= yt .+ ds/2 .* fy

            fy .= fft(gy, 1)

            for n=1:ntau
                fy[n,:] .= fy[n,:]/(1.0+im*ds/2*ltau[n]/ep)
            end

            yf .= ifft(fy, 1)

            for n=1:ntau
                fx[n,1]=(cos(tau[n])*yf[n,1]+sin(tau[n])*yf[n,2])/bx
                fx[n,2]=(cos(tau[n])*yf[n,2]-sin(tau[n])*yf[n,1])/bx
            end

            gx .= xt .+ ds/2. * fx

            fx .= fft(gx, 1)

            for n=1:ntau
                fx[n,:] .= fx[n,:]./(1.0+im*ds/2*ltau[n]/ep)
            end

            xf .= ifft(fx, 1)

            time = time + dt/2

            for n=1:ntau

                et[1,n]=(0.5*cos(real(xf[n,1])/2)*sin(real(xf[n,2])))*(1+0.5*sin(time))
                et[2,n]=(cos(real(xf[n,2]))*sin(real(xf[n,1])/2))*(1+0.5*sin(time))

                interv=(1+0.5*sin(real(xf[n,1]))*sin(real(xf[n,2]))-bx)/bx/ep

                temp[n,1]=(cos(tau[n])*et[1,n]-sin(tau[n])*et[2,n])/bx
                temp[n,2]=(cos(tau[n])*et[2,n]+sin(tau[n])*et[1,n])/bx

                fy[n,1]=temp[n,1]+interv*yf[n,2]
                fy[n,2]=temp[n,2]-interv*yf[n,1]

            end

            gy .= fft(fy, 1)
            yf .= fft(yt, 1)

            for n=1:ntau
                fy[n,:] .= (yf[n,:]*(1.0-im*ds/ep/2.0*ltau[n])+ds*gy[n,:])/(1.0+im*ds/2.0*ltau[n]/ep)
            end

            yf .= yt

            yt .= ifft(fy, 1)

            yf .= (yt .+ yf)./2

            for n=1:ntau
                fx[n,1]=(cos(tau[n])*yf[n,1]+sin(tau[n])*yf[n,2])/bx
                fx[n,2]=(cos(tau[n])*yf[n,2]-sin(tau[n])*yf[n,1])/bx
            end

            gx .= fft(fx,1)
            xf .= fft(xt,1)

            for n=1:ntau
                fx[n,:] .= (xf[n,:]*(1.0-im*ds/ep/2.0*ltau[n])+ds*gx[n,:])/(1.0+im*ds/2.0*ltau[n]/ep)
            end

            xt .= ifft(fx, 1)

            time = time + dt/2.

            #---end imex2 New---

        end

        tilde .= fft(xt, 1)

        temp[1,:] .= 0.

        for n=1:ntau
            temp[1,:]=temp[1,:]+tilde[n,:]*exp(im*ltau[n]*tfinal*bx/ep)
        end

        xxt=real(temp[1,:])

        while ( xxt[1] > xmax )
            xxt[1] = xxt[1] - dimx
        end

        while ( xxt[1] < xmin )
            xxt[1]= xxt[1] + dimx
        end

        while ( xxt[2] > ymax )
            xxt[2]  = xxt[2]  - dimy
        end

        while ( xxt[2]  < ymin )
            xxt[2] = xxt[2]  + dimy
        end

        p.x[1,m] = xxt[1]
        p.x[2,m] = xxt[2]

        tilde .= fft(yt, 1)

        temp[1,:] .= 0.

        for n=1:ntau
            temp[1,:] .= temp[1,:] .+ tilde[n,:]*exp(im*ltau[n]*tfinal*bx/ep)
        end

        p.v[1,m]=real(cos(tfinal*bx/ep)*temp[1,1]+sin(tfinal*bx/ep)*temp[1,2])
        p.v[2,m]=real(cos(tfinal*bx/ep)*temp[1,2]-sin(tfinal*bx/ep)*temp[1,1])

    end

    println( sum(p.v[1,:])+857.95049281063064, " ", sum(p.v[2,:])+593.40700170710875 )

    compute_rho_m6!(f, p)

    #open(unit=851,file='UA3ep0001T1.dat')
    #for i=1,nx
    #for j=1,ny
    #write(851,*)i*dx, j*dy, f.rho(i,j)
    #end
    #write(851,*)
    #end
    #close(851)

end

efd( 16, 1 ) # trigger build

# efd( 16, 204800 )


text:
#
#--------------------------------------#
export gatherScatter
#--------------------------------------#
"""
Q*Q'*u where Q: local -> global operator
"""
function gatherScatter(u,QQtx::AbstractArray,QQty::AbstractArray)

# Qtu = ABu(Qy',Qx',  u); # gather
#QQtu = ABu(Qy ,Qx ,Qtu); # scatter

QQtu = ABu(QQty,QQtx,u); # gather scatter

return QQtu
end
#--------------------------------------#
function gatherScatter(u,msh::Mesh)
@unpack QQtx,QQty = msh
return gatherScatter(u,QQtx,QQty)
end
#--------------------------------------#
#Zygote.@adjoint function gatherScatter(As,Br,u)
#return ABu(As,Br,u),dv->(nothing,nothing,ABu(As',Br',dv));
#end
#--------------------------------------#
function gatherScatter(u,l2g::AbstractArray)

    Gu = NNlib.scatter(+,u,l2g) # gather - nnlib has its terminology reversed
    v  = Gu[l2g]
#   v  = NNlib.gather(Gu,l2g) # scatter

    return v
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#    API implementation for J4 orbit propagator.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

get_epoch(orbp::OrbitPropagatorJ4) = orbp.j4d.epoch

function get_mean_elements(orbp::OrbitPropagatorJ4)
    j4d = orbp.j4d
    orb = KeplerianElements(
        j4d.epoch + j4d.Δt / 86400,
        j4d.al_k * j4d.j4_gc.R0,
        j4d.e_k,
        j4d.i_k,
        j4d.Ω_k,
        j4d.ω_k,
        j4d.f_k
    )

    return orb
end

"""
    init_orbit_propagator(Val(:J4), epoch::Number, a_0::Number, e_0::Number, i_0::Number, Ω_0::Number, ω_0::Number, f_0::Number, dn_o2::Number = 0, ddn_o6::Number = 0; j4_gc::J4_GravCte{T} = j4_gc_egm08) where T
    init_orbit_propagator(Val(:J4), orb_0::Orbit, dn_o2::Number = 0, ddn_o6::Number = 0; j4_gc::J4_GravCte = j4_gc_egm08)

Initialize the J4 orbit propagator.

# Args

- `epoch::Number`: Initial orbit epoch [Julian Day].
- `a_0::Number`: Initial mean semi-major axis [m].
- `e_0::Number`: Initial mean eccentricity.
- `i_0::Number`: Initial mean inclination [rad].
- `Ω_0::Number`: Initial mean right ascension of the ascending node [rad].
- `ω_0::Number`: Initial mean argument of perigee [rad].
- `f_0::Number`: Initial mean true anomaly [rad].
- `dn_o2::Number`: (OPTIONAL) First time derivative of mean motion divided by 2
    \\[rad/s²] (**Default** = 0).
- `ddn_o6::Number`: (OPTIONAL) Second time derivative of mean motion divided by 6
    \\[rad/s³] (**Default** = 0).
- `orb_0::Orbit`: Object of type [`Orbit`](@ref) with the initial mean orbital
    elements [SI].

# Keywords

- `j4_gc::J4_GravCte`: (OPTIONAL) J4 orbit propagator gravitational constants.
    (**Default** = `j4_gc_egm08`).
"""
function init_orbit_propagator(::Val{:J4},
    epoch::Number,
    a_0::Number,
    e_0::Number,
    i_0::Number,
    Ω_0::Number,
    ω_0::Number,
    f_0::Number,
    dn_o2::Number = 0,
    ddn_o6::Number = 0;
    j4_gc::J4_GravCte{T} = j4_gc_egm08
) where T
    # Create the new J4 propagator structure.
    j4d = j4_init(
        epoch,
        a_0,
        e_0,
        i_0,
        Ω_0,
        ω_0,
        f_0,
        dn_o2,
        ddn_o6,
        j4_gc = j4_gc
    )

    # Create and return the orbit propagator structure.
    return OrbitPropagatorJ4(j4d)
end

function init_orbit_propagator(
    ::Val{:J4},
    orb_0::Orbit,
    dn_o2::Number = 0,
    ddn_o6::Number = 0;
    j4_gc::J4_GravCte = j4_gc_egm08
)
    # Convert the orbit representation to Keplerian elements.
    k_0 = convert(KeplerianElements, orb_0)

    return init_orbit_propagator(
        Val(:J4),
        k_0.t,
        k_0.a,
        k_0.e,
        k_0.i,
        k_0.Ω,
        k_0.ω,
        k_0.f,
        dn_o2,
        ddn_o6;
        j4_gc = j4_gc
    )
end

function propagate!(orbp::OrbitPropagatorJ4, t::Number)
    # Auxiliary variables.
    j4d = orbp.j4d

    # Propagate the orbit.
    r_i, v_i = j4!(j4d, t)

    # Return.
    return r_i, v_i
end

function step!(orbp::OrbitPropagatorJ4, Δt::Number)
    # Auxiliary variables.
    j4d = orbp.j4d

    # Propagate the orbit.
    r_i, v_i = j4!(j4d, j4d.Δt + Δt)

    # Return the information about the step.
    return r_i, v_i
end


text:
"""
Abstract replicator. The top node of the Replicator type hierarchy.
"""
abstract type AbstractReplicator end


text:
# Maximum Likelihood estimate with result analysis

using EstimateSDE
using Random
using Distributions
using Printf
using DataStructures
using Plots

## Settings
proc = testprocess(:ICMLA_B)
N = 200
Tav = 1.0
Tkld = 0.5 #0.2
p = 4
nspec = 1000
init_from_actual = false

## Start main
models = OrderedDict()
models["actual"] = proc
printstyled("Actual:\n",color=:blue)
describe_model(proc)

## Data
data = rand(MersenneTwister(5646356),proc,N,Tav)
problem = ProblemLike(data;p=p)

## Maximum Likelihood estimate
tor = MlEstimator()
Ts = [0.5,1,2]*Tkld
rng = MersenneTwister(35723)
init_tors = default_init_tors(rng,Ts)

if init_from_actual
  global proc
  push!(init_tors,InitialModel(proc)); @info "Adding actual as initial model"
end

## Maximum Likelihood estimate
mle,trace = @time estimate(tor,problem;init_tors=init_tors,ret_trace=true,throw_errors=false)

models["ML"] = mle
t_kld = (0:Tkld:N*Tav)
D_ml = kldiv(proc,mle,t_kld)
println("ML estimate: $mle\nKLD : $D_ml")

## Report
# Data
t = [x.t for x in data]; y = [x.y for x in data]
plt_data = plot(t,y,xlabel="t",ylabel="y",markershape=:x,fg_legend=nothing,bg_legend=nothing)
display(plt_data)
# Trace
plts = reporttrace(data,trace,Tkld;proc=proc,t_kld=t_kld)
nothing


text:
using CUDA
using Oceananigans.Architectures: CPU, GPU, AbstractMultiArchitecture

import Base: zeros

zeros(FT, ::CPU, N...) = zeros(FT, N...)
zeros(FT, ::GPU, N...) = CUDA.zeros(FT, N...)

zeros(arch::AbstractArchitecture, grid, N...) = zeros(eltype(grid), arch, N...)


text:
## State proposal Functions

# Jumptypes are:
# bg, move, add,remove, split, merge

function propose_move(rjs::RJStruct, currentstate::BAMFState)
    teststate = StateFlatBg(currentstate)
    # get an emitter
    idx = randID(currentstate.n)
    # move the emitter
    move_emitter!(idx, teststate.x, teststate.y, teststate.photons, rjs)
    return teststate, idx
end
function move_emitter!(ID::Int, x::CuArray, y::CuArray, photons::CuArray, rjs::RJStruct)
    @cuda move_emitter_CUDA!(ID, x, y, photons, rjs.xy_std, rjs.I_std)
end
function move_emitter_CUDA!(ID::Int, x, y, photons, xy_std::Float32, i_std::Float32)
    x[ID] += xy_std * curandn()
    y[ID] += xy_std * curandn()
    photons[ID] += i_std * curandn()
    photons[ID] = max(0, photons[ID])
    return nothing
end
function move_emitter!(ID::Int, x::Vector{Float32}, y::Vector{Float32}, photons::Vector{Float32}, rjs::RJStruct)
    if ID < 1 return nothing end
    x[ID] += rjs.σ_xy * randn()
    y[ID] += rjs.σ_xy * randn()

    sz=rjs.data.sz
    x[ID]=max(0-rjs.bndpixels,min(sz+rjs.bndpixels,x[ID]))
    y[ID]=max(0-rjs.bndpixels,min(sz+rjs.bndpixels,y[ID]))

    photons[ID] += rjs.σ_photons * randn()
    photons[ID] = max(1e-6, photons[ID])
    return nothing
end

function propose_bg(rjs::RJStruct, currentstate::BAMFState)
    teststate = StateFlatBg(currentstate)
    teststate.bg=currentstate.bg+rjs.σ_background*randn()
    return teststate, 0
end

## Add and remove -------------

function propose_add(rjs::RJStruct, currentstate::BAMFState)

    teststate = StateFlatBg(currentstate)

    # calc residum
    residuum = calcresiduum(rjs.modeldata, rjs.data)
    # pick pixel
    ii, jj = arrayrand!(residuum)

    # pick intensity from prior
    photons = rand(rjs.prior_photons)
    if photons < 0
        println(photons)
    end
    # add new emitter
    addemitter!(teststate, Float32(ii), Float32(jj), Float32(photons))

    return teststate, teststate.n
end

function propose_remove(rjs::RJStruct, currentstate::BAMFState)
    teststate = StateFlatBg(currentstate)
    # get an emitter
    idx = randID(currentstate.n)
    if idx > 0
        removeemitter!(teststate, idx)
    end
    return teststate, idx
end

## ---------------------------

function propose_split(rjs::RJStruct, currentstate::BAMFState)

    # println("split")

    if currentstate.n<1
        println("attempting to split 0 sources")
    end

    teststate = StateFlatBg(currentstate)
    # get an emitter
    idx = randID(currentstate.n)

    if idx==0
        return teststate,(idx,teststate.n,0f0,0f0,0f0)

    end

    # these are conserved
    photons_total = teststate.photons[idx]
	mux = teststate.x[idx]
    muy = teststate.y[idx]

    split_std = rjs.σ_split
    u1 = rand(Float32)
	u2 = split_std * randn(Float32)
	u3 = split_std * randn(Float32)

	# Change one emitter
    teststate.photons[idx] = u1 * photons_total
    teststate.x[idx] = mux + u2
    teststate.y[idx] = muy + u3

    # add an emitter
    photons_new = (1 - u1) * photons_total
    x_new = mux - u1 * u2 / (1 - u1)
    y_new = muy - u1 * u3 / (1 - u1)

    addemitter!(teststate, Float32(x_new), Float32(y_new), Float32(photons_new))

    return teststate,(idx,teststate.n,u1,u2,u3)

end

function propose_merge(rjs::RJStruct, currentstate::BAMFState)

    # println("merge")
    teststate = StateFlatBg(currentstate)
    if currentstate.n<2 return teststate,(0,0,0f1,0f1,0f1) end

    # get an emitter
    idx1 = randID(currentstate.n)
    idx2 = findother(currentstate,idx1)

    #keep the lower idx
    if idx1>idx2
        tmp=idx2
        idx2=idx1
        idx1=tmp
    end

    #merged values
    mux = (currentstate.x[idx1]*currentstate.photons[idx1]+
    currentstate.x[idx2]*currentstate.photons[idx2])/
    (currentstate.photons[idx1]+currentstate.photons[idx2])

    muy = (currentstate.y[idx1]*currentstate.photons[idx1]+
    currentstate.y[idx2]*currentstate.photons[idx2])/
    (currentstate.photons[idx1]+currentstate.photons[idx2])

    photons_total = currentstate.photons[idx1]+currentstate.photons[idx2]

    #replace idx1
    teststate.x[idx1]=mux
    teststate.y[idx1]=muy
    teststate.photons[idx1]=photons_total
    #remove second
    removeemitter!(teststate,idx2)

    u1=currentstate.photons[idx1]/photons_total
    u2=currentstate.x[idx1]-currentstate.x[idx2]
    u3=currentstate.y[idx1]-currentstate.y[idx2]

    return teststate,(idx1,idx2,u1,u2,u3)
end


text:
__precompile__()

module ImageFiltering

using Colors, FixedPointNumbers, ImageCore, MappedArrays, FFTViews, OffsetArrays, StaticArrays, ComputationalResources, TiledIteration
using ColorVectorSpace  # for filtering RGB arrays
using Base: Indices, tail, fill_to_length, @pure, depwarn

export Kernel, KernelFactors, Pad, Fill, Inner, NA, NoPad, Algorithm, imfilter, imfilter!, extrema_filter, imgradients, padarray, centered, kernelfactors, reflect

typealias FixedColorant{T<:UFixed} Colorant{T}
typealias StaticOffsetArray{T,N,A<:StaticArray} OffsetArray{T,N,A}

# Needed for type-stability
function Base.transpose{T}(A::StaticOffsetArray{T,2})
    inds1, inds2 = indices(A)
    OffsetArray(transpose(parent(A)), inds2, inds1)
end

module Algorithm
    # deliberately don't export these, but it's expected that they
    # will be used as Algorithm.FFT(), etc.
    abstract Alg
    "Filter using the Fast Fourier Transform" immutable FFT <: Alg end
    "Filter using a direct algorithm" immutable FIR <: Alg end
    "Filter with an Infinite Impulse Response filter" immutable IIR <: Alg end
    "Filter with a cascade of mixed types (IIR, FIR)" immutable Mixed <: Alg end
end
using .Algorithm: Alg, FFT, FIR, IIR, Mixed

Alg{A<:Alg}(r::AbstractResource{A}) = r.settings

include("utils.jl")
include("kernelfactors.jl")
using .KernelFactors: TriggsSdika, IIRFilter, ReshapedOneD, iterdims, kernelfactors

typealias ReshapedVector{T,N,Npre,V<:AbstractVector} ReshapedOneD{T,N,Npre,V}
typealias ArrayType{T} Union{AbstractArray{T}, ReshapedVector{T}}
typealias ReshapedIIR{T,N,Npre,V<:IIRFilter} ReshapedOneD{T,N,Npre,V}
typealias AnyIIR{T} Union{IIRFilter{T}, ReshapedIIR{T}}
typealias ArrayLike{T} Union{ArrayType{T}, AnyIIR{T}}

include("kernel.jl")
using .Kernel
using .Kernel: Laplacian, reflect

typealias NDimKernel{N,K} Union{AbstractArray{K,N},ReshapedOneD{K,N},Laplacian{N}}

include("border.jl")

typealias BorderSpec{T} Union{Pad{0}, Fill{T,0}, Inner{0}}
typealias BorderSpecNoNa{T} Union{Pad{0}, Fill{T,0}, Inner{0}}
typealias BorderSpecAny Union{BorderSpec,NA,NoPad}

include("deprecated.jl")

typealias ProcessedKernel Tuple

include("imfilter.jl")
include("specialty.jl")

include("rank.jl")
using .Rank

function __init__()
    # See ComputationalResources README for explanation
    push!(LOAD_PATH, dirname(@__FILE__))
    # if haveresource(ArrayFireLibs)
    #     @eval using DummyAF
    # end
    pop!(LOAD_PATH)
end

end # module


text:
"""
    functional_cartography(N::T, L::Dict{E,Int64}) where {T<:BinaryNetwork, E}

This function will take the output of a modularity analysis (*i.e.* a network
and a partition), and return a dictionary where every species is associated to
its functional role, as defined in Olesen et al (2005). The first element is the
within-module degree z-score, and the second is the participation coefficient.

#### References

Guimerà, R., Amaral, L.A.N., 2005. Cartography of complex networks: modules and
universal roles. Journal of Statistical Mechanics: Theory and Experiment 2005,
P02001. https://doi.org/10.1088/1742-5468/2005/02/P02001

Guimerà, R., Nunes Amaral, L.A., 2005. Functional cartography of complex
metabolic networks. Nature 433, 895–900. https://doi.org/10.1038/nature03288

Olesen, J.M., Bascompte, J., Dupont, Y.L., Jordano, P., 2007. The modularity of
pollination networks. Proceedings of the National Academy of Sciences 104,
19891–19896. https://doi.org/10.1073/pnas.0706375104
"""
function functional_cartography(N::T, L::Dict{E,Int64}) where {T<:BinaryNetwork, E}
    κ = Dict{last(eltype(N)), Int64}()
    S = Dict{Int64,Tuple{Float64, Float64}}()
    carto = Dict{last(eltype(N)),Tuple{Float64, Float64}}()

    for s in species(N)
        sp_module = L[s]
        k_i = 0
        if s ∈ species(N; dims=2)
            for s2 in N[:,s]
                L[s2] == sp_module && (k_i += 1)
            end
        end
        if s ∈ species(N; dims=1)
            for s2 in N[s,:]
                L[s2] == sp_module && (k_i += 1)
            end
        end
        κ[s] = k_i
    end

    for s in values(L)
        sp = filter(k -> k.second == s, L)
        k_i = [get(κ, x, 0) for x in keys(sp)]
        S[s] = (mean(k_i), std(k_i))
    end

    # Participation coefficient
    for s in species(N)
        ks = zeros(Int64, length(unique(values(L))))
        if s ∈ species(N; dims=2)
            for s2 in N[:,s]
                ks[L[s2]] += 1
            end
        end
        if s ∈ species(N; dims=1)
            for s2 in N[s,:]
                ks[L[s2]] += 1
            end
        end
        carto[s] = ((κ[s] - S[L[s]][1])/S[L[s]][2], 1-sum((ks./degree(N)[s]).^2))
    end

    return carto
end


text:
using Revise
using DataFrames
import CSV
import RetroSignalModel as rs
using Test

@testset "Parameter searching" begin
    include("searching.jl")
end

@testset "Data files" begin
    include("load_data.jl")
end

@testset "Build models" begin
    include("build_model.jl")
end


text:
using BigArrays
using BigArrays.BinDicts
using BigArrays.Infos
using Test
using OffsetArrays

# prepare directory
layerDir = tempname()
@show layerDir
datasetDir = joinpath(layerDir, "6_6_30")
mkdir(layerDir)
mkdir(datasetDir)
infoString = """
{"num_channels": 1, "type": "image", "data_type": "uint8", "scales": [
{"encoding": "gzip", "chunk_sizes": [[100, 100, 5]], "key": "6_6_30", "resolution": [6, 6, 30], "voxel_offset": [-300, -300, -10], "size": [510, 510, 2022]},
{"encoding": "gzip", "chunk_sizes": [[100, 100, 5]], "key": "12_12_30", "resolution": [12, 12, 30], "voxel_offset": [-597, -597, -103], "size": [12286, 11262, 2046]}
]}
"""

write( joinpath(layerDir, "info"), infoString )

@testset "test bigarray construction" begin
    info = Info()
    ba = BigArray(info)
    BigArrays.commit_info(ba)
end

@testset "\ntest BinDict" begin
    h = BinDict(layerDir)
    a = rand(UInt8, 20)
    h["test"] = a
    b = h["test"]
    @test haskey(h, "test")
    @test !haskey(h, "notexist")
    @show a
    @show b
    @test all(a.==b)
end # testset

@testset "\ntest file protocol" begin
    # use sequential mode for debug
    ba = BigArray( ("file:/" * layerDir); mode=:sequential )
    println("info: ", BigArrays.get_info(ba))
    a = rand(UInt8, 200,200,10)
    ba[1:200, 1:200, 1:10] = a
    b = ba[1:200, 1:200, 1:10]
    @test all(a.==parent(b))
end

@testset "\ntest IO of BigArray with backend of BinDict" begin
    ba = BigArray( layerDir )
    a = rand(UInt8, 200,200,10)
    ba[1:200, 1:200, 1:10] = a
    b = ba[1:200, 1:200, 1:10]
    @test all(a.==parent(b))
end # end of testset

@testset "\ntest negative coordinate" begin
    ba = BigArray( BinDict(layerDir) )
    a = rand(UInt8, 200,200,2000)
    ba[-199:0, -99:100, -4:1995] = a
    b = ba[-199:0, -99:100, -4:1995]
    @test all(a.==parent(b))
end # end of testset

#@testset "test sharedarray mode..." begin
#    ba = BigArray( BinDict(layerDir); mode=:sharedarray )
#    a = rand(UInt8, 200,200,2000)
#    ba[-199:0, -99:100, -4:1995] = a
#    b = ba[-199:0, -99:100, -4:1995]
#    @test all(a.==parent(b))
#end # end of testset

@testset "\ntest aligned IO crossing the volume boundary" begin
    ba = BigArray( BinDict(layerDir) )
    a = rand(UInt8, 200,200,10)
    # respect the volume size, the chunk range over volume size will not be written
    ba[101:300, 101:300, 2006:2015] = a
    b = ba[101:300, 101:300, 2006:2015]
    c = zero(a)
    c[1:110, 1:110, 1:7] = a[1:110, 1:110, 1:7]
    @test all(c.==parent(b))
end # end of testset

@testset "\ntest nonaligned IO crossing the volume boundary" begin
    ba = BigArray( BinDict(layerDir) )
    a = rand(UInt8, 190,190,9)
    # respect the volume size, the chunk range over volume size will not be written
    ba[101:290, 101:290, 2006:2014] = a
    b = ba[101:290, 101:290, 2006:2014]
    c = zero(a)
    c[1:110, 1:110, 1:7] = a[1:110, 1:110, 1:7]
    @test all(c.==parent(b))
end # end of testset

infoString = replace(infoString, "gzip" => "zstd")
write( joinpath(layerDir, "info"), infoString )

@testset "test IO of BigArray with backend of BinDict with zstd compression" begin
    ba = BigArray( BinDict(layerDir) )
    a = rand(UInt8, 200,200,10)
    ba[1:200, 1:200, 1:10] = a
    b = ba[1:200, 1:200, 1:10]
    @test all(a.==parent(b))
end # end of testset

infoString = replace(infoString, "zstd"=>"blosclz")
write( joinpath(layerDir, "info"), infoString )

@testset "test IO of BigArray with backend of BinDict with blosclz compression" begin
    ba = BigArray( BinDict(layerDir) )
    a = rand(UInt8, 200,200,10)
    ba[1:200, 1:200, 1:10] = a
    b = ba[1:200, 1:200, 1:10]
    @test all(a.==parent(b))
end # end of testset

#@testset "test merge function with backend of BinDict with blosclz compression" begin
#    ba = BigArray( BinDict(layerDir) )
#    a = rand(UInt8, 200,200,10)
#    @inbounds merge!(ba, OffsetArray(a, 1:200, 1:200, 1:10))
#    @inbounds b = ba[1:200, 1:200, 1:10]
#    @test all(parent(a).==parent(b))
#end # end of testset

infoString = replace(infoString, "blosclz" => "zstd")
write( joinpath(layerDir, "info"), infoString )

@testset "test dataset not aligned starting from 0" begin
    datasetDir = joinpath(layerDir, "12_12_30")
    mkdir(datasetDir)
    ba = BigArray( BinDict(layerDir); mip=2 )
    a = rand(UInt8, 200,200,10)
    ba[204:403, 204:403, 103:112] = a
    b = ba[204:403, 204:403, 103:112]
    @test all(a.==parent(b))
end # end of testset

@testset "test dataset not aligned starting from 0 and negative coordinates" begin
    ba = BigArray( BinDict(layerDir); mip=2 )
    a = rand(UInt8, 200,200,10)
    ba[-96:103, -296:-97, -2:7] = a
    b = ba[-96:103, -296:-97, -2:7]
    @test all(a.==parent(b))
end # end of testset

# clean the temporary directory
rm(layerDir; recursive=true)


text:
using Pseudospectra, Test

@testset "Zoom" begin

    A = Pseudospectra.grcar(40)
    opts = Dict{Symbol,Any}()
    ps_data = new_matrix(A,opts)

    driver!(ps_data,opts,gs)

    zoomin!(gs,ps_data,zkw=1.2+1.8im)
    @test ps_data.zoom_pos == 2

    # this should return to original portrait
    zoomout!(gs,ps_data,zkw=0.0)
    @test length(ps_data.zoom_list) == 2
    @test ps_data.zoom_pos == 1

    # OOB: this should cause a warning:
    @test (@test_logs (:warn,r"^unable to zoom") zoomout!(gs,ps_data,zkw=5.0+3.5im)) == -1

    zoomout!(gs,ps_data,zkw=0.5+2.0im)
    @test ps_data.zoom_pos == 1

    @test length(ps_data.zoom_list) == 3
    @test all(map(iscomputed,ps_data.zoom_list))
end


text:
using Test
using StatsBase

@testset "shape" begin
    s = Kronecker_sampler([1,2,3,4], 2)
    e = rand(s)
    @test typeof(e) == Tuple{Int}
    @test length(e) == 1
    g = rand(s, 7)
    @test eltype(g) == typeof(e)
    @test axes(first(g)) == axes(e)
    @test typeof(g) == Vector{Tuple{Int}}
    @test size(g) == (7,)

    s = Kronecker_sampler(rand(2,3,1), 2)
    e = rand(s)
    @test typeof(e) == NTuple{3, Int}
    @test length(e) == 3
    g = rand(s, 3, 3)
    @test eltype(g) == typeof(e)
    @test axes(first(g)) == axes(e)
    @test typeof(g) == Matrix{NTuple{3, Int}}
    @test size(g) == (3, 3)
end

@testset "hit all" begin
    for (vals, p) in [((4,1,2), 3), ((3,3), 1), ((1,1), 12)]
        n = prod(vals)^p
        m = Int(floor(n*log(n)*3+100))
        for space in [0, m]
            s = Kronecker_sampler(ones(vals...), p; space=space)

            g = rand(s, m)

            @test length(g) == m
            @test length(countmap(g)) == n

            if length(countmap(g)) != n
                println(vals, " ", p, " ", space)
            end
        end
    end
end

max_stars_repo_name: LilithHafner/FBD

max_issues_repo_name: LilithHafner/FBD

max_forks_repo_name: LilithHafner/FBD


text:
getdata(data) = parse.(Int, split(data, ""))
dats = getdata(readline("data/Dec16_data.txt"))

l = length(dats)
for i in eachindex(dats)
pat = repeat([0,1,0,-1], inner = i, outer = 1 + l ÷ 4i)[2:l+1]
res = getindex.(digits.(dats .* pat), Ref(1))
dats .* pat


text:
# TODO:  USER:  Update the paths below with path to the data on your system

# Default
solar_data_path = "/home/eford/Data/SolarSpectra/NEID_solar/"
expres_data_path = "/home/eford/Data/EXPRES/inputs/"
output_dir = "examples/output"

max_stars_repo_name: RvSpectML/EchelleInstruments

max_issues_repo_name: RvSpectML/EchelleInstruments

max_forks_repo_name: RvSpectML/EchelleInstruments


text:
using XLSX, DataFrames

GCPfiledata = Dict(
    2017 => Dict("globalversion" => "v1.3", "landuserange" => "C21:C78", "nationalversion" => "v1.2",
                    "sheet" => "Territorial Emissions GCB", "range" => "A16:IC74"),
    2018 => Dict("globalversion" => "v1.0", "landuserange" => "C21:C79", "nationalversion" => "v1.0",
                    "sheet" => "Territorial Emissions", "range" => "A17:HV76"),
    2019 => Dict("globalversion" => "v1.0", "landuserange" => "D20:D79", "nationalversion" => "v1.0",
                    "sheet" => "Territorial Emissions", "range" => "A17:HW77")
)

RCPregions = Dict(
    # There is already a column "Europe" but it includes Russia and others
    "Europe_alt" => ["EU28", "Iceland", "Norway", "Switzerland", "Turkey", "Macedonia (Republic of)",
                "Albania", "Bosnia and Herzegovina", "Serbia", "Montenegro"],
     # add Europe to OECD in code below
    "OECD" => ["Australia", "Canada", "Japan", "New Zealand", "USA"],
    "REF" => ["Armenia", "Azerbaijan", "Belarus", "Georgia", "Kazakhstan", "Kyrgyzstan", "Moldova",
                "Russian Federation", "Tajikistan", "Turkmenistan", "Ukraine", "Uzbekistan"],
    "ASIA" => ["Afghanistan", "Bangladesh", "Bhutan", "Brunei Darussalam", "Cambodia", "China",
                "Hong Kong", "Macao", "North Korea", "Fiji", "French Polynesia", "India", "Indonesia",
                "Laos", "Malaysia", "Maldives", "Micronesia (Federated States of)", "Mongolia", "Myanmar",
                "Nepal", "New Caledonia", "Pakistan", "Papua New Guinea", "Philippines", "South Korea",
                "Samoa", "Singapore", "Solomon Islands", "Sri Lanka", "Taiwan", "Thailand", "Timor-Leste",
                "Vanuatu", "Viet Nam"],
    "MAF" => ["Algeria", "Angola", "Bahrain", "Benin", "Botswana", "Burkina Faso", "Burundi", "Cameroon",
                "Cape Verde", "Central African Republic", "Chad", "Comoros", "Congo", "Côte d'Ivoire",
                "Democratic Republic of the Congo", "Djibouti", "Egypt", "Equatorial Guinea", "Eritrea",
                "Ethiopia", "Gabon", "Gambia", "Ghana", "Guinea", "Guinea-Bissau", "Iraq", "Iran", "Israel",
                "Hong Kong", "Jordan", "Kenya", "Kuwait", "Lebanon", "Lesotho", "Liberia", "Libya",
                "Madagascar", "Malawi", "Mali", "Mauritania", "Mauritius", "Morocco", "Mozambique",
                "Namibia", "Niger", "Nigeria", "Occupied Palestinian Territory", "Oman", "Qatar", "Rwanda",
                "Réunion", "Saudi Arabia", "Senegal", "Sierra Leone", "Somalia", "South Africa",
                "South Sudan", "Sudan", "Swaziland", "Syria", "Togo", "Tunisia", "Uganda",
                "United Arab Emirates", "Tanzania", "Yemen", "Zambia", "Zimbabwe"],
    "LAM" => ["Argentina", "Aruba", "Bahamas", "Barbados", "Belize", "Bolivia", "Brazil", "Chile",
                "Colombia", "Costa Rica", "Cuba", "Dominican Republic", "Ecuador", "El Salvador",
                "French Guiana", "Grenada", "Guadeloupe", "Guatemala", "Guyana", "Haiti", "Honduras",
                "Jamaica", "Martinique", "Mexico", "Nicaragua", "Panama", "Paraguay", "Peru", "Suriname",
                "Trinidad and Tobago", "Uruguay", "Venezuela"]
)

# Spreadsheet headers not listed above.
# Should be just small islands, EU countries (included in "EU28") and region names.
unlisted = Symbol.(["Andorra", "Anguilla", "Antigua and Barbuda", "Austria", "Belgium", "Bermuda",
            "Bonaire, Saint Eustatius and Saba", "British Virgin Islands", "Bulgaria", "Cayman Islands",
            "Cook Islands", "Croatia", "Curaçao", "Cyprus", "Czech Republic", "Denmark", "Dominica",
            "Estonia", "Faeroe Islands", "Falkland Islands (Malvinas)", "Finland", "France", "Germany",
            "Gibraltar", "Greece", "Greenland", "Hungary", "Ireland", "Italy", "Kiribati", "Kosovo",
            "Latvia", "Liechtenstein", "Lithuania", "Luxembourg", "Malta", "Marshall Islands", "Montserrat",
            "Nauru", "Netherlands", "Niue", "Palau", "Poland", "Portugal", "Romania", "Saint Helena",
            "Saint Lucia", "Sint Maarten (Dutch part)", "Sao Tome and Principe", "Seychelles", "Slovakia",
            "Slovenia", "Spain", "Saint Kitts and Nevis", "Saint Pierre and Miquelon",
            "Saint Vincent and the Grenadines", "Sweden", "Tonga", "Turks and Caicos Islands", "Tuvalu",
            "United Kingdom", "Wallis and Futuna Islands", "KP Annex B", "Non KP Annex B", "OECD", "Non-OECD",
            "Africa", "Asia", "Central America", "Europe", "Middle East", "North America", "Oceania",
            "South America", "Bunkers", "Statistical Difference", "World"])

function importGlobalCarbonProject(year)
    filedata = GCPfiledata[year]
    globalversion, nationalversion = filedata["globalversion"], filedata["nationalversion"]
    globalfile = joinpath(@__DIR__, "..", "SSP", "Global_Carbon_Budget_$year$globalversion.xlsx")
    nationalfile = joinpath(@__DIR__, "..", "SSP", "National_Carbon_Emissions_$year$nationalversion.xlsx")
    allheaders = Symbol.(vcat(values(RCPregions)...))
    # @show allheaders

    # segfault, see https://github.com/felipenoris/XLSX.jl/issues/146
    # XLSX.readtable(nationalfile, filedata["sheet"], filedata["range"])

    data = XLSX.readdata(nationalfile, filedata["sheet"], filedata["range"])
    data = replace(data, "NaN" => 0) # zero OK because only tiny islands & countries are missing data
    # data = replace(data, "NaN" => missing)
    headers = Symbol.(["year"; data[1, 2:end]])
    df = DataFrame(data[2:end,:], headers)

    landuse = XLSX.readdata(globalfile, "Global Carbon Budget", filedata["landuserange"])

    for h in headers[2:end]
        # !(h in allheaders) && print("\"$h\", ")
        !(h in allheaders) && !(h in unlisted) && println("$h missing in RCPregions.")
    end

    for (reg, countries) in RCPregions
        df[reg] = sum(df[c] for c in countries if Symbol(c) in headers)
    end
    df.GLOBALFOSSIL = df.World
    df.OECD = df.OECD + df.Europe_alt
    df.BUNKERS = df.GLOBALFOSSIL - (df.OECD + df.REF + df.ASIA + df.MAF + df.LAM)
    df.Asia = df.ASIA
    df.LANDUSE = vec(landuse) * 1000

    select!(df, [:OECD, :REF, :Asia, :MAF, :LAM, :BUNKERS, :GLOBALFOSSIL, :LANDUSE])
    for col in names(df)
        df[col] = round.(df[col], digits=6)
    end

    preamble = """
    // File automatically generated from spreadsheet data by importGlobalCarbonProject.jl.

    // Data from The Global Carbon Budget $year. Unit [MtC/year].
    // Source: https://www.globalcarbonproject.org/carbonbudget/

    // National (territorial) CO2 emissions 1959-$(year-1), aggregated to RCP/SSP regions.
    // BUNKERS here is sum of actual BUNKERS field and the statistical error (to make regions sum to global total).
    // Asia not capitalized as in the RCP/SSP standard because it is referenced in the web interface.
    // Land use emissions are from the global overview file, the rest are from the national file.

    var CO2emissionhistory = {
    """

    footer = """
    };

    // first year of the above data
    var backgrounddatastart = 1959;
    """

    filename = joinpath(@__DIR__, "UI_backgrounddata.js")
    open(filename, "w") do file
        print(file, preamble)
        for name in names(df)
            name == :GLOBALFOSSIL && continue
            println(file, "    $name: [", join(df[name], ", "), "],")
        end
        print(file, footer)
    end

    println("\nFile $filename written.\n")
    return df
end

# Helper function used to transition from the old file:
#     Global Carbon Budget 2017 - National_Carbon_Emissions_2017v1.1.xlsx
# Reads Excel equation strings and outputs countries, for example:
# ClimateCalculator.getcolumns("=IO17+L17+AI17+CW17+EK17+HD17") =>
#    "IO17", "Australia", "Canada", "Japan", "New Zealand", "USA"
function getcolumns(excelstring)
    year = 2017
    filedata = GCPfiledata[year]
    globalversion, nationalversion = filedata["globalversion"], filedata["nationalversion"]
    globalfile = joinpath(@__DIR__, "..", "SSP", "Global_Carbon_Budget_$year$globalversion.xlsx")
    nationalfile = joinpath(@__DIR__, "..", "SSP", "National_Carbon_Emissions_$year$nationalversion.xlsx")

    data = XLSX.readdata(nationalfile, filedata["sheet"], filedata["range"])
    headers = ["year"; data[1, 2:end]]

    cells = [startswith(c, "IF(") ? split(c, ",")[2] : c for c in split(excelstring[2:end], "+")]
    columnnumbers = XLSX.column_number.(XLSX.CellRef.(cells))
    countries = [cn > length(headers) ? cells[i] : headers[cn] for (i,cn) in enumerate(columnnumbers)]
    for c in countries
        print("\"$c\", ")
    end
end

max_stars_repo_name: niclasmattsson/ClimateCalculator

max_issues_repo_name: niclasmattsson/ClimateCalculator

max_forks_repo_name: niclasmattsson/ClimateCalculator


text:
function compute()::Int
    numerators_product = 1
    denominators_product = 1
    for i ∈ 1:9
        for j ∈ 1:i - 1
            for k ∈ 1:j - 1
                if (k * 10 + i) * j == (i * 10 + j) * k
                    numerators_product *= k
                    denominators_product *= j
                end
            end
        end
    end
    return(denominators_product ÷ gcd(numerators_product, denominators_product))
end

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler


text:
using Acorn
using Test

# Tests for commands
include("cmds/test_bind.jl")
# TODO: include tests for all commands

# TODO: Include additional tests
# rows.jl
# editor.jl
# EditorCOnfig.jl


text:
module GroupedTemporalTerms

include("grouped-lag.jl")

export groupedlag, GroupedLagTerm

end


text:
"""
    start_val(model)

Return a vector of starting values.
Defaults are FABIN 3 starting values for single models and simple starting values for
ensemble models.
"""
function start_val end
# Single Models ----------------------------------------------------------------------------

# splice model and loss functions
start_val(model::AbstractSemSingle; kwargs...) =
    start_val(
        model,
        model.observed,
        model.imply,
        model.optimizer,
        model.loss.functions...;
        kwargs...)

# Fabin 3 starting values for RAM(Symbolic)
start_val(
        model,
        observed,
        imply,
        optimizer,
        args...;
        kwargs...) =
    start_fabin3(model; kwargs...)

# Ensemble Models --------------------------------------------------------------------------
start_val(model::SemEnsemble; kwargs...) = start_simple(model; kwargs...)


text:
const ExprLike = Union{Symbol,Expr}
const Optional{T} = Union{Nothing,T}

is_ptr_to_ptr(ex) = !isnothing(ptrtype(ptrtype(ex)))
is_ptr(ex) = !isnothing(ptr_type(ex))

ptr_type(ex) = @when :(Ptr{$T}) = ex T
ntuple_type(ex) = @when :(NTuple{$N,$T}) = ex T

is_ntuple(ex) = !isnothing(ntuple_type(ex))

is_vulkan_type(name) = name ∈ vcat(spec_structs.name, spec_handles.name)

inner_type(ex) = @when :($T{$(args...)}) = ex map(args) do arg
    @match arg begin
        :(<:$t) => t
        _ => arg
    end
end

function innermost_type(ex::Expr)
    if is_ntuple(ex)
        innermost_type(ntuple_type(ex))
    else
        t = inner_type(ex)
        if !isnothing(t)
            if length(t) > 1
                error("Expected 1 inner type for $ex, found $(length(t)) ($t)")
            else
                innermost_type(first(t))
            end
        else
            ex
        end
    end
end

innermost_type(sym::Symbol) = sym


text:
module TestCheckFoodWeb
    using BioEnergeticFoodWebs
    using Test

    # A correct network will return nothing
    correct_network = [0 0 1 1; 1 1 0 0; 0 1 0 1; 1 1 0 0]
    @test BioEnergeticFoodWebs.check_food_web(correct_network) == nothing

    # A non-square network will fail
    non_square = [0 0 1; 0 1 1; 0 1 0; 1 0 0]
    @test_throws AssertionError BioEnergeticFoodWebs.check_food_web(non_square)

    # A network with antyhing else than 0/1 will fail
    non_binary = [0 0 1 1; 1 1 0 0; 0 1 0 1; 1 1 2 0]
    @test_throws AssertionError BioEnergeticFoodWebs.check_food_web(non_binary)

end

max_stars_repo_name: SimmonsBI/BioEnergeticFoodWebs

max_issues_repo_name: SimmonsBI/BioEnergeticFoodWebs

max_forks_repo_name: SimmonsBI/BioEnergeticFoodWebs


text:
## symbols ##

symbol(s::Symbol) = s
symbol(s::ASCIIString) = symbol(s.data)
symbol(s::UTF8String) = symbol(s.data)
symbol(a::Array{Uint8,1}) =
    ccall(:jl_symbol_n, Any, (Ptr{Uint8}, Int32), a, length(a))::Symbol
symbol(x::Char) = symbol(string(x))

gensym() = ccall(:jl_gensym, Any, ())::Symbol

gensym(s::ASCIIString) = gensym(s.data)
gensym(s::UTF8String) = gensym(s.data)
gensym(a::Array{Uint8,1}) =
    ccall(:jl_tagged_gensym, Any, (Ptr{Uint8}, Int32), a, length(a))::Symbol
gensym(ss::Union(ASCIIString, UTF8String)...) = map(gensym, ss)

macro gensym(names...)
    blk = Expr(:block)
    for name in names
        push!(blk.args, :($(esc(name)) = gensym($(string(name)))))
    end
    push!(blk.args, :nothing)
    return blk
end

## expressions ##

splicedexpr(hd::Symbol, args::Array{Any,1}) = (e=Expr(hd); e.args=args; e)
copy(e::Expr) = (n = Expr(e.head);
                 n.args = astcopy(e.args);
                 n.typ = e.typ;
                 n)
copy(s::SymbolNode) = SymbolNode(s.name, s.typ)
copy(n::GetfieldNode) = GetfieldNode(n.value, n.name, n.typ)

# copy parts of an AST that the compiler mutates
astcopy(x::Union(SymbolNode,GetfieldNode,Expr)) = copy(x)
astcopy(x::Array{Any,1}) = map(astcopy, x)
astcopy(x) = x

isequal(x::Expr, y::Expr) = (is(x.head,y.head) && isequal(x.args,y.args))
isequal(x::QuoteNode, y::QuoteNode) = isequal(x.value, y.value)

function show(io::IO, tv::TypeVar)
    if !is(tv.lb, None)
        show(io, tv.lb)
        print(io, "<:")
    end
    print(io, tv.name)
    if !is(tv.ub, Any)
        print(io, "<:")
        show(io, tv.ub)
    end
end

expand(x) = ccall(:jl_expand, Any, (Any,), x)
macroexpand(x) = ccall(:jl_macroexpand, Any, (Any,), x)

## misc syntax ##

macro eval(x)
    :($(esc(:eval))($(Expr(:quote,x))))
end

## some macro utilities ##

find_vars(e) = find_vars(e, {})
function find_vars(e, lst)
    if isa(e,Symbol)
        if !isdefined(e) || isconst(e)
            # exclude global constants
        else
            push!(lst, e)
        end
    elseif isa(e,Expr) && e.head !== :quote && e.head !== :top
        for x in e.args
            find_vars(x,lst)
        end
    end
    lst
end

# wrap an expression in "let a=a,b=b,..." for each var it references
localize_vars(expr) = localize_vars(expr, true)
function localize_vars(expr, esca)
    v = find_vars(expr)
    # requires a special feature of the front end that knows how to insert
    # the correct variables. the list of free variables cannot be computed
    # from a macro.
    if esca
        v = map(esc,v)
    end
    Expr(:localize, :(()->($expr)), v...)
end

max_stars_repo_name: tobig/julia

max_issues_repo_name: tobig/julia

max_forks_repo_name: tobig/julia


text:
module Turexas

include("tuRLing/env.jl")
include("tuRLing/agent.jl")

end # module


text:
using Distributions

# Checks output from chernikova function for some input matrix
# defines a specified finitely generated cone. That is,
# { x: Ax ≥ 0 } == { B y : y ≥ 0 }
#
# Arguments:
# A::Matrix - constraint matrix (rows correspond to individual constraints)
# B::Matrix - matrix of cones rays (columns correspond to individual rays)
function test_chernikova(A::Matrix{T1}, B::Matrix{T2}) where {T1<:Union{Integer, Rational}, T2<:Integer}
    rays = Cones.norm_cols(B)
    chern_rays = chernikova(A)
    @test Cones.check_columns_same(rays, chern_rays)
end

# Simple two rays example
println("\tExample 1...")
A = [-1//5 4//5;
     4//5 -1//5]
B = [1 4; 4 1]
test_chernikova(A,B)

# Redundant constraints
println("\tExample 2...")
A = [[1 1 0];
     [1 0 1];
     [1 1 1];
     [1 1 1];
     [0 1 1]]
B = [1 0 0;
     0 1 0;
     0 0 1]
test_chernikova(A,B)

# Empty set example
println("\tExample 3...")
A = [-1 -1 -1]
B = Array{Int}(undef, 3, 0)
test_chernikova(A, B)

# Last two examples are broken...

# One ray example
println("\tExample 4...")
A = [[1 -1]; [-1 1]]
B = reshape([1, 1], 2, 1)
test_chernikova(A, B)

# Example from paper
# "Algorithm for finding a general formula for the non-negative solutions of a system of linear inequalities"
# by Chernikova (1964)
println("\tExample 5...")
A = [[3 -4 1 0];
     [2 0 4 -1];
     [-4 -7 2 4];
     [-1 0 20 2];
     [6 -5 -4 2]]

B = [[1, 0, 0, 2] [0, 0, 1, 4] [0, 1, 4, 16] [0, 0, 2, 4]  [0, 32, 128, 336] [6, 0, 0, 6] [21, 12, 0, 42] [52, 40, 4, 120] [20, 0, 32, 4] [24080, 25760, 30800, 53760]]
test_chernikova(A, B)


text:
# To run this script, `cd` to the `./test/fixtures` directory and then, from the Julia terminal, `include("./runner.jl")`.

import JSON

function gen( x, name )
	y = Array( Float64, length( x ) );
	for i in eachindex(x)
		y[i] = erfcinv( x[i] );
	end

	data = Dict([
		("x", x),
		("expected", y)
	]);

	outfile = open( name, "w" );
	write( outfile, JSON.json(data) );
	close( outfile );
end

# 0.5 <= x <= 1.5 (linear regime)
x = linspace( 0.5, 1.5, 3000 );
gen( x, "x_0.5_1.5.json" );

# 0.25 < 2-x < 0.5
x = linspace( 1.5, 1.75, 500 );
gen( x, "x_1.5_1.75.json" );

# 0.25 < x < 0.5
x = linspace( 0.25, 0.5, 500 );
gen( x, "x_0.25_0.5.json" );

# 0.25 < 2-x < 0.00012340980408664937
x = linspace( 1.75, 1.9998, 500 );
gen( x, "x_1.75_1.9998.json" );

# 0.00012340980408664937 < x < 0.25
x = linspace( 0.0002, 0.25, 500 );
gen( x, "x_0.0002_0.25.json" );

# 0.00012340980408664937 < 2-x < 2.220446049250313e-16
x = linspace( 1.9998, 1.9999999999999998, 500 );
gen( x, "x_1.9998_1.9999..8.json" );

# 2.220446049250313e-16 < 2-x < 0
x = linspace( 1.9999999999999998, 2, 500 );
gen( x, "x_1.9999..8_2.json" );

max_stars_repo_name: math-io/erfcinv

max_issues_repo_name: math-io/erfcinv

max_forks_repo_name: math-io/erfcinv


text:
using Catalyst, Test
using LinearAlgebra

include("test_networks.jl")

rn = @reaction_network begin
    (1, 2), A + B <--> C
    (3, 2), D <--> E
    (0.1, 0.2), E <--> F
    6, F --> G
    7, H --> G
    5, 0 --> K
    3, A + B --> Z + C
end

S = netstoichmat(rn)
C = conservationlaws(S)
@test size(C,1) == 3
b = [ 0, 0, 0, 1, 1, 1, 1, 1, 0, 0 ]
@test any(b == C[i,:] for i in 1:size(C,1))

# For the A + B <--> C subsystem one of these must occur
# as a conservation law
D = [ 1 -1 0 0 0 0 0 0 0 0;
      -1 1 0 0 0 0 0 0 0 0
      0  1 1 0 0 0 0 0 0 0 ]
@test any(D[j,:] == C[i,:] for i in 1:size(C,1), j in 1:size(D,1))

C = conservationlaws(rn)
@test size(C,1) == 3
@test Catalyst.get_networkproperties(rn).nullity == 3
@test any(b == C[i,:] for i in 1:size(C,1))
@test any(D[j,:] == C[i,:] for i in 1:size(C,1), j in 1:size(D,1))

Cs_standard = map(conservationlaws, reaction_networks_standard)
@test all(size(C, 1) == 0 for C in Cs_standard)

Cs_hill = map(conservationlaws, reaction_networks_hill)
@test all(size(C, 1) == 0 for C in Cs_hill)

function consequiv(A, B)
    rank([A; B]) == rank(A) == rank(B)
end
Cs_constraint = map(conservationlaws, reaction_networks_constraint)
@test all(consequiv.(Matrix{Int}.(Cs_constraint), reaction_network_constraints))


text:
module HelloWorld

# package code goes here
function hw()
  println("Hello World")
end

end # module