text:
function decode_integer(buf::IOBuffer, prefix_size)
    mask = 0xFF
    if prefix_size != 8
        mask = UInt8(1) << prefix_size - 1
    end
    value = read(buf, UInt8) & mask

    if value < mask
        return value
    end

    # The value does not fit into the prefix bits.
    total = 1
    m = 0
    octet_limit = 5

    while true
        b = read(buf, UInt8)
        total += 1
        value += (UInt(b & 127)) * (1 << m)
        m += 7

        if b & 128 != 128
            return value
        end

        #TODO handle errors here
    end
end

function decode_string(buf::IOBuffer)
    huffman = buf.data[buf.ptr] & 128 == 128
    len = decode_integer(buf, 7)
    str = read(buf, len)
    if huffman
        str = huffman_decode_bytes(str)
    end

    return str
end

function decode_literal(table::DynamicTable, buf::IOBuffer, index::Bool)
    prefix = index ? 6 : 4
    table_index = decode_integer(buf, prefix)

    name =
        if table_index == 0
            decode_string(buf)
        else
            get_header(table, table_index)[1]
        end

    value = decode_string(buf)

    if index
        add_header!(table, (name, value))
    end

    return (name, value)
end

function decode_indexed(table::DynamicTable, buf::IOBuffer)
    index = decode_integer(buf, 7)
    if index == 0
        throw(DecodeError("Index must not be zero."))
    end
    return get_header(table, index)
end

add_hdr(headers::Vector{Tuple{String,String}}, header::HeaderBinary) = push!(headers, (String(copy(header[1])), String(copy(header[2]))))

function decode(table::DynamicTable, buf::IOBuffer)
    headers = Vector{Tuple{String,String}}()

    while !eof(buf)
        initial_octet = buf.data[buf.ptr]

        if initial_octet & 128 == 128 # Indexed
            add_hdr(headers, decode_indexed(table, buf))
        elseif initial_octet & 64 == 64 # Literal with incremental indexing
            add_hdr(headers, decode_literal(table, buf, true))
        elseif initial_octet & 32 == 32 # Size update
            new_size = decode_integer(buf, 5)
            set_max_table_size!(table, new_size)
        elseif initial_octet & 16 == 16 # Literal never indexed
            add_hdr(headers, decode_literal(table, buf, false))
        else # Literal without indexing
            add_hdr(headers, decode_literal(table, buf, false))
        end
    end

    return headers
end

ext: jl


text:
# RUN: julia -e "import Brutus; Brutus.lit(:emit_lowered)" --startup-file=no %s 2>&1 | FileCheck %s

select(c) = 1 + (c ? 2 : 3)
emit(select, Bool)

# CHECK: module  {
# CHECK-NEXT:   func nested @"Tuple{typeof(Main.select), Bool}"(%arg0: !jlir<"typeof(Main.select)">, %arg1: i1) -> i64 attributes {llvm.emit_c_interface} {
# CHECK-NEXT:     %c2_i64 = constant 2 : i64
# CHECK-NEXT:     %c3_i64 = constant 3 : i64
# CHECK-NEXT:     %c1_i64 = constant 1 : i64
# CHECK-NEXT:     %0 = select %arg1, %c2_i64, %c3_i64 : i64
# CHECK-NEXT:     %1 = addi %0, %c1_i64 : i64
# CHECK-NEXT:     return %1 : i64
# CHECK-NEXT:   }
# CHECK-NEXT: }

# CHECK:   llvm.func @"Tuple{typeof(Main.select), Bool}"(%arg0: !llvm.ptr<struct<"struct_jl_value_type", opaque>>, %arg1: i1) -> i64 attributes {llvm.emit_c_interface, sym_visibility = "nested"} {
# CHECK-NEXT:     %0 = llvm.mlir.constant({{[0-9]+}} : i64) : i64
# CHECK-NEXT:     %1 = llvm.mlir.constant({{[0-9]+}} : i64) : i64
# CHECK-NEXT:     %2 = llvm.mlir.constant({{[0-9]+}} : i64) : i64
# CHECK-NEXT:     %3 = llvm.select %arg1, %0, %1 : i1, i64
# CHECK-NEXT:     %4 = llvm.add %3, %2  : i64
# CHECK-NEXT:     llvm.return %4 : i64
# CHECK-NEXT:   }
# CHECK-NEXT:   llvm.func @"_mlir_ciface_Tuple{typeof(Main.select), Bool}"(%arg0: !llvm.ptr<struct<"struct_jl_value_type", opaque>>, %arg1: i1) -> i64 attributes {llvm.emit_c_interface, sym_visibility = "nested"} {
# CHECK-NEXT:     %0 = llvm.call @"Tuple{typeof(Main.select), Bool}"(%arg0, %arg1) : (!llvm.ptr<struct<"struct_jl_value_type", opaque>>, i1) -> i64
# CHECK-NEXT:     llvm.return %0 : i64
# CHECK-NEXT:   }
# CHECK-NEXT: }

max_stars_repo_name: femtomc/brutus

max_issues_repo_name: femtomc/brutus

max_forks_repo_name: femtomc/brutus


text:
# ---
# title: 102. Binary Tree Level Order Traversal
# id: problem102
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Tree, Breadth-first Search
# link: <https://leetcode.com/problems/binary-tree-level-order-traversal/description/>
# hidden: true
# ---
#
# Given a binary tree, return the _level order_ traversal of its nodes' values.
# (ie, from left to right, level by level).
#
# For example:
# Given binary tree `[3,9,20,null,null,15,7]`,
#
#
#
#         3
#        / \
#       9  20
#         /  \
#        15   7
#
#
# return its level order traversal as:
#
#
#
#     [
#       [3],
#       [9,20],
#       [15,7]
#     ]
#
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using Test
using TreeTools
using TreeKnit, TreeKnit.SplitGraph

println("##### main #####")

nwk1 = "((1:1,2:1):1,((3:1,4:1):1,5:1,6:1):1)"
nwk2 = "((6:1,2:1):1,((3:1,4:1):1,5:1,1:1):1)"

t1 = node2tree(TreeTools.parse_newick(nwk1))
t2 = node2tree(TreeTools.parse_newick(nwk2))

## Testing runopt
ct1 = deepcopy(t1)
ct2 = deepcopy(t2)
out = TreeKnit.runopt(t1, t2; likelihood_sort=false, verbose=true)
##
trees = Dict(1=>deepcopy(t1), 2=>deepcopy(t2));
@testset "Pre-resolve" begin
	MCCs = TreeKnit.computeMCCs!(
		trees, OptArgs(likelihood_sort=false, γ=3);
		preresolve = true,
	)
	@test MCCs[1,2] == [["1"], ["2"], ["6"], ["3", "4", "5"]]
end

trees2 = Dict(1=>deepcopy(t1), 2=>deepcopy(t2));
solutions = [
	[["1"], ["2"], ["3", "4", "5", "6"]],
	[["2"], ["6"], ["1", "3", "4", "5"]],
	[["1"], ["6"], ["2", "3", "4", "5"]],
	[["1"], ["2", "6"], ["3", "4", "5"]],
	[["6"], ["1", "2"], ["3", "4", "5"]],
]
@testset "Dyn-resolve" begin
	MCCs = TreeKnit.computeMCCs!(
		trees2, OptArgs(likelihood_sort=false, γ=3);
		preresolve=false,
	)
	@test in(MCCs[1,2], solutions)
end

max_stars_repo_name: PierreBarrat/OldTreeKnit

max_issues_repo_name: PierreBarrat/OldTreeKnit

max_forks_repo_name: PierreBarrat/OldTreeKnit


text:
export ScalePointCloud,
    RotatePointCloud,
    ReAlignPointCloud,
    NormalizePointCloud,
    ScaleTriMesh,
    RotateTriMesh,
    ReAlignTriMesh,
    NormalizeTriMesh,
    TranslateTriMesh,
    OffsetTriMesh,
    TriMeshToVoxelGrid,
    PointCloudToVoxelGrid,
    VoxelGridToTriMesh,
    PointCloudToTriMesh,
    TriMeshToPointCloud,
    VoxelGridToPointCloud

"""
    ScalePointCloud(factor::Number; inplace::Bool=true)

Scale PointCloud with a given scaling factor `factor`.

`factor` should be strictly greater than `0.0` for obvious reason.
`inplace` is optional keyword argument, to make transformation in-place.
If `inplace` is set to `false`, it will create deepcopy of PointCloud.
Given `factor`, this transform scale each point in PointCloud, ie. `point = point * factor`

See also: [`scale`](@ref), [`scale!`](@ref)
"""
struct ScalePointCloud <: AbstractTransform
    factor::Float32
    inplace::Bool
end

function ScalePointCloud(factor::Number; inplace::Bool = true)
    factor > 0.0 || error("factor must be greater than 0.0")
    ScalePointCloud(Float32(factor), inplace)
end

@functor ScalePointCloud

function (t::ScalePointCloud)(pcloud::PointCloud)
    t.inplace || (pcloud = deepcopy(pcloud))
    scale!(pcloud, t.factor)
    return pcloud
end

Base.show(io::IO, t::ScalePointCloud) =
    print(io, "$(typeof(t))(factor=$(t.factor); inplace=$(t.inplace))")

"""
    RotatePointCloud(rotmat::AbstractArray{<:Number,2}; inplace::Bool=true)

Rotate PointCloud with a given rotation matrix `rotmat`.

`rotmat` must be `AbstractArray{<:Number,2}` of size `(3,3)`.
`inplace` is optional keyword argument, to make transformation in-place
If `inplace` is set to `false`, it will create deepcopy of PointCloud.
Given `rotmat`, this transform will rotate each point coordinates (ie. x,y,z) in PointCloud.

See also: [`rotate`](@ref), [`rotate!`](@ref)
"""
struct RotatePointCloud <: AbstractTransform
    rotmat::AbstractArray{Float32,2}
    inplace::Bool
end

function RotatePointCloud(rotmat::AbstractArray{<:Number,2}; inplace::Bool = true)
    size(rotmat) == (3, 3) ||
        error("rotmat must be (3,3) array, but instead got $(size(rotmat)) array")
    return RotatePointCloud(Float32.(rotmat), inplace)
end

@functor RotatePointCloud

function (t::RotatePointCloud)(pcloud::PointCloud)
    t.inplace || (pcloud = deepcopy(pcloud))
    rotate!(pcloud, t.rotmat)
    return pcloud
end

Base.show(io::IO, t::RotatePointCloud) =
    print(io, "$(typeof(t))(rotmat; inplace=$(t.inplace))")

"""
    ReAlignPointCloud(target::PointCloud; inplace::Bool=true)
    ReAlignPointCloud(target::AbstractArray{<:Number, 2}; inplace::Bool=true)

Re-Align PointCloud to axis aligned bounding box of `target` PointCloud.

`input` PointCloud and `target` PointCloud should be of same size.
`inplace` is optional keyword argument, to make transformation in-place
If `inplace` is set to `false`, it will create deepcopy of PointCloud.

See also: [`realign`](@ref), [`realign!`](@ref)
"""
struct ReAlignPointCloud <: AbstractTransform
    t_min::AbstractArray{Float32,2}
    t_max::AbstractArray{Float32,2}
    inplace::Bool
end

function ReAlignPointCloud(target::PointCloud, index::Number = 1; inplace::Bool = true)
    points = target[index]
    t_min = minimum(points, dims = 2)
    t_max = maximum(points, dims = 2)
    ReAlignPointCloud(t_min, t_max, inplace)
end

ReAlignPointCloud(target::AbstractArray{<:Number}; inplace::Bool = true) =
    ReAlignPointCloud(PointCloud(target), inplace = inplace)

@functor ReAlignPointCloud

function (t::ReAlignPointCloud)(pcloud::PointCloud)
    t.inplace || (pcloud = deepcopy(pcloud))
    realign!(pcloud, t.t_min, t.t_max)
    return pcloud
end

Base.show(io::IO, t::ReAlignPointCloud) =
    print(io, "$(typeof(t))(target=PointCloud(...); inplace=$(t.inplace))")

"""
    NormalizePointCloud(; inplace::Bool=true)

Normalize PointCloud with mean at origin and unit standard deviation.

`inplace` is optional keyword argument, to make transformation in-place
If `inplace` is set to `false`, it will create deepcopy of PointCloud.

See also: [`normalize`](@ref), [`normalize!`](@ref)
"""
struct NormalizePointCloud <: AbstractTransform
    inplace::Bool
end

NormalizePointCloud(; inplace::Bool = true) = NormalizePointCloud(inplace)

@functor NormalizePointCloud

function (t::NormalizePointCloud)(pcloud::PointCloud)
    t.inplace || (pcloud = deepcopy(pcloud))
    normalize!(pcloud)
    return pcloud
end

Base.show(io::IO, t::NormalizePointCloud) = print(io, "$(typeof(t))(;inplace=$(t.inplace))")

#TODO: Add support for multidimension factor
"""
    ScaleTriMesh(factor::Number; inplace::Bool=true)

Scale TriMesh with a given scaling factor `factor`.

`factor` should be strictly greater than `0.0` for obvious reason.
`inplace` is optional keyword argument, to make transformation in-place.
If `inplace` is set to `false`, it will create deepcopy of TriMesh.
Given `factor`, this transform scale each vertices in TriMesh, ie. `point = point * factor`

See also: [`scale`](@ref), [`scale!`](@ref)
"""
struct ScaleTriMesh <: AbstractTransform
    factor::Float32
    inplace::Bool
end

function ScaleTriMesh(factor::Number; inplace::Bool = true)
    factor > 0.0 || error("factor must be greater than 0.0")
    ScaleTriMesh(Float32(factor), inplace)
end

@functor ScaleTriMesh

function (t::ScaleTriMesh)(m::TriMesh)
    t.inplace || (m = deepcopy(m))
    scale!(m, t.factor)
    return m
end

Base.show(io::IO, t::ScaleTriMesh) =
    print(io, "$(typeof(t))(factor=$(t.factor); inplace=$(t.inplace))")

"""
    RotateTriMesh(rotmat::AbstractArray{<:Number,2}; inplace::Bool=true)

Rotate vertices of TriMesh with a given rotation matrix `rotmat`.

`rotmat` must be `AbstractArray{<:Number,2}` of size `(3,3)`.
`inplace` is optional keyword argument, to make transformation in-place
If `inplace` is set to `false`, it will create deepcopy of TriMesh.
Given `rotmat`, this transform will rotate each vertices coordinates (ie. x,y,z) in TriMesh.

See also: [`rotate`](@ref), [`rotate!`](@ref)
"""
struct RotateTriMesh <: AbstractTransform
    rotmat::AbstractArray{Float32,2}
    inplace::Bool
end

function RotateTriMesh(rotmat::AbstractArray{<:Number,2}; inplace::Bool = true)
    size(rotmat) == (3, 3) ||
        error("rotmat must be (3,3) array, but instead got $(size(rotmat)) array")
    return RotateTriMesh(Float32.(rotmat), inplace)
end

@functor RotateTriMesh

function (t::RotateTriMesh)(m::TriMesh)
    t.inplace || (m = deepcopy(m))
    rotate!(m, t.rotmat)
    return m
end

Base.show(io::IO, t::RotateTriMesh) =
    print(io, "$(typeof(t))(rotmat; inplace=$(t.inplace))")

"""
    ReAlignTriMesh(target::TriMesh; inplace::Bool=true)

Re-Align TriMesh to axis aligned bounding box of mesh at `index` in TriMesh `target`.

`inplace` is optional keyword argument, to make transformation in-place
If `inplace` is set to `false`, it will create deepcopy of TriMesh.

See also: [`realign`](@ref), [`realign!`](@ref)
"""
struct ReAlignTriMesh <: AbstractTransform
    t_min::AbstractArray{Float32,2}
    t_max::AbstractArray{Float32,2}
    inplace::Bool
end

function ReAlignTriMesh(target::TriMesh, index::Integer = 1; inplace::Bool = true)
    verts = get_verts_list(target)[index]
    t_min = minimum(verts, dims = 2)
    t_max = maximum(verts, dims = 2)
    ReAlignTriMesh(t_min, t_max, inplace)
end

@functor ReAlignTriMesh

function (t::ReAlignTriMesh)(m::TriMesh)
    t.inplace || (m = deepcopy(m))
    realign!(m, t.t_min, t.t_max)
    return m
end

Base.show(io::IO, t::ReAlignTriMesh) =
    print(io, "$(typeof(t))(target=TriMesh(...); inplace=$(t.inplace))")

"""
    NormalizeTriMesh(; inplace::Bool=true)

Normalize TriMesh with mean at origin and unit standard deviation.

`inplace` is optional keyword argument, to make transformation in-place
If `inplace` is set to `false`, it will create deepcopy of TriMesh.

See also: [`normalize`](@ref), [`normalize!`](@ref)
"""
struct NormalizeTriMesh <: AbstractTransform
    inplace::Bool
end

NormalizeTriMesh(; inplace::Bool = true) = NormalizeTriMesh(inplace)

@functor NormalizeTriMesh

function (t::NormalizeTriMesh)(m::TriMesh)
    t.inplace || (m = deepcopy(m))
    normalize!(m)
    return m
end

Base.show(io::IO, t::NormalizeTriMesh) = print(io, "$(typeof(t))(;inplace=$(t.inplace))")

"""
    TranslateTriMesh(vector::AbstractArray{<:Number}; inplace::Bool=true)

Translate TriMesh with given translation `vector`.

`inplace` is optional keyword argument, to make transformation in-place
If `inplace` is set to `false`, it will create deepcopy of TriMesh.

See also: [`translate`](@ref), [`translate!`](@ref)
"""
struct TranslateTriMesh <: AbstractTransform
    vector::AbstractArray{Float32,1}
    inplace::Bool
end

function TranslateTriMesh(vector::AbstractArray{<:Number,1}; inplace::Bool = true)
    (size(vector) == (3,)) ||
        error("vector must be (3, ), but instead got $(size(vector)) array")
    TranslateTriMesh(Float32.(vector), inplace)
end

TranslateTriMesh(vector::Number; inplace::Bool = true) =
    TranslateTriMesh(fill(vector, (3,)); inplace = inplace)

@functor TranslateTriMesh

function (t::TranslateTriMesh)(m::TriMesh)
    t.inplace || (m = deepcopy(m))
    translate!(m, t.vector)
    return m
end

Base.show(io::IO, t::TranslateTriMesh) =
    print(io, "$(typeof(t))(vector=$(t.vector);inplace=$(t.inplace))")

"""
    OffsetTriMesh(offset_verts::AbstractArray{<:Number,2}; inplace::Bool=true)

Add offset to the TriMesh by given offset vertices `offset_verts`

`inplace` is optional keyword argument, to make transformation in-place
If `inplace` is set to `false`, it will create deepcopy of TriMesh.

See also: [`offset`](@ref), [`offset!`](@ref)
"""
struct OffsetTriMesh <: AbstractTransform
    offset_verts::AbstractArray{Float32,2}
    inplace::Bool
end

OffsetTriMesh(offset_verts::AbstractArray{<:Number,2}; inplace::Bool = true) =
    OffsetTriMesh(Float32.(offset_verts), inplace)

@functor OffsetTriMesh

function (t::OffsetTriMesh)(m::TriMesh)
    t.inplace || (m = deepcopy(m))
    offset!(m, t.offset_verts)
    return m
end

Base.show(io::IO, t::OffsetTriMesh) =
    print(io, "$(typeof(t))(offset_verts; inplace=$(t.inplace))")

"""
    TriMeshToVoxelGrid(resolution::Int=32)

Converts a TriMesh to VoxelGrid having specified `resolution`.

See also: [`TriMesh`](@ref), [`VoxelGrid`](@ref)
"""
struct TriMeshToVoxelGrid <: AbstractTransform
    resolution::Int
    TriMeshToVoxelGrid(res::Int = 32) = new(res)
end

(t::TriMeshToVoxelGrid)(m::TriMesh) = VoxelGrid(m, t.resolution)

Base.show(io::IO, t::TriMeshToVoxelGrid) =
    print(io, "$(typeof(t))(resolution=$(t.resolution))")

"""
    PointCloudToVoxelGrid(resolution::Int=32)

Converts a PointCloud to VoxelGrid having specified `resolution`.

See also: [`PointCloud`](@ref), [`VoxelGrid`](@ref)
"""
struct PointCloudToVoxelGrid <: AbstractTransform
    resolution::Int
    PointCloudToVoxelGrid(res::Int = 32) = new(res)
end

(t::PointCloudToVoxelGrid)(p::PointCloud) = VoxelGrid(p, t.resolution)

Base.show(io::IO, t::PointCloudToVoxelGrid) =
    print(io, "$(typeof(t))(resolution=$(t.resolution))")

"""
    VoxelGridToTriMesh(; threshold=0.5, algo=:Exact)

Converts a VoxelGrid to TriMesh.

`threshold` is the threshold from which to make binary voxels, and `algo`
is the mode to be used to convert binary voxels. Available `algo` are
[:Exact, :MarchingCubes, :MarchingTetrahedra, :NaiveSurfaceNets].

See also: [`TriMesh`](@ref), [`VoxelGrid`](@ref)
"""
struct VoxelGridToTriMesh <: AbstractTransform
    threshold::Float32
    algo::Symbol
end

function VoxelGridToTriMesh(; thresh::Number = 0.5f0, algo = :MarchingCubes)
    algo in [:Exact, :MarchingCubes, :MarchingTetrahedra, :NaiveSurfaceNets] || error(
        "given algo=$(algo) is not supported. Accepted algos are {:Exact,:MarchingCubes, :MarchingTetrahedra, :NaiveSurfaceNets}.",
    )
    (0 <= thresh <= 1) || error("given threshold=$(thresh) is not between [0,1]")
    return VoxelGridToTriMesh(Float32(thresh), algo)
end

(t::VoxelGridToTriMesh)(v::VoxelGrid) = TriMesh(v; thresh = t.threshold, algo = t.algo)

Base.show(io::IO, t::VoxelGridToTriMesh) =
    print(io, "$(typeof(t))((threshold=$(t.threshold), algo=$(t.algo))")

"""
    PointCloudToTriMesh(resolution::Int=32)

Converts a PointCloud to TriMesh having specified `resolution`.

See also: [`PointCloud`](@ref), [`TriMesh`](@ref)
"""
struct PointCloudToTriMesh <: AbstractTransform
    resolution::Int
    algo::Symbol
    function PointCloudToTriMesh(res::Int = 32; algo = :MarchingCubes)
        algo in [:Exact, :MarchingCubes, :MarchingTetrahedra, :NaiveSurfaceNets] || error(
            "given algo=$(algo) is not supported. Accepted algos are {:Exact,:MarchingCubes, :MarchingTetrahedra, :NaiveSurfaceNets}.",
        )
        return new(res, algo)
    end
end

(t::PointCloudToTriMesh)(p::PointCloud) = TriMesh(p, t.resolution)

Base.show(io::IO, t::PointCloudToTriMesh) =
    print(io, "$(typeof(t))(resolution=$(t.resolution))")

"""
    TriMeshToPointCloud(npoints::Int=1000)

Converts a TriMesh to PointCloud having `npoints`.

See also: [`TriMesh`](@ref), [`PointCloud`](@ref)
"""
struct TriMeshToPointCloud <: AbstractTransform
    npoints::Int
    function TriMeshToPointCloud(npoints::Int = 1024)
        npoints >= 0 || error("npoints cannot be less than 0")
        return new(npoints)
    end
end

(t::TriMeshToPointCloud)(m::TriMesh) = PointCloud(m, t.npoints)

Base.show(io::IO, t::TriMeshToPointCloud) = print(io, "$(typeof(t))(npoints=$(t.npoints))")

"""
    VoxelGridToPointCloud(npoints::Int=1000; threshold=0.5, algo=:Exact)

Converts a VoxelGrid to PointCloud having `npoints`.

`threshold` is the threshold from which to make binary voxels, and `algo`
is the mode to be used to convert binary voxels. Available `algo` are
[:Exact, :MarchingCubes, :MarchingTetrahedra, :NaiveSurfaceNets].

See also: [`PointCloud`](@ref), [`VoxelGrid`](@ref)
"""
struct VoxelGridToPointCloud <: AbstractTransform
    npoints::Int
    threshold::Float32
    algo::Symbol
end

function VoxelGridToPointCloud(
    npoints::Int = 1024;
    thresh::Number = 0.5f0,
    algo::Symbol = :MarchingCubes,
)
    npoints >= 0 || error("npoints cannot be less than 0")
    algo in [:Exact, :MarchingCubes, :MarchingTetrahedra, :NaiveSurfaceNets] || error(
        "given algo=$(algo) is not supported. Accepted algos are {:Exact,:MarchingCubes, :MarchingTetrahedra, :NaiveSurfaceNets}.",
    )
    (0 <= thresh <= 1) || error("given threshold=$(thresh) is not between [0,1]")
    return VoxelGridToPointCloud(npoints, Float32(thresh), algo)
end

(t::VoxelGridToPointCloud)(v::VoxelGrid) =
    PointCloud(v, t.npoints; thresh = t.threshold, algo = t.algo)

Base.show(io::IO, t::VoxelGridToPointCloud) = print(
    io,
    "$(typeof(t))(npoints=$(t.npoints), threshold=$(t.threshold), algo=$(t.algo))",
)


text:
using DataFrames

abstract type NetworkNode{A} end

@def network_node begin
    node_id::String
    area::Float64  # area in km^2
end

# function get_climate_data(node::NetworkNode, climate_data::DataFrame)
#     tgt::String = node.node_id
#     rain_prefix = "pr_"
#     et_prefix = "wvap_"

#     rain = climate_data[Symbol(rain_prefix*tgt)]
#     et = climate_data[Symbol(et_prefix*tgt)]

#     return rain, et
# end

max_stars_repo_name: ConnectedSystems/camgraph

max_issues_repo_name: ConnectedSystems/camgraph

max_forks_repo_name: ConnectedSystems/camgraph


text:
using Manifolds, Documenter

makedocs(
    # for development, we disable prettyurls
    format = Documenter.HTML(prettyurls = false),
    modules = [Manifolds],
    sitename = "Manifolds",
    pages = [
        "Home" => "index.md",
        "Manifolds" => [
            "Basic manifolds" => [
                "Euclidean" => "manifolds/euclidean.md",
                "Rotations" => "manifolds/rotations.md",
                "Sphere" => "manifolds/sphere.md"
            ],
            "Combined manifolds" => [
                "Product manifold" => "manifolds/product.md"
            ],
            "Manifold decorators" => [
                "Array manifold" => "manifolds/array.md",
                "Metric manifold" => "manifolds/metric.md"
            ]
        ],
        "Distributions" => "distributions.md",
        "Library" => [
            "Internals" => "lib/internals.md"
        ]
    ]
)


text:
f(x, y) = 2x + y           # call-like form on left-hand side
g = x -> x * x             # arrow notation
@assert g(f(2, 9)) == 169

max_stars_repo_name: idrougge/ple

max_issues_repo_name: idrougge/ple

max_forks_repo_name: idrougge/ple


text:
function StructDefinition{false}(spec::SpecStruct)
    p = Dict(
        :category => :struct,
        :decl => :($(struct_name(spec)) <: VulkanStruct{$(needs_deps(spec))}),
    )
    p[:fields] = [:(vks::$(spec.name))]
    needs_deps(spec) && push!(p[:fields], :(deps::Vector{Any}))
    for member in parent_handles(spec)
        handle_type = remove_vk_prefix(member.type)
        name = wrap_identifier(member)
        field_type = is_optional(member) ? :(OptionalPtr{$handle_type}) : handle_type
        push!(p[:fields], :($name::$field_type))
    end
    StructDefinition{false}(spec, p)
end

function Constructor(def::StructDefinition{false})
    spec = def.spec
    cconverted_members = getindex(children(spec), findall(is_semantic_ptr, children(spec).type))
    cconverted_ids = map(wrap_identifier, cconverted_members)
    p = Dict(
        :category => :function,
        :name => name(def),
        :relax_signature => true,
    )
    if needs_deps(spec)
        p[:body] = quote
            $(
                (
                    :($id = cconvert($(m.type), $id)) for
                    (m, id) ∈ zip(cconverted_members, cconverted_ids)
                )...
            )
            deps = Any[$((cconverted_ids)...)]
            vks = $(spec.name)($(map(vk_call, children(spec))...))
            $(name(def))(vks, deps, $(wrap_identifier.(parent_handles(spec))...))
        end
    else
        p[:body] = :($(name(def))($(spec.name)($(map(vk_call, children(spec))...)), $(wrap_identifier.(parent_handles(spec))...)))
    end
    potential_args = filter(x -> x.type ≠ :VkStructureType, children(spec))
    add_func_args!(p, spec, potential_args)
    Constructor(p, def, def.spec)
end

StructDefinition{true}(spec::Spec) = StructDefinition{true}(StructDefinition{false}(spec))

function StructDefinition{true}(def::StructDefinition{false})
    (; spec) = def
    p = Dict(
        :category => :struct,
        :decl => :($(struct_name(spec, true)) <: HighLevelStruct),
        :fields => hl_fields(spec),
    )
    StructDefinition{true}(spec, p)
end

function hl_fields(spec::Union{SpecStruct,SpecUnion})
    fields = Expr[]
    for member in filter(!drop_field, children(spec))
        T = hl_type(member)
        if hl_is_optional(member)
            T = :(OptionalPtr{$T})
        end
        push!(fields, :($(wrap_identifier(member))::$T))
    end
    fields
end

drop_field(x::Spec) = drop_arg(x) || x.name == :sType

function Constructor(T::StructDefinition{false}, x::StructDefinition{true})
    (; spec) = T
    p = Dict(
        :category => :function,
        :name => name(T),
        :args => [:(x::$(name(x)))],
        :short => true,
    )
    args = Expr[]
    kwargs = Expr[]
    for member in filter(!drop_field, children(spec))
        id = wrap_identifier(member)
        id_deref = :(x.$id)
        Tsym = hl_type(member)
        _Tsym = @match innermost_type(Tsym) begin
            GuardBy(is_hl) => @match Tsym begin
                :(Vector{$T}) => :(Vector{$(Symbol(:_, T))})
                :(NTuple{$N,$T}) => :(NTuple{$N,$(Symbol(:_, T))})
                ::Symbol => Symbol(:_, Tsym)
                _ => error("Unhandled conversion for type $Tsym")
            end
            _ => nothing
        end
        ex = isnothing(_Tsym) ? id_deref : :(convert_nonnull($_Tsym, $id_deref))
        if is_optional(member)
            push!(kwargs, ex == id_deref ? ex : Expr(:kw, id, ex))
        else
            push!(args, ex)
        end
    end
    p[:body] = reconstruct_call(Dict(:name => p[:name], :args => args, :kwargs => kwargs))
    Constructor(p, T, x)
end

function Constructor(T::StructDefinition{true}, x::SpecStruct)
    p = Dict(
        :category => :function,
        :name => name(T),
        :args => [:(x::$(x.name))],
        :body => :($(name(T))($(filter(!isnothing, from_vk_call.(filter(!drop_field, children(x))))...))),
        :short => true,
    )
    :pNext in x.members.name && push!(p[:args], :(next_types::Type...))
    Constructor(p, T, x)
end

function Constructor(T::StructDefinition{true}, x::StructDefinition{false})
    (; spec) = x
    p = Dict(
        :category => :function,
        :name => name(T),
        :args => [:(x::$(name(x)))],
        :short => !needs_deps(spec),
    )
    :pNext in spec.members.name && push!(p[:args], :(next_types::Type...))
    p[:body] = if needs_deps(spec)
        quote
            (; deps) = x
            GC.@preserve deps $(name(T))(x.vks, next_types...)
        end
    else
        Expr(:call, name(T), :(x.vks))
    end
    Constructor(p, T, x)
end

function Convert(T::StructDefinition{false}, x::StructDefinition{true})
    p = Dict(
        :category => :function,
        :name => :convert,
        :args => [:(T::Type{$(name(T))}), :(x::$(name(x)))],
        :body => :(T(x)),
        :short => true,
    )
    Convert(T, x, p)
end

function Constructor(def::StructDefinition{true})
    spec = def.spec
    p = Dict(
        :category => :function,
        :name => name(def),
        :args => [],
        :kwargs => [],
        :short => true,
        :relax_signature => true,
    )
    args = []
    for member in filter(!drop_field, children(spec))
        id = wrap_identifier(member)
        if is_optional(member)
            push!(p[:kwargs], Expr(:kw, id, hl_default(member)))
        else
            push!(p[:args], :($id::$(signature_type(hl_type(member)))))
        end
        push!(args, id)
    end
    p[:body] = reconstruct_call(Dict(:name => p[:name], :args => args))
    Constructor(p, def, def.spec)
end

function hl_default(member::SpecStructMember)
    !embeds_sentinel(member) && return :C_NULL
    @match hl_type(member) begin
        :String => ""
        :(Vector{$t}) => :($t[])
        _ => default(member)
    end
end

hl_is_optional(member::SpecStructMember) = is_optional(member) && !embeds_sentinel(member) && member.name ≠ :pNext

function embeds_sentinel(member::SpecStructMember)
    @match hl_type(member) begin
        :(Vector{$_}) => return !is_optional(member) # optional means C_NULL can be provided, and may have different semantics than for an empty vector
        :String || :UInt32 || :UInt64 || :Int32 || :Int => return true
        _ => nothing
    end

    @match member.type begin
        GuardBy(in(spec_enums.name)) || GuardBy(in(spec_bitmasks.name)) || GuardBy(is_flag) => true
        _ => false
    end
end

function StructureType(spec::SpecStruct)
    stype = structure_types[spec.name]
    types = [spec.name, struct_name(spec.name), struct_name(spec.name, true)]
    ex = :(structure_type(@nospecialize(_::Union{$((:(Type{$T}) for T in types)...)})) = $stype)
    StructureType(spec, ex)
end

max_stars_repo_name: JuliaGPU/VulkanAbstraction

max_issues_repo_name: JuliaGPU/VulkanAbstraction

max_forks_repo_name: JuliaGPU/VulkanAbstraction


text:
using ApproximateBayesianComputing
const ABC = ApproximateBayesianComputing
using Distributions, Random

import ApproximateBayesianComputing.CompositeDistributions.CompositeDist
import ApproximateBayesianComputing.TransformedBetaDistributions.LinearTransformedBeta

Random.seed!(1234)

# Function to adjust originally proposed model parameters, so that they will be valid
function normalize_theta!(theta::Array)
  theta[2:end]./= sum(theta[2:end])
  return theta
end

# Set Prior for Population Parameters
max_rate = 3.0
theta_true = vcat(1.0,[1.0,2.0,3.0,4.0,5.0])
normalize_theta!(theta_true)
nbins = length(theta_true)-1

#param_prior = CompositeDist(vcat(Uniform(0.0,max_rate),ContinuousDistribution[Uniform(0.0,1.0) for i in 1:nbins]))
param_prior = CompositeDist([Uniform(0.0,max_rate),Dirichlet(ones(nbins))])

# Function to test if the proposed model parameter are valid
function is_valid(theta::Array)
    # Test total rate parameter from uniform distribution
    if !(minimum(param_prior.dist[1])<=theta[1]<=maximum(param_prior.dist[1]))
        return false
    end
    # Test parametrs of categorical distribution
    if !Distributions.isprobvec(theta[2:end])
        return false
    end
    return true
end

# Code to generate simulated data given array of model parameters
num_data_default = 1000
function draw_dirchlet_multinomial_with_poisson_rate(theta::Array)
   counts = zeros(Int64,length(theta)-1)
   n = rand(Poisson(theta[1]))
   d_cat = Categorical(theta[2:end])
   for i in 1:n
     counts[rand(d_cat)] += 1
   end
   return counts
end

function gen_data(theta::Array, n::Integer = num_data_default)
   data = Array{Int64}( undef, (length(theta)-1, n) )
   for i in 1:n
      data[:,i] .= draw_dirchlet_multinomial_with_poisson_rate(theta)
   end
   return data
end

function calc_mean_per_bin(data::Array{Int64,2})
   vec(mean(data,dims=2))
end

calc_dist_l1(x::Array{Float64},y::Array{Float64}) = sum(abs.(x.-y))/length(x)
calc_dist_l2(x::Array{Float64},y::Array{Float64}) = sum(abs2.(x.-y))/length(x)

#using Distributions
using SpecialFunctions
using Statistics

# https://en.wikipedia.org/wiki/Trigamma_function
function trigamma_x_gr_4(x::T) where T<: Real
   1/x + 0.5/x^2 + 1/(6*x^3) - 1/(30*x^5) + 1/(42*x^7) - 1/(30*x^9) + 5/(66*x^11) - 691/(2730*x^13) + 7/(6*x^15)
end

function trigamma_x_lt_4(x::T) where T<: Real
  n = floor(Int64,5-x)
  z = x+n
  val = trigamma_x_gr_4(z)
  for i in 1:n
    z -= 1
    val += 1/z^2
  end
  val
end

function trigamma(x::T) where T<: Real
   x >= 4 ? trigamma_x_gr_4(x) : trigamma_x_lt_4(x)
end

function var_weighted(x::AbstractArray{Float64,1}, w::AbstractArray{Float64,1} )
  #println("# size(x) = ",size(x), " size(w) = ", size(w)); flush(stdout)
  @assert(length(x)==length(w) )
  sumw = sum(w)
  @assert( sumw > 0. )
  if(sumw!= 1.0)
     w /= sum(w)
     sumw = 1.0
  end
  sumw2 = sum(w.*w)
  xbar =  sum(x.*w)
  covar = sum((x.-xbar).*(x.-xbar) .* w) * sumw/(sumw*sumw-sumw2)
end

function make_proposal_dist_multidim_beta(theta::AbstractArray{Float64,2}, weights::AbstractArray{Float64,1},  tau_factor::Float64; verbose::Bool = false)

  function mom_alpha(x_bar::T, v_bar::T) where T<: Real
    x_bar * (((x_bar * (1 - x_bar)) / v_bar) - 1)
  end
  function mom_beta(x_bar::T, v_bar::T) where T<: Real
    (1 - x_bar) * (((x_bar * (1 - x_bar)) / v_bar) - 1)
  end
  # For algorithm, see https://scholarsarchive.byu.edu/cgi/viewcontent.cgi?article=2613&context=etd
  function fit_beta_mle(x::AbstractArray{T,1}; tol::T = 1e-6, max_it::Int64 = 10, init_guess::AbstractArray{T,1} = Array{T}(undef,0), w::AbstractArray{T,1} = Array{T}(undef,0), verbose::Bool = false ) where T<: Real
    lnxbar =   length(w)>1 ? Statistics.mean(log.(x),AnalyticWeights(w)) : Statistics.mean(log.(x))
    ln1mxbar = length(w)>1 ? Statistics.mean(log.(1.0.-x),AnalyticWeights(w)) : Statistics.mean(log.(1.0.-x))

    function itterate( mle_guess::Vector{T} ) where T<:Real
       (alpha, beta) = (mle_guess[1], mle_guess[2])
       dgab = digamma(alpha+beta)
       g1 = dgab - digamma(alpha) + lnxbar
       g2 = dgab - digamma(beta) + ln1mxbar
       tgab = trigamma(alpha+beta)
       G = [dgab-trigamma(alpha) tgab; tgab tgab-trigamma(beta)]
       mle_guess -= G \ [g1, g2]
    end

    local mle_new
    if length(init_guess) != 2
       xbar = length(w)>1 ? Statistics.mean(x,AnalyticWeights(w)) : Statistics.mean(x)
       vbar = length(w)>1 ? Statistics.varm(x,xbar,AnalyticWeights(w)) : Statistics.varm(x,xbar)
       mle_new = (vbar < xbar*(1.0-xbar)) ? [mom_alpha(xbar, vbar), mom_beta(xbar,vbar)] : ones(T,2)
    else
       mle_new = init_guess
    end
    if verbose
       println("it = 0: ", mle_new)
    end
    if any(mle_new.<=zero(T))
       println("# Warning: mean= ", xbar, " var= ",var," (alpha,beta)_init= ",mle_new," invalid, reinitializing to (1,1)")
       verbose = true
       mle_new = ones(T,2)
    end
    for i in 1:max_it
       mle_old = mle_new
       mle_new = itterate( mle_old )
       epsilon = max(abs.(mle_old.-mle_new))
       if verbose
          println("# it = ", i, ": ", mle_new, " max(Delta alpha, Delta beta)= ", epsilon)
       end
       if epsilon < tol
          break
       end
    end
    return mle_new
  end
  function make_beta(x::AbstractArray{T,1}, w::AbstractArray{T,1};
              xmean::T = Statistics.mean(x,AnalyticWeights(w)),
              xvar::T = Statistics.varm(x,xmean,AnalyticWeights(w)), tau_factor::T=one(T) ) where T<:Real
       alpha_beta = (xvar < xmean*(1.0-xmean)) ? [mom_alpha(xmean, xvar), mom_beta(xmean,xvar)] : ones(T,2)
       if any(alpha_beta.<=zero(T))
          alpha_beta = fit_beta_mle(x, w=w, init_guess=alpha_beta, verbose=true)
       end
       if any(alpha_beta.<=zero(T))
          alpha_beta = ones(T,2)
       else
          if minimum(alpha_beta)>1.5*tau_factor && sum(alpha_beta)>=20.0*tau_factor
             alpha_beta ./= tau_factor
          end
       end
       Beta(alpha_beta[1], alpha_beta[2])
  end
  function make_beta_transformed(x::AbstractArray{T,1}, w::AbstractArray{T,1}; xmin::T=zero(T), xmax::T=one(T),
              xmean::T = Statistics.mean(x,AnalyticWeights(w)),
              xvar::T = Statistics.varm(x,xmean,AnalyticWeights(w)), tau_factor::T=one(T) ) where T<:Real
       alpha_beta = (xvar < xmean*(1.0-xmean)) ? [mom_alpha(xmean, xvar), mom_beta(xmean,xvar)] : ones(T,2)
       if any(alpha_beta.<=zero(T))
          alpha_beta = fit_beta_mle(x, w=w, init_guess=alpha_beta, verbose=true)
       end
       if any(alpha_beta.<=zero(T))
          alpha_beta = ones(T,2)
       else
          if minimum(alpha_beta)>1.5*tau_factor && sum(alpha_beta)>=20.0*tau_factor
             alpha_beta ./= tau_factor
          end
       end
       LinearTransformedBeta(alpha_beta[1], alpha_beta[2], xmin=xmin, xmax=xmax)
  end

  theta_mean =  sum(theta.*weights',dims=2) # weighted mean for parameters
  theta_var = ABC.var_weighted(theta'.-theta_mean',weights)  # scaled, weighted covar for parameters
  tau_factor_indiv = fill(tau_factor,length(theta_var))
  if verbose
     println("total: ",theta_mean[1]," ",theta_var[1])
  end
  for i in 2:size(theta,1)
    mean_ratio = sum(theta[1,:].*theta[i,:].*weights) /(theta_mean[1]*theta_mean[i]) # weighted mean for parameters
    var_ratio = var_weighted(vec(theta[1,:].*theta[i,:]).-(theta_mean[1]*theta_mean[i]),weights)/(2 * theta_mean[1] * theta_var[i]) # scaled, weighted covar for parameters
    if verbose
       println("i=",i,": ",theta_mean[i]," ",theta_var[i]," ratios: ",mean_ratio, " ",var_ratio)
    end
    var_ratio  = var_ratio  >= one(var_ratio)  ? var_ratio  : one(var_ratio)
    tau_factor_indiv[i] = tau_factor*var_ratio
  end
  if verbose
     flush(Base.stdout)
  end
  #=
  println("mean= ",theta_mean)
  println("var= ",theta_var)
  for i in 1:length(theta_mean)
     println("a= ",alpha(theta_mean[i],tau_factor*theta_var[i]), "  b= ",beta(theta_mean[i],tau_factor*theta_var[i]))
  end
  =#

  dist = ApproximateBayesianComputing.CompositeDistributions.CompositeDist( vcat(
         make_beta_transformed(theta[1,:]./max_rate, weights, xmin=0.0, xmax=max_rate, xmean=theta_mean[1]/max_rate, xvar=theta_var[1]/max_rate^2, tau_factor=tau_factor_indiv[1]),
         ContinuousDistribution[ make_beta(theta[i,:], weights, xmean=theta_mean[i], xvar=theta_var[i], tau_factor=tau_factor_indiv[i]) for i in 2:size(theta,1) ] ))

end

function make_proposal_dist_multidim_beta(pop::abc_population_type, tau_factor::Float64; verbose::Bool = false)
	make_proposal_dist_multidim_beta(pop.theta, pop.weights, tau_factor, verbose=verbose)
end

# Tell ABC what it needs to know for a simulation
abc_plan = abc_pmc_plan_type(gen_data,calc_mean_per_bin,calc_dist_l1, param_prior;
                is_valid=is_valid,normalize=normalize_theta!,
                make_proposal_dist=make_proposal_dist_multidim_beta,
                epsilon_init=0.6, epsilon_reduction_factor=0.75001,
                tau_factor=1.2,target_epsilon=0.00001*nbins,
                num_part=100,num_max_attempt=100
                );

# Generate "true/observed data" and summary statistics
data_true = abc_plan.gen_data(theta_true)
ss_true = abc_plan.calc_summary_stats(data_true)
println("theta_true = ",theta_true," ss_true = ",ss_true, " d= 0")
dist_draws_from_true = [abc_plan.calc_dist(ss_true,abc_plan.calc_summary_stats(abc_plan.gen_data(theta_true))) for i in 1:100]
println("dist from true param: mean = ", mean(dist_draws_from_true), " std = ", std(dist_draws_from_true))

abc_plan.target_epsilon = mean(dist_draws_from_true)-1*std(dist_draws_from_true)

# Run ABC simulation
@time pop_out = run_abc(abc_plan,ss_true;verbose=true);


text:
abstract type AbstractBregmanPotential end

# Usually, Tv is Float64
struct EuclideanBregmanPotential{Tweights<:AbstractArray}<:AbstractBregmanPotential
    # the potential is given by x ↦ x' diagm(weights) x / 2
    weights::Tweights
end

function EuclideanBregmanPotential(weights)
    return EuclideanBregmanPotential{typeof(weights)}(weights)
end

function evaluate_potential(potential::EuclideanBregmanPotential, x)
    return #value of potential
end

function evaluate_gradient(potential::EuclideanBregmanPotential, x)
    return #handwritten gradient
end

function evaluate_gradient(potential::AbstractBregmanPotential, x)
    return #using autograd and evaluate_potential
end


text:
@testset "974.subarray-sums-divisible-by-k.jl" begin
    @test subarrays_div_by_k([4,5,0,-2,-3,1], 5) == 7
end


text:
using GtkTextUtils, Gtk
import GtkTextUtils: offset

@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

import Gtk: GtkTextIter, GLib.mutable, hasselection

w = GtkWindow()
b = GtkTextBuffer()
b.text[String] = "test"
v = GtkTextView(b)

push!(w,v)
showall(w)

begin
    it = GtkTextIter(b)
    get_gtk_property(it,:line,Int64)
end

word, its, ite = get_current_line_text(b)
@test word == "test"
@test offset(its) == 0
@test offset(ite) == 4

word,its,ite = select_word(GtkTextIter(b),b,false)
str="test2"

function replace_text2(b::GtkTextBuffer,its,ite,str::AbstractString)
    pos = offset(its)+1
    splice!(b,its:ite)
    insert!(b,GtkTextIter(b,pos),str)
end

replace_text2(b,its,ite,"test2")

word, its, ite = get_current_line_text(b)
@test word == "test2"

destroy(w)


text:
### Computing fields and other physical quantities ###
using LinearAlgebra
import LinearAlgebra: transpose!
import CartesianGrids: convective_derivative!

"""
    vorticity(sol::ODESolution,sys::NavierStokes,t)

Return the vorticity field associated with solution vector `sol` on the grid in `sys`,
at time(s) `t`.
"""
@inline vorticity(w::Nodes{Dual,NX,NY},sys::NavierStokes{NX,NY},t::Real) where {NX,NY} =
            vorticity!(sys.velocity_cache.Wn,w,sys,t)

function vorticity!(w::Nodes{Dual,NX,NY},w_in::Nodes{Dual,NX,NY},sys::NavierStokes{NX,NY},t::Real) where {NX,NY}
  w .= 0.0
  _vorticity!(w,w_in,sys,t)
  _vorticity_sheet!(w,sys,t)
  return w
end

@inline _vorticity!(w::Nodes{Dual,NX,NY},win::Nodes{Dual,NX,NY},
                    sys::NavierStokes{NX,NY},t::Real) where {NX,NY} = w .+= win/cellsize(sys)

_vorticity_sheet!(w::Nodes{Dual,NX,NY},sys::NavierStokes{NX,NY,0},t::Real) where {NX,NY} = w

_vorticity_sheet!(w::Nodes{Dual,NX,NY},sys::NavierStokes{NX,NY,N,MT,FS,
                                              ExternalInternalFlow},t::Real) where {NX,NY,N,MT,FS} = w

function _vorticity_sheet!(w::Nodes{Dual,NX,NY},sys::NavierStokes{NX,NY,N,MT,FS,SD},t::Real) where {NX,NY,N,MT,FS,SD}
  @unpack velocity_cache = sys
  @unpack Sb, Δus = velocity_cache
  #if isnothing(sys.Rn)
  regop = _regularization(sys)
  Rn = RegularizationMatrix(regop,Sb,w)
  #else
  #  Rn = sys.Rn
  #end
  surface_velocity_jump!(Δus,sys,t)
  cross!(Sb,Δus,normals(sys))
  w .+= Rn*Sb
  return w
end

function velocity!(u::Edges{Primal,NX,NY},w::Nodes{Dual,NX,NY},sys::NavierStokes{NX,NY,N},t::Real) where {NX,NY,N}
    u .= 0.0
    _velocity_vorticity!(u,w,sys)
    _velocity_single_layer!(u,sys,t)
    _velocity_freestream!(u,sys,t)
    return u
end

function _velocity_vorticity!(u::Edges{Primal,NX,NY},w::Nodes{Dual,NX,NY},sys::NavierStokes{NX,NY,N}) where {NX,NY,N}
  @unpack velocity_cache = sys
  @unpack Wn, Vf = velocity_cache
  Wn .= 0.0
  _unscaled_streamfunction_vorticity!(Wn,w,sys)
  Vf .= 0.0
  curl!(Vf,Wn)
  u .+= Vf
  return u
end

@inline _velocity_single_layer!(u::Edges{Primal,NX,NY},sys::NavierStokes{NX,NY,N,MT,FS,
                                      ExternalInternalFlow},t::Real) where {NX,NY,N,MT,FS} = u

@inline _velocity_single_layer!(u::Edges{Primal,NX,NY},sys::NavierStokes{NX,NY,0},t::Real) where {NX,NY} = u

function _velocity_single_layer!(u::Edges{Primal,NX,NY},sys::NavierStokes{NX,NY,N,MT,FS,SD},t::Real) where {NX,NY,N,MT,FS,SD}
    @unpack velocity_cache = sys
    @unpack Vf, Sc, Sb, Δus = velocity_cache
    surface_velocity_jump!(Δus,sys,t)
    pointwise_dot!(Sb,Δus,normals(sys))
    Sb .*= cellsize(sys)
    sys.slc(Sc,Sb)
    _unscaled_scalarpotential!(Sc,Sc,sys)
    fill!(Vf,0.0)
    grad!(Vf,Sc)
    u .+= Vf
    return u
end

function _velocity_freestream!(u::Edges{Primal,NX,NY},sys::NavierStokes{NX,NY,N,MT,FS,SD},t::Real) where {NX,NY,N,MT,FS,SD}
  U∞, V∞ = freestream(t,sys)
  u.u .+= U∞
  u.v .+= V∞
  return u
end

"""
    velocity(sol::ODESolution,sys::NavierStokes,t)

Return the velocity field associated with solution vector `sol` on the grid in `sys`,
at time(s) `t`.
"""
velocity(w::Nodes{Dual,NX,NY},a...) where {NX,NY} = velocity!(Edges(Primal,(NX,NY)),w,a...)

function streamfunction!(ψ::Nodes{Dual,NX,NY},w::Nodes{Dual,NX,NY},sys::NavierStokes{NX,NY},t::Real) where {NX,NY}
  ψ .= 0.0
  _unscaled_streamfunction_vorticity!(ψ,w,sys)
  ψ .*= cellsize(sys)
  _streamfunction_freestream!(ψ,sys,t)
  return ψ
end

function _unscaled_streamfunction_vorticity!(ψ::Nodes{Dual,NX,NY},w::Nodes{Dual,NX,NY},sys::NavierStokes{NX,NY}) where {NX,NY}
  @unpack streamfunction_cache, L = sys
  @unpack Sn = streamfunction_cache
  Sn .= 0.0
  ldiv!(Sn,L,w)
  ψ .-= Sn
  return ψ
end

function _streamfunction_freestream!(ψ::Nodes{Dual,NX,NY},sys::NavierStokes{NX,NY},t::Real) where {NX,NY}
  xg, yg = coordinates(ψ,sys.grid)
  U∞, V∞ = freestream(t,sys)
  ψ .+= U∞*yg' .- V∞*xg
  return ψ
end

"""
    streamfunction(sol::ODESolution,sys::NavierStokes,t)

Return the streamfunction field associated with solution vector `sol` on the grid in `sys`,
at time(s) `t`.
"""
streamfunction(w::Nodes{Dual,NX,NY},sys::NavierStokes{NX,NY},t::Real) where {NX,NY} = streamfunction!(Nodes(Dual,(NX,NY)),w,sys,t)

function scalarpotential!(ϕ::Nodes{Primal,NX,NY},dil::Nodes{Primal,NX,NY},sys::NavierStokes{NX,NY},t::Real) where {NX,NY}
  _unscaled_scalarpotential!(ϕ,dil,sys)
  # This assumes that dil is the physical dilation times the grid spacing
  ϕ .*= cellsize(sys)
  return ϕ
end

function _unscaled_scalarpotential!(ϕ::Nodes{Primal,NX,NY},dil::Nodes{Primal,NX,NY},sys::NavierStokes{NX,NY}) where {NX,NY}
  ldiv!(ϕ,sys.L,dil)
  return ϕ
end

"""
    scalarpotential(sol::ODESolution,sys::NavierStokes,t)

Return the scalar potential field associated with solution vector `sol` on the grid in `sys`,
at time(s) `t`.
"""
scalarpotential(dil::Nodes{Primal,NX,NY},sys::NavierStokes{NX,NY},t::Real) where {NX,NY} = scalarpotential!(Nodes(Primal,(NX,NY)),dil,sys,t)

function convective_derivative!(out::Edges{Primal,NX,NY},u::Edges{Primal,NX,NY},
                        sys::NavierStokes{NX,NY},t::Real) where {NX,NY}
  out .= u
  _unscaled_convective_derivative!(out,sys)
  Δx⁻¹ = 1/cellsize(sys)
  out .*= Δx⁻¹
  return out
end

"""
    convective_derivative(sol::ODESolution,sys::NavierStokes,t)

Return the convective derivative associated with solution vector `sol` on the grid in `sys`,
at time(s) `t`.
"""
convective_derivative(u::Edges{Primal,NX,NY},sys::NavierStokes{NX,NY},t::Real) where {NX,NY} =
      convective_derivative!(Edges(Primal,(NX,NY)),u,sys,t)

"""
    pressure(sol::ODESolution,sys::NavierStokes,t)

Return the pressure field associated with solution vector `sol` on the grid in `sys`,
at time(s) `t`.
"""
function pressure(u::ConstrainedSystems.ArrayPartition,sys::NavierStokes{NX,NY},t) where {NX,NY}
    @unpack pressure_cache = sys
    @unpack dvp, dvf, lp = pressure_cache
    w, τ = state(u), constraint(u)
    ns_rhs_velocity!(dvp,w,sys,t)
    fill!(dvf,0.0)
    _vel_ns_op_constraint_force!(dvf,τ,sys)
    dvp .-= dvf
    fill!(lp,0.0)
    divergence!(lp,dvp)
    press = zero(lp)
    ldiv!(press,sys.L,lp)
    press .*= cellsize(sys)
    return press
end

"""
    vorticity(integrator)

Return the vorticity field associated with `integrator` at its current state.
""" function vorticity end

"""
    velocity(integrator)

Return the velocity field associated with `integrator` at its current state.
""" function velocity end

"""
    streamfunction(integrator)

Return the streamfunction field associated with `integrator` at its current state.
""" function streamfunction end

"""
    scalarpotential(integrator)

Return the scalar potential field associated with `integrator` at its current state.
""" function scalarpotential end

"""
    convective_derivative(integrator)

Return the convective derivative associated with `integrator` at its current state.
""" function convective_derivative end

"""
    pressure(integrator)

Return the pressure field associated with `integrator` at its current state.
""" function pressure end

for fcn in (:vorticity,:velocity,:streamfunction,:scalarpotential,:convective_derivative)

  @eval $fcn(s::ConstrainedSystems.ArrayPartition,sys::NavierStokes,t) = $fcn(state(s),sys,t)

end

for fcn in (:vorticity,:velocity,:streamfunction,:scalarpotential,:convective_derivative,:pressure)
  @eval $fcn(integ::ConstrainedSystems.OrdinaryDiffEq.ODEIntegrator) = $fcn(integ.u,integ.p,integ.t)

  @eval $fcn(sol::ConstrainedSystems.OrdinaryDiffEq.ODESolution,sys::NavierStokes,t) = $fcn(sol(t),sys,t)

  @eval $fcn(sol::ConstrainedSystems.OrdinaryDiffEq.ODESolution,sys::NavierStokes,t::AbstractArray) =
      map(ti -> $fcn(sol(ti),sys,ti),t)

end

# Surface field quantities

"""
    traction(integrator)

Return the traction at all of the surface points associated with `integrator` at
its current state.
""" function traction end

"""
    pressurejump(integrator)

Return the pressure jump at all of the surface points associated with `integrator` at
its current state.
""" function pressurejump end

for fcn in (:traction,:pressurejump)
  @eval $fcn(s::ConstrainedSystems.ArrayPartition,a...;kwargs...) = $fcn(constraint(s),a...;kwargs...)
  @eval $fcn(integ::ConstrainedSystems.OrdinaryDiffEq.ODEIntegrator) = $fcn(integ.u,integ.p,integ.t)
  @eval $fcn(sol::ODESolution,sys::NavierStokes,t) = $fcn(sol(t),sys,t)

  @eval $fcn(sol::ODESolution,sys::NavierStokes,t::AbstractArray) =
        map(ti -> $fcn(sol(ti),sys,ti),t)

end

"""
    traction(sol::ODESolution,sys::NavierStokes,t)

Return the traction at all of the surface points associated with solution `sol` on
grid in `sys` at time(s) t.
"""
traction(τ,sys::NavierStokes{NX,NY,0,MT,FS,SD},t) where {NX,NY,MT,FS,SD} = τ

traction(τ::VectorData{N},sys::NavierStokes{NX,NY,N,MT,FS,ExternalInternalFlow},t) where {NX,NY,N,MT,FS} = τ

function traction(τ::VectorData{N},sys::NavierStokes{NX,NY,N,MT,FS,SD},t) where {NX,NY,N,MT,FS,SD}
    @unpack surfacetraction_cache, bodies = sys
    @unpack Vb, Δus, Sb = surfacetraction_cache
    nrm = normals(bodies)
    relative_surface_velocity!(Δus,sys,t)
    pointwise_dot!(Sb,nrm,Δus)
    surface_velocity_jump!(Δus,sys,t)
    product!(Vb,Δus,Sb)
    return τ + Vb
end

"""
    pressurejump(sol::ODESolution,sys::NavierStokes,t)

Return the pressure jump at all of the surface points associated with solution `sol` on
grid in `sys` at time(s) t.
"""
function pressurejump(τ::VectorData{N},sys::NavierStokes{NX,NY,N},t) where {NX,NY,N}
    @unpack surfacetraction_cache, bodies = sys
    @unpack Sb = surfacetraction_cache

    #_hasfilter(sys) ? (τf = similar(τ); τf .= sys.Cf*force(τ,sys)) : τf = deepcopy(force(τ,sys))
    #τf = deepcopy(τ)
    #return nrm.u∘τf.u + nrm.v∘τf.v # This might need some modification

    press = zero(Sb)
    nrm = normals(bodies)
    pointwise_dot!(press,nrm,traction(τ,sys,t))
    return -press
end

## Total quantities

force(τ,sys::NavierStokes{NX,NY,0},t,bodyi::Int) where {NX,NY} = Vector{Float64}(), Vector{Float64}()

function force(τ::VectorData{N},sys::NavierStokes,t,bodyi::Int) where {N}
    @unpack surfacetraction_cache, bodies = sys
    @unpack Vb = surfacetraction_cache
    product!(Vb,traction(τ,sys,t),areas(bodies))
    fxvec = sum(Vb.u,bodies,bodyi)
    fyvec = sum(Vb.v,bodies,bodyi)
    return fxvec, fyvec
end

force(s::ConstrainedSystems.ArrayPartition,sys,t,bodyi) = force(constraint(s),sys,t,bodyi)

"""
    force(integ,bodyindex)

Given the state of the system in integrator `integ`, return the current force
on the body with index `bodyindex`.
"""
force(integ::ConstrainedSystems.OrdinaryDiffEq.ODEIntegrator,bodyi) = force(integ.u,integ.p,integ.t,bodyi)

"""
    force(sol,sys::NavierStokes,bodyindex) -> Tuple{Vector,Vector}

Given the solution history vector `sol` and the system `sys`, return the force
history on the body with index `bodyindex` as a tuple of vectors, one for
each component.
"""
function force(sol::ConstrainedSystems.OrdinaryDiffEq.ODESolution,sys::NavierStokes,bodyi::Int)
    fx = map((u,t) -> force(u,sys,t,bodyi)[1],sol.u,sol.t)
    fy = map((u,t) -> force(u,sys,t,bodyi)[2],sol.u,sol.t)
    fx, fy
end

function moment(τ::VectorData{N},bodies::BodyList,sys::NavierStokes,t,bodyi::Int;center=(0.0,0.0)) where {N}
    @unpack surfacetraction_cache = sys
    @unpack Vb = surfacetraction_cache
    xc, yc = center
    x, y = collect(bodies)
    dX = VectorData(x .- xc,  y .- yc)
    product!(Vb,traction(τ,sys,t),areas(bodies))
    return sum(cross(dX,Vb),bodies,bodyi)
end

function moment(s::ConstrainedSystems.ArrayPartition,sys::NavierStokes{NX,NY,N,MovingPoints},t,bodyi;center=(0.0,0.0)) where {NX,NY,N}
    bodies = deepcopy(sys.bodies)
    tl! = RigidTransformList(aux_state(s))
    tl!(bodies)
    moment(constraint(s),bodies,sys,t,bodyi;center=center)
end

moment(s::ConstrainedSystems.ArrayPartition,sys::NavierStokes{NX,NY,N,StaticPoints},t,bodyi;center=(0.0,0.0)) where {NX,NY,N} =
        moment(constraint(s),sys.bodies,sys,t,bodyi;center=center)

"""
    moment(integ,bodyindex[,center=(0,0))

Given the state of the system in integrator `integ`, return the current moment
on the body with index `bodyindex`.  The center of the
moment can be passed as an optional tuple argument.
"""
moment(integ::ConstrainedSystems.OrdinaryDiffEq.ODEIntegrator,bodyi;center=(0.0,0.0)) =
      moment(integ.u,integ.p,integ.t,bodyi,center=center)

"""
    moment(sol,sys::NavierStokes,bodyindex[,center=(0,0)]) -> Vector

Given the solution history vector `sol` and the system `sys`, return the moment
history on the body with index `bodyindex` as a vector. The center of the
moment can be passed as an optional tuple argument.
"""
function moment(sol::ConstrainedSystems.OrdinaryDiffEq.ODESolution,sys::NavierStokes,bodyi::Int;center=(0.0,0.0))
    mom = map((u,t) -> moment(u,sys,t,bodyi,center=center),sol.u,sol.t)
end


text:
mutable struct StateSpace
    A::Matrix{Float64}
    B::Matrix{Float64}
    C::Matrix{Float64}
    D::Matrix{Float64}
    ns::Integer
    nu::Integer
    ny::Integer

    function StateSpace(A::Array, B::Array, C::Array, D::Array=Array{Float64,2}(undef,0,2))::StateSpace
        this = new();

        this.ns = size(A)[1];
        this.nu = size(B)[2];
        this.ny = size(C)[1];

        # Check dimensions
        if (size(A)[1] != size(A)[2]) error("Matrix A must be square."); end
        if (size(B)[1] != size(A)[1]) error("Matrices A and B must have same number of rows."); end
        if (size(C)[2] != size(A)[2]) error("Matrices A and C must have same number of cols."); end

        this.A = float(A);
        this.B = float(B);
        this.C = float(C);

        if isempty(D)
            this.D = zeros(this.ny,this.nu);
        else
            if (size(D)[1] != size(C)[1]) error("Matrices C and D must have same number of rows."); end
            if (size(D)[2] != size(B)[2]) error("Matrices B and D must have same number of cols."); end
            this.D = float(D);
        end
        return(this);
    end
end

function CheckIODimensions(sys::StateSpace,no::Integer,ni::Integer)
	if (ni<1 || ni>sys.nu) error("Invalid input specified.") end
    if (no<1 || no>sys.ny) error("Invalid output specified.") end
end

function Poles(S::StateSpace)
    lam = eigvals(S.A);
    return lam;
end

function NaturalFrequencies(S::StateSpace)::Vector{Float64}
    lam = Poles(S);
    ri = real.(lam);
    ii = imag.(lam);
    om = sqrt.(ri.^2 + ii.^2);
    return(om);
end

function Damping(S::StateSpace)::Vector{Float64}
    lam = Poles(S);
    ri = real.(lam);
    ii = imag.(lam);
    om = sqrt.(ri.^2 + ii.^2);
    d = abs.(ri./om);
    return(d);
end

function rifd(S::StateSpace,printFlag::Bool=false)
    lam = Poles(S);
    rr = Vector{Float64}(undef,length(lam));
    ii = Vector{Float64}(undef,length(lam));
    om = Vector{Float64}(undef,length(lam));
    d = Vector{Float64}(undef,length(lam));
    if printFlag
        @printf("   Real \t   Imag \t  Omega \t  Damping\n")
    end
    for (i,l) in enumerate(lam)
        rr[i] = real(l)
        ii[i] = imag(l)
        om[i] = sqrt(rr[i]^2+ii[i]^2)
        d[i]= abs(rr[i]/om[i])
        if printFlag
            @printf("%+.6f \t %+.6f \t %.6f \t %.6f\n",rr[i],ii[i],om[i],d[i])
        end
    end
    return(rr,ii,om,d)
end

function TransmissionZeros(S::StateSpace,no::Integer,ni::Integer)::Vector{Float64}

    CheckIODimensions(S,no,ni);

    B = S.B[:,ni];
    C = reshape(S.C[no,:],1,length(S.C[no,:])); # Must ensure row vector
    D = S.D[no,ni];
    L = [S.A B; C D];
    II = Matrix{Float64}(I, size(S.A)[1], size(S.A)[2]);
    M = [II 0*B;0*C 0*D];
    e,v = eigen(L,M);
    z = filter(x-> x!=Inf,e)
    return z;
end

function DCGain(S::StateSpace)::Matrix{Float64}
    lam = Poles(S);
    PolesAtOrigin = filter(x->x==eps(Float64),abs.(lam));
    if isempty(PolesAtOrigin)
        return(S.D - S.C*inv(S.A)*S.B);
    else
        return(ones(S.ny,S.nu)*Inf);
    end
end

function ControllabilityMatrix(S::StateSpace)::Matrix{Float64}
    R = Matrix{Float64}(undef,S.ns,S.ns*S.nu);
    for i=0:(S.ns-1)
        i1 = i*S.nu+1;
        i2 = (i+1)*S.nu;
        R[:,i1:i2] = S.A^i*S.B;
    end
    return(R);
end

function ObservabilityMatrix(S::StateSpace)::Matrix{Float64}
    R = Matrix{Float64}(undef,S.ns*S.ny,S.ns);
    for i=0:(S.ns-1)
        i1 = i*S.ny+1;
        i2 = (i+1)*S.ny;
        R[i1:i2,:] = S.C*S.A^i;
    end
    return(R);
end

function IsControllable(S::StateSpace)::Bool
    if (rank(ControllabilityMatrix(S))==S.ns)
        return true;
    else
        return false;
    end
end

function IsObservable(S::StateSpace)::Bool
    if (rank(ObservabilityMatrix(S))==S.ns)
        return true;
    else
        return false;
    end
end

max_stars_repo_name: isrlab/RobustControl

max_issues_repo_name: isrlab/RobustControl

max_forks_repo_name: isrlab/RobustControl


text:
function init_model(swap_axes; data = nothing , model_params = nothing)
    if random_seed != nothing
        @eval @everywhere Random.seed!($random_seed)
    end
    global lm
    if data == nothing
        points_dict = load_data(data_path, groups_count, prefix = data_prefix)
        if swap_axes != nothing
            points_dict = axes_swapper(points_dict,swap_axes)
        end
        preprocessing!(points_dict,local_dim,global_preprocessing,local_preprocessing)
    else
        points_dict = data
    end

    if model_params == nothing
        model_hyperparams = model_hyper_params(global_hyper_params,local_hyper_params,α,γ,η,global_weight,local_weight,total_dim,local_dim)
    else
        model_hyperparams = model_params
    end
    blocal_hyper_params = model_hyperparams.local_hyper_params
    groups_dict = Dict()
    for (k,v) in points_dict
        labels = rand(1:initial_local_clusters,(size(v,2),1))
        labels_subcluster = rand(1:4,(size(v,2),1))
        weights = ones(initial_local_clusters) * (1/initial_local_clusters)
        local_clusters = local_cluster[]
        if isa(blocal_hyper_params,Array)
            localised_params = model_hyper_params(model_hyperparams.global_hyper_params
            ,blocal_hyper_params[k],model_hyperparams.α,
            model_hyperparams.γ,model_hyperparams.η,model_hyperparams.global_weight,
            model_hyperparams.local_weight,model_hyperparams.total_dim,model_hyperparams.local_dim)
        else
            localised_params = model_hyperparams
        end
        groups_dict[k] = local_group(localised_params,v,labels,labels_subcluster,local_clusters,Float64[],k)
    end
    if isa(model_hyperparams.global_hyper_params,topic_modeling_hyper)
        global is_tp = true
    else
        global is_tp = false
    end

    @eval global groups_dict = $groups_dict
    if mp
        num_of_workers = nworkers()
        for w in workers()
            @spawnat w global groups_dict = Dict()
        end

        @sync for (index,group) in groups_dict
            @spawnat ((index % num_of_workers)+2) set_group(index,group)
        end
    end
    return hdp_shared_features(model_hyperparams,groups_dict,global_cluster[],Float64[])
end

function init_first_clusters!(hdp_model::hdp_shared_features)
    for (k,v) in hdp_model.groups_dict
        v.local_clusters = []
        for i=1:initial_local_clusters
            push!(v.local_clusters, create_first_local_cluster(v, initial_global_clusters))
        end
    end
    global_c = []
    for i=1:initial_global_clusters
        push!(global_c,create_first_global_cluster(hdp_model.model_hyperparams, hdp_model.groups_dict, i))
    end
    global global_clusters_vector  = global_c
    for w in workers()
        @eval @spawnat $w global global_clusters_vector = $global_c
    end
    # @eval @everywhere global global_clusters_vector = $global_c
end

function hdp_shared_features(model_params, swap_axes = nothing;multiprocess = false)
    cur_dir = pwd()
    include(model_params)
    cd(cur_dir)
    hdp_model = init_model(swap_axes)
    init_first_clusters!(hdp_model)
    groups_stats = Dict()
    global num_of_workers = nworkers()
    if @isdefined split_delays
        @everywhere global split_delays = split_delays
    else
        @everywhere global split_delays = true
    end
    @everywhere global hard_clustering = hard_clustering
    global posterior_history = []
    global word_ll_history = []
    global topic_count = []
    global mp = multiprocess
    for i=1:iterations
        println("Iteration: " * string(i))
        println("Global Counts: " * string([x.clusters_count for x in global_clusters_vector]))
        final = false
        no_more_splits = false
        if i > iterations - argmax_sample_stop #We assume the cluters k has been setteled by now, and a low probability random split can do dmg
            final = true
        end
        if i >= iterations - split_stop
            no_more_splits = true
        end
        model_iteration(hdp_model,final,no_more_splits)

    end
    hdp_model.global_clusters = global_clusters_vector
    return hdp_model, posterior_history, word_ll_history, topic_count
end

function model_iteration(hdp_model,final, no_more_splits,burnout = 5)
    groups_stats = Vector{local_group_stats}(undef,length(groups_dict))
    @everywhere global burnout_period = 5
    sample_global_clusters_params!(hdp_model)
    global global_clusters_vector = global_clusters_vector
    refs= Dict()
    if mp
        for w in workers()
            refs[w] = remotecall(set_global_clusters_vector, w, global_clusters_vector)
        end
        for w in workers()
            fetch(refs[w])
        end
    end
    begin
        if mp
            @sync for (index,group) in hdp_model.groups_dict
                lc = group.local_clusters
                groups_stats[index] = @spawnat ((index % num_of_workers)+2) group_step(index,lc, final)
            end
        else
            Threads.@threads for index=1:length(groups_dict)
                group=groups_dict[index]
                lc = group.local_clusters
                groups_stats[index] = group_step(index,lc, final)
            end
        end
        for index=1:length(groups_dict)
            update_group_from_stats!(hdp_model.groups_dict[index], fetch(groups_stats[index]))
        end
    end
    sample_clusters_labels!(hdp_model, (hard_clustering ? true : final))
    remove_empty_clusters!(hdp_model)
    update_suff_stats_posterior!(hdp_model,collect(1:length(global_clusters_vector)))
    hdp_model.global_clusters = global_clusters_vector
    push!(posterior_history,calc_global_posterior(hdp_model))
    # if isa(hdp_model.model_hyperparams.global_hyper_params, topic_modeling_hyper)
    #     word_ll = calc_avg_word(hdp_model)
    #     println("Per Word LL:" * string(word_ll))
    #     push!(word_ll_history,word_ll)
    #     push!(topic_count,length(global_clusters_vector))
    # end
    if no_more_splits == false
        # println(length((global_clusters_vector)))
        indices = check_and_split!(hdp_model, final)
        i = 1
        while i < length(indices)
            for (index,group) in hdp_model.groups_dict
                split_global_cluster!(group,indices[i],indices[i+1])
            end
            i+= 2
        end
        # println(length((global_clusters_vector)))
        indices = check_and_merge!(hdp_model, final)
        if length(indices) > 0
            for (index,group) in hdp_model.groups_dict
                merge_global_cluster!(group,indices[1],indices[2])
            end
        end
        remove_empty_clusters!(hdp_model)
        if length(indices) > 0
            println("merged:" * string(indices))
        end
    end
end

function create_default_priors(gdim,ldim,prior_type::Symbol)
    if prior_type == :niw
        g_prior = niw_hyperparams(1.0,
            zeros(gdim),
            gdim+3,
            Matrix{Float64}(I, gdim, gdim)*1)
        l_prior = niw_hyperparams(1.0,
            zeros(ldim),
            ldim+3,
            Matrix{Float64}(I, ldim, ldim)*1)
    else
        g_prior = multinomial_hyper(ones(gdim)*500.0)
        l_prior = multinomial_hyper(ones(ldim)*500.0)
    end
    return g_prior, l_prior
end

@everywhere function swap_axes_worker(swap_vec)
    for (k,v) in groups_dict
        v.points = v.points[swap_vec,:]
    end
end

function create_swap_vec(dim,glob_mapping, index)
    swap_vec = zeros(dim)
    reverse_swap_vec = zeros(dim)
    front_index = 1
    back_index = dim
    for (k,v) in enumerate(glob_mapping)
        if k == index || v == 1
            swap_vec[front_index] = k
            reverse_swap_vec[k] = front_index
            front_index+=1
        else
            swap_vec[back_index] = k
            reverse_swap_vec[k] = back_index
            back_index -= 1
        end
    end
    return Int.(swap_vec), Int.(reverse_swap_vec)
end

function calc_global_posterior(hdp_model::hdp_shared_features, ismnm = false)
    pts_count = 0.0
    log_posterior = log(hdp_model.model_hyperparams.γ)
    for (k,group) in hdp_model.groups_dict
        pts_count += size(group.points,2)
        # if ismnm
            log_posterior+= calc_group_posterior(group)
        # end
    end
    log_posterior-= logabsgamma(pts_count)[1]
    for cluster in hdp_model.global_clusters
        if cluster.cluster_params.cluster_params.suff_statistics.N == 0
            continue
        end
        #posterior_param = update_posterior_evidence(model.hyper, model.clusters_params.sufficient_stats, index)
        log_posterior += log_marginal_likelihood(cluster.cluster_params.cluster_params.hyperparams,
            cluster.cluster_params.cluster_params.posterior_hyperparams,
            cluster.cluster_params.cluster_params.suff_statistics)
        log_posterior += log(hdp_model.model_hyperparams.γ) + logabsgamma(cluster.cluster_params.cluster_params.suff_statistics.N)[1]
        # println(cluster.cluster_params.cluster_params.suff_statistics)
    end
    return log_posterior
end

function calc_avg_word(hdp_model::hdp_shared_features)
    # return 0
    global_preds = get_model_global_pred(hdp_model)
    group_mixtures = Dict([k => [x/ sum(counts(v,length(hdp_model.global_clusters))) for x in counts(v,length(hdp_model.global_clusters))] for (k,v) in global_preds])
    # word_count = Dict([k => counts(Int.(v.points[:]), length(hdp_model.model_hyperparams.global_hyper_params.α)) for (k,v) in hdp_model.groups_dict])
    cluster_dists = [x.cluster_params.cluster_params.distribution.α for x in hdp_model.global_clusters]
    total_likelihood = 0.0
    total_points = 0
    for (k,v) in group_mixtures
        parr = zeros(length(hdp_model.model_hyperparams.global_hyper_params.α),length(v))
        for (index,part) in enumerate(v)
            parr[:,index] = cluster_dists[index] .+ log(part)
        end
        parr = exp.(parr)
        parr[isnan.(parr)] .= 0
        wordp = sum(parr,dims = 2)
        wordp = log.(wordp)
        wordp[isnan.(wordp)] .= 0
        rel_pts = hdp_model.groups_dict[k].points
        word_counts = counts(Int.(rel_pts), length(hdp_model.model_hyperparams.global_hyper_params.α))
        cluster_ll = wordp .* word_counts
        cluster_ll[isnan.(cluster_ll)] .= 0
        total_points += sum(word_counts)
        # println(any(isnan.(cluster_ll)))
        total_likelihood += sum(cluster_ll)
    end
    # println("ll: " * string(total_likelihood) * "   pts: " * string(total_points))
    return total_likelihood / total_points
end

function calc_group_posterior(group::local_group)
    log_posterior = log(group.model_hyperparams.α) - logabsgamma(size(group.points,2))[1]
    for cluster in group.local_clusters
        if cluster.cluster_params.cluster_params.suff_statistics.N == 0
            continue
        end
        log_posterior += log_marginal_likelihood(cluster.cluster_params.cluster_params.hyperparams,
            cluster.cluster_params.cluster_params.posterior_hyperparams,
            cluster.cluster_params.cluster_params.suff_statistics)
    end
    return log_posterior
end

function k_mean_likelihood(likehood_rating,k)
    min_likelihood = minimum(likehood_rating)
    max_likelihood = maximum(likehood_rating)
    k_means = zeros(k)
    k_means[1] = min_likelihood
    k_means[k] = 0
    k_interval = abs(max_likelihood - min_likelihood) / k
    for i=2:k-1
        k_means[i] = k_means[i-1] + k_interval
    end
    centers = reshape(k_means,1,k)
    R = kmeans!(reshape(likehood_rating,1,:), centers; maxiter=100)
    return R.centers, assignments(R)
end

function hdp_fit(data, α,γ,prior,iters, initial_custers = 1,burnout = 5;multiprocess=false)
    dim = size(data[1],1)
    gdim = dim
    gprior,lprior = create_default_priors(gdim,dim-gdim,:niw)
    return vhdp_fit(data,gdim, α,γ,α,prior,lprior,iters,initial_custers,burnout,multiprocess=multiprocess)
end

function vhdp_fit(data,gdim, α,γ,η,prior::Symbol,iters, initial_custers = 1,burnout = 5;multiprocess=false)
    dim = size(data[1],1)
    gprior,lprior = create_default_priors(gdim,dim-gdim,prior)
    return vhdp_fit(data,gdim, α,γ,η,gprior,lprior,iters, initial_custers,burnout,multiprocess=multiprocess)
end

function vhdp_fit(data,gdim, α,γ,η,gprior::distribution_hyper_params,lprior,iters, initial_custers = 1,burnout = 5;multiprocess=false)
    global random_seed = nothing
    global initial_local_clusters = initial_custers
    global initial_global_clusters = initial_custers
    global mp = multiprocess
    dim = size(data[1],1)
    model_hyperparams = model_hyper_params(gprior,lprior,α,γ,η,1.0,1.0,dim,gdim + 1)
    model = init_model(nothing; data = data , model_params = model_hyperparams)
    global posterior_history = []
    global word_ll_history = []
    global topic_count = []
    @everywhere global split_delays = true
    global burnout_period = burnout
    if mp
        for w in workers()
            @spawnat w set_burnout(burnout)
        end
    end
    global num_of_workers = nworkers()
    iter = 1
    total_time = 0

    init_first_clusters!(model)
    for i=1:iters
        tic = time()
        model_iteration(model,false,false,burnout)
        toc = time() -tic
        println("Iteration: " * string(i) * "|| Global Counts: " * string([x.clusters_count for x in global_clusters_vector]) * "|| iter time: " * string(toc))
        total_time+= toc

    end
    model.global_clusters = global_clusters_vector
    return model, total_time, posterior_history,word_ll_history,topic_count
end

function remove_all_zeros_from_data(data)
    concatenated_vals = reduce(hcat,values(data))
    all_zeros = []
    for i=1:size(concatenated_vals,1)
        if all(concatenated_vals[i,:] .== 0)
            push!(all_zeros,i)
        end
    end
    allvec = [i for i=1:size(concatenated_vals,1) if !(i in(all_zeros))]
    new_data = deepcopy(data)
    for (k,v) in data
        vlen = size(v,2)
        new_data[k] = v[allvec,:]
    end
    return new_data
end

function get_model_global_pred(model)
    return Dict([k=>create_global_labels(v) for (k,v) in model.groups_dict])
end

function results_stats(pred_dict, gt_dict)
    avg_nmi = 0
    for i=1:length(pred_dict)
        nmi = mutualinfo(pred_dict[i],gt_dict[i])
        avg_nmi += nmi
    end
    return avg_nmi / length(pred_dict)
end


text:
@naiveVariationalRule(:node_type     => NAutoRegressiveMovingAverageX,
                      :outbound_type => Message{GaussianMeanPrecision},
                      :inbound_types => (Nothing, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution),
                      :name          => VariationalNARMAXOutNPPPPPP)

@naiveVariationalRule(:node_type     => NAutoRegressiveMovingAverageX,
                      :outbound_type => Message{GaussianWeightedMeanPrecision},
                      :inbound_types => (ProbabilityDistribution, Nothing, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution),
                      :name          => VariationalNARMAXIn1PNPPPPP)

@naiveVariationalRule(:node_type     => NAutoRegressiveMovingAverageX,
                      :outbound_type => Message{GaussianWeightedMeanPrecision},
                      :inbound_types => (ProbabilityDistribution, ProbabilityDistribution, Nothing, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution),
                      :name          => VariationalNARMAXIn2PPNPPPP)

@naiveVariationalRule(:node_type     => NAutoRegressiveMovingAverageX,
                      :outbound_type => Message{GaussianWeightedMeanPrecision},
                      :inbound_types => (ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, Nothing, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution),
                      :name          => VariationalNARMAXIn3PPPNPPP)

@naiveVariationalRule(:node_type     => NAutoRegressiveMovingAverageX,
                      :outbound_type => Message{GaussianWeightedMeanPrecision},
                      :inbound_types => (ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, Nothing, ProbabilityDistribution, ProbabilityDistribution),
                      :name          => VariationalNARMAXIn4PPPPNPP)

@naiveVariationalRule(:node_type     => NAutoRegressiveMovingAverageX,
                      :outbound_type => Message{GaussianWeightedMeanPrecision},
                      :inbound_types => (ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, Nothing, ProbabilityDistribution),
                      :name          => VariationalNARMAXIn5PPPPPNP)

@naiveVariationalRule(:node_type     => NAutoRegressiveMovingAverageX,
                      :outbound_type => Message{Gamma},
                      :inbound_types => (ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, Nothing),
                      :name          => VariationalNARMAXIn6PPPPPPN)

max_stars_repo_name: wmkouw/NARMAX

max_issues_repo_name: wmkouw/NARMAX

max_forks_repo_name: wmkouw/NARMAX


text:
module My_package

# Write your package code here.
export test_fun

"""
    test_fun(n)

Returns sum of 1:n
"""
function test_fun(n)
return 1:n|> sum
end
end


text:
@testset "Ramification" begin
  K, a = Hecke.rationals_as_number_field()
  Q = Hecke.quaternion_algebra2(K, -1, -1)
  @test (@inferred Hecke.ramified_infinite_places(Q)) == real_places(K)

  Q = Hecke.quaternion_algebra2(QQ, -1, -1)
  QQQ = direct_product(Q, Q)[1]
  C, = center(QQQ)
  flds = first.(Hecke.components(Hecke.Field, C))
  @test (@inferred Hecke.ramified_infinite_places_of_center(QQQ)) == real_places.(flds)

  Q = matrix_algebra(QQ, 3)
  @test (@inferred schur_index(Q, inf)) == 1
  @test iseichler(Q)

  # Test (-1, -1/QQ)
  Q = Hecke.AlgQuat(QQ, QQ(-1), QQ(-1))
  @test !issplit(Q, 2)
  @test schur_index(Q, 2) == 2
  @test schur_index(Q, ZZ(2)) == 2
  @test schur_index(Q, inf) == 2
  @test !issplit(Q, inf)
  @test schur_index(Q, 3) == 1
  @test schur_index(Q, ZZ(3)) == 1
  @test issplit(Q, ZZ(3))
  @test schur_index(Q) == 2
  @test !issplit(Q)
  @test !iseichler(Q)

  Q, = AlgAss(Q)
  @test !issplit(Q, 2)
  @test schur_index(Q, 2) == 2
  @test schur_index(Q, ZZ(2)) == 2
  @test schur_index(Q, inf) == 2
  @test !issplit(Q, inf)
  @test schur_index(Q, 3) == 1
  @test schur_index(Q, ZZ(3)) == 1
  @test issplit(Q, ZZ(3))
  @test schur_index(Q) == 2
  @test !issplit(Q)
  @test !iseichler(Q)

  # Test Mat_2((-1, -1)/QQ)
  Q = Hecke.AlgQuat(QQ, QQ(-1), QQ(-1))
  A = matrix_algebra(QQ, Q, 2)
  @test !issplit(A, 2)
  @test schur_index(A, 2) == 2
  @test schur_index(A, ZZ(2)) == 2
  @test schur_index(A, inf) == 2
  @test !issplit(A, inf)
  @test schur_index(A, 3) == 1
  @test schur_index(A, ZZ(3)) == 1
  @test issplit(A, ZZ(3))
  @test schur_index(A) == 2
  @test !issplit(A)
  @test iseichler(A)

  K, = Hecke.rationals_as_number_field()
  OK = maximal_order(K)
  # Test (-1, -1; QQ)
  Q = Hecke.AlgQuat(K, K(-1), K(-1))
  p2 = prime_decomposition(OK, 2)[1][1]
  p3 = prime_decomposition(OK, 3)[1][1]
  rlp = real_places(K)[1]
  @test !issplit(Q, p2)
  @test schur_index(Q, p2) == 2
  @test schur_index(Q, rlp) == 2
  @test !issplit(Q, rlp)
  @test schur_index(Q, p3) == 1
  @test issplit(Q, p3)
  @test schur_index(Q) == 2
  @test !issplit(Q)
  @test !iseichler(Q)

  Q, = AlgAss(Q)
  @test !issplit(Q, p2)
  @test schur_index(Q, p2) == 2
  @test schur_index(Q, rlp) == 2
  @test !issplit(Q, rlp)
  @test schur_index(Q, p3) == 1
  @test issplit(Q, p3)
  @test schur_index(Q) == 2
  @test !issplit(Q)
  @test !iseichler(Q)

  A = Hecke.AlgQuat(QQ, QQ(1), QQ(1))
  @test schur_index(A) == 1
  @test issplit(A)
  @test iseichler(A)
  A = Hecke.AlgQuat(K, K(1), K(1))
  @test schur_index(A) == 1
  @test issplit(A)
  @test iseichler(A)

  K, = quadratic_field(5)
  A = Hecke.AlgQuat(K, K(1), K(1))
  @test schur_index(A) == 1
  @test issplit(A)
  @test iseichler(A)

  K, = quadratic_field(5)
  A = Hecke.AlgQuat(K, K(-1), K(-1))
  @test schur_index(A) == 2
  @test !issplit(A)
  @test !iseichler(A)

  A = matrix_algebra(K, 3)
  @test iseichler(A)

  Qx, x = QQ["x"]
  K, a = number_field(x^3 - 2, "a")
  A = Hecke.AlgQuat(K, K(-1), K(-1))
  @test schur_index(A) == 2
  @test !issplit(A)
  @test iseichler(A)
  p = complex_places(K)
  @test issplit(A, p[1])
end


text:
module Sophus

include("so2.jl")

end # module


text:
"""
Constrained-optimisation-based control allocator.
"""
struct ConstrainedAllocator <: StaticAllocator
    B
    u_min
    u_max
    u
    function ConstrainedAllocator(B, u_min, u_max; p=Inf)
        @assert size(u_min) == size(u_max)
        dim_input = length(u_min)
        u = Convex.Variable(dim_input)
        new(B, u_min, u_max, u)
    end
end

"""
# Variables
ν: virtual input
# Notes
ν = B*u where u: control input
"""
function (allocator::ConstrainedAllocator)(ν, Λ=Diagonal(ones(size(ν)));
                                           p=Inf, silent_solver=true)
    @unpack u_min, u_max, u, B = allocator
    prob = minimize(
                    norm(u, p)
                    # sum(u)
                    + 1e5*norm(ν - B*Λ*u, 1)  # exact penalty method
                   )
    prob.constraints += [
                         u .>= u_min;  # min
                         u .<= u_max;  # max
                         # ν == B*Λ*u  # equality; replaced by exact penalty method
                        ]
    Convex.solve!(prob,
                  # SCS.Optimizer();
                  Mosek.Optimizer();  # best
                  silent_solver=silent_solver,
                  verbose=false)
    u.value[:]  # dim: n×1 -> n
end


text:
using ClimateBase, Test, Dates
using Statistics
using ClimateBase.Interpolations
import StatsBase
Time = ClimateBase.Time
cd(@__DIR__)

# Create the artificial dimensional arrays
function monthly_insolation(t::TimeType, args...)
    d = ClimateBase.monthspan(t)
    mean(insolation(τ, args...) for τ in d)
end

lats = -86:4:86
lons = collect(0.5:10:360)
t = Date(2000, 3, 15):Month(1):Date(2020, 2, 15)
tdaily = Date(2000, 3, 15):Day(1):Date(2020, 3, 14)
thourly = DateTime(2000, 3, 15):Hour(1):DateTime(2000, 6, 14)

d = (Lon(lons), Lat(lats), Time(t))
A = zeros([length(x) for x in (lons, lats, t)]...)
B = copy(A)

# generate solar rad
for i in 1:length(lats)
    θ = lats[i]
    s = monthly_insolation.(t, θ)
    for j in 1:length(lons)
        c = cosd(lons[j])
        A[j, i, :] .= c.*s
        B[j, i, :] .= s
    end
end

A = ClimArray(A, d; name = "insolation")
B = ClimArray(B, d; attrib = Dict("a" => 2)) # on purpose without name

# Create similar arrays with Coord dimension
reduced_points = [round(Int, length(lons)*cosd(θ)) for θ in lats]
coords = ClimateBase.reduced_grid_to_points(lats, reduced_points)
sort!(coords; by = reverse) # critical!

coord_dim = Coord(coords, (Lon, Lat))
C = zeros(length(coords), length(t))
for (i, (lon, θ)) in enumerate(coords)
    C[i, :] .= monthly_insolation.(t, θ) .* cosd(lon)
end

C = ClimArray(C, (coord_dim, Time(t)); name = "has_coords", attrib = Dict("a" => 2))

# %%
@testset "ClimateBase.jl tests" begin
include("general_tests.jl")
include("temporal_tests.jl")
include("space_tests.jl")
include("space_coord_tests.jl")
include("io_tests.jl")
include("advanced_tests.jl")
end


text:
module TestExamples2DEulerMulticomponent

using Test
using Trixi

include("test_trixi.jl")

# pathof(Trixi) returns /path/to/Trixi/src/Trixi.jl, dirname gives the parent directory
EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "tree_2d_dgsem")

@testset "Compressible Euler Multicomponent" begin
  # NOTE: Some of the L2/Linf errors are comparably large. This is due to the fact that some of the
  #       simulations are set up with dimensional states. For example, the reference pressure in SI
  #       units is 101325 Pa, i.e., pressure has values of O(10^5)

  @trixi_testset "elixir_eulermulti_shock_bubble.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_shock_bubble.jl"),
      l2   = [7.37942029e+01, 8.93064845e-01, 5.79494667e+04, 1.82993538e-01, 1.11657182e-02],
      linf = [1.96909734e+02, 7.50913818e+00, 1.58891250e+05, 8.18676275e-01, 8.01984099e-02],
      tspan = (0.0, 0.001))
  end

  @trixi_testset "elixir_eulermulti_ec.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_ec.jl"),
      l2   = [0.050182236154087095, 0.050189894464434635, 0.2258715597305131, 0.06175171559771687],
      linf = [0.3108124923284472, 0.3107380389947733, 1.054035804988521, 0.29347582879608936])
  end

  @trixi_testset "elixir_eulermulti_es.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_es.jl"),
      l2   = [4.96549498e-02, 4.96558239e-02, 2.24252879e-01, 4.08716843e-03, 8.17433686e-03, 1.63486737e-02, 3.26973475e-02],
      linf = [2.48866593e-01, 2.48362494e-01, 9.31003945e-01, 1.74513568e-02, 3.49027135e-02, 6.98054270e-02, 1.39610854e-01])
  end

  @trixi_testset "elixir_eulermulti_convergence_ec.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_convergence_ec.jl"),
      l2   = [1.22902255e-04, 1.22902255e-04, 1.88673979e-04, 4.85423218e-05, 9.70846435e-05],
      linf = [6.72281924e-04, 6.72281924e-04, 1.26622928e-03, 2.84384418e-04, 5.68768837e-04])
  end

  @trixi_testset "elixir_eulermulti_convergence_es.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_convergence_es.jl"),
      l2   = [2.1855868685251012e-6, 2.1855868686504143e-6, 6.339327237279491e-6, 8.886020745630879e-7, 1.7772041491261757e-6],
      linf = [1.3465800098977354e-5, 1.3465800097645086e-5, 4.7306482458431276e-5, 5.182977349305062e-6, 1.0365954698610125e-5])
  end

  @trixi_testset "elixir_eulermulti_convergence_es.jl with flux_chandrashekar" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_convergence_es.jl"),
      l2   = [1.74779183e-06, 1.74779183e-06, 5.62550396e-06, 6.15682712e-07, 1.23136542e-06],
      linf = [1.51966473e-05, 1.51966473e-05, 5.51065526e-05, 5.58096521e-06, 1.11619304e-05],
      volume_flux = flux_chandrashekar)
  end
end

end # module


text:
################################################################################
## Copyright 2016. Juri Hinz and Jeremy Yee.
## Bellman recursion
################################################################################

## Bellman type definition

immutable Bellman

    ## Members
    value::Array{Float64, 4} # subgradient rep of value functions
    evalue::Array{Float64, 4} # subgradient rep of expected value functions
    policy::Array{Int64, 3} # prescribed policy
    method::String # slow, neighbours or fast method for continuation value
    grid::Array{Float64, 2} # the grid used
    control::Array{Float64} # prob for position transitions

    ## Parameters
    gnum::Int64 # grid size
    snum::Int64 # state size
    dnum::Int64 # disturbance sampling size
    pnum::Int64 # number of positions
    anum::Int64 # number of actions
    tnum::Int64 # number of time points
    fullcontrol::Bool # is the positions fully controllable

    ## Method using only row rearrangement

    function Bellman(css::CSS)

        method = "row-rearrange"

        ## Extract info from CSS
        gnum = css.gnum
        snum = css.snum
        dnum = css.dnum
        pnum = css.pnum
        anum = css.anum
        tnum = css.tnum
        fullcontrol = css.fullcontrol
        grid = css.grid
        control = css.control

        ## Containers
        value = zero(Array{Float64}(gnum, snum, pnum, tnum))
        evalue = zero(Array{Float64}(gnum, snum, pnum, tnum - 1))
        policy = zero(Array{Int64}(gnum, pnum, tnum - 1))
        container = zero(Matrix{Float64}(anum, snum))

        ## Backwards induction
        print("Bellman: ", tnum, ".")
        for i = 1:gnum
            for p = 1:pnum
                value[i, :, p, tnum] = css.Scrap(grid[i, :], p)
            end
        end
        t = tnum - 1
        while (t > 0)
            print(t, ".")
            for p = 1:pnum
                evalue[:, :, p, t] = expected(value[:, :, p, t + 1], grid,
                                              css.disturb, css.weight)
            end
            for i = 1:gnum
                for p = 1:pnum
                    container = zero(container)
                    if css.fullcontrol # full control of positions
                        for a = 1:anum
                            container[a, :] = evalue[i, :, control[p, a], t]
                            container[a, :] += css.Reward(t, grid[i, :], p, a)
                        end
                    else # partial control of positions
                        for a = 1:anum
                            for pp = 1:pnum
                                container[a, :] += control[p, a, pp] * evalue[i, :, pp, t]
                            end
                            container[a, :] += css.Reward(t, grid[i, :], p, a)
                        end
                    end
                    policy[i, p, t] = indmax(container * grid[i, :])
                    value[i, :, p, t] = container[policy[i, p, t], :]
                end
            end
            t = t - 1
        end

        print("End.\n")
        return new(value, evalue, policy, method, grid, control, gnum, snum,
                   dnum, pnum, anum, tnum, fullcontrol)

    end

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

    ## Method using nearest neighbours

    function Bellman(css::CSS, search::Function, k::Int64)

        method = "nearest-neighbors"

        ## Extract info from CSS
        gnum = css.gnum
        snum = css.snum
        dnum = css.dnum
        pnum = css.pnum
        anum = css.anum
        tnum = css.tnum
        fullcontrol = css.fullcontrol
        grid = css.grid
        control = css.control

        ## Containers
        value = zero(Array{Float64}(gnum, snum, pnum, tnum))
        evalue = zero(Array{Float64}(gnum, snum, pnum, tnum - 1))
        policy = zero(Array{Int64}(gnum, pnum, tnum - 1))
        container = zero(Matrix{Float64}(anum, snum))

        ## Batch processing of neighbours for faster speed
        indices = Array{Array{Int64, 1}}(gnum * dnum)
        temp_indices = Array{Array{Int64, 1}}(gnum)
        tgrid = transpose(grid)
        tree = search(tgrid)
        for i = 1:dnum
            temp_indices = knn(tree, css.disturb[:, :, i] * tgrid, k, true)[1]
            indices[((i - 1) * gnum + 1):(i * gnum)] = temp_indices
        end

        ## Backwards induction
        print("Bellman: ", tnum, ".")
        for i = 1:gnum
            for p = 1:pnum
                value[i, :, p, tnum] = css.Scrap(grid[i, :], p)
            end
        end
        t = tnum - 1
        while (t > 0)
            print(t, ".")
            for p = 1:pnum
                evalue[:, :, p, t] = expected_host(value[:, :, p, t + 1], grid,
                                                   css.disturb, css.weight, indices)
            end
            for i = 1:gnum
                for p = 1:pnum
                    container = zero(container)
                    if css.fullcontrol # full control of positions
                        for a = 1:anum
                            container[a, :] = evalue[i, :, control[p, a], t]
                            container[a, :] += css.Reward(t, grid[i, :], p, a)
                        end
                    else # partial control of positions
                        for a = 1:anum
                            for pp = 1:pnum
                                container[a, :] += control[p, a, pp] * evalue[i, :, pp, t]
                            end
                            container[a, :] += css.Reward(t, grid[i, :], p, a)
                        end
                    end
                    policy[i, p] = indmax(container * grid[i, :])
                    value[i, :, p, t] = container[policy[i, p], :]
                end
            end
            t = t - 1
        end

        print("End.\n")
        return new(value, evalue, policy, method, grid, control, gnum, snum,
                   dnum, pnum, anum, tnum, fullcontrol)

    end

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

    ## Method using the permutation matrices

    function Bellman(css::CSS, X::PermMat)

        method = "fast"

        ## Extract info from CSS
        gnum = css.gnum
        snum = css.snum
        dnum = css.dnum
        pnum = css.pnum
        anum = css.anum
        tnum = css.tnum
        fullcontrol = css.fullcontrol
        grid = css.grid
        control = css.control

        ## Containers
        value = zero(Array{Float64}(gnum, snum, pnum, tnum))
        evalue = zero(Array{Float64}(gnum, snum, pnum, tnum - 1))
        policy = zero(Array{Int64}(gnum, pnum, tnum - 1))
        container = zero(Matrix{Float64}(anum, snum))

        ## Backwards induction
        print("Bellman: ", tnum, ".")
        for i = 1:gnum
            for p = 1:pnum
                value[i, :, p, tnum] = css.Scrap(grid[i, :], p)
            end
        end
        t = tnum - 1
        while (t > 0)
            print(t, ".")
            for p = 1:pnum
                evalue[:, :, p, t] = expected_fast(value[:, :, p, t + 1], X)
            end
            for i = 1:gnum
                for p = 1:pnum
                    container = zero(container)
                    if css.fullcontrol # full control of positions
                        for a = 1:anum
                            container[a, :] = evalue[i, :, control[p, a], t]
                            container[a, :] += css.Reward(t, grid[i, :], p, a)
                        end
                    else # partial control of positions
                        for a = 1:anum
                            for pp = 1:pnum
                                container[a, :] += control[p, a, pp] * evalue[i, :, pp, t]
                            end
                            container[a, :] += css.Reward(t, grid[i, :], p, a)
                        end
                    end
                    policy[i, p] = indmax(container * grid[i, :])
                    value[i, :, p, t] = container[policy[i, p], :]
                end
            end
            t = t - 1
        end

        print("End.\n")
        return new(value, evalue, policy, method, grid, control, gnum, snum,
                   dnum, pnum, anum, tnum, fullcontrol)

    end

end


text:
function is_symbol_and_type(ex)
	isa(ex, Expr) && ex.head == :(::) && length(ex.args) == 2 && isa(ex.args[1],Symbol)
end

function parse_argument_data(exprs::Expr)
	vars = Any([nothing, nothing, []])
	for ex in exprs.args[1:end-1]

		## x::Mytype
		if is_symbol_and_type(ex)
			sym, typ = ex.args
			vars = hcat(vars, Any([nothing, nothing, []]))

			vars[1, end] = sym
			vars[2, end] = typ

		## a < x < b or a > x > b
		elseif isa(ex, Expr) && length(ex.args) == 5
			if is_symbol_and_type(ex.args[3])
				sym, typ = ex.args[3].args
			elseif isa(ex.args[3], Symbol)  # Default to Float64 here
				sym, typ = ex.args[3], Float64
			else
				error("Could not argument data $ex.")
			end
			vars = hcat(vars, Any([nothing, nothing, []]))
			vars[1, end] = sym
			vars[2, end] = typ

			if ex.args[2] == :< && ex.args[4] == :<     ## a<x<b
				push!(vars[3, end], ex.args[1], ex.args[5], true, true)
			elseif ex.args[2] == :> && ex.args[4] == :> ## a>x>b
				push!(vars[3, end], ex.args[5], ex.args[1], true, true)
			else
				error("Could not parse argument data $ex")
			end

		elseif isa(ex, Expr) && length(ex.args) == 3 ## a<<x<b, a>>x>b, a<x<<b or a>x>>b
			if ex.args[1] == :<

				## a<<x<b
				if isa(ex.args[2], Expr) &&	length(ex.args[2].args) == 3 && ex.args[2].args[1] == :<<
					if is_symbol_and_type(ex.args[2].args[3])
						sym, typ = ex.args[2].args[3].args
					elseif isa(ex.args[2].args[3], Symbol)  # Default to Float64 here
						sym, typ = ex.args[2].args[3], Float64
					else
						error("Could not argument data $ex.")
					end
					vars = hcat(vars, Any([nothing, nothing, []]))
					vars[1, end] = sym
					vars[2, end] = typ
					push!(vars[3, end], ex.args[2].args[2], ex.args[3], false, true)

				## a<x<<b
				elseif isa(ex.args[3], Expr) && length(ex.args[3].args) == 3 && ex.args[3].args[1] == :<<
					if is_symbol_and_type(ex.args[3].args[2])
						sym, typ = ex.args[3].args[2].args
					elseif isa(ex.args[3].args[2], Symbol)  # Default to Float64 here
						sym, typ = ex.args[3].args[2], Float64
					else
						error("Could not argument data $ex.")
					end
					vars = hcat(vars, Any([nothing, nothing, []]))
					vars[1, end] = sym
					vars[2, end] = typ
					push!(vars[3, end], ex.args[2], ex.args[3].args[3], true, false)
				else
					error("Could not argument data $ex.")
				end

			elseif ex.args[1] == :>

				## a>>x>b
				if isa(ex.args[2], Expr) && length(ex.args[2].args) == 3 && ex.args[2].args[1] == :>>
					if is_symbol_and_type(ex.args[2].args[3])
						sym, typ = ex.args[2].args[3].args
					elseif isa(ex.args[2].args[3], Symbol)  # Default to Float64 here
						sym, typ = ex.args[2].args[3], Float64
					else
						error("Could not argument data $ex.")
					end
					vars = hcat(vars, Any([nothing, nothing, []]))
					vars[1, end] = sym
					vars[2, end] = typ
					push!(vars[3, end], ex.args[3], ex.args[2].args[2], true, false)

				## a>x>>b
				elseif isa(ex.args[3], Expr) && length(ex.args[3].args) == 3 && ex.args[3].args[1] == :>>
					if is_symbol_and_type(ex.args[3].args[2])
						sym, typ = ex.args[3].args[2].args
					elseif isa(ex.args[3].args[2], Symbol)  # Default to Float64 here
						sym, typ = ex.args[3].args[2], Float64
					else
						error("Could not argument data $ex.")
					end
					vars = hcat(vars, Any([nothing, nothing, []]))
					vars[1, end] = sym
					vars[2, end] = typ
					push!(vars[end, 3], ex.args[3].args[3], ex.args[2], false, true)
				else
					error("Could not parse argument data $ex")
				end
			else
				error("Could not parse argument data $ex")
			end
		# push arguments to the last column
		elseif isa(vars[1, end], Symbol)
			push!(vars[3, end], ex)
		else
			error("Could not parse argument data $ex.")
		end
	end
	return vars[:,2:end]
end