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

"""
    ValueShapes.default_datatype(T::Type)

Return a default specific type U that is more specific than T, with U <: T.

e.g.

    ValueShapes.default_datatype(Real) == Float64
    ValueShapes.default_datatype(Complex) == Complex{Float64}
"""
function default_datatype end

@inline default_datatype(::Type{>:Int}) = Int
@inline default_datatype(::Type{>:Float64}) = Float64
@inline default_datatype(::Type{>:Real}) = Float64
@inline default_datatype(::Type{>:Complex{Float64}}) = Complex{Float64}
@inline default_datatype(T::Type) = T

"""
    abstract type AbstractValueShape

An `AbstractValueShape` combines type and size information.

Subtypes are defined for shapes of scalars (see [`ScalarShape`](@ref)),
arrays (see [`ArrayShape`](@ref)), constant values
(see [`ConstValueShape`](@ref)) and `NamedTuple`s (see
[`NamedTupleShape`](@ref)).

Subtypes of `AbstractValueShape` must support `eltype`, `size` and
[`totalndof`](@ref).

Value shapes can be used as constructors to generate values of the given
shape with undefined content. If the element type of the shape is an abstract
or union type, a suitable concrete type will be chosen automatically, if
possible (see [`ValueShapes.default_datatype`](@ref)):

```julia
shape = ArrayShape{Real}(2,3)
A = shape(undef)
typeof(A) == Array{Float64,2}
size(A) == (2, 3)
valshape(A) == ArrayShape{Float64}(2,3)
```

Use

    (shape::AbstractValueShape)(data::AbstractVector{<:Real})::eltype(shape)

to view a flat vector of anonymous real values
as a value of the given shape:

```julia
data = [1, 2, 3, 4, 5, 6]
shape(data) == [1 3 5; 2 4 6]
```

In return,

    Base.Vector{<:Real}(undef, shape::AbstractValueShape)

will create a suitable uninitialized vector of the right length to hold such
flat data for the given shape. If no type `T` is given, a suitable data
type will be chosen automatically.

When dealing with multiple vectors of flattened data, use

    shape.(data::ArraysOfArrays.AbstractVectorOfSimilarVectors)

ValueShapes supports this via specialized broadcasting.

In return,

    ArraysOfArrays.VectorOfSimilarVectors{<:Real}(shape::AbstractValueShape)

will create a suitable vector (of length zero) of vectors that can hold
flattened data for the given shape. The result will be a
`VectorOfSimilarVectors` wrapped around a 2-dimensional `ElasticArray`.
This way, all data is stored in a single contiguous chunk of memory.
"""
abstract type AbstractValueShape end

export AbstractValueShape

# Reserve broadcasting semantics for value shapes:
@inline Base.Broadcast.broadcastable(shape::AbstractValueShape) =
    throw(ArgumentError("broadcasting over `AbstractValueShape`s is reserved"))

function _valshapeoftype end

"""
    ValueShapes.default_unshaped_eltype(shape::AbstractValueShape)

Returns the default real array element type to use for unshaped
representations of data with shape `shape`.

Subtypes of `AbstractValueShape` must implemenent
`ValueShapes.default_unshaped_eltype`.
"""
function default_unshaped_eltype end

"""
    ValueShapes.shaped_type(shape::AbstractValueShape, ::Type{T}) where {T<:Real}
    ValueShapes.shaped_type(shape::AbstractValueShape)

Returns the type the will result from reshaping a real-valued vector (of
element type `T`, if specified) with `shape`.

Subtypes of `AbstractValueShape` must implement

    ValueShapes.shaped_type(shape::AbstractValueShape, ::Type{T}) where {T<:Real}
"""
function shaped_type end

shaped_type(shape::AbstractValueShape) = shaped_type(shape, default_unshaped_eltype(shape))

"""
    valshape(x)::AbstractValueShape
    valshape(acc::ValueAccessor)::AbstractValueShape

Get the value shape of an arbitrary value, resp. the shape a `ValueAccessor`
is based on, or the shape of the variates for a `Distribution`.
"""
function valshape end
export valshape

@inline valshape(x::T) where T = _valshapeoftype(T)

"""
    elshape(x)::AbstractValueShape

Get the shape of the elements of x
"""
function elshape end
export elshape

@inline elshape(x::T) where T = _valshapeoftype(eltype(T))

@inline elshape(A::AbstractArray{<:AbstractArray}) = ArrayShape{eltype(eltype(A))}(innersize(A)...)

"""
    totalndof(shape::AbstractValueShape)

Get the total number of degrees of freedom of values of the given shape.

Equivalent to the length of a vector that would result from flattening the
data into a sequence of real numbers, excluding any constant values.
"""
function totalndof end
export totalndof

"""
    unshaped(x)::AbstractVector{<:Real}

Retrieve the unshaped underlying data of x, assuming x is a structured view
(based on some [`AbstractValueShape`](@ref)) of a flat/unstructured
real-valued data vector.

Example:

```julia
shape = NamedTupleShape(
    a = ScalarShape{Real}(),
    b = ArrayShape{Real}(2, 3)
)
data = [1, 2, 3, 4, 5, 6, 7]
x = shape(data)
@assert unshaped(x) == data
@assert unshaped(x.a) == view(data, 1:1)
@assert unshaped(x.b) == view(data, 2:7)
```
"""
function unshaped end
export unshaped

unshaped(x::SubArray{T,0}) where {T<:Real} = view(parent(x), x.indices[1]:x.indices[1])
unshaped(x::AbstractArray{T,1}) where {T<:Real} = x
unshaped(x::Base.ReshapedArray{T,N,<:AbstractArray{T,1}}) where {T<:Real,N} = parent(x)

"""
    stripscalar(x)

Dereference value `x`.

If x is a scalar-like object, like a 0-dimensional array or a `Ref`,
`stripscalar` returns it's inner value. Otherwise, `x` is returned unchanged.

Useful to strip shaped scalar-like views of their 0-dim array semantics
(if present), but leave array-like views unchanged.

Example:

```julia
data = [1, 2, 3]
shape1 = NamedTupleShape(a = ScalarShape{Real}(), b = ArrayShape{Real}(2))
x1 = shape1(data)
@assert x1 isa AbstractArray{<:NamedTuple,0}
@assert stripscalar(x) isa NamedTuple

shape2 = ArrayShape{Real}(3)
x2 = shape2(data)
@assert x2 isa AbstractArray{Int,1}
@assert ref(x2) isa AbstractArray{Int,1}
```
"""
function stripscalar end
export stripscalar

stripscalar(x::Any) = x
stripscalar(x::Ref) = x[]
stripscalar(x::AbstractArray{T,0}) where T = x[]

function _checkcompat(shape::AbstractValueShape, data::AbstractVector{<:Real})
    n_shape = totalndof(shape)
    n_data = length(eachindex(data))
    if n_shape != n_data
        throw(ArgumentError("Data vector of length $(n_data) incompatible with value shape with $(n_shape) degrees of freedom"))
    end
    nothing
end

function _checkcompat_inner(shape::AbstractValueShape, data::AbstractVector{<:AbstractVector{<:Real}})
    n_shape = totalndof(shape)
    n_data = prod(innersize(data))
    if n_shape != n_data
        throw(ArgumentError("Data vector of length $(n_data) incompatible with value shape with $(n_shape) degrees of freedom"))
    end
    nothing
end

@inline function _apply_shape_to_data(shape::AbstractValueShape, data::AbstractVector{<:Real})
    @boundscheck _checkcompat(shape, data)
    view(data, ValueAccessor(shape, 0))
end

@static if VERSION >= v"1.3"
    @inline function (shape::AbstractValueShape)(data::AbstractVector{<:Real})
        _apply_shape_to_data(shape, data)
    end
else
    # With Julia < v1.3, need to define (shape::SomeShape)(...) explicitly for
    # each shape type, see source code for the respective shapes.
end

Base.Vector{T}(::UndefInitializer, shape::AbstractValueShape) where {T <: Real} =
    Vector{T}(undef, totalndof(shape))

Base.Vector{<:Real}(::UndefInitializer, shape::AbstractValueShape) =
    Vector{default_unshaped_eltype(shape)}(undef, shape)

ArraysOfArrays.VectorOfSimilarVectors{T}(shape::AbstractValueShape) where {T<:Real} =
    VectorOfSimilarVectors(ElasticArray{T}(undef, totalndof(shape), 0))

const VSBroadcasted1{F,T} = Base.Broadcast.Broadcasted{
    <:Base.Broadcast.AbstractArrayStyle{1},
    <:Any,
    F,
    <:Tuple{T}
}

const VSBroadcasted2{F,T1,T2} = Base.Broadcast.Broadcasted{
    <:Base.Broadcast.AbstractArrayStyle{1},
    <:Any,
    F,
    <:Tuple{T1,T2}
}

# Specialize (::AbstractValueShape).(::AbstractVector{<:AbstractVector{<:Real}}):
Base.copy(instance::VSBroadcasted1{<:AbstractValueShape,AbstractVector{<:AbstractVector{<:Real}}}) =
    broadcast(view, instance.args[1], Ref(ValueAccessor(instance.f, 0)))

function _zerodim_array(x::T) where T
    A = Array{T,0}(undef)
    A[] = x
end


text:
using PowerSimulations
using Test

@testset "Device Constructors" begin
    include("thermalgen_testing.jl")
    include("renewables_testing.jl")
    include("load_testing.jl")
    #include("hydro_testing.jl")
    #include("storage_testing.jl")
end

@testset "Network Constructors" begin
    include("powermodels_testing.jl")
    include("network_testing.jl")
end

@testset "Services Constructors" begin
    include("service_testing.jl")
end

@testset "Model Constructors" begin
    include("model_testing.jl")
    #include("model_solve_testing.jl")
    #include("buildED_CN_testing.jl")
    #include("buildED_NB_testing.jl")
end

@testset "Simulation routines" begin
    include("simulations_testing.jl")
end

ext: jl


text:
module AWSMetadata

using DataStructures: OrderedDict
using HTTP
using JSON

"""
    parse_aws_metadata()

Update the AWS service metadata definitions in `metadata.json` by parsing the `AWS SDK JS`.

Check to see if the AWS service exists in the metadata, if not then a new service has been
created and we should generate low and high level wrappers for it. Otherwise, check to see
if the service SHA hash matches what we keep on file, if they differ then the API has been
updated and we need to re-generate low and high level wrappers for the service.
"""
function parse_aws_metadata()
    """

    """
    function _process_service(file, version)
        data_changed = true
        push!(metadata, file["name"] => Dict("version" => version, "sha" => file["sha"]))
        push!(services_modified, file)
    end

    """
    Only get the latest version of a Service, we don't care for supporting older ones
    """
    function _filter_latest_service_version(services)
        seen_services = String[]
        new_list = OrderedDict{String, Any}[]

        for service in services
            service_name, version = _get_service_info(service)

            if !(service_name in seen_services)
                push!(seen_services, service_name)
                push!(new_list, service)
            end
        end

        return new_list
    end

    """
    Get the service name and its version
    """
    function _get_service_info(service)
        filename = join(split(service["name"], '.')[1:end-2],'.')
        filename = split(filename, '-')
        service_name = join(filename[1:end-3], '-')
        version = join(filename[end-2:end], '-')

        return (service_name, version)
    end

    """
    Get a list of all full Service files from the AWS SDK JS GitHub repository
    """
    function _get_aws_sdk_js_files()
        headers = ["User-Agent" => "JuliaCloud/AWSCore.jl"]
        url = "https://api.github.com/repos/aws/aws-sdk-js/contents/apis"
        req = HTTP.get(url, headers)
        files = JSON.parse(String(req.body), dicttype=OrderedDict)
        filter!(f -> occursin(r".normal.json$", f["name"]), files)  # Only get ${Service}.normal.json files
        files = _filter_latest_service_version(files)

        return files
    end

    metadata_path = joinpath(@__DIR__, "metadata.json")
    metadata = JSON.parsefile(metadata_path, dicttype=OrderedDict)

    files = _get_aws_sdk_js_files()

    data_changed = false
    services_modified = OrderedDict{String, Any}[]

    # Loop through all Service files, if its a new service or the SHA hashes do not match add them to a Dictionary to have their definitions (re)generated
    for file in files
        service_name, version = _get_service_info(file)
        filename = file["name"]

        if !haskey(metadata, filename)
            println(service_name, " does not exist in metadata.")
            _process_service(file, version)
        else
            if metadata[filename]["sha"] != file["sha"]
                println(service_name, " sha hashes do not match, updating.")
                _process_service(file, version)
            end
        end
    end

    # Regenerate all low level definitions, and only modified/new high level definitions
    if data_changed
        _generate_low_level_wrapper(files)
        _generate_high_level_wrapper(services_modified)
        open(metadata_path, "w") do f
            print(f, json(OrderedDict(metadata), 2))
        end
    end
end

function _generate_low_level_wrapper(services::Array{OrderedDict{String, Any}})
    service_definitions = _generate_service_definitions(services)

    template = """
    module Services

    using AWSCorePrototype

    $(join(service_definitions, "\n"))

    end
    """

    services_path = joinpath(@__DIR__, "AWSCorePrototypeServices.jl")
    open(services_path, "w") do f
        print(f, template)
    end
end

function _generate_service_definitions(services::Array{OrderedDict{String, Any}})
    service_definitions = String[]

    for service in services
        service_name = service["name"]
        println("Generating low level wrapper for $service_name")
        request = HTTP.get(service["download_url"])

        service_metadata = JSON.parse(String(request.body))["metadata"]

        definition = _generate_service_definition(service_metadata)
        if definition != nothing
            push!(service_definitions, definition)
        end
    end

    return service_definitions
end

function _generate_service_definition(service::Dict{String, Any})
    println("Generating Metadata definitions for ", service["serviceId"])
    request_protocol = service["protocol"]
    service_name = service["endpointPrefix"]
    service_id = replace(lowercase(service["serviceId"]), ' ' => '_')
    api_version = service["apiVersion"]

    if request_protocol == "rest-xml"
        return "const $service_id = AWSCorePrototype.RestXMLService(\"$service_name\", \"$api_version\")"
    elseif request_protocol in ["ec2", "query"]
        return "const $service_id = AWSCorePrototype.QueryService(\"$service_name\", \"$api_version\")"
    elseif request_protocol == "rest-json"
        return "const $service_id = AWSCorePrototype.RestJSONService(\"$service_name\", \"$api_version\")"
    elseif request_protocol == "json"
        json_version = service["jsonVersion"]
        target = service["targetPrefix"]
        return "const $service_id = AWSCorePrototype.JSONService(\"$service_name\", \"$api_version\", \"$json_version\", \"$target\")"
    else
        # This should ping an alarm, we need to update this source code as its something we haven't seen before
        println(service_name, " uses a new protocol ", request_protocol)
    end
end

function _generate_high_level_wrapper(services::Array{OrderedDict{String, Any}})
    for service in services
        service_name = service["name"]
        println("Generating high level wrapper for $service_name")
        service = JSON.parse(String(HTTP.get(service["download_url"]).body))
        service_name = lowercase(service["metadata"]["serviceId"])
        service_name = replace(service_name, ' ' => '_')
        operations = service["operations"]
        shapes = service["shapes"]

        protocol = service["metadata"]["protocol"]

        operations = _generate_wrapper(service_name, protocol, operations, shapes)

        service_path = joinpath(@__DIR__, "services/$service_name.jl")
        open(service_path, "w") do f
            println(f, "include(\"../AWSCorePrototypeServices.jl\")")
            println(f, "using .Services: $service_name\n")
            print(f, join(operations, "\n"))
        end
    end
end

function _generate_wrapper(service_name::String, protocol::String, operations::Dict{String, Any}, shapes::Dict{String, Any})
    """
    Clean up the documentation by removing any HTML tags, and Julia escaping characters
    """
    function _documentation_cleaning(documentation)
        documentation = replace(documentation, r"\<.*?\>" => "")
        documentation = replace(documentation, '$' => ' ')
        documentation = replace(documentation, "\\" => ' ')

        return documentation
    end

    """
    Get the required and optional parameters for a given operation
    """
    function _get_parameters(input, shapes)
        required_parameters = Dict()
        optional_parameters = Dict()

        input_shape = shapes[input["shape"]]

        if haskey(input_shape, "required")
            for parameter in input_shape["required"]
                required_parameters[parameter] = _documentation_cleaning(get(input_shape["members"][parameter], "documentation", ""))
            end
        end

        if haskey(input_shape, "members")
            for parameter in input_shape["members"]
                parameter_name = parameter[1]

                if !haskey(required_parameters, parameter_name)
                    optional_parameters[parameter_name] = _documentation_cleaning(get(parameter[2], "documentation", ""))
                end
            end
        end

        return (required_parameters, optional_parameters)
    end

    """
    Generate the high level definition for a given operation
    """
    function _generate_operation_definition(name, protocol, method, request_uri, required_parameters, optional_parameters, documentation)
        required_param_keys = collect(keys(required_parameters))

        operation_definition = """
        \"\"\"
            $name()

        $documentation
        """

        # Add in the required parameters if applicable
        if !isempty(required_parameters)
            operation_definition = operation_definition * """

            Required Parameters
            $(json(required_parameters, 2))"""
        end

        # Add in the optional parameters if applicable
        if !isempty(optional_parameters)
            operation_definition = operation_definition * """

            Optional Parameters
            $(json(optional_parameters, 2))"""
        end

        operation_definition = operation_definition * "\"\"\""

        # Depending on the protocol type of the operation we need to generate a different definition
        if protocol in ["json", "query", "ec2"]
            if !isempty(required_parameters)
                operation_definition = operation_definition * "\n$name(args) = $service_name(\"$name\", args)\n"
            else
                operation_definition = operation_definition * """\n
                    $name() = $service_name(\"$name\")
                    $name(args) = $service_name(\"$name\", args)
                    """
            end
        elseif protocol == "rest-json"
            if !isempty(required_parameters)
                operation_definition = operation_definition * "\n$name(args) = $service_name(\"$method\", \"$request_uri\", args)\n"
            else
                operation_definition = operation_definition * """\n
                $name() = $service_name(\"$method\", \"$request_uri\")
                $name(args) = $service_name(\"$method\", \"$request_uri\", args)
                """
            end
        elseif protocol == "rest-xml"
            if !isempty(required_parameters)
                operation_definition = operation_definition * """
                $name($(join(required_param_keys, ", "))) = $service_name(\"$method\", \"$request_uri\")
                $name($(join(required_param_keys, ", ")), args) = $service_name(\"$method\", \"$request_uri\", args)
                $name(a...; b...) = $name(a..., b)
                """
            else
                operation_definition = operation_definition * """
                $name() = $service_name(\"$method\", \"$request_uri\")
                $name(args) $service_name(\"$method\", \"$request_uri\", args)
                $name(a...; b...) = $name(a..., b)
                """
            end
        end

        return operation_definition
    end

    operation_definitions = []

    for operation in operations
        operation = operation[2]
        name = operation["name"]
        method = operation["http"]["method"]
        request_uri = operation["http"]["requestUri"]

        documentation = ""

        if haskey(operation, "documentation")
            documentation = _documentation_cleaning(operation["documentation"])
        end

        required_parameters = Dict()
        optional_parameters = Dict()

        if haskey(operation, "input")
            required_parameters, optional_parameters = _get_parameters(operation["input"], shapes)
        end

        operation_definition = _generate_operation_definition(name, protocol, method, request_uri, required_parameters, optional_parameters, documentation)

        push!(operation_definitions, operation_definition)
    end

    return operation_definitions
end

end

max_stars_repo_name: mattBrzezinski/Julia-AWS-Status

max_issues_repo_name: mattBrzezinski/Julia-AWS-Status

max_forks_repo_name: mattBrzezinski/Julia-AWS-Status


text:
function localstiffness(p::CG2D{1},c,m::Mesh,cache)
    D = [(c.p3[1] - c.p2[1]) (c.p1[1] - c.p3[1]) (c.p2[1] - c.p1[1]);
          (c.p3[2] - c.p2[2]) (c.p1[2] - c.p3[2]) (c.p2[2] - c.p1[2])]
    D'D/(4*cache.elemJ[cid(c)])
end

function localstiffness(p::CG2D,c,m::Mesh,cache)
    pdphi = AC274.dphi(p)
    Al = zeros(Float64,𝖓(p),𝖓(p))
    invA = inv(AC274.Ak(c)')
    for i = 1:𝖓(p), j = 1:𝖓(p)
        Al[i,j] =
        cache.elemJ[cid(c)]*AC274.do_quad_ref(x->dot(invA*pdphi[i](x),invA*pdphi[j](x)),p)
    end
    Al
end

𝒩(p::CG2D{1}) = length(mesh(p).vertices)
𝒩(p::CG2D{2}) = length(mesh(p).vertices) + length(p.dualmesh.vertices)

# Assemble stiffness matrix
function stiffness{pp}(p::CG2D{pp},is∂D,cache=AC274.generateMatrices(p))
    m = mesh(p)

    A = spzeros(𝒩(p),𝒩(p))
    for c in mesh(p)
        Al = localstiffness(p,c,m,cache)
        for i = 1:𝖓(p), j = 1:𝖓(p)
            A[ℳ(p,c,pp,i),ℳ(p,c,pp,j)] += Al[i,j]
        end
    end
    #applyA∂D!(A,calciᴰ(p,is∂D,cache.neighbors))
    A
end

function applyA∂D!(A,iᴰ)
    for i in iᴰ
        A[i,1:end] = spzeros(1,size(A,2))
        A[i,i] = 1.0
    end
    A
end
applyA∂D!(A,is∂D,cache,m::Mesh) = applyA∂D!(A,calciᴰ(m,is∂D,cache.neighbors))

function loadV(p::CG2D, is∂D, f, gN, gD, cache=AC274.generateMatrices(p))
    iᴰ = calciᴰ(p,is∂D,cache.neighbors)
    F = loadVbase(p, f, cache)
    apply∂N!(p, F, is∂D, gN, cache)
    apply∂D!(p, F, gD, iᴰ)
    F
end

function loadVbase(p::CG2D, f, cache)
    basis, dphi, ns = cache.basis, cache.dphi, cache.nodes
    m = mesh(p)
    F = zeros(𝒩(p))

    for c in mesh(p)
        for i = 1:nbf(p)
            F[ℳ(p,c,porder(p),i)] += cache.elemJ[cid(c)]*do_quad_ref(x->(f(x)*basis[i](x)),p)
        end
    end

    F
end

function apply∂N!(p, F, is∂D, gN, cache)
    m = mesh(p)
    for c in m
        for face in c
            if !has_neighbor(m, c, face, cache) && !is∂D(face.p1)
                vRHS = zeros(Float64,nbf(p))
                gx, gw = cache.Δ1quadpoints
                # Manually inline quadrature
                factor = sqrt(dot(face.p2-face.p1,face.p2-face.p1))/2
                for i in 1:length(gw)
                    point⁻, w = gx[i], gw[i]
                    cpoint⁻ = ∂(c,face,point⁻)
                    AC274.evaluate_basis2d!(vRHS,p,-gN(face.p1)*factor*w,cpoint⁻)
                end
                for i = 1:nbf(p)
                    F[ℳ(p,c,porder(p),i)] += vRHS[i]
                end
            end
        end
    end
end

function apply∂D!(p, F, gD, iᴰ)
    idxs = Int64[i for i in iᴰ]
    𝒩 = length(p.mesh.vertices) + length(p.dualmesh.vertices)
    for i in idxs
        i > 𝒩 && (i=mod1(i,𝒩))
        v = i > length(p.mesh.vertices) ?
            p.dualmesh.vertices[i-length(p.mesh.vertices)] :
            p.mesh.vertices[i]
        F[i] = gD(v)
    end
end

# Calculate indicies for dirichlet nodes indicated by a given filter function is∂D
function calciᴰ{pp}(p::CG2D{pp},is∂D,neighbors = computeNeighbors(mesh(p)))
    iᴰ = Set{Int64}()
    for c in mesh(p)
        for edge in c
            if !has_neighbor(mesh(p), c, edge, neighbors) && is∂D(edge.p1) && is∂D(edge.p2)
                c = mesh(p)[cid(edge)]
                id = AC274.edgeid(c,edge)
                # Add both to the dirichlet boundary
                push!(iᴰ,ℳ(p,c,pp,id))
                push!(iᴰ,ℳ(p,c,pp,mod1(id+1,3)))
                pp == 2 && push!(iᴰ,ℳ(p,c,pp,3+id))
            end
        end
    end
    iᴰ
end


text:
# Julia wrapper for header: /usr/include/libavutil/buffer.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    av_buffer_alloc,
    av_buffer_allocz,
    av_buffer_create,
    av_buffer_default_free,
    av_buffer_ref,
    av_buffer_unref,
    av_buffer_is_writable,
    av_buffer_get_opaque,
    av_buffer_get_ref_count,
    av_buffer_make_writable,
    av_buffer_realloc,
    av_buffer_pool_init,
    av_buffer_pool_init2,
    av_buffer_pool_uninit,
    av_buffer_pool_get

function av_buffer_alloc(size::Integer)
    ccall((:av_buffer_alloc, libavutil), Ptr{AVBufferRef}, (Cint,), size)
end

function av_buffer_allocz(size::Integer)
    ccall((:av_buffer_allocz, libavutil), Ptr{AVBufferRef}, (Cint,), size)
end

function av_buffer_create(data, size::Integer, free, opaque, flags::Integer)
    ccall((:av_buffer_create, libavutil), Ptr{AVBufferRef}, (Ptr{UInt8}, Cint, Ptr{Cvoid}, Ptr{Cvoid}, Cint), data, size, free, opaque, flags)
end

function av_buffer_default_free(opaque, data)
    ccall((:av_buffer_default_free, libavutil), Cvoid, (Ptr{Cvoid}, Ptr{UInt8}), opaque, data)
end

function av_buffer_ref(buf)
    ccall((:av_buffer_ref, libavutil), Ptr{AVBufferRef}, (Ptr{AVBufferRef},), buf)
end

function av_buffer_unref(buf)
    ccall((:av_buffer_unref, libavutil), Cvoid, (Ptr{Ptr{AVBufferRef}},), buf)
end

function av_buffer_is_writable(buf)
    ccall((:av_buffer_is_writable, libavutil), Cint, (Ptr{AVBufferRef},), buf)
end

function av_buffer_get_opaque(buf)
    ccall((:av_buffer_get_opaque, libavutil), Ptr{Cvoid}, (Ptr{AVBufferRef},), buf)
end

function av_buffer_get_ref_count(buf)
    ccall((:av_buffer_get_ref_count, libavutil), Cint, (Ptr{AVBufferRef},), buf)
end

function av_buffer_make_writable(buf)
    ccall((:av_buffer_make_writable, libavutil), Cint, (Ptr{Ptr{AVBufferRef}},), buf)
end

function av_buffer_realloc(buf, size::Integer)
    ccall((:av_buffer_realloc, libavutil), Cint, (Ptr{Ptr{AVBufferRef}}, Cint), buf, size)
end

function av_buffer_pool_init(size::Integer, alloc)
    ccall((:av_buffer_pool_init, libavutil), Ptr{AVBufferPool}, (Cint, Ptr{Cvoid}), size, alloc)
end

function av_buffer_pool_init2(size::Integer, opaque, alloc, pool_free)
    ccall((:av_buffer_pool_init2, libavutil), Ptr{AVBufferPool}, (Cint, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}), size, opaque, alloc, pool_free)
end

function av_buffer_pool_uninit(pool)
    ccall((:av_buffer_pool_uninit, libavutil), Cvoid, (Ptr{Ptr{AVBufferPool}},), pool)
end

function av_buffer_pool_get(pool)
    ccall((:av_buffer_pool_get, libavutil), Ptr{AVBufferRef}, (Ptr{AVBufferPool},), pool)
end


text:
#################
#     TESTS
#################

using Nemo, FastArithmetic, Base.Test

function testMonomialDual()
  print("monomialToDual(), dualToMonomial()... ")

  k,u=FiniteField(5,1,"u")
  K,t=PolynomialRing(k,"t")
  P=t^2+t+1
  a0=[k(1),k(0)]
  a1=[k(0),k(1)]
  a2=a0+a1
  c=[k(4),k(2)]

  @test monomialToDual(a0,P)+monomialToDual(a1,P)==monomialToDual(a2,P)

  @test a0==dualToMonomial(monomialToDual(a0,P),P)

  @test a1==dualToMonomial(monomialToDual(a1,P),P)

  @test c==dualToMonomial(monomialToDual(c,P),P)

  @test dualToMonomial(a0,P)+dualToMonomial(a1,P)==dualToMonomial(a2,P)

  @test a0==monomialToDual(dualToMonomial(a0,P),P)

  @test a1==monomialToDual(dualToMonomial(a1,P),P)

  @test c==monomialToDual(dualToMonomial(c,P),P)
  println("PASS")
end

function testMulRem()
  print("mulT(), mulTmid(), remT()... ")

  k,u=FiniteField(5,1,"u")
  K,t=PolynomialRing(k,"t")
  P=t^2-t+2
  Q=2*t-1
  a=[k(2),k(1),k(1),k(0),k(3)]
  b=[k(-1),k(0),k(1),k(2),k(3),k(-1)]

  @test mulT(a,P,2)==t+4

  @test mulTmid(a,P,2)==t+4

  @test mulT(b,Q,4)==4*t^3+3*t^2+2*t+1

  @test mulTmid(b,Q,4)==4*t^3+3*t^2+2*t+1

  l,v=FiniteField(1993,1,"v")
  S,x=PolynomialRing(l,"x")
  fibo=[l(1),l(1)]
  W=x^2-x-1
  Y=t+1

  @test remT(fibo,W,7)==[l(1),l(1),l(2),l(3),l(5),l(8),l(13)]

  @test remT([k(1)],Y,3)==[k(1),k(4),k(1)]

  println("PASS")
end

function testRemTnaif()
  print("remTnaif()... ")
  l,v=FiniteField(1993,1,"v")
  S,x=PolynomialRing(l,"x")
  fibo=[l(1),l(1)]
  W=x^2-x-1

  @test remTnaif(fibo,W,8)==[l(1),l(1),l(2),l(3),l(5),l(8),l(13),l(21)]

  println("PASS")
end

function testBerlekampMassey()
  print("berlekampMassey()... ")

  l,v=FiniteField(1993,1,"v")
  S,x=PolynomialRing(l,"x")

  @test berlekampMassey([l(1),l(1),l(2),l(3),l(5),l(8),l(13),l(21)],2,S)==x^2-x-1

  @test berlekampMassey([l(1),l(1),l(2),l(3)],2,S)==x^2-x-1

  println("PASS")
end

function testEmbedProject()
  print("embed(), project(), computeR()... ")

  k,u=FiniteField(5,1,"u")
  K,t=PolynomialRing(k,"t")
  P=t^2+t+1
  Q=t^3+t+1
  R=computeR(P,Q)
  F=ResidueRing(ZZ,5)
  S,x=PolynomialRing(F,"x")

  @test parent(R)==parent(P)

  @test degree(R)==degree(P)*degree(Q)

  @test length(factor(R))==1

  up=monomialToDual([k(1)],P)
  uq=monomialToDual([k(1)],Q)
  a0=[k(1),k(0)]
  a1=[k(0),k(1)]
  a2=a0+a1
  b0=[k(1),k(0),k(0)]
  b2=[k(0),k(0),k(1)]
  b3=b0+b2

  @test embed(a0,P,uq,Q)+embed(a1,P,uq,Q)==embed(a2,P,uq,Q)

  @test embed(up,P,b0,Q)+embed(up,P,b2,Q)==embed(up,P,b3,Q)

  @test project(dualToMonomial(embed(monomialToDual(a0,P),P,uq,Q),R),P,Q)==a0

  @test project(dualToMonomial(embed(monomialToDual(a1,P),P,uq,Q),R),P,Q)==a1

  @test project(dualToMonomial(embed(up,P,monomialToDual(b0,Q),Q),R),Q,P)==b0

  @test project(dualToMonomial(embed(up,P,monomialToDual(b2,Q),Q),R),Q,P)==b2

  println("PASS")
end

function testPhi1()
  print("phi1(), inversePhi1()... ")

  k,u=FiniteField(5,1,"u")
  K,t=PolynomialRing(k,"t")
  P=t^2+t+1
  Q=t^3+t+1
  b=fq_nmod[k(1) k(2);k(4) k(2); k(1) k(3)]
  a=phi1(b,P,Q)
  R=computeR(P,Q)
  c=dualToMonomial(a,R)
  B=inversePhi1(c,P,Q)

  B=transpose(B)

  for j in 1:3
    B[:,j]=dualToMonomial(B[:,j],P)
  end

  B=transpose(B)

  for j in 1:2
    B[:,j]=monomialToDual(B[:,j],Q)
  end

  @test b==B

  up=monomialToDual([k(1)],P)
  uq=monomialToDual([k(1)],Q)
  r=embed(up,P,uq,Q)
  ur=fq_nmod[k(1) k(0);k(0) k(0); k(0) k(0)]

  for i in 1:2
    ur[:,i]=monomialToDual(ur[:,i],Q)
  end

  d=phi1(ur,P,Q)

  @test r==d

  J=transpose(inversePhi1([k(0),k(1),k(0),k(0),k(0),k(0)],P,Q))

  for i in 1:3
    J[:,i]=dualToMonomial(J[:,i],P)
  end

  @test J==[k(0) k(0) k(0);k(0) k(1) k(0)]

  println("PASS")
end

function testPhi2()
  print("phi2(), inversePhi2()... ")

  k,u=FiniteField(5,1,"u")
  K,t=PolynomialRing(k,"t")
  P=t^2+t+1
  Q=t^3+t+1
  R=computeR(P,Q)
  un=fq_nmod[k(1) k(0) k(0); k(0) k(0) k(0)]
  xy=fq_nmod[k(0) k(0) k(0); k(0) k(1) k(0)]
  b1=fq_nmod[k(2) k(1) k(4); k(3) k(1) k(4)]
  b2=fq_nmod[k(3) k(2) k(3); k(1) k(2) k(0)]

  @test phi2(un,P,Q,R)==[k(1),k(0),k(0),k(0),k(0),k(0)]

  @test phi2(xy,P,Q,R)==[k(0),k(1),k(0),k(0),k(0),k(0)]

  @test phi2(b1+b2,P,Q,R)==phi2(b1,P,Q,R)+phi2(b2,P,Q,R)

  b2p=transpose(copy(b2))

  for i in 1:2
    b2p[:,i]=monomialToDual(b2p[:,i],Q)
  end

  @test dualToMonomial(phi1(b2p,P,Q),R)==phi2(b2,P,Q,R)

  c=inversePhi2(monomialToDual([k(0),k(1),k(0),k(0),k(0),k(0)],R),P,Q,R)

  for i in 1:3
    c[:,i]=dualToMonomial(c[:,i],P)
  end

  c=transpose(c)

  for i in 1:2
    c[:,i]=dualToMonomial(c[:,i],Q)
  end

  c=transpose(c)

  @test c==xy

  a1=[k(1),k(3),k(2),k(2),k(0),k(1)]
  a2=[k(3),k(3),k(0),k(4),k(1),k(2)]

  @test inversePhi2(a1+a2,P,Q,R)==inversePhi2(a1,P,Q,R)+inversePhi2(a2,P,Q,R)

  q=inversePhi2(a1,P,Q,R)

  for i in 1:3
    q[:,i]=dualToMonomial(q[:,i],P)
  end

  q=transpose(q)

  for i in 1:2
    q[:,i]=dualToMonomial(q[:,i],Q)
  end

  q=transpose(q)

  @test monomialToDual(phi2(q,P,Q,R),computeR(P,Q))==a1

  aa=inversePhi2(monomialToDual(phi2(b1,P,Q,R),R),P,Q,R)

  for i in 1:3
    aa[:,i]=dualToMonomial(aa[:,i],P)
  end

  aa=transpose(aa)

  for i in 1:2
    aa[:,i]=dualToMonomial(aa[:,i],Q)
  end

  aa=transpose(aa)

  @test aa==b1

  g=inversePhi1(dualToMonomial(a2,R),P,Q)

  for i in 1:2
    g[:,i]=monomialToDual(g[:,i],Q)
  end

  @test transpose(g)==inversePhi2(a2,P,Q,R)

  println("PASS")
end

function testAll()

  testMonomialDual()
  testMulRem()
  testRemTnaif()
  testBerlekampMassey()
  testEmbedProject()
  testPhi1()
  testPhi2()

  println("\nAll tests passed successfully.\n")

end

testAll()

max_stars_repo_name: artichaut/fastarithmetic

max_issues_repo_name: artichaut/fastarithmetic

max_forks_repo_name: artichaut/fastarithmetic


text:
const to = TimerOutput()

function prettytime(t::Int)
    if t < 1e3
        value, units = t, "ns"
    elseif t < 1e6
        value, units = t / 1e3, "μs"
    elseif t < 1e9
        value, units = t / 1e6, "ms"
    else
        if t / 1e9 < 60
          value, units = t / 1e9, "s"
        elseif t / 1e9 < 3600
          value, units = t / (1e9*60), "m"
        else
          value, units = t / (1e9*60*60), "h"
        end
    end

    if round(value) >= 100
        str = string(@sprintf("%.0f", value), units)
    elseif round(value * 10) >= 100
        str = string(@sprintf("%.1f", value), units)
    else
        str = string(@sprintf("%.2f", value), units)
    end
    return lpad(str, 6, " ")
end
function prettytime(t::AbstractFloat)
    if t < 1e3
        value, units = t, "ns"
    elseif t < 1e6
        value, units = t / 1e3, "μs"
    elseif t < 1e9
        value, units = t / 1e6, "ms"
    else
        if t / 1e9 < 60
          value, units = t / 1e9, "s"
        elseif t / 1e9 < 3600
          value, units = t / (1e9*60), "m"
        else
          value, units = t / (1e9*60*60), "h"
        end
    end

    if round(value) >= 100
        str = string(@sprintf("%.0f", value), units)
    elseif round(value * 10) >= 100
        str = string(@sprintf("%.1f", value), units)
    else
        str = string(@sprintf("%.2f", value), units)
    end
    return lpad(str, 6, " ")
end


text:
"""
    to_dict(object) -> Dict

Serialize `object` into a dictionary of type `Dict{Symbol,Any}`.
Values stored into this dictionary must be of any native
Julia data type (Real, Integer, String, etc.)
Dictionary keys are `Symbols`.
"""
function to_dict(object)
  @abstractmethod
end

"""
    check_dict(::Type{T},dict::Dict) where T

Check validity of a dictionary `dict` for an object of type `T`.
It runs succesfully if the dictionary is valid for a particular
type or throws an error in any other case.
"""
function check_dict(::Type{T},dict::Dict) where T
  @abstractmethod
end

"""
    from_dict(::Type{T},dict::Dict) where T

De-serialize an object of type `T` from the dictionary `dict`.
Values stored into this dictionary must be of any native
Julia data type (Real, Integer, String, etc.)
Dictionary keys are `Symbols`.
"""
function from_dict(::Type{T},dict::Dict) where T
  @abstractmethod
end


text:
module Openscad

abstract type Op end
abstract type Group <: Op end
abstract type Terminal <: Op end

export Circle, Square, End,
    Unify, Difference,
    Translate,
    toString, editor

include("./3dObjects.jl")

include("./2dObjects.jl")

include("./Booleans.jl")

include("./Transforms.jl")

struct End <: Terminal end
toString(s::End) = ";"

function editor(tree::Op, show::Bool = true)
    tempFile = toString(tree)
    path = "/tmp/$(hash(tree)).scad"
    open(path, "w") do f
        write(f, tempFile)
    end
    if show
        sleep(0.5)
        run(`openscad $path`)
    end
    tempFile
end

end


text:
#****************************************************************************
#                     Molecular Dynamics Potentials (MDP)
#                            CESMIX-MIT Project
#
#  Contributing authors: Ngoc-Cuong Nguyen (cuongng@mit.edu, exapde@gmail.com)
#****************************************************************************

function cppfiles(foldername,filename, stropu, strcpu, strgpu, potnum)

#foldername = "app";
opufile = "opu" * filename;
cpufile = "cpu" * filename;
gpufile = "gpu" * filename;

if potnum==1
    stropu = replace(stropu, "int ng)" => "int ng, int potnum)");
    stropu = replace(stropu, "int);" => "int, int);");
    strcpu = replace(strcpu, "int ng)" => "int ng, int potnum)");
    strcpu = replace(strcpu, "int);" => "int, int);");
    strgpu = replace(strgpu, "int ng)" => "int ng, int potnum)");
    strgpu = replace(strgpu, "int);" => "int, int);");
end

#display(foldername * "/" * opufile * ".cpp")

ioopu = open(foldername * "/" * opufile * ".cpp", "w");
write(ioopu, stropu);
close(ioopu);

iocpu = open(foldername * "/" * cpufile * ".cpp", "w");
write(iocpu, strcpu);
close(iocpu);

iogpu = open(foldername * "/" * gpufile * ".cu", "w");
write(iogpu, strgpu);
close(iogpu);

end

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP


text:
#! /usr/bin/julia

# Rosetta Code, Priority queue

using Base.Collections

test = ["Clear drains" 3;
        "Feed cat" 4;
        "Make tea" 5;
        "Solve RC tasks" 1;
        "Tax return" 2]

task = PriorityQueue(Base.Order.Reverse)
for i in 1:size(test)[1]
    enqueue!(task, test[i,1], test[i,2])
end

println("Tasks, completed according to priority:")
while !isempty(task)
    (t, p) = peek(task)
    dequeue!(task)
    println("    \"", t, "\" has priority ", p)
end

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
using MCMCBenchmarks, Test

@testset "Plot Tests " begin
    cd( @__DIR__)
    results = CSV.read("Plot_Data.csv")

    #Plot parameter recovery
    recoveryPlots = plotrecovery(results, (mu=0, sigma=1), (:sampler, :Nd))
    @test typeof(recoveryPlots) ==  Array{Plots.Plot, 1}
    @test length(recoveryPlots) == 6

    #Plot mean run time as a function of number of data points (Nd) for each sampler
    meantimePlot = plotsummary(results, :Nd, :time, (:sampler, );yscale=:log10)
    @test typeof(meantimePlot) ==  Array{Plots.Plot, 1}
    @test length(meantimePlot) == 2

    #Plot mean allocations as a function of number of data points (Nd) for each sampler
    meanallocPlot = plotsummary(results, :Nd, :allocations, (:sampler, );yscale=:log10,
      ylabel="Allocations (log scale)")
    @test typeof(meanallocPlot) ==  Array{Plots.Plot, 1}
    @test length(meanallocPlot) == 1

    #Plot mean ess per second of number of data points (Nd) for each sampler
    essPS = plotsummary(results, :Nd, :ess_ps, (:sampler, ))
    @test typeof(essPS) ==  Array{Plots.Plot, 1}
    @test length(essPS) == 2

    #Plot density of effective sample size as function of number of data points (Nd) for each sampler
    essPlots = plotdensity(results, :ess, (:sampler, :Nd))
    @test typeof(essPlots) ==  Array{Plots.Plot, 1}
    @test length(essPlots) == 4

    #Plot density of rhat as function of number of data points (Nd) for each sampler
    rhatPlots = plotdensity(results, :r_hat, (:sampler, :Nd))
    @test typeof(rhatPlots) ==  Array{Plots.Plot, 1}
    @test length(rhatPlots) == 2

    #Plot density of time as function of number of data points (Nd) for each sampler
    timePlots = plotdensity(results, :time, (:sampler, :Nd))
    @test typeof(timePlots) ==  Array{Plots.Plot, 1}
    @test length(timePlots) == 2

    #Plot density of gc time percent as function of number of data points (Nd) for each sampler
    gcPlots = plotdensity(results, :gcpercent, (:sampler, :Nd))
    @test typeof(gcPlots) ==  Array{Plots.Plot, 1}
    @test length(gcPlots) == 1

    #Plot density of memory allocations as function of number of data points (Nd) for each sampler
    memPlots = plotdensity(results, :allocations, (:sampler, :Nd);xscale=:log10,
      xlabel="Allocations (log scale)")
    @test typeof(memPlots) ==  Array{Plots.Plot, 1}
    @test length(memPlots) == 1

    #Plot density of megabytes allocated as function of number of data points (Nd) for each sampler
    megPlots = plotdensity(results, :megabytes, (:sampler, :Nd))
    @test typeof(megPlots) ==  Array{Plots.Plot, 1}
    @test length(megPlots) == 1

    #Scatter plot of epsilon and effective sample size as function of number of data points (Nd) for each sampler
    scatterPlots = plotscatter(results, :epsilon, :ess, (:sampler, :Nd))
    @test typeof(scatterPlots) ==  Array{Plots.Plot, 1}
    @test length(scatterPlots) == 4
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

struct AppsEventAuthorizationsApi <: SwaggerApi
    client::Swagger.Client
end

function _swaggerinternal_appsEventAuthorizationsList(_api::AppsEventAuthorizationsApi, event_context::String; cursor=nothing, limit=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", Dict{String, Any}, "/apps.event.authorizations.list", ["slackAuth"])
    Swagger.set_param(_ctx.query, "event_context", event_context)  # type String
    Swagger.set_param(_ctx.query, "cursor", cursor)  # type String
    Swagger.set_param(_ctx.query, "limit", limit)  # type Int32
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded", "application/json"] : [_mediaType])
    return _ctx
end

"""

Get a list of authorizations for the given event context. Each authorization represents an app installation that the event is visible to.
Param: event_context::String (required)
Param: cursor::String
Param: limit::Int32
Return: Dict{String, Any}
"""
function appsEventAuthorizationsList(_api::AppsEventAuthorizationsApi, event_context::String; cursor=nothing, limit=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_appsEventAuthorizationsList(_api, event_context; cursor=cursor, limit=limit, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function appsEventAuthorizationsList(_api::AppsEventAuthorizationsApi, response_stream::Channel, event_context::String; cursor=nothing, limit=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_appsEventAuthorizationsList(_api, event_context; cursor=cursor, limit=limit, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

export appsEventAuthorizationsList


text:
"""
    collisionBGK(mom, tempdist, input)

Computes the BGK collision operation for a `D1Q3` and a `D2Q9` lattice Boltzmann algorithm.

The collision operator is the central object in kinetic theory.
However there is no general solution for this operator.
One work around is to make assumptions and try to model them.
The BGK operator is one of the fairly successful operators.
It stems from the idea that the distributions are close to their equilibria and that all of them *relax* with a single relaxation time `τ`.

# Math
The collision operator is usually the right hand side of the kinetic equation

`` \\partial_t f + \\frac{\\mathbf{p}}{m}\\cdot\\nabla f + \\mathbf{F}\\cdot \\partial_{\\mathbf{p}} f = \\Omega . ``

After tedious mathematics and expansions one can obtain the simplest collision operator the BGK as

`` \\Omega_{BGK} = -\\frac{f - f^{eq}}{\\tau}, ``

with `τ` being the relaxation parameter.
Such the collision process in the lattice Boltzmann method tries to equilibrate the system.

If there is an additional force acting on the fluid, i.e. `gravity` than this part is usually included in the numerical implementation of the collision operator.
The forcing vector of source term is often called `` \\mathbf{S}_{\\alpha} `` and given by

`` \\mathbf{S}_{\\alpha} = \\Delta t \\frac{3 w_{\\alpha}}{\\mathvbf{c}_{\\alpha}^2}\\mathbf{c}_{\\alpha,i} \\mathbf{F}_i. ``

# Example

## One spatial dimension
```jldoctest
julia> using JuSwalbe

julia> input, mom, force, dist = minimalsetup1d(10)
(JuSwalbe.Inputconstants
  lx: Int64 10
  ly: Int64 512
  maxruntime: Int64 100000
  dumping: Int64 1000
  τ: Float64 1.0
  gravity: Float64 0.0
  γ: Float64 0.01
  δ: Float64 1.0
  μ: Float64 0.16666666666666666
  kbt: Float64 0.0
, JuSwalbe.Macroquant{Array{Float64,1},Array{Float64,1}}
  height: Array{Float64}((10,)) [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]
  velocity: Array{Float64}((10,)) [0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1]
  pressure: Array{Float64}((10,)) [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
  energy: Array{Float64}((10,)) [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
, JuSwalbe.Forces{Array{Float64,1}}
  slip: Array{Float64}((10,)) [0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1]
  h∇p: Array{Float64}((10,)) [-0.1, -0.1, -0.1, -0.1, -0.1, -0.1, -0.1, -0.1, -0.1, -0.1]
  bathymetry: Array{Float64}((10,)) [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
  thermal: Array{Float64}((10,)) [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
, JuSwalbe.DistributionD1Q3{Array{Float64,1}}
  f0: Array{Float64}((10,)) [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]
  f1: Array{Float64}((10,)) [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
  f2: Array{Float64}((10,)) [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
)

julia> force.slip = zeros(10)
10-element Array{Float64,1}:
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0

julia> force.h∇p = zeros(10)
10-element Array{Float64,1}:
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0

julia> newdist = collisionBGK(mom, force, dist, input)
JuSwalbe.DistributionD1Q3{Array{Float64,1}}
  f0: Array{Float64}((10,)) [0.99, 0.99, 0.99, 0.99, 0.99, 0.99, 0.99, 0.99, 0.99, 0.99]
  f1: Array{Float64}((10,)) [0.05500000000000001, 0.05500000000000001, 0.05500000000000001, 0.05500000000000001, 0.05500000000000001, 0.05500000000000001, 0.05500000000000001, 0.05500000000000001, 0.05500000000000001, 0.05500000000000001]
  f2: Array{Float64}((10,)) [-0.045, -0.045, -0.045, -0.045, -0.045, -0.045, -0.045, -0.045, -0.045, -0.045]

julia> equi = calc_equilibrium_distribution(mom, gravity=input.gravity)
JuSwalbe.DistributionD1Q3{Array{Float64,1}}
  f0: Array{Float64}((10,)) [0.99, 0.99, 0.99, 0.99, 0.99, 0.99, 0.99, 0.99, 0.99, 0.99]
  f1: Array{Float64}((10,)) [0.05500000000000001, 0.05500000000000001, 0.05500000000000001, 0.05500000000000001, 0.05500000000000001, 0.05500000000000001, 0.05500000000000001, 0.05500000000000001, 0.05500000000000001, 0.05500000000000001]
  f2: Array{Float64}((10,)) [-0.045, -0.045, -0.045, -0.045, -0.045, -0.045, -0.045, -0.045, -0.045, -0.045]

julia> using Test; @test equi.f0 == newdist.f0 # without forces and τ = 1 they have to be equal
Test Passed
```

## Two spatial dimensions
```jldoctest
julia> using JuSwalbe, Test

julia> input, mom, force, dist = minimalsetup2d(10,5)
(JuSwalbe.Inputconstants
  lx: Int64 10
  ly: Int64 5
  maxruntime: Int64 100000
  dumping: Int64 1000
  τ: Float64 1.0
  gravity: Float64 0.0
  γ: Float64 0.01
  δ: Float64 1.0
  μ: Float64 0.16666666666666666
  kbt: Float64 0.0
, JuSwalbe.Macroquant{Array{Float64,2},JuSwalbe.Twovector{Array{Float64,2}}}
  height: Array{Float64}((10, 5)) [1.0 1.0 … 1.0 1.0; 1.0 1.0 … 1.0 1.0; … ; 1.0 1.0 … 1.0 1.0; 1.0 1.0 … 1.0 1.0]
  velocity: JuSwalbe.Twovector{Array{Float64,2}}
  pressure: Array{Float64}((10, 5)) [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]
  energy: Array{Float64}((10, 5)) [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]
, JuSwalbe.Forces{JuSwalbe.Twovector{Array{Float64,2}}}
  slip: JuSwalbe.Twovector{Array{Float64,2}}
  h∇p: JuSwalbe.Twovector{Array{Float64,2}}
  bathymetry: JuSwalbe.Twovector{Array{Float64,2}}
  thermal: JuSwalbe.Twovector{Array{Float64,2}}
, JuSwalbe.DistributionD2Q9{Array{Float64,2}}
  f0: Array{Float64}((10, 5)) [1.0 1.0 … 1.0 1.0; 1.0 1.0 … 1.0 1.0; … ; 1.0 1.0 … 1.0 1.0; 1.0 1.0 … 1.0 1.0]
  f1: Array{Float64}((10, 5)) [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]
  f2: Array{Float64}((10, 5)) [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]
  f3: Array{Float64}((10, 5)) [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]
  f4: Array{Float64}((10, 5)) [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]
  f5: Array{Float64}((10, 5)) [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]
  f6: Array{Float64}((10, 5)) [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]
  f7: Array{Float64}((10, 5)) [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]
  f8: Array{Float64}((10, 5)) [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]
)

julia> dist
JuSwalbe.DistributionD2Q9{Array{Float64,2}}
  f0: Array{Float64}((10, 5)) [1.0 1.0 … 1.0 1.0; 1.0 1.0 … 1.0 1.0; … ; 1.0 1.0 … 1.0 1.0; 1.0 1.0 … 1.0 1.0]
  f1: Array{Float64}((10, 5)) [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]
  f2: Array{Float64}((10, 5)) [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]
  f3: Array{Float64}((10, 5)) [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]
  f4: Array{Float64}((10, 5)) [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]
  f5: Array{Float64}((10, 5)) [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]
  f6: Array{Float64}((10, 5)) [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]
  f7: Array{Float64}((10, 5)) [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]
  f8: Array{Float64}((10, 5)) [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]

julia> equi = calc_equilibrium_distribution(mom, gravity=input.gravity)
JuSwalbe.DistributionD2Q9{Array{Float64,2}}
  f0: Array{Float64}((10, 5)) [0.9994666666666666 0.9994666666666666 … 0.9994666666666666 0.9994666666666666; 0.9994666666666666 0.9994666666666666 … 0.9994666666666666 0.9994666666666666; … ; 0.9994666666666666 0.9994666666666666 … 0.9994666666666666 0.9994666666666666; 0.9994666666666666 0.9994666666666666 … 0.9994666666666666 0.9994666666666666]
  f1: Array{Float64}((10, 5)) [0.006733333333333333 0.006733333333333333 … 0.006733333333333333 0.006733333333333333; 0.006733333333333333 0.006733333333333333 … 0.006733333333333333 0.006733333333333333; … ; 0.006733333333333333 0.006733333333333333 … 0.006733333333333333 0.006733333333333333; 0.006733333333333333 0.006733333333333333 … 0.006733333333333333 0.006733333333333333]
  f2: Array{Float64}((10, 5)) [-0.006599999999999999 -0.006599999999999999 … -0.006599999999999999 -0.006599999999999999; -0.006599999999999999 -0.006599999999999999 … -0.006599999999999999 -0.006599999999999999; … ; -0.006599999999999999 -0.006599999999999999 … -0.006599999999999999 -0.006599999999999999; -0.006599999999999999 -0.006599999999999999 … -0.006599999999999999 -0.006599999999999999]
  f3: Array{Float64}((10, 5)) [-0.006599999999999999 -0.006599999999999999 … -0.006599999999999999 -0.006599999999999999; -0.006599999999999999 -0.006599999999999999 … -0.006599999999999999 -0.006599999999999999; … ; -0.006599999999999999 -0.006599999999999999 … -0.006599999999999999 -0.006599999999999999; -0.006599999999999999 -0.006599999999999999 … -0.006599999999999999 -0.006599999999999999]
  f4: Array{Float64}((10, 5)) [0.006733333333333333 0.006733333333333333 … 0.006733333333333333 0.006733333333333333; 0.006733333333333333 0.006733333333333333 … 0.006733333333333333 0.006733333333333333; … ; 0.006733333333333333 0.006733333333333333 … 0.006733333333333333 0.006733333333333333; 0.006733333333333333 0.006733333333333333 … 0.006733333333333333 0.006733333333333333]
  f5: Array{Float64}((10, 5)) [-3.3333333333333335e-5 -3.3333333333333335e-5 … -3.3333333333333335e-5 -3.3333333333333335e-5; -3.3333333333333335e-5 -3.3333333333333335e-5 … -3.3333333333333335e-5 -3.3333333333333335e-5; … ; -3.3333333333333335e-5 -3.3333333333333335e-5 … -3.3333333333333335e-5 -3.3333333333333335e-5; -3.3333333333333335e-5 -3.3333333333333335e-5 … -3.3333333333333335e-5 -3.3333333333333335e-5]
  f6: Array{Float64}((10, 5)) [-0.0031666666666666666 -0.0031666666666666666 … -0.0031666666666666666 -0.0031666666666666666; -0.0031666666666666666 -0.0031666666666666666 … -0.0031666666666666666 -0.0031666666666666666; … ; -0.0031666666666666666 -0.0031666666666666666 … -0.0031666666666666666 -0.0031666666666666666; -0.0031666666666666666 -0.0031666666666666666 … -0.0031666666666666666 -0.0031666666666666666]
  f7: Array{Float64}((10, 5)) [-3.3333333333333335e-5 -3.3333333333333335e-5 … -3.3333333333333335e-5 -3.3333333333333335e-5; -3.3333333333333335e-5 -3.3333333333333335e-5 … -3.3333333333333335e-5 -3.3333333333333335e-5; … ; -3.3333333333333335e-5 -3.3333333333333335e-5 … -3.3333333333333335e-5 -3.3333333333333335e-5; -3.3333333333333335e-5 -3.3333333333333335e-5 … -3.3333333333333335e-5 -3.3333333333333335e-5]
  f8: Array{Float64}((10, 5)) [0.0034999999999999996 0.0034999999999999996 … 0.0034999999999999996 0.0034999999999999996; 0.0034999999999999996 0.0034999999999999996 … 0.0034999999999999996 0.0034999999999999996; … ; 0.0034999999999999996 0.0034999999999999996 … 0.0034999999999999996 0.0034999999999999996; 0.0034999999999999996 0.0034999999999999996 … 0.0034999999999999996 0.0034999999999999996]

julia> force.slip = JuSwalbe.Twovector(x=zeros(10,5), y=zeros(10,5))
JuSwalbe.Twovector{Array{Float64,2}}
  x: Array{Float64}((10, 5)) [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]
  y: Array{Float64}((10, 5)) [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]

julia> force.h∇p = JuSwalbe.Twovector(x=zeros(10,5), y=zeros(10,5))
JuSwalbe.Twovector{Array{Float64,2}}
  x: Array{Float64}((10, 5)) [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]
  y: Array{Float64}((10, 5)) [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]

julia> newdist = collisionBGK(mom, force, dist, input)
JuSwalbe.DistributionD2Q9{Array{Float64,2}}
  f0: Array{Float64}((10, 5)) [0.9994666666666666 0.9994666666666666 … 0.9994666666666666 0.9994666666666666; 0.9994666666666666 0.9994666666666666 … 0.9994666666666666 0.9994666666666666; … ; 0.9994666666666666 0.9994666666666666 … 0.9994666666666666 0.9994666666666666; 0.9994666666666666 0.9994666666666666 … 0.9994666666666666 0.9994666666666666]
  f1: Array{Float64}((10, 5)) [0.006733333333333333 0.006733333333333333 … 0.006733333333333333 0.006733333333333333; 0.006733333333333333 0.006733333333333333 … 0.006733333333333333 0.006733333333333333; … ; 0.006733333333333333 0.006733333333333333 … 0.006733333333333333 0.006733333333333333; 0.006733333333333333 0.006733333333333333 … 0.006733333333333333 0.006733333333333333]
  f2: Array{Float64}((10, 5)) [-0.006599999999999999 -0.006599999999999999 … -0.006599999999999999 -0.006599999999999999; -0.006599999999999999 -0.006599999999999999 … -0.006599999999999999 -0.006599999999999999; … ; -0.006599999999999999 -0.006599999999999999 … -0.006599999999999999 -0.006599999999999999; -0.006599999999999999 -0.006599999999999999 … -0.006599999999999999 -0.006599999999999999]
  f3: Array{Float64}((10, 5)) [-0.006599999999999999 -0.006599999999999999 … -0.006599999999999999 -0.006599999999999999; -0.006599999999999999 -0.006599999999999999 … -0.006599999999999999 -0.006599999999999999; … ; -0.006599999999999999 -0.006599999999999999 … -0.006599999999999999 -0.006599999999999999; -0.006599999999999999 -0.006599999999999999 … -0.006599999999999999 -0.006599999999999999]
  f4: Array{Float64}((10, 5)) [0.006733333333333333 0.006733333333333333 … 0.006733333333333333 0.006733333333333333; 0.006733333333333333 0.006733333333333333 … 0.006733333333333333 0.006733333333333333; … ; 0.006733333333333333 0.006733333333333333 … 0.006733333333333333 0.006733333333333333; 0.006733333333333333 0.006733333333333333 … 0.006733333333333333 0.006733333333333333]
  f5: Array{Float64}((10, 5)) [-3.3333333333333335e-5 -3.3333333333333335e-5 … -3.3333333333333335e-5 -3.3333333333333335e-5; -3.3333333333333335e-5 -3.3333333333333335e-5 … -3.3333333333333335e-5 -3.3333333333333335e-5; … ; -3.3333333333333335e-5 -3.3333333333333335e-5 … -3.3333333333333335e-5 -3.3333333333333335e-5; -3.3333333333333335e-5 -3.3333333333333335e-5 … -3.3333333333333335e-5 -3.3333333333333335e-5]
  f6: Array{Float64}((10, 5)) [-0.0031666666666666666 -0.0031666666666666666 … -0.0031666666666666666 -0.0031666666666666666; -0.0031666666666666666 -0.0031666666666666666 … -0.0031666666666666666 -0.0031666666666666666; … ; -0.0031666666666666666 -0.0031666666666666666 … -0.0031666666666666666 -0.0031666666666666666; -0.0031666666666666666 -0.0031666666666666666 … -0.0031666666666666666 -0.0031666666666666666]
  f7: Array{Float64}((10, 5)) [-3.3333333333333335e-5 -3.3333333333333335e-5 … -3.3333333333333335e-5 -3.3333333333333335e-5; -3.3333333333333335e-5 -3.3333333333333335e-5 … -3.3333333333333335e-5 -3.3333333333333335e-5; … ; -3.3333333333333335e-5 -3.3333333333333335e-5 … -3.3333333333333335e-5 -3.3333333333333335e-5; -3.3333333333333335e-5 -3.3333333333333335e-5 … -3.3333333333333335e-5 -3.3333333333333335e-5]
  f8: Array{Float64}((10, 5)) [0.0034999999999999996 0.0034999999999999996 … 0.0034999999999999996 0.0034999999999999996; 0.0034999999999999996 0.0034999999999999996 … 0.0034999999999999996 0.0034999999999999996; … ; 0.0034999999999999996 0.0034999999999999996 … 0.0034999999999999996 0.0034999999999999996; 0.0034999999999999996 0.0034999999999999996 … 0.0034999999999999996 0.0034999999999999996]

julia> @test equi.f0 == newdist.f0 # No force and τ=1 means fnew has to be feq!
Test Passed

julia> @test equi.f2 == newdist.f2
Test Passed
```

# References
In general a very good reference concerning the method at all is the book written by Timm Krüger
-[The Lattice Boltzmann Method: Principles and Practice](https://www.springer.com/gp/book/9783319446479)

## One spatial dimension
Here is the work by Chopard et al. really worth reading
- [Asymmetric lattice Boltzmann model for shallow water flows](https://www.sciencedirect.com/science/article/abs/pii/S0045793013003599)

## Two spatial dimensions
In two dimensions there is the problematic of how to add forces.
One way to justify is to do a Chapman Enskog expansion and match the shallow water equations with a source term.
Apparently there is not a unique solution but a class of possible forcings.
The first approach by Salmon does not include the lattice weights, which results in isotropy problems.
On the other hand Zhous approach is not as acqurate as he says, which was pointed out by Chopard.
Nevertheless here are some references educate yourself which is the fitting one
-[The lattice Boltzmann method as a basis for ocean circulation modeling, by Salmon](https://www.ingentaconnect.com/content/jmr/jmr/1999/00000057/00000003/art00005)
-[Lattice Boltzmann Methods for Shallow Water Flows, by Zhou](https://www.springer.com/gp/book/9783540407461)
-[An evaluation of force terms in the lattice Boltzmann models in simulating shallow water flows over complex topography, by Peng](https://onlinelibrary.wiley.com/doi/pdf/10.1002/fld.4726)

"""
function collisionBGK(mom::JuSwalbe.Macroquant{Vector{T},Vector{T}}, forces::JuSwalbe.Forces{Vector{T}}, tempdist::JuSwalbe.DistributionD1Q3{Vector{T}}, input::JuSwalbe.Inputconstants) where {T<:Number}
    # Get the size, type and allocate result array
    len = length(mom.height)
    weights = [T(2/3) T(1/6) T(1/6)]
    clat = [T(0) T(1) T(-1)]
    τ = input.τ
    gravity = input.gravity
    csquared = T(3)
    allforces = zeros(T, len)
    allforces = sum(hcat(forces.slip, forces.h∇p, forces.thermal, forces.bathymetry), dims=2)[:,1]

    newdist = JuSwalbe.DistributionD1Q3(f0=zeros(T, len),
                                        f1=zeros(T, len),
                                        f2=zeros(T, len))

    eqdist = calc_equilibrium_distribution(mom; gravity=T(gravity))

    newdist.f0 .= T(1-1/τ) * tempdist.f0 .+ T(1/τ) * eqdist.f0
    newdist.f1 .= T(1-1/τ) * tempdist.f1 .+ T(1/τ) * eqdist.f1 .+ weights[2] * csquared * clat[2] * allforces
    newdist.f2 .= T(1-1/τ) * tempdist.f2 .+ T(1/τ) * eqdist.f2 .+ weights[3] * csquared * clat[3] * allforces

    return newdist
end

function collisionBGK(mom::JuSwalbe.Macroquant{Matrix{T},JuSwalbe.Twovector{Matrix{T}}}, forces::JuSwalbe.Forces{JuSwalbe.Twovector{Matrix{T}}}, tempdist::JuSwalbe.DistributionD2Q9{Matrix{T}}, input::JuSwalbe.Inputconstants) where {T<:Number}
  # Get the size, type and allocate result array
  width, thick = size(mom.height)
  weights = [T(4/9) T(1/9) T(1/9) T(1/9) T(1/9) T(1/36) T(1/36) T(1/36) T(1/36)]
  clat = [T(0) T(0); T(1) T(0); T(0) T(1); T(-1) T(0); T(0) T(-1); T(1) T(1); T(-1) T(1); T(-1) T(-1); T(1) T(-1)]
  clatsquare = [T(0) T(1) T(1) T(1) T(1) T(2) T(2) T(2) T(2)]
  τ = input.τ
  gravity = input.gravity
  csquared = T(3)
  allforces = dist2array(forces)
  # println(allforces)

  newdist = JuSwalbe.DistributionD2Q9(f0=zeros(T, (width, thick)),
                                      f1=zeros(T, (width, thick)),
                                      f2=zeros(T, (width, thick)),
                                      f3=zeros(T, (width, thick)),
                                      f4=zeros(T, (width, thick)),
                                      f5=zeros(T, (width, thick)),
                                      f6=zeros(T, (width, thick)),
                                      f7=zeros(T, (width, thick)),
                                      f8=zeros(T, (width, thick)))

  eqdist = calc_equilibrium_distribution(mom; gravity=T(gravity))
  # Collision operation for a D2Q9 shallow water lattice Boltzmann.
  newdist.f0 .= T(1-1/τ) * tempdist.f0 .+ T(1/τ) * eqdist.f0

  newdist.f1 .= (T(1-1/τ) * tempdist.f1 .+ T(1/τ) * eqdist.f1
              .+ (weights[2] * csquared)/clatsquare[2] * (clat[2,1] * sum(allforces.x, dims=3)[:, :, 1] .+ clat[2,2] * sum(allforces.y, dims=3)[:, :, 1]))
  newdist.f2 .= (T(1-1/τ) * tempdist.f2 .+ T(1/τ) * eqdist.f2
              .+ (weights[3] * csquared)/clatsquare[3] * (clat[3,1] * sum(allforces.x, dims=3)[:, :, 1] .+ clat[3,2] * sum(allforces.y, dims=3)[:, :, 1]))
  newdist.f3 .= (T(1-1/τ) * tempdist.f3 .+ T(1/τ) * eqdist.f3
              .+ (weights[4] * csquared)/clatsquare[4] * (clat[4,1] * sum(allforces.x, dims=3)[:, :, 1] .+ clat[4,2] * sum(allforces.y, dims=3)[:, :, 1]))
  newdist.f4 .= (T(1-1/τ) * tempdist.f4 .+ T(1/τ) * eqdist.f4
              .+ (weights[5] * csquared)/clatsquare[5] * (clat[5,1] * sum(allforces.x, dims=3)[:, :, 1] .+ clat[5,2] * sum(allforces.y, dims=3)[:, :, 1]))
  newdist.f5 .= (T(1-1/τ) * tempdist.f5 .+ T(1/τ) * eqdist.f5
              .+ (weights[6] * csquared)/clatsquare[6] * (clat[6,1] * sum(allforces.x, dims=3)[:, :, 1] .+ clat[6,2] * sum(allforces.y, dims=3)[:, :, 1]))
  newdist.f6 .= (T(1-1/τ) * tempdist.f6 .+ T(1/τ) * eqdist.f6
              .+ (weights[7] * csquared)/clatsquare[7] * (clat[7,1] * sum(allforces.x, dims=3)[:, :, 1] .+ clat[7,2] * sum(allforces.y, dims=3)[:, :, 1]))
  newdist.f7 .= (T(1-1/τ) * tempdist.f7 .+ T(1/τ) * eqdist.f7
              .+ (weights[8] * csquared)/clatsquare[8] * (clat[8,1] * sum(allforces.x, dims=3)[:, :, 1] .+ clat[8,2] * sum(allforces.y, dims=3)[:, :, 1]))
  newdist.f8 .= (T(1-1/τ) * tempdist.f8 .+ T(1/τ) * eqdist.f8
              .+ (weights[9] * csquared)/clatsquare[9] * (clat[9,1] * sum(allforces.x, dims=3)[:, :, 1] .+ clat[9,2] * sum(allforces.y, dims=3)[:, :, 1]))

  return newdist
end

"""
    streamdistperiodic!(dist)

Streams the lattice Boltzmann distribution accoriding to their velocity vector.

The so called streaming process transports the distribution functions according to thier lattice velocity.
In one dimension this means that the zero velocity distribution stays at its lattice point while the positive velocity advances by one lattice node.
Of course the negative velocity component `.f2` decreases its lattice index by one.

In two dimensions there are nine of these operations according to the `D2Q9` lattice vectors.

# Math
After the collision process the newly calculated distribution functions need to be streamed.
They performe one interation or streaming per time update.
Such the distributions after collision `` f^{\\ast}_i `` are transported to their respective new position

`` f^{\\ast}_i(\\mathbf{x} + \\mathbf{c}_i \\Delta t)``

One has to be careful in the proximity of boundaries, because they could destroy mass which is unphysical.

# Example
```jldoctest
julia> using JuSwalbe

julia> dist = JuSwalbe.DistributionD2Q9(f0=ones(5,5), f1=fill(0.1, (5,5)), f2=fill(0.1, (5,5)),
                                        f3=fill(0.1, (5,5)), f4=fill(0.1, (5,5)), f5=zeros(5,5),
                                        f6=zeros(5,5), f7=zeros(5,5), f8=zeros(5,5))
JuSwalbe.DistributionD2Q9{Array{Float64,2}}
  f0: Array{Float64}((5, 5)) [1.0 1.0 … 1.0 1.0; 1.0 1.0 … 1.0 1.0; … ; 1.0 1.0 … 1.0 1.0; 1.0 1.0 … 1.0 1.0]
  f1: Array{Float64}((5, 5)) [0.1 0.1 … 0.1 0.1; 0.1 0.1 … 0.1 0.1; … ; 0.1 0.1 … 0.1 0.1; 0.1 0.1 … 0.1 0.1]
  f2: Array{Float64}((5, 5)) [0.1 0.1 … 0.1 0.1; 0.1 0.1 … 0.1 0.1; … ; 0.1 0.1 … 0.1 0.1; 0.1 0.1 … 0.1 0.1]
  f3: Array{Float64}((5, 5)) [0.1 0.1 … 0.1 0.1; 0.1 0.1 … 0.1 0.1; … ; 0.1 0.1 … 0.1 0.1; 0.1 0.1 … 0.1 0.1]
  f4: Array{Float64}((5, 5)) [0.1 0.1 … 0.1 0.1; 0.1 0.1 … 0.1 0.1; … ; 0.1 0.1 … 0.1 0.1; 0.1 0.1 … 0.1 0.1]
  f5: Array{Float64}((5, 5)) [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]
  f6: Array{Float64}((5, 5)) [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]
  f7: Array{Float64}((5, 5)) [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]
  f8: Array{Float64}((5, 5)) [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]

julia> dist.f1[2,2] = 1; dist.f2[2,2] = 1; dist.f5[2,2] = 1; dist.f3[2,2] = 1; dist.f8[2,2] = 1
1

julia> newdist = streamdistperiodic!(dist)
JuSwalbe.DistributionD2Q9{Array{Float64,2}}
  f0: Array{Float64}((5, 5)) [1.0 1.0 … 1.0 1.0; 1.0 1.0 … 1.0 1.0; … ; 1.0 1.0 … 1.0 1.0; 1.0 1.0 … 1.0 1.0]
  f1: Array{Float64}((5, 5)) [0.1 0.1 … 0.1 0.1; 0.1 0.1 … 0.1 0.1; … ; 0.1 0.1 … 0.1 0.1; 0.1 0.1 … 0.1 0.1]
  f2: Array{Float64}((5, 5)) [0.1 0.1 … 0.1 0.1; 0.1 0.1 … 0.1 0.1; … ; 0.1 0.1 … 0.1 0.1; 0.1 0.1 … 0.1 0.1]
  f3: Array{Float64}((5, 5)) [0.1 1.0 … 0.1 0.1; 0.1 0.1 … 0.1 0.1; … ; 0.1 0.1 … 0.1 0.1; 0.1 0.1 … 0.1 0.1]
  f4: Array{Float64}((5, 5)) [0.1 0.1 … 0.1 0.1; 0.1 0.1 … 0.1 0.1; … ; 0.1 0.1 … 0.1 0.1; 0.1 0.1 … 0.1 0.1]
  f5: Array{Float64}((5, 5)) [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]
  f6: Array{Float64}((5, 5)) [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]
  f7: Array{Float64}((5, 5)) [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]
  f8: Array{Float64}((5, 5)) [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]

julia> dist.f1
5×5 Array{Float64,2}:
 0.1  0.1  0.1  0.1  0.1
 0.1  1.0  0.1  0.1  0.1
 0.1  0.1  0.1  0.1  0.1
 0.1  0.1  0.1  0.1  0.1
 0.1  0.1  0.1  0.1  0.1

julia> newdist.f1
5×5 Array{Float64,2}:
 0.1  0.1  0.1  0.1  0.1
 0.1  0.1  0.1  0.1  0.1
 0.1  1.0  0.1  0.1  0.1
 0.1  0.1  0.1  0.1  0.1
 0.1  0.1  0.1  0.1  0.1

julia> dist.f8
5×5 Array{Float64,2}:
 0.0  0.0  0.0  0.0  0.0
 0.0  1.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0

julia> newdist.f8
5×5 Array{Float64,2}:
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 1.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0

```

# References
Any book on the lattice Boltzmann method should do.
- [The Lattice Boltzmann Method: Principles and Practice](https://www.springer.com/gp/book/9783319446479)

See also [`collisionBGK`](@ref)
"""
function streamdistperiodic!(dist::JuSwalbe.DistributionD1Q3{Vector{T}}) where {T<:Number}
    len = length(dist.f0)
    f1dummy = zeros(T, len)
    f2dummy = zeros(T, len)
    f1dummy .= circshift(dist.f1, 1)
    f2dummy .= circshift(dist.f2, -1)

    dist.f1 .= f1dummy
    dist.f2 .= f2dummy

    return dist
end

function streamdistperiodic!(dist::JuSwalbe.DistributionD2Q9{Matrix{T}}) where {T<:Number}
  width, thick = size(dist.f0)
  newdist = zeros(T, (width, thick, 9))
  # Shifting the distributions by the nine lattice velocities
  distarray = dist2array(dist)
  latvel = [(0, 0); (1, 0); (0, 1); (-1, 0); (0, -1); (1, 1); (-1, 1); (-1, -1); (1, -1)]
  for i in 2:9
    newdist[:, :, i] .= circshift(distarray[i, :, :], latvel[i])
  end

  dist = JuSwalbe.DistributionD2Q9(f0=dist.f0, f1=newdist[:, :, 2], f2=newdist[:, :, 3],
                                   f3=newdist[:, :, 4], f4=newdist[:, :, 5], f5=newdist[:, :, 6],
                                   f6=newdist[:, :, 7], f7=newdist[:, :, 8], f8=newdist[:, :, 9])

  return dist
end

"""
    streamdistperiodic(dist)

Streams the lattice Boltzmann distribution accoriding to their velocity vector.

The so called streaming process transports the distribution functions according to thier lattice velocity.
In one dimension this means that the zero velocity distribution stays at its lattice point while the positive velocity advances by one lattice node.
Of course the negative velocity component `.f2` decreases its lattice index by one.

In two dimensions there are nine of these operations according to the `D2Q9` lattice vectors.

# Math
After the collision process the newly calculated distribution functions need to be streamed.
They performe one interation or streaming per time update.
Such the distributions after collision `` f^{\\ast}_i `` are transported to their respective new position

`` f^{\\ast}_i(\\mathbf{x} + \\mathbf{c}_i \\Delta t)``

One has to be careful in the proximity of boundaries, because they could destroy mass which is unphysical.

# Example
```jldoctest
julia> using JuSwalbe

julia> dist = JuSwalbe.DistributionD2Q9(f0=ones(5,5), f1=fill(0.1, (5,5)), f2=fill(0.1, (5,5)),
                                        f3=fill(0.1, (5,5)), f4=fill(0.1, (5,5)), f5=zeros(5,5),
                                        f6=zeros(5,5), f7=zeros(5,5), f8=zeros(5,5))
JuSwalbe.DistributionD2Q9{Array{Float64,2}}
  f0: Array{Float64}((5, 5)) [1.0 1.0 … 1.0 1.0; 1.0 1.0 … 1.0 1.0; … ; 1.0 1.0 … 1.0 1.0; 1.0 1.0 … 1.0 1.0]
  f1: Array{Float64}((5, 5)) [0.1 0.1 … 0.1 0.1; 0.1 0.1 … 0.1 0.1; … ; 0.1 0.1 … 0.1 0.1; 0.1 0.1 … 0.1 0.1]
  f2: Array{Float64}((5, 5)) [0.1 0.1 … 0.1 0.1; 0.1 0.1 … 0.1 0.1; … ; 0.1 0.1 … 0.1 0.1; 0.1 0.1 … 0.1 0.1]
  f3: Array{Float64}((5, 5)) [0.1 0.1 … 0.1 0.1; 0.1 0.1 … 0.1 0.1; … ; 0.1 0.1 … 0.1 0.1; 0.1 0.1 … 0.1 0.1]
  f4: Array{Float64}((5, 5)) [0.1 0.1 … 0.1 0.1; 0.1 0.1 … 0.1 0.1; … ; 0.1 0.1 … 0.1 0.1; 0.1 0.1 … 0.1 0.1]
  f5: Array{Float64}((5, 5)) [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]
  f6: Array{Float64}((5, 5)) [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]
  f7: Array{Float64}((5, 5)) [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]
  f8: Array{Float64}((5, 5)) [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]

julia> dist.f1[2,2] = 1; dist.f2[2,2] = 1; dist.f5[2,2] = 1; dist.f3[2,2] = 1; dist.f8[2,2] = 1
1

julia> newdist = streamdistperiodic(dist)
JuSwalbe.DistributionD2Q9{Array{Float64,2}}
  f0: Array{Float64}((5, 5)) [1.0 1.0 … 1.0 1.0; 1.0 1.0 … 1.0 1.0; … ; 1.0 1.0 … 1.0 1.0; 1.0 1.0 … 1.0 1.0]
  f1: Array{Float64}((5, 5)) [0.1 0.1 … 0.1 0.1; 0.1 0.1 … 0.1 0.1; … ; 0.1 0.1 … 0.1 0.1; 0.1 0.1 … 0.1 0.1]
  f2: Array{Float64}((5, 5)) [0.1 0.1 … 0.1 0.1; 0.1 0.1 … 0.1 0.1; … ; 0.1 0.1 … 0.1 0.1; 0.1 0.1 … 0.1 0.1]
  f3: Array{Float64}((5, 5)) [0.1 1.0 … 0.1 0.1; 0.1 0.1 … 0.1 0.1; … ; 0.1 0.1 … 0.1 0.1; 0.1 0.1 … 0.1 0.1]
  f4: Array{Float64}((5, 5)) [0.1 0.1 … 0.1 0.1; 0.1 0.1 … 0.1 0.1; … ; 0.1 0.1 … 0.1 0.1; 0.1 0.1 … 0.1 0.1]
  f5: Array{Float64}((5, 5)) [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]
  f6: Array{Float64}((5, 5)) [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]
  f7: Array{Float64}((5, 5)) [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]
  f8: Array{Float64}((5, 5)) [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]

julia> dist.f1
5×5 Array{Float64,2}:
 0.1  0.1  0.1  0.1  0.1
 0.1  1.0  0.1  0.1  0.1
 0.1  0.1  0.1  0.1  0.1
 0.1  0.1  0.1  0.1  0.1
 0.1  0.1  0.1  0.1  0.1

julia> newdist.f1
5×5 Array{Float64,2}:
 0.1  0.1  0.1  0.1  0.1
 0.1  0.1  0.1  0.1  0.1
 0.1  1.0  0.1  0.1  0.1
 0.1  0.1  0.1  0.1  0.1
 0.1  0.1  0.1  0.1  0.1

julia> dist.f8
5×5 Array{Float64,2}:
 0.0  0.0  0.0  0.0  0.0
 0.0  1.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0

julia> newdist.f8
5×5 Array{Float64,2}:
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 1.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0

```

# References
Any book on the lattice Boltzmann method should do.
- [The Lattice Boltzmann Method: Principles and Practice](https://www.springer.com/gp/book/9783319446479)

See also [`collisionBGK`](@ref)
"""
function streamdistperiodic(dist::JuSwalbe.DistributionD1Q3{Vector{T}}) where {T<:Number}
  len = length(dist.f0)
  new_dist = JuSwalbe.DistributionD1Q3(f0 = dist.f0,
                                       f1 = zeros(len),
                                       f2 = zeros(len))
  f1dummy = zeros(T, len)
  f2dummy = zeros(T, len)
  f1dummy .= circshift(dist.f1, 1)
  f2dummy .= circshift(dist.f2, -1)

  new_dist.f1 .= f1dummy
  new_dist.f2 .= f2dummy

  return new_dist
end

function streamdistperiodic(dist::JuSwalbe.DistributionD2Q9{Matrix{T}}) where {T<:Number}
  width, thick = size(dist.f0)
  arr = zeros(width, thick, 9)

  distarray = dist2array(dist)
  latvel = [(0, 0); (1, 0); (0, 1); (-1, 0); (0, -1); (1, 1); (-1, 1); (-1, -1); (1, -1)]
  for i in 2:9
    arr[:, :, i] .= circshift(distarray[i, :, :], latvel[i])
  end
  new_dist = JuSwalbe.DistributionD2Q9(f0 = dist.f0,
                                       f1 = arr[:, :, 2],
                                       f2 = arr[:, :, 3],
                                       f3 = arr[:, :, 4],
                                       f4 = arr[:, :, 5],
                                       f5 = arr[:, :, 6],
                                       f6 = arr[:, :, 7],
                                       f7 = arr[:, :, 8],
                                       f8 = arr[:, :, 9])

  return new_dist
end

max_stars_repo_name: Zitzeronion/JuSwalbe

max_issues_repo_name: Zitzeronion/JuSwalbe


text:
using DecomposedMDPSolver
using Flux
using Test

Na = 5
Np = 10
solutions = [(x) -> rand(Na) for i=1:Np ]

## Weights network
base = Chain(Dense(4, 32, relu), Dense(32, Na))
attn = Chain(Dense(4, 32, relu), Dense(32, Np+1), softmax)
a2t_model = A2TNetwork(base, attn, solutions)
@test size(a2t_model(rand(4))) == (Na, 1)

a2t_model(rand(4))
input = rand(4)
i = 1
s = 10
a2t_model.solutions[s](input[:, i])

## Constant Weights
base = Chain(Dense(4, 32, relu), Dense(32, Na))
attn = Chain(ConstantLayer(Np+1), softmax)
a2t_model = A2TNetwork(base, attn, solutions)
@test size(a2t_model(rand(4))) == (Na, 1)

## Version that was failing -- local approx policy eval

sols = rand(100)
function val_fn(s)
    sols[1] = rand()
    return sols[1]
end
base = Chain(Dense(2, 32, relu), Dense(32, 1, σ))
attn = Chain(Dense(2, 32, relu), Dense(32, 2, exp))
solutions = [val_fn]

model = A2TNetwork(base, attn, solutions)

for l in model
    @test !isnothing(l)
end

S, G = rand(2, 100), rand(1,100)
data = Flux.Data.DataLoader((S, G), batchsize=32, shuffle = true)
opt = ADAM()
Flux.train!((x, y) -> Flux.mse(model(x), y), Flux.params(model), data, opt)

## Testing the state transformation network
base = Chain(Dense(2, 32, relu), Dense(32, 1, σ))
attn = Chain(Dense(2, 32, relu), Dense(32, 3, exp))
strans = [Chain(Dense(2, 32, relu), Dense(32, 2, relu)), Chain(Dense(2, 32, relu), Dense(32, 2, relu))]
solutions = [Chain(Dense(2, 32, relu), Dense(32, 1, σ)), Chain(Dense(2, 32, relu), Dense(32, 1, σ))]

model = A2TSTNetwork(base, attn, strans, solutions)
@test size(strans[1](rand(2))) == (2,)
@test size(strans[1](rand(2, 100))) == (2,100)
@test size(model(rand(2, 100))) == (1,100)

i = rand(2, 100)
v = [s(i) for s in model.solutions]
out = model(i)

S, G = rand(2, 100), rand(1,100)
data = Flux.Data.DataLoader((S, G), batchsize=32, shuffle = true)
opt = ADAM()
Flux.train!((x, y) -> Flux.mse(model(x), y), Flux.params(model), data, opt)

v2 = [s(i) for s in model.solutions]
out2 = model(i)

@test v2 == v
@test all(out2 .!= out)
for l in model
    @test !isnothing(l)
end

## Test the fine-tune A2T network

base = Chain(Dense(4, 32, relu), Dense(32, 9, σ))
attn = Chain(Dense(4, 32, relu), Dense(32, 3), softmax)
solutions = [Chain(Dense(4, 32, relu), Dense(32, 9, σ)), Chain(Dense(4, 32, relu), Dense(32, 9, σ))]
finetune = [Chain(Dense(13, 9)), Chain(Dense(13, 9))]

model = A2TFTNetwork(base, attn, solutions, finetune)

model(rand(4))
model(rand(4, 100))

S, G = rand(4, 100), rand(9,100)

v = [s(S) for s in model.solutions]
out = model(S)
val = Flux.mse(model(S), G)

data = Flux.Data.DataLoader((S, G), batchsize=32, shuffle = true)
opt = ADAM()
Flux.train!((x, y) -> Flux.mse(model(x), y), Flux.params(model), data, opt)

v2 = [s(S) for s in model.solutions]
out2 = model(S)
val2 = Flux.mse(model(S), G)

@test val2 < val
@test all(out[:] != out2[:])
@test all( v .== v2)

for l in model
    @test !isnothing(l)
end

@test deepcopy(model) isa A2TFTNetwork

## Test the fine-tune network
base = Chain(Dense(2, 2, sigmoid), Dense(2, 2, sigmoid), Dense(2,2))
mynet = FTNetwork(base, [1,2])

params_to_change = deepcopy(params(mynet.net))[1:4]
params_to_stay_same = deepcopy(params(mynet.net))[5:6]

@test length(Flux.trainable(mynet)) == 2

S = rand(2,100)
G = rand(2, 100)

data = Flux.Data.DataLoader((S, G), batchsize=32, shuffle = true)
opt = ADAM()
Flux.train!((x, y) -> Flux.mse(mynet(x), y), Flux.params(mynet), data, opt)

@test all(params_to_change .!= params(mynet.net)[1:4])
@test all(params_to_stay_same .== params(mynet.net)[5:6])

## Full test with A2T network + finetuning
base = Chain(Dense(4, 32, sigmoid), Dense(32, 9, σ))
attn = Chain(Dense(4, 32, sigmoid), Dense(32, 3), softmax)
solutions = [Chain(Dense(4, 32, sigmoid), Dense(32, 9, σ)), Chain(Dense(4, 32, sigmoid), Dense(32, 9, σ))]
mysols = [FTNetwork(net, [2]) for net in solutions]
model = A2TNetwork(base, attn, mysols, false)

base_params = deepcopy(params(model.base))
@test all(base_params .== params(model.base))
attn_params = deepcopy(params(model.attn))
@test all(attn_params .== params(model.attn))
params_to_change1 = deepcopy(params(model.solutions[1].net))[3:4]
@test all(params_to_change1 .== params(model.solutions[1].net)[3:4])
params_to_change2 = deepcopy(params(model.solutions[2].net))[3:4]
@test all(params_to_change2 .== params(model.solutions[2].net)[3:4])
params_to_stay_same1 = deepcopy(params(model.solutions[1].net))[1:2]
@test all(params_to_stay_same1 .== params(model.solutions[1].net)[1:2])
params_to_stay_same2 = deepcopy(params(model.solutions[2].net))[1:2]
@test all(params_to_stay_same2 .== params(model.solutions[2].net)[1:2])

@test length(Flux.trainable(model)) == 4
S = rand(4, 100)
G = rand(9, 100)

@test length(Flux.params(model)) == 12

data = Flux.Data.DataLoader((S, G), batchsize=32, shuffle = true)
opt = ADAM()
Flux.train!((x, y) -> Flux.mse(model(x), y), Flux.params(model), data, opt)

@test all(base_params .!= params(model.base))
@test all(attn_params .!= params(model.attn))
@test all(params_to_change1 .!= params(model.solutions[1].net)[3:4])
@test all(params_to_change2 .!= params(model.solutions[2].net)[3:4])
@test all(params_to_stay_same1 .== params(model.solutions[1].net)[1:2])
@test all(params_to_stay_same2 .== params(model.solutions[2].net)[1:2])


text:
using LinearAlgebra
using Random, Distributions
using Plots
Plots.scalefontsizes(1.5)

# Global Terms
Random.seed!(1234)

# Logistic Regression

function f(β)
    #=
    Objective Function
    x: n x d
    y: d x 1
    β: d x 1
    λ: penalty
    =#
    neg_yxβ = - X * β .* y
    log_term = log.(1 .+ exp.(neg_yxβ))
    # log_term = min.(log_term, 10000)  # Prevent Overflow
    summation = sum(log_term) / n
    return summation + λ/2 * norm(β,2)^2
end

function ∇f(β, σ)
    #=
    This is already noisy gradient.
    β: 100 x 1
    σ: noise
    =#
    neg_yxβ =  - X * β .* y
    weight = exp.(neg_yxβ) ./ (1 .+ exp.(neg_yxβ))

    # Replace nan in weight
    for idx in findall(isnan, weight)
        weight[idx] = 0
    end

    sums = sum(weight .* (-X .* y),dims=1) / n
    # @show sums
    # @show size(sums)
    gradient = sums' + λ .* β
    return gradient + rand(Normal(0,σ^2), d)
end
# MASG

function MASG(f, ∇f, nk, α, μ, σ)
    #=
        f:  objective functin
        ∇f: gradient of f
        X:  design matrix
        y:  {-1,1} class
        a_n:an array containing each step each iteration
        K:  total number of iterations
        α:  an array containing step sizes
        μ:  μ-strongly convex
        λ:  penalty
    =#
    x_0 = x_1 = zeros(d)

    # Record Histories
    x_history = zeros(sum(nk))
    f_history = zeros(sum(nk))
    ymk_history = zeros(sum(nk))
    grad_history = zeros(sum(nk))
    γ = 1

    for k in 1:K
        x_0 = x_1
        for m in 1:nk[k]
            w = (1 - sqrt(μ * α[k])) / (1 + sqrt(μ * α[k]))
            @show w
            x_updated = (1 + w) .* x_1 - w .* x_0
            @show x_updated
            fx_updated = f(x_updated)
            ymk_history[γ] = fx_updated
            x_0 = x_1
            x_1 = x_updated .- α[k] .* ∇f(x_updated, σ)
            grad_history[γ] = norm(∇f(x_updated, σ),2)
            # Record History
            x_history[γ] = norm(x_1)
            f_history[γ] = f(x_1)
            γ = γ + 1
        end
    end
    return x_history, f_history, ymk_history, grad_history
end

function compute_nk(κ, p, Δ, σ, label)
    C = floor(sqrt(κ) * log(2^(p+2))) + 1
    grid = [1:K;]
    nk = 2 .^ grid
    if label == 3.7
        nk[1] = floor(sqrt(κ) * log(κ)) + 1
    elseif label == 4
        nk[1] = floor(sqrt(κ) * log(2 * L * Δ / (σ^2 * sqrt(κ)))) + 1
    end
    return nk
end

function compute_MASG_α(L)
    grid = [1:K;]
    α = 1 ./ (2 .^ (2 .* [1:K;]) .* L)
    α[1] = 1 / L
    return α
end

# Vanilla Gradient Descent
function GD(f, ∇f, σ)
    x = zeros(d)
    f_history = zeros(η)
    x_history = zeros(η)
    for i in 1:η
        x = x - 1 / L * ∇f(x, σ)
        f_history[i] = f(x)
        x_history[i] = norm(x,2)
    end
    return x_history, f_history
end

# Accelerated Gradient Descent
function AGD(f, ∇f, σ)
    x_0 = x_1 = zeros(d)
    f_history = zeros(η)
    x_history = zeros(η)
    for i in 1:η
        y = x_1 + (i-1)/(i+2) *(x_1-x_0)
        x_0 = x_1
        x_1 = y - 1/L * ∇f(y, σ)

        f_history[i] = f(x_1)
        x_history[i] = norm(x_1,2)
    end
    return x_history, f_history
end

# AC-SA
## Bregman Distance
function D(β1, β2)
    #=
    Prox function, Bregman Distance
    =#
    return f(β1) - f(β2) - ∇f(β2) * (β1 - β2)
end

# μAGD

function compute_μAGD_a()
    grid = [1:η;]
    return 1 ./ (κ .* sqrt.(grid))
end

function compute_μAGD_A(a)
    A = zeros(η)
    A[1] = a[1]
    for i in 2:η
        A[i] = A[i-1] + a[i]
    end
    return A
end

function μAGD(f, ∇f, A, a, σ, mode)
    #=
    Mode is for selecting restart + slow down algorithm
    if 0, use vanilla
    if 2, use restart and slowdown-2
    =#
    f_history = zeros(η)
    x_history = zeros(η)

    xk = x0 = zeros(d)
    summation = a[1] .* ∇f(xk,σ)
    yk = x0 .- summation / λ
    f_history[1] = f(xk)
    a_square_sum = a[1]^2
    z0 = λ .* xk
    changed_flag = false

    for k in 2:η
        xk = A[k-1] / A[k] .* yk + a[k] / A[k] .* (x0 .- summation ./ λ)
        summation = summation .+ a[k] .* ∇f(xk, σ)
        yk = A[k-1] / A[k] .* yk + a[k] ./ A[k] .* (x0 .- summation ./ λ)
        zk = - summation .+ z0
        a_square_sum += a[k]^2

        if mode == 2 && norm(zk,2) <= (a_square_sum * d) && changed_flag == false
            a = compute_μAGD_a()
            A = compute_μAGD_A(a)
            changed_flag == true
        end

        # Record history
        f_history[k] = f(xk)
        x_history[k] = norm(xk, 2)
    end
    return x_history, f_history
end

## Data setup
global n = 1000    # number of datapoints
global d = 100    # number of dimensions
global λ = 0.01   # Penalty for beta term
global X = rand(Normal(0,1), (n,d))
β = rand(Normal(0, 1), (d, 1))
global y = sign.(X * β)

# Experiments
global K = 12      # number of iterations for MASG
global η = 8000   # Number of iterations for other methods
global κ = 1000   # condition number from the paper
global L = maximum(eigvals(X' * X))
global μ = L / κ

σ = sqrt(1e-1)

# MASG
α_MASG = compute_MASG_α(L)
## MASG Thm 3.7
nk_7 = compute_nk(κ, 1, norm(β,2), σ, 3.7)
x_MASG_7, f_MASG_7, ymk_MASG_7, grad_MASG_7 = MASG(f, ∇f, nk_7, α_MASG, μ, σ)

cumulative = zeros(12)
cumulative[1] = 219
for i in 2:12
    cumulative[i] = cumulative[i-1] + nk_7[i]
end
cumulative

plot(x_MASG_7,xscale=:log,yscale=:log,label="x",legend=:bottomleft,w=2)
plot!(ymk_MASG_7,xscale=:log,yscale=:log,label="intermediate variable",w=2)
plot!(grad_MASG_7, label="gradient",w=.5,alpha=0.5,color=:black)
plot!(f_MASG_7,xscale=:log,yscale=:log, label="Objective",w=2)
vline!(cumulative, label="stage",w=1,color=:grey)
xlabel!("Iteration Count")
ylabel!("log(f)")
savefig("./figures/illustration.png")
plot(α_MASG,label="step size")
xlabel!("stage")
savefig("./figures/stepsize.png")
## MASG Thm 4
nk_4 = compute_nk(κ, 1, norm(β,2), σ, 4.0)
x_MASG_4, f_MASG_4, ymk_MASG_4 = MASG(f, ∇f, nk_4, α_MASG, μ, σ)

# Vanilla Gradient Descent
x_GD, f_GD = GD(f, ∇f, σ)

# Vanilla Accelerated Gradient Descent
x_AGD, f_AGD = AGD(f, ∇f, σ)

# μAGD+
μAGD_a = compute_μAGD_a()
μAGD_A = compute_μAGD_A([1:η;] ./ L)
x_μAGD_0, f_μAGD_0 = μAGD(f,∇f,μAGD_A,[1:η;] ./ L,σ,0)
x_μAGD_2, f_μAGD_2 = μAGD(f,∇f,μAGD_A,[1:η;] ./ L,σ,2)

plot(f_GD, label="GD",legend=:topleft, xscale = :log, yscale = :log,w=2,linestyle=:solid,color=:red)
plot!(f_AGD, label="AGD", w=2, linestyle=:solid,color=:purple)
plot!(f_μAGD_0, label= "AGD+", xscale = :log, yscale = :log, color=:blue,w=2,linestyle=:dashdot)
plot!(f_μAGD_2, label= "AGD+v.2", xscale = :log, yscale = :log, color=:grey,w=2,linestyle=:dashdot)
plot!(f_MASG_7, label="MASG",w=2, linestyle=:dash, color=:green)
plot!(f_MASG_4, label="MASG*",w=2, linestyle=:dash, color=:black)
xlabel!("Iteration")
ylabel!("f")
savefig("./figures/10000_4.png")

max_stars_repo_name: chenyangzhu/noisy-gradients

max_issues_repo_name: chenyangzhu/noisy-gradients

max_forks_repo_name: chenyangzhu/noisy-gradients


text:
module OneClassSampling

using Memento

const LOGGER = getlogger(@__MODULE__)

function __init__()
    Memento.register(LOGGER)
    Memento.config!(LOGGER, "warn"; fmt="[{level} | {name}]: {msg}")
end

include("sampler_base.jl")

export
    Sampler,
    sample,
    SamplingException,

    RAPID,

    HDS,
    RandomRatioSampler, RandomNSampler,
    CSSVDD, BPS, DAEDS, DBSRSVDD, FBPE, HSC, HSR, IESRSVDD, KFNCBD, KMSVDD, NDPSR, RSSVDD,

    PreFilteringWrapper,

    KDECache,
    kde,

    ThresholdStrategy,
    FixedThresholdStrategy, NumLabelTresholdStrategy, MaxDensityThresholdStrategy,
    OutlierPercentageThresholdStrategy, GroundTruthThresholdStrategy,
    calculate_threshold, calculate_threshold!, split_masks, predict,
    PWC

end


text:
@testset "Grade projections" begin
    @testset "Projection to grade $r" for r ∈ 0:4 begin
            Pᵣ = x -> grade_projection(x, Val(r))
            Aᵣ, Bᵣ = As[r+1], Bs[r+1]
            @test Pᵣ(SA) ≈ As[r+1]
            @test Pᵣ(SB) ≈ Bs[r+1]
            @test Pᵣ(SC) ≈ Cs[r+1]
            @test Pᵣ(A * B) ≈ (-1) ^ (r * (r - 1) ÷ 2) * Pᵣ(reverse(A * B))
            @test Pᵣ(A * Bᵣ * C) ≈ Pᵣ(reverse(C) * Bᵣ * reverse(A))
            for s ∈ 0:4
                Bₛ = Bs[s+1]
                @test Pᵣ(Aᵣ * Bₛ) ≈ Pᵣ(reverse(Bₛ) * Aᵣ) ≈ (-1) ^ (s * (s - 1) ÷ 2) * Pᵣ(Bₛ * Aᵣ)
            end
        end
    end
end


text:
function bob(stimulus::AbstractString)
    stimulus = strip(stimulus)

    if issilence(stimulus)
        return "Fine. Be that way!"
    elseif isshouting(stimulus)
        return "Whoa, chill out!"
    elseif isquestion(stimulus)
        return "Sure."
    else
        return "Whatever."
    end
end

issilence(stimulus::AbstractString) = isempty(stimulus)
isquestion(stimulus::AbstractString) = endswith(stimulus, '?')

function isshouting(stimulus::AbstractString)
    all(isupper, stimulus) && return true
    !any(isalpha, stimulus) && return false

    for c in stimulus
        # ignore all non-letter chars
        if isalpha(c) && !isupper(c)
            return false
        end
    end

    return true
end