text:
""" Take one time step with three forcing functions on u, v, w. """
function time_step_with_forcing_functions(arch)
    @inline Fu(x, y, z, t) = exp(π * z)
    @inline Fv(x, y, z, t) = cos(42 * x)
    @inline Fw(x, y, z, t) = 1.0

    grid = RegularCartesianGrid(size=(1, 1, 1), extent=(1, 1, 1))
    model = IncompressibleModel(grid=grid, architecture=arch, forcing=(u=Fu, v=Fv, w=Fw))
    time_step!(model, 1, euler=true)

    return true
end

@inline Fu_discrete_func(i, j, k, grid, clock, model_fields) = @inbounds -model_fields.u[i, j, k]
@inline Fv_discrete_func(i, j, k, grid, clock, model_fields, params) = @inbounds - model_fields.v[i, j, k] / params.τ
@inline Fw_discrete_func(i, j, k, grid, clock, model_fields, params) = @inbounds - model_fields.w[i, j, k]^2 / params.τ

""" Take one time step with a DiscreteForcing function. """
function time_step_with_discrete_forcing(arch)

    Fu = Forcing(Fu_discrete_func, discrete_form=true)

    grid = RegularCartesianGrid(size=(1, 1, 1), extent=(1, 1, 1))
    model = IncompressibleModel(grid=grid, architecture=arch, forcing=(u=Fu,))
    time_step!(model, 1, euler=true)

    return true
end

""" Take one time step with ParameterizedForcing forcing functions. """
function time_step_with_parameterized_discrete_forcing(arch)

    Fv = Forcing(Fv_discrete_func, parameters=(τ=60,), discrete_form=true)
    Fw = Forcing(Fw_discrete_func, parameters=(τ=60,), discrete_form=true)

    grid = RegularCartesianGrid(size=(1, 1, 1), extent=(1, 1, 1))
    model = IncompressibleModel(grid=grid, architecture=arch, forcing=(v=Fv, w=Fw))
    time_step!(model, 1, euler=true)

    return true
end

""" Take one time step with a Forcing forcing function with parameters. """
function time_step_with_parameterized_continuous_forcing(arch)

    u_forcing = Forcing((x, y, z, t, ω) -> sin(ω * x), parameters=π)

    grid = RegularCartesianGrid(size=(1, 1, 1), extent=(1, 1, 1))
    model = IncompressibleModel(grid=grid, architecture=arch, forcing=(u=u_forcing,))
    time_step!(model, 1, euler=true)

    return true
end

""" Take one time step with a Forcing forcing function with parameters. """
function time_step_with_single_field_dependent_forcing(arch, fld)

    forcing = NamedTuple{(fld,)}((Forcing((x, y, z, t, u) -> -u, field_dependencies=:u),))

    grid = RegularCartesianGrid(size=(1, 1, 1), extent=(1, 1, 1))
    model = IncompressibleModel(grid=grid, architecture=arch, forcing=forcing)
    time_step!(model, 1, euler=true)

    return true
end

""" Take one time step with a Forcing forcing function with parameters. """
function time_step_with_multiple_field_dependent_forcing(arch)

    u_forcing = Forcing((x, y, z, t, v, w, T) -> sin(v) * exp(w) * T, field_dependencies=(:v, :w, :T))

    grid = RegularCartesianGrid(size=(1, 1, 1), extent=(1, 1, 1))
    model = IncompressibleModel(grid=grid, architecture=arch, forcing=(u=u_forcing,))
    time_step!(model, 1, euler=true)

    return true
end

""" Take one time step with a Forcing forcing function with parameters. """
function time_step_with_parameterized_field_dependent_forcing(arch)

    u_forcing = Forcing((x, y, z, t, u, p) -> sin(p.ω * x) * u, parameters=(ω=π,), field_dependencies=:u)

    grid = RegularCartesianGrid(size=(1, 1, 1), extent=(1, 1, 1))
    model = IncompressibleModel(grid=grid, architecture=arch, forcing=(u=u_forcing,))
    time_step!(model, 1, euler=true)

    return true
end

function relaxed_time_stepping(arch)
    x_relax = Relaxation(rate = 1/60,   mask = GaussianMask{:x}(center=0.5, width=0.1),
                                      target = LinearTarget{:x}(intercept=π, gradient=ℯ))

    y_relax = Relaxation(rate = 1/60,   mask = GaussianMask{:y}(center=0.5, width=0.1),
                                      target = LinearTarget{:y}(intercept=π, gradient=ℯ))

    z_relax = Relaxation(rate = 1/60,   mask = GaussianMask{:z}(center=0.5, width=0.1),
                                      target = π)

    grid = RegularCartesianGrid(size=(1, 1, 1), extent=(1, 1, 1))
    model = IncompressibleModel(grid=grid, architecture=arch, forcing=(u=x_relax, v=y_relax, w=z_relax))
    time_step!(model, 1, euler=true)

    return true
end

@testset "Forcings" begin
    @info "Testing forcings..."

    for arch in archs
        @testset "Forcing function time stepping [$(typeof(arch))]" begin
            @info "  Testing forcing function time stepping [$(typeof(arch))]..."

            @testset "Non-parameterized forcing functions [$(typeof(arch))]" begin
                @info "      Testing non-parameterized forcing functions [$(typeof(arch))]..."
                @test time_step_with_forcing_functions(arch)
                @test time_step_with_discrete_forcing(arch)
            end

            @testset "Parameterized forcing functions [$(typeof(arch))]" begin
                @info "      Testing parameterized forcing functions [$(typeof(arch))]..."
                @test time_step_with_parameterized_continuous_forcing(arch)
                @test time_step_with_parameterized_discrete_forcing(arch)
            end

            @testset "Field-dependent forcing functions [$(typeof(arch))]" begin
                @info "      Testing field-dependent forcing functions [$(typeof(arch))]..."

                for fld in (:u, :v, :w, :T)
                    @test time_step_with_single_field_dependent_forcing(arch, fld)
                end

                @test time_step_with_multiple_field_dependent_forcing(arch)
                @test time_step_with_parameterized_field_dependent_forcing(arch)
            end

            @testset "Relaxation forcing functions [$(typeof(arch))]" begin
                @info "      Testing relaxation forcing functions [$(typeof(arch))]..."
                @test relaxed_time_stepping(arch)
            end
        end
    end
end


text:
import Deformula: _deint, deformulaZeroToInf, deformulaMinusOneToOne, deint
using Test

@testset "Deformula.jl" begin
    result = _deint(deformulaZeroToInf) do x
        0.5*x^-0.5*exp(-x^0.5)
    end
    @test result.s ≈ 1.0
    result = _deint(deformulaMinusOneToOne) do x
        y = x + 1
        0.5*y^-0.5*exp(-y^0.5)
    end
    @test result.s ≈ 0.7568830631028768
end

@testset "Deformula2" begin
    result = _deint(deformulaZeroToInf) do x
        0.5*x^-0.5*exp(-x^0.5)
    end
    @test sum(result.w) * result.h ≈ 1.0
end

@testset "Deformula3" begin
    result = deint(0.0, Inf64) do x
        0.5*x^-0.5*exp(-x^0.5)
    end
    @test sum(result.w) * result.h ≈ 1.0
end

@testset "Deformula4" begin
    result = deint(0.0, 10.0) do x
        if 0 <= x <= 10
            1
        else
            0
        end
    end
    println(result.x)
    @test result.s ≈ 10.0
end

@testset "Deformula5" begin
    result = deint(1.0, Inf64) do x
        1.0/sqrt(2.0*pi) * exp(-(x - 1.0)^2/2.0)
    end
    println(result.x)
    @test result.s ≈ 0.5
end

@testset "Deformula6" begin
    dfm(x; m, lambda) = lambda * x^(m-1)* exp(-lambda * x^m) + lambda * m * x^(m-1) * log(x) * exp(-lambda * x^m) + lambda * m * x^(m-1) * exp(-lambda * x^m) * (-lambda * x^m * log(x))
    result = deint(0.0, Inf64) do x
        dfm(x, m=2.0, lambda=1.0)
    end
    println(length(result.x))
    println(result.s)
    # @test result.s ≈ 0.5
end

ext: jl


text:
# """
# Keep enough history for the given lag indices.  Intended for time-series.
# Given lag list [0, 2, 10] for series {yₜ}, we care about storing/indexing [yₜ, yₜ₋₂, yₜ₋₁₀]
# Note that we must keep all data points from (t) --> (t - iₙ), so data storage scales
# with the largest lag.
# """
# type Window{T <: Real, VECI <: AVec{Int}} <: OnlineStat
#     buf::CircularBuffer{T}
#     lagIndices::VECI
#     n::Int
# end
#
# function Window{T<:Real}(::Type{T}, lagIndices::AVec{Int})
#     @assert minimum(lagIndices) >= 0
#     Window(CircularBuffer(T, maximum(lagIndices)+1), lagIndices, 0)
# end
#
# # replicates a typical circular buffer
# Window{T<:Real}(::Type{T}, sz::Int) = Window(T, sz-1:-1:0)
#
# # ------------------------------
#
# statenames(o::Window) = [:lags, :nobs]
# state(o::Window) = Any[lags(o), nobs(o)]
#
# lags{T}(o::Window{T}) = T[x for x in o]
# Base.convert(::Type{Array}, o::Window) = lags(o)
#
# function update!{T<:Real}(o::Window{T}, x::Real)
#     push!(o.buf, convert(T, x))
#     o.n += 1
#     return
# end
#
# function Base.empty!{T<:Real}(o::Window{T})
#     o.buf = CircularBuffer(T, maximum(o.lagIndices)+1)
#     o.n = 0
#     return
# end
#
# # ------------------------------
#
# _bufferIndex(o::Window, i::Int) = length(o.buf) - o.lagIndices[i]
# bufferIndex(o::Window, i::Int) = (i < 1 || i > length(o)) ? error("Idx $i out of range. ", o) : _bufferIndex(o, i)
#
#
# Base.getindex(o::Window, i::Int) = o.buf[bufferIndex(o, i)]
# Base.unsafe_getindex(o::Window, i::Int) = o.buf[_bufferIndex(o, i)]
#
# function Base.setindex!{T}(o::Window, data::T, i::Int)
#   o.buf[bufferIndex(o, i)] = data
#   nothing
# end
# function Base.unsafe_setindex!{T}(o::Window, data::T, i::Int)
#   o.buf[_bufferIndex(o, i)] = data
#   nothing
# end
#
# # iteration
# Base.start(o::Window) = 1
# Base.done(o::Window, state::Int) = state > length(o)
# Base.next(o::Window, state::Int) = (o[state], state+1)
#
# # note: length is 0 until the underlying buffer is full
# Base.length(o::Window) = (isfull(o) ? capacity(o) : 0)
# Base.size(o::Window) = (length(o),)
#
# QuickStructs.capacity(o::Window) = length(o.lagIndices)
# QuickStructs.isfull(o::Window) = isfull(o.buf)
#
# # Base.push!{T}(o::Window{T}, data::T) = push!(o.buf, data)
# # Base.append!{T}(o::Window{T}, datavec::AbstractVector{T}) = append!(o.buf, datavec)


text:
"""
Uses `bwbench()` to measure the memory bandwidth for an increasing number of threads (`1:max_nthreads`).
"""
function bwscaling(; max_nthreads = Threads.nthreads(), kwargs...)
    max_nthreads ≤ Threads.nthreads() || throw(ArgumentError("max_nthreads must be ≤ Threads.nthreads()."))
    sdaxpy_results = zeros(max_nthreads)
    for n in 1:max_nthreads
        df = bwbench(; nthreads=n, kwargs...)
        sdaxpy_results[n] = last(df.var"Rate (MB/s)")
        println("$n Thread(s): SDaxpy Bandwidth (MB/s) is ", round(sdaxpy_results[n], digits=2))
        flush(stdout)
    end

    # print results
    println("\n\nScaling results:")
    data = hcat(1:max_nthreads, sdaxpy_results)
    pretty_table(data; header=["# Threads", "SDaxpy Bandwidth (MB/s)"])
    return data
end

function flopsscaling(; max_nthreads=Threads.nthreads())
    triad_results = zeros(max_nthreads)
    for n in 1:max_nthreads
        # bench
        N = 120_000_000
        alignment = 64
        a = allocate(Float64, N, alignment)
        b = allocate(Float64, N, alignment)
        c = allocate(Float64, N, alignment)
        scalar = 3.0

        Nperthread = floor(Int, N / n)
        rest = rem(N, n)
        thread_indices = collect(Iterators.partition(1:N, Nperthread))

        times = zeros(10)
        for k in 1:10
            times[k] = @elapsed triad_kernel(a, b, c, scalar; nthreads=n, thread_indices=thread_indices)
        end

        mintime = minimum(times)
        flops = 2 * N
        flop_rate = 1.0e-06 * flops / mintime
        # bench end

        triad_results[n] = flop_rate
        println("$n Thread(s): Triad Performance (MFlop/s) is ", round(triad_results[n], digits=2))
        flush(stdout)
    end

    # print results
    println("\n\nScaling results:")
    data = hcat(1:max_nthreads, triad_results)
    pretty_table(data; header=["# Threads", "Triad Performance (MFlop/s)"])
    return data
end

max_stars_repo_name: carstenbauer/BandwidthBenchmark

max_issues_repo_name: carstenbauer/BandwidthBenchmark

max_forks_repo_name: carstenbauer/BandwidthBenchmark


text:
# Tests for the Transaction system.
# Author: Noah C. Benson <n@nben.net>

@testset "TX" begin

    @testset "tx1" begin
        vals = Int[]
        vols = Vector{Volatile{Int}}[]
        log = Actor{typeof(stdout)}(stdout)
        worker() = begin
            global vals
            global vols
            while vols[end][1][] < 0
                # Pick a random layer and cell to calculate
                l = rand((2,3,4,5))
                layer = vols[l]
                k = rand(1:length(layer))
                # synconize for the rest.
                tx() do
                    #println("$l $k $(current_tx() === nothing) $(getfield(_current_tx, :values)))")
                    res = layer[k]
                    if res[] < 0
                        up = vols[l - 1]
                        a = up[k*2 - 1]
                        b = up[k*2]
                        if a[] >= 0 && b[] >= 0
                            res[] = a[] + b[]
                            send(log) do log
                                println("Layer $l, item $k complete $(objectid(current_task())).")
                                log
                            end
                        end
                    end
                end
            end
            send(log) do log
                println("Ending Woker $(objectid(current_task())).")
                log
            end
        end
        tx1_run(n::Int) = begin
            global vals
            global vols
            vals = rand(0:1000, 16)
            vols = Vector{Volatile{Int}}[
                [Volatile{Int}(u)   for u in vals],
                [Volatile{Int}(-1) for _ in 1:8],
                [Volatile{Int}(-1) for _ in 1:4],
                [Volatile{Int}(-1) for _ in 1:2],
                [Volatile{Int}(-1)]]
            # Spawn n threads.
            threads = [(Threads.@spawn worker()) for _ in 1:n]
            for th in threads
                Threads.wait(th)
            end
            println("  ---->  $(vols[end][1][]) == $(sum(vals))")
            return (vols[end][1][] == sum(vals))
        end
        @test tx1_run(1)
        @test tx1_run(2)
        @test tx1_run(3)
        @test tx1_run(4)
        @test tx1_run(5)
        @test tx1_run(6)
        @test tx1_run(7)
        @test tx1_run(8)
        @test tx1_run(9)
        @test tx1_run(10)
        @test tx1_run(11)
    end

end


text:
module NDSpline
import CxxWrap

CxxWrap.wrap_module(joinpath(dirname(@__FILE__),"../deps/src/libspline"))

function (fun::Array{Spline})(x::Array)
    g(f_i) = f_i(x)[1]
  return g.(fun)
end

function (fun::Array{Spline})(x::Real)
    g(f_i) = f_i(x)
  return g.(fun)
end

function fit(s::Spline,X::Array{Float64,2},y::Vector{Float64},k::Vector{Int},ordered::Bool)
  fit(s,X,size(X)[2],y,k,ordered)
end

export Spline,fit

end


text:
using Base.Test
using JuMP
using MIPVerify
using MIPVerify: check_size, increment!
isdefined(:TestHelpers) || include("../../TestHelpers.jl")
using TestHelpers: get_new_model

@testset "conv2d.jl" begin

    @testset "Conv2d" begin
        @testset "Base.show" begin
            filter = rand(3, 3, 2, 5)
            p = Conv2d(filter)
            io = IOBuffer()
            Base.show(io, p)
            @test String(take!(io)) == "Conv2d(2, 5, kernel_size=(3, 3), stride=(1, 1), padding=same)"
        end
        @testset "With Bias" begin
            @testset "Matched Size" begin
                out_channels = 5
                filter = rand(3, 3, 2, out_channels)
                bias = rand(out_channels)
                p = Conv2d(filter, bias)
                @test p.filter == filter
                @test p.bias == bias
            end
            @testset "Unmatched Size" begin
                filter_out_channels = 4
                bias_out_channels = 5
                filter = rand(3, 3, 2, filter_out_channels)
                bias = rand(bias_out_channels)
                @test_throws AssertionError Conv2d(filter, bias)
            end
        end
        @testset "No Bias" begin
            filter = rand(3, 3, 2, 5)
            p = Conv2d(filter)
            @test p.filter == filter
        end
        @testset "JuMP Variables" begin
            m = Model()
            filter_size = (3, 3, 2, 5)
            filter = map(_ -> @variable(m), CartesianRange(filter_size))
            p = Conv2d(filter)
            @test p.filter == filter
        end
        @testset "check_size" begin
            filter = rand(3, 3, 2, 5)
            p = Conv2d(filter)
            @test check_size(p, (3, 3, 2, 5)) === nothing
            @test_throws AssertionError check_size(p, (3, 3, 2, 4))
        end
    end

    @testset "increment!" begin
        @testset "Real * Real" begin
            @test 7 == increment!(1, 2, 3)
        end
        @testset "JuMP.AffExpr * Real" begin
            m = get_new_model()
            x = @variable(m, start=100)
            y = @variable(m, start=1)
            s = 5*x+3*y
            t = 3*x+2*y
            increment!(s, 2, t)
            @test getvalue(s)==1107
            increment!(s, t, -1)
            @test getvalue(s)==805
            increment!(s, x, 3)
            @test getvalue(s)==1105
            increment!(s, y, 2)
            @test getvalue(s)==1107
        end
    end

    @testset "conv2d" begin
        srand(1)
        input_size = (1, 4, 4, 2)
        input = rand(0:5, input_size)
        filter_size = (3, 3, 2, 1)
        filter = rand(0:5, filter_size)
        bias_size = (1, )
        bias = rand(0:5, bias_size)
        true_output_raw = [
            49 74 90 56;
            67 118 140 83;
            66 121 134 80;
            56 109 119 62
        ]
        true_output = reshape(transpose(true_output_raw), (1, 4, 4, 1))
        p = Conv2d(filter, bias)
        @testset "Numerical Input, Numerical Layer Parameters" begin
            evaluated_output = MIPVerify.conv2d(input, p)
            @test evaluated_output == true_output
        end
        @testset "Numerical Input, Variable Layer Parameters" begin
            m = get_new_model()
            filter_v = map(_ -> @variable(m), CartesianRange(filter_size))
            bias_v = map(_ -> @variable(m), CartesianRange(bias_size))
            p_v = Conv2d(filter_v, bias_v)
            output_v = MIPVerify.conv2d(input, p_v)
            @constraint(m, output_v .== true_output)
            solve(m)

            p_solve = MIPVerify.Conv2d(getvalue(filter_v), getvalue(bias_v))
            solve_output = MIPVerify.conv2d(input, p_solve)
            @test solve_output≈true_output
        end
        @testset "Variable Input, Numerical Layer Parameters" begin
            m = get_new_model()
            input_v = map(_ -> @variable(m), CartesianRange(input_size))
            output_v = MIPVerify.conv2d(input_v, p)
            @constraint(m, output_v .== true_output)
            solve(m)

            solve_output = MIPVerify.conv2d(getvalue(input_v), p)
            @test solve_output≈true_output
        end
    end
end


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

"""
    abstract type DifferentiationAlgorithm

*Experimental feature, not part of stable public API.*

Abstract type for integrated autocorrelation length estimation algorithms.
"""
abstract type DifferentiationAlgorithm end
export DifferentiationAlgorithm

MaybeDiffAlg = Union{DifferentiationAlgorithm,Missing}

function vjp_algorithm end
vjp_algorithm(::Any) = missing

function promote_vjp_algorithm end
promote_vjp_algorithm() = missing
promote_vjp_algorithm(a::MaybeDiffAlg) = a
promote_vjp_algorithm(a::DifferentiationAlgorithm, b::DifferentiationAlgorithm) = a
promote_vjp_algorithm(a::DifferentiationAlgorithm, b::Missing) = a
promote_vjp_algorithm(a::Missing, b::DifferentiationAlgorithm) = b
promote_vjp_algorithm(a::Missing, b::Missing) = b
function promote_vjp_algorithm(a::MaybeDiffAlg, b::MaybeDiffAlg, cs::Vararg{MaybeDiffAlg,N}) where {N}
    promote_vjp_algorithm(promote_vjp_algorithm(a, b), cs...)
end

generic_vjp_algorithm(obj::T) where {T} = promote_vjp_algorithm(map(vjp_algorithm, object_contents(obj))...)
vjp_algorithm(f::F) where {F<:Function} = generic_vjp_algorithm(f)
vjp_algorithm(d::AbstractDensity) = generic_vjp_algorithm(d)

# ToDo: Switch default from ForwardDiff to Zygote:
default_vjp_algorithm() = ForwardDiffAD()
#default_vjp_algorithm() = ZygoteAD()

@inline function some_vjp_algorithm(x)
    alg = vjp_algorithm(x)
    ismissing(alg) ? default_vjp_algorithm() : alg
end

function jvp_algorithm end
jvp_algorithm(::Any) = missing

function promote_jvp_algorithm end
promote_jvp_algorithm() = missing
promote_jvp_algorithm(a::MaybeDiffAlg) = a
promote_jvp_algorithm(a::DifferentiationAlgorithm, b::DifferentiationAlgorithm) = a
promote_jvp_algorithm(a::DifferentiationAlgorithm, b::Missing) = a
promote_jvp_algorithm(a::Missing, b::DifferentiationAlgorithm) = b
promote_jvp_algorithm(a::Missing, b::Missing) = b
function promote_jvp_algorithm(a::MaybeDiffAlg, b::MaybeDiffAlg, cs::Vararg{MaybeDiffAlg,N}) where {N}
    promote_jvp_algorithm(promote_jvp_algorithm(a, b), cs...)
end

generic_jvp_algorithm(obj::T) where {T} = promote_jvp_algorithm(map(jvp_algorithm, object_contents(obj))...)
jvp_algorithm(f::F) where {F<:Function} = generic_jvp_algorithm(f)

default_jvp_algorithm() = ForwardDiffAD()

function some_jvp_algorithm(x)
    alg = jvp_algorithm(x)
    ismissing(alg) ? default_jvp_algorithm() : alg
end

"""
    valgradof(f::Function, algorithm::DifferentiationAlgorithm  = vjp_algorithm(f))

*Experimental feature, not part of stable public API.*

Generates a function that calculates both value and gradient of `f` at given
points. The differentiation algorithm used depends on `f`.

The function `f` must support `ValueShapes.varshape(f)` and
`ValueShapes.unshaped(f)`.

Returns function `valgrad_f` with

```julia
f_at_x, grad_of_f_at_x = valgrad_f(x)

grad_of_f_at_x = zero(x)
f_at_x = valgrad_f(!, grad_of_f_at_x, x)
```
"""
function valgradof end
export valgradof

abstract type GradientFunction <: Function end


text:
@dict_readable struct ParameterInformationCapabilities <: Outbound
    labelOffsetSupport::Union{Bool,Missing}
end

@dict_readable struct SignatureInformationCapabilities <: Outbound
    documentationFormat::Union{Vector{String},Missing}
    parameterInformation::Union{ParameterInformationCapabilities,Missing}
end

@dict_readable struct SignatureHelpClientCapabilities <: Outbound
    dynamicRegistration::Union{Bool,Missing}
    signatureInformation::Union{SignatureInformationCapabilities,Missing}
    contextSupport::Union{Bool,Missing}
end

const SignatureHelpTriggerKind = Int
const SignatureHelpTriggerKinds = (Invoked = 1,
                                   TriggerCharacter = 2,
                                   ContentChange = 3)

struct SignatureHelpOptions <: Outbound
    triggerCharacters::Union{Vector{String},Missing}
    retriggerCharacters::Union{Vector{String},Missing}
end

struct SignatureHelpRegistrationOptions <: Outbound
    documentSelector::Union{DocumentSelector,Nothing}
    triggerCharacters::Union{Vector{String},Missing}
    retriggerCharacters::Union{Vector{String},Missing}
end

struct ParameterInformation <: Outbound
    label::Union{String,Tuple{Int,Int}}
    documentation::Union{String,MarkupContent,Missing}
end

struct SignatureInformation <: Outbound
    label::String
    documentation::Union{String,MarkedString,Missing}
    parameters::Union{Vector{ParameterInformation},Missing}
end

struct SignatureHelp <: Outbound
    signatures::Vector{SignatureInformation}
    activeSignature::Union{Int,Missing}
    activeParameter::Union{Int,Missing}
end

@dict_readable struct SignatureHelpContext <: Outbound
    triggerKind::SignatureHelpTriggerKind
    triggerCharacter::Union{String,Missing}
    isRetrigger::Bool
    activeSignatureHelp::Union{SignatureHelp,Missing}
end

@dict_readable struct SignatureHelpParams <: Outbound
    textDocument::TextDocumentIdentifier
    position::Position
    context::Union{SignatureHelpContext,Missing}
    workDoneToken::Union{Int,String,Missing} # ProgressToken
end


text:
function find_vortices3D(psi, X)
    x = X[1]; y = X[2];
    Nx = length(x); Ny = length(y);
    if length(X) == 3
        z = X[3];
        Nz = length(z);

        for i = 1:Nz
            psi_z = Torus(psi[:,:,i], x, y);

            if i == 1
                vort_array = vortex_array(findvortices(psi_z));
                v = z[i]*ones(length(vort_array[:, 1]));
                vort_array = hcat(vort_array, v);
                vort_array = view(vort_array, :, [1,2,4,3]);
            else
                temp = vortex_array(findvortices(psi_z));
                v = z[i]*ones(length(temp[:, 1]));
                temp = hcat(temp, v);
                temp = view(temp, :, [1, 2, 4, 3]);
                vort_array = vcat(vort_array, temp);
            end
        end
    else

        vort_array = vortex_array(findvortices(psi));
    end
    temp = zeros(length(vort_array[:, 1]))
    vort_array = hcat(vort_array, temp);
    return vort_array;
end

function find_vortices3D_v2(psi, X)
    x = X[1]; y = X[2];
    z = X[3];
    vorts_3d = Array{Float64, 2}[]
    for zidx in 1:length(z)
        vorts = vortex_array(findvortices(Torus(psi[:, :, zidx], x, y)));
        if (length(vorts) == 0)
            push!(vorts_3d, vorts)
        else
            temp = zeros(length(vorts[:, 1]))
            zval = z[zidx].*ones(length(vorts[:, 1]))
            vorts = hcat(vorts, temp);
            vorts = hcat(vorts, zval);
            vorts =view(vorts, :, [1, 2, 5, 3, 4])
            push!(vorts_3d, vorts);
        end
    end
    return vorts_3d;
end

function find_vortices3D_v2(psi)
    @unpack ψ,x,y,z = psi
    vorts_3d = Array{Float64, 2}[]
    for zidx in 1:length(z)
        vorts = vortex_array(findvortices(Torus(ψ[:, :, zidx], x, y)));
        if (length(vorts) == 0)
            # push!(vorts_3d, vorts)
        else
            temp = zeros(length(vorts[:, 1]))
            zval = z[zidx].*ones(length(vorts[:, 1]))
            vorts = hcat(vorts, temp);
            vorts = hcat(vorts, zval);
            vorts =view(vorts, :, [1, 2, 5, 3, 4])
            push!(vorts_3d, vorts);
        end
    end
    return vorts_3d;
end

# Need to pass X and define x = X[1] etc
function vortex_boundary_bottom(v,X, dx, dy, dz)
    x = X[1]; y = X[2]; z = X[3];
    at_boundary = true;

    if ((v[1, 1] - x[1] > dx) && (x[end] - v[1, 1] > dx)) # Then not at x Boundary
        if ((v[1, 2] - y[1] > dy) && (y[end] - v[1, 2] > dy)) # " y boundary
            if ((v[1, 3] - z[1] > dz) && (z[end] - v[1, 3] > dz))
                at_boundary = false;
            end
        end
    end
    return at_boundary;
end

function vortex_boundary_top(v, X, dx, dy, dz)
    x = X[1]; y = X[2]; z = X[3];
    at_boundary = true;
    if ((v[end, 1] - x[1] > dx) && (x[end] - v[end, 1] > dx)) # Then not at x Boundary
        if ((v[end, 2] - y[1] > dy) && (y[end] - v[end, 2] > dy)) # " y boundary
            if ((v[end, 3] - z[1] > dz) && (z[end] - v[end, 3] > dz))
                at_boundary = false;
            end
        end
    end
    return at_boundary;
end

function vortex_link_bottom(vidx, vorts_label, X)
    x = X[1]; y = X[2]; z = X[3];
    dx = x[2]-x[1]; dy = y[2]-y[1]; dz = z[2]-z[1];
    v = vorts_label[vidx];
    for i in 1:length(vorts_label)
        vi = vorts_label[i];
        if vi[1,5] != v[1,5] # different label
            # Find matching z level
            index = 0;
            for j in 1:length(vi[:, 1])
                if vi[j, 3] == v[1, 3]
                    index = j;
                end
            end
            # If matching z level exists
            if index != 0
                # If z level is at bottom of vi
                if vi[index, 3] == vi[1, 3]
                    vorts_label[i][:, 5] .= vorts_label[vidx][1, 5]
                    if !vortex_boundary_top(vi, X, dx, dy, dz)
                        vortex_link_top(i, vorts_label, X);
                    end
                elseif vi[index, 3] == vi[end, 3]
                    vorts_label[i][:, 5] .= vorts_label[vidx][1, 5]
                    if !vortex_boundary_bottom(vi,X, dx, dy, dz)
                        vortex_link_bottom(i, vorts_label, X);
                    end
                end
            end
        end
    end
end

function vortex_link_top(vidx, vorts_label, X)
    x = X[1]; y = X[2]; z = X[3];
    dx = x[2]-x[1]; dy = y[2]-y[1]; dz = z[2]-z[1];
    v = vorts_label[vidx];
    for i in 1:length(vorts_label)
        vi = vorts_label[i];
        if vi[1,5] != v[1,5] # different label
            # Find matching z level
            index = 0;
            for j in 1:length(vi[:, 1])
                if vi[j, 3] == v[end, 3]
                    index = j;
                end
            end
            # If matching z level exists
            if index != 0
                # If z level is at bottom of vi
                if vi[index, 3] == vi[1, 3]
                    vorts_label[i][:, 5] .= vorts_label[vidx][1, 5]
                    if !vortex_boundary_top(vi,X, dx, dy, dz)
                        vortex_link_top(i, vorts_label, X);
                    end
                elseif vi[index, 3] == vi[end, 3]
                    vorts_label[i][:, 5] .= vorts_label[vidx][1, 5]

                    if !vortex_boundary_bottom(vi,X, dx, dy, dz)
                        vortex_link_bottom(i, vorts_label, X);
                    end
                end
            end
        end
    end
end

function sort_vorts_label(vorts_label)
    NUM_VORTS = length(vorts_label)
    temp = Array{Float64, 2}[]
    i = 1;
    for i in 1:length(vorts_label)
        for j in 1:NUM_VORTS

            if (vorts_label[i][1, 5] == j)
               #println(j)
                push!(temp, vorts_label[i]);
            end
        end
    end
    return temp;
end

function psi_itp_2D(psi_itp, zval, X)
    x = X[1]; y = X[2];
    Nx = length(x); Ny = length(y);
    psi_itp2D = im.*zeros(Nx, Ny);
    for i in 1:Nx
        for j in 1:Ny
                psi_itp2D[i, j] = psi_itp(x[i], y[j], zval);
        end
    end
    return psi_itp2D;
end


text:
using Compat, FastGroupBy
import Base: Forward
import SortingAlgorithms: RADIX_MASK, RADIX_SIZE, uint_mapping
using BenchmarkTools

using SortingAlgorithms, Base.Order, Compat, IndexedTables, DataFrames
import Base: ht_keyindex, rehash!, _setindex!
import SortingAlgorithms: uint_mapping, RADIX_SIZE, RADIX_MASK
import DataFrames: DataFrame, AbstractDataFrame
import IndexedTables: IndexedTable, column
import PooledArrays.PooledArray
import CategoricalArrays.CategoricalArray
import SplitApplyCombine.groupreduce

using StatsBase

const N = 1_000_000
const K = 100
by = rand(Int32(1):Int32(N/K), N);
val = similar(by);
T = Int

using DataStructures

# define a heap that invovles the byval and S
struct ByVal{T<:Integer, S <: Number}
    index::Integer
    by::T
    val::S
end

import Base.isless
isless(x::ByVal,y::ByVal) = x.index < y.index

"WAAAAAAAYYY too slow uses less memory by building a heap"
function sumby_lessmem_heap{T, S<:Number}(by::AbstractVector{T},  val::AbstractVector{S})::Dict{T,S}
    o = Forward
    lo = 1
    hi = length(by)

    if lo >= hi;  return by;  end

    # Make sure we're sorting a bits type
    if !isbits(T)
      error("the by vector must be of bits types (e.g. Float64, Integer) got $T")
    end

    # Init
    iters = ceil(Integer, sizeof(T)*8/RADIX_SIZE)
    bin = zeros(UInt32, 2^RADIX_SIZE, iters)

    # Histogram for each element, radix
    for i = lo:hi
        v = uint_mapping(o, by[i])
        for j = 1:iters
            idx = @compat(Int((v >> (j-1)*RADIX_SIZE) & RADIX_MASK)) + 1
            @inbounds bin[idx,j] += 1
        end
    end

    # Sort!
    len = hi-lo+1
    for j = 1:iters
        # Unroll first data iteration, check for degenerate case
        v = uint_mapping(o, by[hi])
        idx = @compat(Int((v >> (j-1)*RADIX_SIZE) & RADIX_MASK)) + 1

        # are all values the same at this radix?
        if bin[idx,j] == len;  continue;  end

        cbin = cumsum(bin[:,j])
        ci = cbin[idx]

        # ByValHeap -- keeps the data in a heap
        bvheap = binary_maxheap(ByVal)

        # if ci == hi then do nothing
        if ci < hi
          push!(bvheap, ByVal(ci, by[ci], val[ci]))
          by[ci] = by[hi]
          val[ci] = val[hi]
        end
        cbin[idx] -= 1

        loopto = lo-1
        i = hi - 1 # i is used to keep track of where we are up to

        # Finish the loop...
        #   @inbounds for i in hi-1:-1:lo
        while i > 0
          if isempty(bvheap)
              loopto = lo-1
          else
              loopto = top(bvheap).index
          end
          while i > loopto
              v = uint_mapping(o, by[i])
              idx = @compat(Int((v >> (j-1)*RADIX_SIZE) & RADIX_MASK)) + 1
              ci = cbin[idx]
              if ci > i # it's ok to just put it in place
                  by[ci] = by[i]
                  val[ci] = val[i]
              elseif ci < i # keep the value in a heap for retrieval
                  push!(bvheap, ByVal(ci, by[ci], val[ci]))
                  by[ci] = by[i]
                  val[ci] = val[i]
                  # make the loop end sooner
                  if ci > loopto; loopto = ci; end
              end

              cbin[idx] -= 1
              i -= 1
          end
          # if it got here without being 0 then it must be up to a heap value
          if i != 0
              tmp = pop!(bvheap)
              v = uint_mapping(o, tmp.by)
              idx = @compat(Int((v >> (j-1)*RADIX_SIZE) & RADIX_MASK)) + 1
              ci = cbin[idx]
              if ci > i # it's ok to just put it in place
                  by[ci] = tmp.by
                  val[ci] = tmp.val
              elseif ci < i # keep the value in a heap for retrieval
                  push!(bvheap, ByVal(ci, by[ci], val[ci]))
                  by[ci] = tmp.by
                  val[ci] = tmp.val
              end

              cbin[idx] -= 1
              i -= 1
          end
        end
    end

    sumby_contiguous(by, val)
end

sumby_lessmem_heap([1:10...],[1:10]

function abc()
  by1, val1 = copy(by), copy(val)
  @elapsed sumby_lessmem_heap(by1, val1)
end

function def()
  by1, val1 = copy(by), copy(val)
  @elapsed sumby(by1, val1)
end

a1 = mean([abc() for i = 1:5][2:end])
b1 = mean([def() for i = 1:5][2:end])
1 - a1/b1

function ghi()
  by1, val1 = copy(by), copy(val)
  @elapsed sumby_lessmem_chain(by1, val1)
end
c1 = mean([ghi() for i = 1:5][2:end])
1 - c1/b1


text:
using CUDNN: cudnnConvolutionBackwardBias, cudnnConvolutionBackwardFilter, cudnnConvolutionBackwardData, cudnnConvolutionForward
using CUDNN: CD, CUDNN_CROSS_CORRELATION, CUDNN_CONVOLUTION
using CUDArt

# x = reshape(Float64[1:9],3,3,1,1); tx = CudaArray(x)
# w = reshape(Float64[1:4],2,2,1,1); tw = CudaArray(w)
pads = (0,0); strides = (2,2)
corrd = CD(padding=pads, stride=strides, mode=CUDNN_CROSS_CORRELATION)
x = rand(Float64[1:10;],5,5,1,1); tx = CudaArray(x)
w = rand(Float64[1:10;],2,2,1,1); tw = CudaArray(w)
@show size(x)
@show size(w)
@show pads
@show strides
ty = cudnnConvolutionForward(tx, tw;convDesc=corrd); y = to_host(ty)
@show size(y)
dy = rand(Float64[1:10;],size(ty)); tdy = CudaArray(dy)
# dy = reshape(Float64[5 6 7 8],2,2,1,1); tdy = CudaArray(dy) #2x2
println(x)
println(w)
println("dy:")
println(dy)
dx = zeros(x); tdx = CudaArray(dx)
dw = zeros(w); tdw = CudaArray(dw)

println("corr dw:")
cudnnConvolutionBackwardFilter(tx, tdy, tdw; convDesc=corrd)
println(to_host(tdw))
println("corr dx:")
cudnnConvolutionBackwardData(tw, tdy, tdx; convDesc=corrd)
println(to_host(tdx))

println("testing pooling..")
using CUDNN: cudnnPoolingForward, cudnnPoolingBackward, cudnnGetPoolingNdForwardOutputDim, PD, CUDNN_POOLING_MAX
# x = reshape(Float64[5 2 3 4 5 6 7 7 3 1 4 5], 3, 4, 1, 1); tx = CudaArray(x)
x = rand(Float64[1:10;],10,10,1,1); tx = CudaArray(x)
println("x:"); println(x)
pd2 = PD(2,3,0,1, CUDNN_POOLING_MAX)
ydims = cudnnGetPoolingNdForwardOutputDim(pd2, tx)
@show ydims


text:
#=
solves:
- Julia version:
- Author: qd4314
- Date: 2019-03-21
=#

function linFunction(phiVector::AbstractVector, KS::KineticSolver)

    # reshape vector input to tensor
    phi = reshape(phiVector, (KS.DecompNumberElements, KS.nx, KS.ny))

    OSigma = KS.OSigma
    M = KS.DecompM

    # compute sweeping rhs which is O*Sigma*phi
    OSigmaPhi = zeros(KS.nquadpoints, KS.nx, KS.ny)
    @einsum OSigmaPhi[q, i, j] = OSigma[q, moments] * phi[moments, i, j]

    # multiply with sigmaS
    @inbounds for j in KS.rangex, i in KS.rangey
        OSigmaPhi[:, i, j] = OSigmaPhi[:, i, j] .* KS.SigmaS[i, j]
    end

    # perform sweeping
    Y = sweeping(KS, deepcopy(KS.psiC), OSigmaPhi, KS.ScatteringKernelConvolution)

    #Y= sweeping(KS,0.0*deepcopy(KS.psiC),OSigmaPhi,KS.ScatteringKernelConvolution);

    KS.psiC .= Y

    # multiply with M
    MY = zeros(KS.DecompNumberElements, KS.nx, KS.ny)
    @einsum MY[m, i, j] = M[m, q] * Y[q, i, j]

    # substract from phi to get A*phi output
    Aphi = phi .- MY

    # reshape for vector output
    AphiVector = reshape(Aphi, prod((KS.DecompNumberElements, KS.nx, KS.ny)))

    return AphiVector

end

function corysolve!(KS::KineticSolver)

    ### DEBUGGING
    # turn off scattering and absorption
    #KS.SigmaS = zeros(size(KS.SigmaS));
    #KS.SigmaT = zeros(size(KS.SigmaT));
    #KS.convolutionmagnitude = 0.0;
    close("all")
    io = open("mytestfile.txt", "a")

    writedlm(io, "start!!!!!!!!")
    close(io)
    # set convolution magnitude

    # set max number of sweeping steps
    KS.nSweepingSteps = 1000 # 1 is one-shot
    KS.epsSweeping = 1e-3#1e-7;

    # update SigmaT if convolution is active

    if KS.ictype == "linesource"
        println("Problem Type is ", KS.ictype)
        sigmaA = 0.0 # hard coded for linesource
    else
        sigmaA = KS.SigmaT - KS.SigmaS
    end
    KS.SigmaT =
        KS.SigmaS +
        sigmaA .* ones(size(KS.SigmaT)) +
        KS.convolutionmagnitude * ones(size(KS.SigmaT))

    # set estimate of contraction number
    L = KS.convolutionmagnitude / (1 / KS.dt + minimum(KS.SigmaT) + sqrt(2) / KS.dx)
    KS.contractionEstimate = L / (1 - L)

    # set moment order
    order = 0
    KS.DecompNumberElements = (order + 1) * (order + 1)

    # define total number of unknowns
    Nentries = KS.DecompNumberElements * KS.ny * KS.nx

    # decomposition order of kernel is Norder = (order+1)*(order+1)+1
    KS.DecompO, KS.DecompSigma, KS.DecompM = kerneldecomposition(KS, order) # Henyey Grenstein
    # precompute convolution matrix if needed
    if KS.convolutionmagnitude > 0
        KS.ScatteringKernelConvolution = computeScatteringKernelConvolution(KS)
    else
        KS.ScatteringKernelConvolution = zeros(KS.nquadpoints, KS.nquadpoints)
    end# precompute O*Sigma
    KS.OSigma = KS.DecompO * KS.DecompSigma

    # write Decomposition matrices for physical scattering
    M = KS.DecompM
    Sigma = KS.DecompSigma
    O = KS.DecompO
    OSigma = KS.OSigma

    # Set initial condition
    psi = deepcopy(KS.psi0)
    psi0 = deepcopy(KS.psi0)
    KS.lastpsi = deepcopy(KS.psi0)
    KS.psiC = deepcopy(KS.psi0)

    phiNew = zeros(KS.DecompNumberElements, KS.ny, KS.nx)
    @einsum phiNew[m, i, j] = M[m, q] * psi[q, i, j]

    rho = computeMoments(KS, psi)
    println("Initial mass = $(sum(rho[:])).")

    # define scattering kernel for source iteration
    S = deepcopy(KS.ScatteringKernelConvolution)
    for q = 1:size(S)[1]
        S[q, :] .*= KS.Q.weights[:]
    end

    KS.kernelTimesWeightsSparse = sparse(S)

    @showprogress 0.01 "Computing... " for k = 1:KS.nt

        # adapt dt if k*dt would be larger than tEnd
        dt = ifelse(k * KS.dt < KS.tEnd, KS.dt, KS.tEnd - k * KS.dt)

        # compute rhs
        rhsquads = sweeping(
            KS,
            deepcopy(KS.lastpsi),
            KS.lastpsi ./ KS.dt + KS.Source,
            KS.ScatteringKernelConvolution,
        )

        rhsmoms = zeros(KS.DecompNumberElements, KS.ny, KS.nx)
        @einsum rhsmoms[m, i, j] = M[m, q] * rhsquads[q, i, j]

        ## reshape to vector
        rrhsmoms = reshape(rhsmoms, Nentries)

        # define linear Map for Krylov solver
        D = LinearMap((x) -> linFunction(x, KS), Nentries, Nentries; ismutating = false) # muss in jeder Zeit iteration neu definiert werden, da KS sich aendert!

        # initial guess for Krylov solver
        #x0 = zeros(KS.DecompNumberElements,KS.nx,KS.ny)
        #@einsum x0[m,i,j] = M[m,q]*psi0[q,i,j]
        #x0 = reshape(x0,Nentries)

        # solve system with GMRES
        #println("############ starting krylov")
        println(sqrt(eps(real(eltype(rrhsmoms)))))
        phiNew = IterativeSolvers.gmres!(rrhsmoms, D, rrhsmoms, verbose = false)
        #println("done krylov")

        #phiNew = IterativeSolvers.gmres!(x0,D,rrhsmoms,verbose=false)
        phiNew = reshape(phiNew, (KS.DecompNumberElements, KS.nx, KS.ny))

        #visMoments(KS,phiNew,"OutputGMRES")

        # determine PsiNew
        OSigmaPhi = zeros(KS.nquadpoints, KS.nx, KS.ny)
        @einsum OSigmaPhi[q, i, j] = OSigma[q, moments] * phiNew[moments, i, j]

        # multiply with sigmaS
        @inbounds for j in KS.rangex, i in KS.rangey
            OSigmaPhi[:, i, j] = OSigmaPhi[:, i, j] .* KS.SigmaS[i, j]
        end

        # perform sweeping to update psi
        psiNew = sweeping(
            KS,
            deepcopy(KS.lastpsi),
            KS.Source .+ KS.lastpsi / KS.dt + OSigmaPhi,
            KS.ScatteringKernelConvolution,
        )
        KS.lastpsi .= psiNew
    end

    rho = computeMoments(KS, KS.lastpsi)
    println("Final mass = $(sum(rho[:])).")

    vis(KS, KS.lastpsi)
    return rho
end

function solveWithSweeping!(KS::KineticSolver)

    # set max number of sweeping steps
    KS.nSweepingSteps = 1000
    KS.epsSweeping = 1e-3#1e-7;

    # update SigmaT if convolution is active
    if KS.ictype == "linesource"
        println("Problem Type is ", KS.ictype)
        sigmaA = 0.0 # hard coded for linesource
    else
        sigmaA = KS.SigmaT - KS.SigmaS
    end
    KS.SigmaT =
        KS.SigmaS +
        sigmaA .* ones(size(KS.SigmaT)) +
        KS.convolutionmagnitude * ones(size(KS.SigmaT))

    # set estimate of contraction number
    L =
        (KS.convolutionmagnitude + maximum(KS.SigmaS)) /
        (1 / KS.dt + minimum(KS.SigmaT) + sqrt(2) / KS.dx)
    KS.contractionEstimate = L / (1 - L)

    # compute artificial scattering kernel
    KS.ScatteringKernelConvolution = computeScatteringKernelConvolution(KS) # Convo

    # Set initial condition
    psi = zeros(KS.nquadpoints, KS.ny, KS.nx)
    psi = deepcopy(KS.psi0)
    KS.lastpsi = deepcopy(KS.psi0)
    rho = computeMoments(KS, psi)
    println("Initial mass = $(sum(rho[:])).")

    nicewrite(
        string(KS.outputfolderprefix, "data/source.txt"),
        computeMoments(KS, KS.Source),
    )

    # define scattering kernel for source iteration
    asKernel = deepcopy(KS.ScatteringKernelConvolution) # artificial inscattering
    physKernel = ones(size(asKernel)) ./ 4.0 / pi # add physical inscattering
    g = 0.01
    HG = computeScatteringKernelHenyeyGreenstein(KS, g)
    # decomposition order of kernel is Norder = (order+1)*(order+1)+1
    KS.DecompO, KS.DecompSigma, KS.DecompM = kerneldecomposition(KS, 1) # Henyey Grenstein

    KS.kernelTimesWeightsSparse = sparse(asKernel)

    #println("HG = ",HG)
    #println("OSigmaM = ",KS.DecompO*KS.DecompSigma*KS.DecompM)
    #println("wQ = ",KS.Q.weights[:])
    for q = 1:size(physKernel)[1]
        asKernel[q, :] .*= KS.Q.weights[:]
        physKernel[q, :] .*= KS.Q.weights[:]
    end
    physKernel = KS.DecompO * KS.DecompSigma * KS.DecompM

    @showprogress 0.01 "Computing... " for k = 1:KS.nt

        # adapt dt if k*dt would be larger than tEnd
        dt = ifelse(k * KS.dt < KS.tEnd, KS.dt, KS.tEnd - k * KS.dt)

        # first euler step
        psiNew = sweepingWithKernels(
            KS,
            deepcopy(KS.lastpsi),
            KS.Source .+ KS.lastpsi / KS.dt,
            asKernel,
            physKernel,
        )
        KS.lastpsi .= psiNew
        #vis(KS,psi)
        #return;
    end
    vis(KS, KS.lastpsi) # do all post computational plotting
    rho = computeMoments(KS, KS.lastpsi)
    println("End mass = $(sum(rho[:])).")
end # solve

function solveRSNWithSweeping!(KS::KineticSolver)

    # set max number of sweeping steps
    KS.nSweepingSteps = 1000
    KS.epsSweeping = 1e-3#1e-7;

    # update SigmaT if convolution is active
    if KS.ictype == "linesource"
        println("Problem Type is ", KS.ictype)
        sigmaA = 0.0 # hard coded for linesource
    else
        sigmaA = KS.SigmaT - KS.SigmaS
    end
    KS.SigmaT =
        KS.SigmaS +
        sigmaA .* ones(size(KS.SigmaT)) +
        KS.convolutionmagnitude * ones(size(KS.SigmaT))

    # set estimate of contraction number
    L =
        (KS.convolutionmagnitude + maximum(KS.SigmaS)) /
        (1 / KS.dt + minimum(KS.SigmaT) + sqrt(2) / KS.dx)
    KS.contractionEstimate = L / (1 - L)

    # compute artificial scattering kernel
    KS.ScatteringKernelConvolution = computeScatteringKernelConvolution(KS) # Convo

    # Set initial condition
    psi = zeros(KS.nquadpoints, KS.ny, KS.nx)
    psi = deepcopy(KS.psi0)
    KS.lastpsi = deepcopy(KS.psi0)
    rho = computeMoments(KS, psi)
    println("Initial mass = $(sum(rho[:])).")

    nicewrite(
        string(KS.outputfolderprefix, "data/source.txt"),
        computeMoments(KS, KS.Source),
    )

    # define scattering kernel for source iteration
    asKernel = deepcopy(KS.ScatteringKernelConvolution) # artificial inscattering
    physKernel = ones(size(asKernel)) ./ 4.0 / pi # add physical inscattering
    g = 0.01
    HG = computeScatteringKernelHenyeyGreenstein(KS, g)
    # decomposition order of kernel is Norder = (order+1)*(order+1)+1
    KS.DecompO, KS.DecompSigma, KS.DecompM = kerneldecomposition(KS, 1) # Henyey Grenstein

    KS.kernelTimesWeightsSparse = sparse(asKernel)

    #println("HG = ",HG)
    #println("OSigmaM = ",KS.DecompO*KS.DecompSigma*KS.DecompM)
    #println("wQ = ",KS.Q.weights[:])
    for q = 1:size(physKernel)[1]
        asKernel[q, :] .*= KS.Q.weights[:]
        physKernel[q, :] .*= KS.Q.weights[:]
    end
    physKernel = KS.DecompO * KS.DecompSigma * KS.DecompM

    @showprogress 0.01 "Computing... " for k = 1:KS.nt

        # adapt dt if k*dt would be larger than tEnd
        dt = ifelse(k * KS.dt < KS.tEnd, KS.dt, KS.tEnd - k * KS.dt)

        # first euler step
        psiNew = sweepingWithKernels(
            KS,
            deepcopy(KS.lastpsi),
            KS.Source .+ KS.lastpsi / KS.dt,
            asKernel,
            physKernel,
        )
        KS.lastpsi .= psiNew

        # perform rotation and interpolation step, or convolution step, or other modifications
        methods!(KS, KS.lastpsi, psiNew, k)
        #vis(KS,psi)
        #return;
    end
    vis(KS, KS.lastpsi) # do all post computational plotting
    rho = computeMoments(KS, KS.lastpsi)
    println("End mass = $(sum(rho[:])).")
end # solve

function solve!(KS::KineticSolver)
    #close("all")
    psi = zeros(KS.nquadpoints, KS.ny, KS.nx)
    psiLastTimeStep = zeros(KS.nquadpoints, KS.ny, KS.nx)
    flux = zeros(KS.nquadpoints, KS.ny, KS.nx)
    psi = deepcopy(KS.psi0) # Set initial condition
    rho = computeMoments(KS, psi)
    println("Initial mass = $(sum(rho[:])).")

    #nicewrite(string(KS.outputfolderprefix,"data/rho0.txt"),computeMoments(KS,psi))
    nicewrite(
        string(KS.outputfolderprefix, "data/source.txt"),
        computeMoments(KS, KS.Source),
    )

    summarizeAsPlot(KS, psi) #close("all")

    # precompute convolution matrix if needed
    if KS.convolutionmagnitude > 0
        KS.ScatteringKernelConvolution = computeScatteringKernelConvolution(KS)
    else
        KS.ScatteringKernelConvolution = zeros(KS.nquadpoints, KS.nquadpoints)
    end

    # low rank computation?
    if KS.rank > 0
        dA = zeros(KS.nquadpoints, (KS.ny) * (KS.nx))
        psireshaped = reshape(psi, KS.nquadpoints, KS.ny * KS.nx)
        s0 = svd(psireshaped)
        U0 = s0.U[1:end, 1:KS.rank]
        S0 = diagm(0 => s0.S[1:KS.rank])
        V0 = s0.V[1:end, 1:KS.rank]
        U1 = zeros(size(U0))
        V1 = zeros(size(V0))
        S1 = zeros(size(S0))
        K1 = zeros(size(U0 * S0))
        L1 = zeros(size(V0 * S1'))
        Y0 = zeros(KS.nquadpoints, (KS.ny) * (KS.nx))
    end

    @showprogress 0.01 "Computing... " for k = 1:KS.nt

        # force periodicity (if wanted)
        forcePeriodicity!(KS, psi)

        # adapt dt if k*dt would be larger than tEnd
        dt = ifelse(k * KS.dt < KS.tEnd, KS.dt, KS.tEnd - k * KS.dt)

        # store old psi to average later on
        #@inbounds for q=1:length(psi);	psiLastTimeStep[q] = psi[q]	end
        psiLastTimeStep .= psi

        # first euler step
        solveFlux!(KS, psi, flux)
        Euler!(KS, psi, flux, dt)

        # second euler step
        solveFlux!(KS, psi, flux)
        Euler!(KS, psi, flux, dt)

        #average
        psi .= (psi .+ psiLastTimeStep) ./ 2.0

        # add source contribution
        psi .+= KS.Source .* KS.dt

        # Note: So far, low rank and rSN DO NOT WORK TOGETHER. Either lowrank or rSN.
        if KS.rank > 0
            @inbounds for q = 1:length(psi)
                dA[q] = psi[q] - psiLastTimeStep[q]
            end
            K1 = U0 * S0 + dA * V0
            U1, S1 = qr(K1)
            U1 = Matrix(U1)
            S1 = S1 - U1' * dA * V0
            L1 = V0 * S1' + dA' * U1
            V1, S1 = qr(L1)
            V1 = Matrix(V1)
            S1 = S1'
            Y1 = U1 * S1 * V1'
            for i = 1:length(Y0), Y0[i] in Y1[i]
            end
            for i = 1:length(U0), U0[i] in U1[i]
            end
            for i = 1:length(V0), V0[i] in V1[i]
            end
            for i = 1:length(S0), S0[i] in S1[i]
            end
            for i = 1:length(psi), psi[i] in Y0[i]
            end
        else
            # perform rotation and interpolation step, or convolution step, or other modifications
            methods!(KS, psi, psiLastTimeStep, k)
        end
    end
    vis(KS, psi) # do all post computational plotting
end # solve

max_stars_repo_name: vavrines/CSN

max_issues_repo_name: vavrines/CSN

max_forks_repo_name: vavrines/CSN


text:
# Author: Youngjun Kim, youngjun@stanford.edu
# Date: 06/16/2014

module SimpleTCASSensorImpl

export
    initialize,
    update,

    updateSensor,

    SimpleTCASSensor,
    SimpleTCASSRState,
    SimpleTCASSensorInput,
    SimpleTCASSensorOutput

using AbstractSensorImpl
using AbstractSensorInterfaces
using CommonInterfaces

using Base.Test

import CommonInterfaces.initialize
import CommonInterfaces.update
import AbstractSensorInterfaces.updateSensor

type SimpleTCASSRState

    t::Float64
    x::Float64
    y::Float64
    h::Float64
    vx::Float64
    vy::Float64
    vh::Float64
end

type SimpleTCASSensorInput

    states::Vector{SimpleTCASSRState}
end

type SimpleTCASSensorOutput

    t::Float64

    r::Float64
    r_d::Float64

    a::Float64
    a_d::Float64

    h::Vector{Float64}
    h_d::Vector{Float64}
end

type SimpleTCASSensor <: AbstractSensor

  aircraft_number::Int
  states::Vector{SimpleTCASSRState}
  output::SimpleTCASSensorOutput

  function SimpleTCASSensor(aircraft_number::Int)
    obj = new()

    obj.aircraft_number = aircraft_number
    obj.output = SimpleTCASSensorOutput(0.0,0.0,0.0,0.0,0.0,Float64[],Float64[])

    return obj
  end
end

function updateSensor(sr::SimpleTCASSensor, input::SimpleTCASSensorInput)

  sr.states = input.states

  if sr.aircraft_number == 1
    intruder_number = 2
  else
    intruder_number = 1
  end

  t1, x1, y1, h1, vx1, vy1, vh1 = sr.states[sr.aircraft_number].t, sr.states[sr.aircraft_number].x, sr.states[sr.aircraft_number].y, sr.states[sr.aircraft_number].h, sr.states[sr.aircraft_number].vx, sr.states[sr.aircraft_number].vy, sr.states[sr.aircraft_number].vh

  t2, x2, y2, h2, vx2, vy2, vh2 = sr.states[intruder_number].t, sr.states[intruder_number].x, sr.states[intruder_number].y, sr.states[intruder_number].h, sr.states[intruder_number].vx, sr.states[intruder_number].vy, sr.states[intruder_number].vh

  #@test t1 == t2

  t = t1

  dxy = [(x2 - x1), (y2 - y1)]
  dvxy = [(vx2 - vx1), (vy2 - vy1)]
  r = norm(dxy)
  r_d = dot(dxy,dvxy) / norm(dxy)

  a = abs(h1 - h2)
  a_d = sign(h2 - h1) * (vh2 - vh1)

  h = [h1, h2]
  h_d = [vh1, vh2]

  sr.output.t = t
  sr.output.r = r
  sr.output.r_d = r_d
  sr.output.a = a
  sr.output.a_d = a_d
  sr.output.h = h
  sr.output.h_d = h_d

  return sr.output
end

update(sr::SimpleTCASSensor, input) = update(sr,convert(SimpleTCASSensorInput, input))

update(sr::SimpleTCASSensor, input::SimpleTCASSensorInput) = updateSensor(sr, input)

function initialize(sr::SimpleTCASSensor)

end

end

max_stars_repo_name: sisl/RLES-SISLES

max_issues_repo_name: sisl/RLES-SISLES

max_forks_repo_name: sisl/RLES-SISLES


text:
module Position
using JSON3, HTTP

export listPositions, listOpenPositions, getPosition, closePosition, closePositionFull

"""
Individual data for positions (A Long or Short)

# Fields
- 'pl': Profit / Loss
- 'resettablePL': Profit / Loss since last reset
- 'units': Number of units in the trade
- 'unrealizedPL': Unrealized profit / loss in the posistion
"""
mutable struct posData
    pl # The profit / Loss
    resettablePL # Profit / Loss since last reset
    units # Number of units in the trade
    unrealizedPL # Unrealized profit / loss of position

    posData() = new()
end

"""
Detailed Position struct from Oanda

# Fields
- 'instrument': The instrument of the position
- 'long::posData': Data for longs on the position
- 'short::posData': Data for shorts on the position
- 'pl': Overall profit / loss for the position
- 'resettablePL': Profit / Loss since last reset
- 'unrealizedPL': Unrealized profit / loss in the posistion
"""
mutable struct position
    instrument # The instrument
    # I have both of the following to enable people with hedging enabled
    long::posData # Data for a Long trade on the position
    short::posData # Data for a Short trade on the position
    pl # The profit / loss for this position
    resettablePL # The profit / loss since last reset for this position
    unrealizedPL # The total unrealised profit / loss for this position

    position() = new()
end

"Top layer for positions endpoint"
mutable struct positionTopLayer
    positions::Vector{position}

    positionTopLayer() = new()
end

"Top layer for positions endpoint"
mutable struct positionTopLayerSingle
    position::position

    positionTopLayerSingle() = new()
end

"Coerce a given Position into its proper types (Used internally)"
function coercePos(pos::position)
    pos.long.pl = parse(Float32, pos.long.pl)
    pos.long.resettablePL = parse(Float32, pos.long.resettablePL)
    pos.long.units = parse(Int32, pos.long.units)
    pos.long.unrealizedPL = parse(Float32, pos.long.unrealizedPL)

    pos.short.pl = parse(Float32, pos.short.pl)
    pos.short.resettablePL = parse(Float32, pos.short.resettablePL)
    pos.short.units = parse(Int32, pos.short.units)
    pos.short.unrealizedPL = parse(Float32, pos.short.unrealizedPL)

    pos.pl = parse(Float32, pos.pl)
    pos.resettablePL = parse(Float32, pos.resettablePL)
    pos.unrealizedPL = parse(Float32, pos.unrealizedPL)

    return pos
end

# Declaring JSON3 struct types
JSON3.StructType(::Type{position}) = JSON3.Mutable()
JSON3.StructType(::Type{positionTopLayer}) = JSON3.Mutable()
JSON3.StructType(::Type{positionTopLayerSingle}) = JSON3.Mutable()
JSON3.StructType(::Type{posData}) = JSON3.Mutable()

"""
    listPositions(config)

Returns a list of current positions

Returns an object of type 'Vector{position}'

# Arguments
- 'config::config': A valid config object
"""
function listPositions(config)
    r = HTTP.request(
        "GET",
        string(
            "https://",
            config.hostname,
            "/v3/accounts/",
            config.account,
            "/positions",
        ),
        [
         "Authorization" => string("Bearer ", config.token),
         "Accept-Datetime-Format" => config.datetime,
        ],
    )
    data = JSON3.read(r.body, positionTopLayer)
    data = data.positions

    temp = Vector{position}()
    for pos in data
        pos = coercePos(pos)
        push!(temp, pos)
    end

    return temp
end

"""
    listOpenPositions(config)

Returns a list of current positions that have an open trade

Returns an object of 'Vector{position}'

# Arguments
- 'config::config': A valid config object
"""
function listOpenPositions(config)
    r = HTTP.request(
        "GET",
        string(
            "https://",
            config.hostname,
            "/v3/accounts/",
            config.account,
            "/openPositions",
        ),
        [
         "Authorization" => string("Bearer ", config.token),
         "Accept-Datetime-Format" => config.datetime,
        ],
    )
    data = JSON3.read(r.body, positionTopLayer)
    data = data.positions

    temp = Vector{position}()
    for pos in data
        pos = coercePos(pos)
        push!(temp, pos)
    end

    return temp
end

"""
    getPosition(config, instrument)

Returns position data for a specified instrument

Returns an object of type 'position'

# Arguments
- 'config::config': A valid config object
- 'instrument': The instrument to get info from
"""
function getPosition(config, instrument)
    r = HTTP.request(
        "GET",
        string(
            "https://",
            config.hostname,
            "/v3/accounts/",
            config.account,
            "/positions/",
            instrument,
        ),
        [
         "Authorization" => string("Bearer ", config.token),
         "Accept-Datetime-Format" => config.datetime,
        ],
    )
    data = JSON3.read(r.body, positionTopLayerSingle)

    data = coercePos(data.position)

    return data
end

"""
    closePositionFull(config, instrument, long=true)

Closes a position completely

Returns true on success

# Arguments
- 'config::config': A valid config object
- 'instrument': instrument to close
- 'long::Bool': True to close long false to close short
"""
function closePositionFull(config, instrument, long::Bool=true)
    data = ""
    if long
        data = "{\"longUnits\": \"ALL\"}"
    else
        data = "{\"shortUnits\": \"ALL\"}"
    end
    r = HTTP.request(
        "PUT",
        string(
            "https://",
            config.hostname,
            "/v3/accounts/",
            config.account,
            "/positions/",
            instrument,
            "/close"
        ),
        [
         "Authorization" => string("Bearer ", config.token),
         "Accept-Datetime-Format" => config.datetime,
          "Content-Type" => "application/json",
        ],
        data,
    )

    return true
end

"""
    closePosition(config, instrument, LongUnits=NONE, ShortUnits=NONE)

Closes a positions units based on input

Returns true on success

# Arguments
- 'config::config': A valid config object
- 'instrument': The instrument to act on
- 'longUnits': The number of long units to close
- 'shortUnits': The number of short units to close
"""
function closePosition(config, instrument, longUnits="NONE", shortUnits="NONE")
    r = HTTP.request(
        "PUT",
        string(
            "https://",
            config.hostname,
            "/v3/accounts/",
            config.account,
            "/positions/",
            instrument,
            "/close"
        ),
        [
         "Authorization" => string("Bearer ", config.token),
         "Accept-Datetime-Format" => config.datetime,
          "Content-Type" => "application/json",
        ],
        string("{\"longUnits\": \"", longUnits,"\",\n\"shortUnits\": \"", shortUnits,"\"}"),
    )

    return true
end

end

max_stars_repo_name: CalebDepatie/Julianda

max_forks_repo_name: CalebDepatie/Julianda


text:
isdefined(Base, :__precompile__) && __precompile__()

module CUDArt

using CUDAdrv

using Compat
import Compat: UTF8String, ASCIIString

export
    # pointer symbols
    CudaPtr, rawpointer, CUDA_NULL,
    # array symbols
    AbstractCudaArray, AbstractCudaVector, AbstractCudaMatrix,
    CudaArray, CudaVector, CudaMatrix, CudaVecOrMat,
    CudaPitchedArray, HostArray, to_host,
    # other symbols
    device, devices, device_reset, attribute, capability,
    driver_version, runtime_version,
    pitchel, pitchbytes,
    device_synchronize, synchronize,
    Stream, null_stream, cudasleep,
    destroy, free, cudafinalizer, wait_free

import Base: ==, -, +, getindex, setindex!
import Base: length, size, ndims, eltype, similar, pointer, stride,
    copy, convert, reinterpret, show, summary,
    copy!, get!, fill!, launch, wait, unsafe_convert, vec

# Prepare the CUDA runtime API bindings
include("libcudart-6.5.jl")
import .CUDArt_gen
const rt = CUDArt_gen

# To load PTX code, we also need access to the driver API module utilities
if is_windows()
    const libcuda = Libdl.find_library(["nvcuda"], [""])
else # linux or mac
    const libcuda = Libdl.find_library(["libcuda"], ["/usr/lib/", "/usr/local/cuda/lib"])
end

if isempty(libcuda)
    error("CUDA driver API library cannot be found")
end

include("version.jl")
include("types.jl")
include("device.jl")
include("stream.jl")
#include("event.jl")
include("pointer.jl")
include("arrays.jl")
include("execute.jl")

if isdefined(Base, :__precompile__)
    include("precompile.jl")
    _precompile_()
end

function __init__()
    c_async_send_cudastream[] = cfunction(async_send_cudastream, Void, (rt.cudaStream_t, rt.cudaError_t, Ptr{Void}))
end

end


text:
"""
    newdims = Dims_create(nnodes::Integer, dims)

A convenience function for selecting a balanced Cartesian grid of a total of `nnodes`
nodes, for example to use with [`MPI.Cart_create`](@ref).

`dims` is an array or tuple of integers specifying the number of nodes in each dimension.
The function returns an array `newdims` of the same length, such that if `newdims[i] =
dims[i]` if `dims[i]` is non-zero, and `prod(newdims) == nnodes`, and values `newdims` are
as close to each other as possible.

`nnodes` should be divisible by the product of the non-zero entries of `dims`.

# External links
$(_doc_external("MPI_Dims_create"))
"""
function Dims_create(nnodes::Integer, dims)
    dims = Cint[dim for dim in dims]
    @mpichk ccall((:MPI_Dims_create, libmpi), Cint,
                  (Cint, Cint, Ptr{Cint}), nnodes, length(dims), dims)
    return dims
end

"""
    comm_cart = Cart_create(comm::Comm, dims; periodic=map(_->false, dims), reorder=false)

Create new MPI communicator with Cartesian topology information attached.

`dims` is an array or tuple of integers specifying the number of MPI processes in each
coordinate direction, and `periodic` is an array or tuple of `Bool`s indicating the
periodicity of each coordinate. `prod(dims)` must be less than or equal to the size of
`comm`; if it is smaller than some processes are returned a null communicator.

If `reorder == false` then the rank of each process in the new group is identical to its
rank in the old group, otherwise the function may reorder the processes.

See also [`MPI.Dims_create`](@ref).

# External links
$(_doc_external("MPI_Cart_create"))
"""
function Cart_create(comm::Comm, dims; periodic = map(_->false, dims), reorder=false)
    if !(dims isa Array{Cint})
        dims = Cint[dim for dim in dims]
    end
    if !(periodic isa Array{Cint})
        periodic = Cint[flag for flag in periodic]
    end
    length(dims) == length(periodic) || error("dims and periodic arguments are not the same length")
    comm_cart = Comm()
    # int MPI_Cart_create(MPI_Comm comm_old, int ndims, const int dims[],
    #                     const int periods[], int reorder, MPI_Comm *comm_cart)
    @mpichk ccall((:MPI_Cart_create, libmpi), Cint,
                  (MPI_Comm, Cint, Ptr{Cint}, Ptr{Cint}, Cint, Ptr{MPI_Comm}),
                  comm, length(dims), dims, periodic, reorder, comm_cart)
    if comm_cart != COMM_NULL
        finalizer(free, comm_cart)
    end
    comm_cart
end

"""
    rank = Cart_rank(comm::Comm, coords)

Determine process rank in communicator `comm` with Cartesian structure.  The `coords`
array specifies the 0-based Cartesian coordinates of the process. This is the inverse of [`MPI.Cart_coords`](@ref)

# External links
$(_doc_external("MPI_Cart_rank"))
"""
function Cart_rank(comm::Comm, coords)
    if !(coords isa Vector{Cint})
        coords = Cint[coord for coord in coords]
    end
    rank = Ref{Cint}()
    # int MPI_Cart_rank(MPI_Comm comm, const int coords[], int *rank)
    @mpichk ccall((:MPI_Cart_rank, libmpi), Cint,
                  (MPI_Comm, Ptr{Cint}, Ptr{Cint}),
                  comm, coords, rank)
    Int(rank[])
end

"""
    dims, periods, coords = Cart_get(comm::Comm)

Obtain information on the Cartesian topology of dimension `N` underlying the
communicator `comm`. This is specified by two `Cint` arrays of `N` elements
for the number of processes and periodicity properties along each Cartesian dimension.
A third `Cint` array is returned, containing the Cartesian coordinates of the calling process.

# External links
$(_doc_external("MPI_Cart_get"))
"""
function Cart_get(comm::Comm)
    maxdims = Cartdim_get(comm)
    # preallocate with nontrivial values
    dims    = Cint[-1 for i = 1:maxdims]
    periods = Cint[-1 for i = 1:maxdims]
    coords  = Cint[-1 for i = 1:maxdims]
    # int MPI_Cart_get(MPI_Comm comm, int maxdims, int dims[], int periods[], int coords[])
    @mpichk ccall((:MPI_Cart_get, libmpi), Cint,
                  (MPI_Comm, Cint, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}),
                  comm, maxdims, dims, periods, coords)
    return dims, periods, coords
end

"""
    ndims = Cartdim_get(comm::Comm)

Return number of dimensions of the Cartesian topology associated with the communicator `comm`.

# External links
$(_doc_external("MPI_Cartdim_get"))
"""
function Cartdim_get(comm::Comm)
    dims = Ref{Cint}()
    # int MPI_Cartdim_get(MPI_Comm comm, int *ndims)
    @mpichk ccall((:MPI_Cartdim_get, libmpi), Cint,
                  (MPI_Comm, Ptr{Cint}),
                  comm, dims)
    return Int(dims[])
end

"""
    coords = Cart_coords(comm::Comm, rank::Integer=Comm_rank(comm))

Determine coordinates of a process with rank `rank` in the Cartesian communicator
`comm`. If no `rank` is provided, it returns the coordinates of the current process.

Returns an integer array of the 0-based coordinates. The inverse of [`Cart_rank`](@ref).

# External links
$(_doc_external("MPI_Cart_coords"))
"""
function Cart_coords(comm::Comm, rank::Integer=Comm_rank(comm))
    maxdims = Cartdim_get(comm)
    coords = Vector{Cint}(undef, maxdims)
    # int MPI_Cart_coords(MPI_Comm comm, int rank, int maxdims, int coords[])
    @mpichk ccall((:MPI_Cart_coords, libmpi), Cint,
                  (MPI_Comm, Cint, Cint, Ptr{Cint}),
                  comm, rank, maxdims, coords)
    return Int.(coords)
end

"""
    rank_source, rank_dest = Cart_shift(comm::Comm, direction::Integer, disp::Integer)

Return the source and destination ranks associated to a shift along a given
direction.

# External links
$(_doc_external("MPI_Cart_shift"))
"""
function Cart_shift(comm::Comm, direction::Integer, disp::Integer)
    rank_source = Ref{Cint}()
    rank_dest   = Ref{Cint}()
    # int MPI_Cart_shift(MPI_Comm comm, int direction, int disp,
    #                    int *rank_source, int *rank_dest)
    @mpichk ccall((:MPI_Cart_shift, libmpi), Cint,
                  (MPI_Comm, Cint, Cint, Ptr{Cint}, Ptr{Cint}),
                  comm, direction, disp, rank_source, rank_dest)
    Int(rank_source[]), Int(rank_dest[])
end

"""
    comm_sub = Cart_sub(comm::Comm, remain_dims)

Create lower-dimensional Cartesian communicator from existent Cartesian
topology.

`remain_dims` should be a boolean vector specifying the dimensions that should
be kept in the generated subgrid.

# External links
$(_doc_external("MPI_Cart_sub"))
"""
function Cart_sub(comm::Comm, remain_dims)
    if !(remain_dims isa Array{Cint})
        remain_dims = Cint[dim for dim in remain_dims]
    end
    comm_sub = Comm()
    # int MPI_Cart_sub(MPI_Comm comm, const int remain_dims[], MPI_Comm *comm_new)
    @mpichk ccall((:MPI_Cart_sub, libmpi), Cint,
                  (MPI_Comm, Ptr{Cint}, Ptr{MPI_Comm}),
                  comm, remain_dims, comm_sub)
    if comm_sub != COMM_NULL
        finalizer(free, comm_sub)
    end
    comm_sub
end

struct Unweighted
end
Base.cconvert(::Type{Ptr{Cint}}, ::Unweighted) = Consts.MPI_UNWEIGHTED[]
const UNWEIGHTED = Unweighted()

struct WeightsEmpty
end
Base.cconvert(::Type{Ptr{Cint}}, ::WeightsEmpty) = Consts.MPI_WEIGHTS_EMPTY[]
const WEIGHTS_EMPTY = WeightsEmpty()

"""
    graph_comm = Dist_graph_create_adjacent(comm::Comm, sources::Vector{Cint}, destinations::Vector{Cint}; source_weights::Union{Vector{Cint}, Unweighted, WeightsEmpty}=UNWEIGHTED, destination_weights::Union{Vector{Cint}, Unweighted, WeightsEmpty}=UNWEIGHTED, reorder=false, infokws...)

Create a new communicator from a given directed graph topology, described by local incoming and outgoing edges on an existing communicator.

# Arguments
- `comm::Comm`: The communicator on which the distributed graph topology should be induced.
- `sources::Vector{Cint}`: The local, incoming edges on the rank of the calling process.
- `destinations::Vector{Cint}`: The local, outgoing edges on the rank of the calling process.
- `source_weights::Union{Vector{Cint}, Unweighted, WeightsEmpty}=MPI.UNWEIGHTED`: The edge weights of the local, incoming edges.
- `destinations_weights::Union{Vector{Cint}, Unweighted, WeightsEmpty}=MPI.UNWEIGHTED`: The edge weights of the local, outgoing edges.
- `reorder::Bool=false`: If set true, then the MPI implementation can reorder the source and destination indices.

# Example
We can generate a ring graph `1 --> 2 --> ... --> N --> 1`, where N is the number of ranks in the communicator, as follows
```julia
julia> rank = MPI.Comm_rank(comm);
julia> N = MPI.Comm_size(comm);
julia> sources = Cint[mod(rank-1, N)];
julia> destinations = Cint[mod(rank+1, N)];
julia> graph_comm = Dist_graph_create_adjacent(comm, sources, destinations);
```

# External links
$(_doc_external("MPI_Dist_graph_create_adjacent"))
"""
function Dist_graph_create_adjacent(comm::Comm, sources::Vector{Cint}, destinations::Vector{Cint}; source_weights::Union{Vector{Cint}, Unweighted, WeightsEmpty}=UNWEIGHTED, destination_weights::Union{Vector{Cint}, Unweighted, WeightsEmpty}=UNWEIGHTED, reorder=false, infokws...)
    graph_comm = Comm()
    # int MPI_Dist_graph_create_adjacent(MPI_Comm comm_old,
    #       int indegree, const int sources[],
    #       const int sourceweights[],
    #       int outdegree, const int destinations[],
    #       const int destweights[],
    #       MPI_Info info, int reorder, MPI_Comm *comm_dist_graph)
    @mpichk ccall((:MPI_Dist_graph_create_adjacent, libmpi), Cint,
        (MPI_Comm, Cint, Ptr{Cint}, Ptr{Cint}, Cint, Ptr{Cint}, Ptr{Cint}, MPI_Info, Cint, Ptr{MPI_Comm}),
        comm,length(sources),sources,source_weights,length(destinations),destinations,destination_weights,Info(infokws...),reorder,graph_comm) v"2.2"

    return graph_comm
end

"""
    graph_comm = Dist_graph_create(comm::Comm, sources::Vector{Cint}, degrees::Vector{Cint}, destinations::Vector{Cint}; weights::Union{Vector{Cint}, Unweighted, WeightsEmpty}=UNWEIGHTED, reorder=false, infokws...)

Create a new communicator from a given directed graph topology, described by incoming and outgoing edges on an existing communicator.

# Arguments
- `comm::Comm`: The communicator on which the distributed graph topology should be induced.
- `sources::Vector{Cint}`: An array with the ranks for which this call will specify outgoing edges.
- `degrees::Vector{Cint}`: An array with the number of outgoing edges for each entry in the sources array.
- `destinations::Vector{Cint}`: An array containing with lenght of the sum of the entries in the degrees array
                                describing the ranks towards the edges point.
- `weights::Union{Vector{Cint}, Unweighted, WeightsEmpty}=MPI.UNWEIGHTED`: The edge weights of the specified edges.
- `reorder::Bool=false`: If set true, then the MPI implementation can reorder the source and destination indices.

# Example
We can generate a ring graph `1 --> 2 --> ... --> N --> 1`, where N is the number of ranks in the communicator, as follows
```julia
julia> rank = MPI.Comm_rank(comm);
julia> N = MPI.Comm_size(comm);
julia> sources = Cint[rank];
julia> degrees = Cint[1];
julia> destinations = Cint[mod(rank-1, N)];
julia> graph_comm = Dist_graph_create(comm, sources, degrees, destinations)
```

# External links
$(_doc_external("MPI_Dist_graph_create"))
"""
function Dist_graph_create(comm::Comm, sources::Vector{Cint}, degrees::Vector{Cint}, destinations::Vector{Cint}; weights::Union{Vector{Cint}, Unweighted, WeightsEmpty}=UNWEIGHTED, reorder=false, infokws...)
    graph_comm = Comm()
    # int MPI_Dist_graph_create(MPI_Comm comm_old, int n, const int sources[],
    #       const int degrees[], const int destinations[],
    #       const int weights[],
    #       MPI_Info info, int reorder, MPI_Comm *comm_dist_graph)
    @mpichk ccall((:MPI_Dist_graph_create, libmpi), Cint,
        (MPI_Comm, Cint, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, MPI_Info, Cint, Ptr{MPI_Comm}),
        comm,length(sources),sources,degrees,destinations,weights,Info(infokws...),reorder,graph_comm) v"2.2"

    return graph_comm
end

"""
    indegree, outdegree, weighted = Dist_graph_neighbors_count(graph_comm::Comm)

Return the number of in and out edges for the calling processes in a distributed graph topology and a flag indicating whether the distributed graph is weighted.

# Arguments
- `graph_comm::Comm`: The communicator of the distributed graph topology.

# Example
Let us assume the following graph `0 <--> 1 --> 2`, which has no weights on its edges, then the
process with rank 1 will obtain the following result from calling the function

```julia-repl
julia> Dist_graph_neighbors_count(graph_comm)
(1,2,false)
```

# External links
$(_doc_external("MPI_Dist_graph_neighbors_count"))
"""
function Dist_graph_neighbors_count(graph_comm::Comm)
    indegree = Ref{Cint}()
    outdegree = Ref{Cint}()
    weighted = Ref{Cint}()
    # int MPI_Dist_graph_neighbors_count(MPI_Comm comm, int *indegree, int *outdegree, int *weighted)
    @mpichk ccall((:MPI_Dist_graph_neighbors_count, libmpi), Cint,
        (MPI_Comm, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}),
        graph_comm,indegree,outdegree,weighted) v"2.2"

    return (indegree[], outdegree[], weighted[] != 0)
end

"""
    Dist_graph_neighbors!(graph_comm::Comm, sources::Vector{Cint}, source_weights::Vector{Cint}, destinations::Vector{Cint}, destination_weights::Vector{Cint})

Return the neighbors and edge weights of the calling process in a distributed graph topology.

# Arguments
- `graph_comm::Comm`: The communicator of the distributed graph topology.
- `sources::Vector{Cint}`: A preallocated vector, which will be filled with the ranks
                        of the processes whose edges pointing towards the calling process. The
                        length is exactly the indegree returned by [`MPI.Dist_graph_neighbors_count`](@ref).
- `source_weights::Vector{Cint}`: A preallocated vector, which will be filled with the weights
                        associated to the edges pointing towards the calling process. The
                        length is exactly the indegree returned by [`MPI.Dist_graph_neighbors_count`](@ref).
- `destinations::Vector{Cint}`: A preallocated vector, which will be filled with the ranks
                        of the processes towards which the edges of the calling process point. The
                        length is exactly the outdegree returned by [`MPI.Dist_graph_neighbors_count`](@ref).
- `destination_weights::Vector{Cint}`: A preallocated vector, which will be filled with the weights
                        associated to the edges of the outgoing edges of the calling process point. The
                        length is exactly the outdegree returned by [`MPI.Dist_graph_neighbors_count`](@ref).

# Example
Let us assume the following graph `0 <-3-> 1 -4-> 2`, then the process with rank 1 will require to
preallocate a sources vector of length 1 and a destination vector of length 2. The call will fill
the vectors as follows:

```julia-repl
julia> Dist_graph_neighbors!(graph_comm, sources, source_weights, destinations, destination_weights);
julia> sources
[0]
julia> source_weights
[3]
julia> destinations
[0,2]
julia> destination_weights
[3,4]
```

Note that the edge between ranks 0 and 1 can have a different weight depending on wether it is the
incoming edge "`(0,1)"` or the outgoing one "`(1,0)"`.

# External links
$(_doc_external("MPI_Dist_graph_neighbors"))
"""
function Dist_graph_neighbors!(graph_comm::Comm, sources::Vector{Cint}, source_weights::Vector{Cint}, destinations::Vector{Cint}, destination_weights::Vector{Cint})
    # int MPI_Dist_graph_neighbors(MPI_Comm comm,
    #       int maxindegree, int sources[], int sourceweights[],
    #       int maxoutdegree, int destinations[], int destweights[])
    @mpichk ccall((:MPI_Dist_graph_neighbors, libmpi), Cint,
        (MPI_Comm, Cint, Ptr{Cint}, Ptr{Cint}, Cint, Ptr{Cint}, Ptr{Cint}),
        graph_comm,length(sources),sources,source_weights,length(destinations),destinations,destination_weights) v"2.2"
end

"""
    Dist_graph_neighbors!(graph_comm::Comm, sources::Vector{Cint}, destinations::Vector{Cint})

Return the neighbors of the calling process in a distributed graph topology without edge weights.

# Arguments
- `graph_comm::Comm`: The communicator of the distributed graph topology.
- `sources::Vector{Cint}`: A preallocated vector, which will be filled with the ranks of the
                        processes whose edges pointing towards the calling process. The
                        length is exactly the indegree returned by [`MPI.Dist_graph_neighbors_count`](@ref).
- `destinations::Vector{Cint}`: A preallocated vector, which will be filled with the ranks
                        of the processes towards which the edges of the calling process point. The
                        length is exactly the outdegree returned by [`MPI.Dist_graph_neighbors_count`](@ref).

# Example
Let us assume the following graph `0 <--> 1 --> 2`, then the process with rank 1 will require to
preallocate a sources vector of length 1 and a destination vector of length 2. The call will fill
the vectors as follows:

```julia-repl
julia> Dist_graph_neighbors!(graph_comm, sources, destinations);
julia> sources
[0]
julia> destinations
[0,2]
```

# External links
$(_doc_external("MPI_Dist_graph_neighbors"))
"""
function Dist_graph_neighbors!(graph_comm::Comm, sources::Vector{Cint}, destinations::Vector{Cint})
    source_weights = Array{Cint}(undef,0)
    destination_weights = Array{Cint}(undef,0)
    Dist_graph_neighbors!(graph_comm, sources::Vector{Cint}, source_weights, destinations::Vector{Cint}, destination_weights)
end