text:
#  MIT License
#
#  Copyright (c) 2020 Michael J Simms. All rights reserved.
#
#  Permission is hereby granted, free of charge, to any person obtaining a copy
#  of this software and associated documentation files (the "Software"), to deal
#  in the Software without restriction, including without limitation the rights
#  to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
#  copies of the Software, and to permit persons to whom the Software is
#  furnished to do so, subject to the following conditions:
#
#  The above copyright notice and this permission notice shall be included in all
#  copies or substantial portions of the Software.
#
#  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
#  IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
#  FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
#  AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
#  LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
#  OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
#  SOFTWARE.

module Powers

export nearest_power_of_2

function nearest_power_of_2(num)
	n = num > 0 ? num - 1 : 0

	n |= n >> 1
	n |= n >> 2
	n |= n >> 4
	n |= n >> 8
	n |= n >> 16
	n = n + 1

	n
end

end

max_stars_repo_name: msimms/LibMath

max_issues_repo_name: msimms/LibMath

max_forks_repo_name: msimms/LibMath


text:
function invoke_tfunc(interp::CustomInterpreter, @nospecialize(ft), @nospecialize(types), @nospecialize(argtype), sv::InferenceState)
    ret = @invoke invoke_tfunc(interp::AbstractInterpreter, ft, types, argtype, sv::InferenceState)

    return ret
end

function builtin_tfunction(interp::CustomInterpreter, @nospecialize(f), argtypes::Array{Any,1},
                           sv::Union{InferenceState,Nothing})
    ret = @invoke builtin_tfunction(interp::AbstractInterpreter, f, argtypes::Array{Any,1},
                                    sv::Union{InferenceState,Nothing})

    return ret
end

function return_type_tfunc(interp::CustomInterpreter, argtypes::Vector{Any}, sv::InferenceState)
    ret = @invoke return_type_tfunc(interp::AbstractInterpreter, argtypes::Vector{Any}, sv::InferenceState)

    return ret
end

ext: jl

max_stars_repo_name: aviatesk/AbstractInterpreterPlayground

max_issues_repo_name: aviatesk/AbstractInterpreterPlayground

max_forks_repo_name: aviatesk/AbstractInterpreterPlayground


text:
# test_direct_discretization.jl

# direct discretization for elliptic problem in 2D/3D
function get_problem_2d(n)
    A = MultigridMethod((n,m) -> p₂₂(n), (n, n), V(3, 2))
    b = fill(1,(n-1)*(n-1))
    return (A,b)
end

@testset "Direct discretization of elliptic PDE, 2D" begin
    for n in [16 32 64 128 256 512 1024]
        (A, b) = get_problem_2d(n)
        A.grids[1].b .= b # copy rhs
        for item in Base.Iterators.take(A,15) end # iterate
        @test A.resnorm[end] < 1/n^2
        log(A,2)
    end
end

function get_problem_3d(n)
    A = MultigridMethod((n,m,k) -> p₂₃(n), (n, n, n), V(3, 2))
    b = fill(1,(n-1)*(n-1)*(n-1))
    return (A,b)
end

@testset "Direct discretization of elliptic PDE, 2D" begin
    for n in [4 8 16 32 64 128]
        (A, b) = get_problem_3d(n)
        A.grids[1].b .= b # copy rhs
        for item in Base.Iterators.take(A,15) end # iterate
        @test A.resnorm[end] < 1/n^3
        log(A,3)
    end
end


text:
using CorpusLoaders
using Test
using Base.Iterators
using MultiResolutionIterators

@testset "Basic use" begin
    train = load(CoNLL2000(), "train")
    test = load(CoNLL2000(), "test")

    @test length(train) > 0
    @test length(test) > 0
end

@testset "Using flatten_levels" begin
    train = load(CoNLL2000())
    num_sents = 10
    dataset = train[1:10]
    words = full_consolidate(flatten_levels(dataset, (!lvls)(CoNLL2000, :word)))
    @test length(words) > num_sents
    @test length(words) == sum(length.(dataset))
    @test typeof(words) == Vector{CorpusLoaders.POSTaggedWord}

    plain_words = word.(words)
    @test typeof(plain_words) <: Vector{String}
    pos_tags = part_of_speech.(words)
    @test typeof(pos_tags) <: Vector{String}
end


text:
import LinearAlgebra: det

#CirclePointRadius methods

CirclePointRadius(x::Int, y::Int, args...; kwargs...) = CirclePointRadius(Point(x,y), args...; kwargs...)
CirclePointRadius(p::CartesianIndex{2}, args...; kwargs...) = CirclePointRadius(Point(p), args...; kwargs...)

draw!(img::AbstractArray{T, 2}, circle::CirclePointRadius, color::T) where {T<:Colorant} = draw!(img, Ellipse(circle), color)

#CircleThreePoints methods

CircleThreePoints(x1::Int, y1::Int, x2::Int, y2::Int, x3::Int, y3::Int) =
    CircleThreePoints(Point(x1,y1), Point(x2,y2), Point(x3,y3))
CircleThreePoints(p1::CartesianIndex{2}, p2::CartesianIndex{2}, p3::CartesianIndex{2}) =
    CircleThreePoints(Point(p1), Point(p2), Point(p3))

function draw!(img::AbstractArray{T, 2}, circle::CircleThreePoints, color::T) where T<:Colorant
    ind = axes(img)
    x1 = circle.p1.x; y1 = circle.p1.y
    x2 = circle.p2.x; y2 = circle.p2.y
    x3 = circle.p3.x; y3 = circle.p3.y
    D = det([[x1, x2, x3] [y1, y2, y3] [1, 1, 1]])
    ! isapprox(D, 0, rtol = 0, atol = maximum(abs, (x1, x2, x3)) * maximum(abs, (y1, y2, y3)) * eps(Float64)) || error("Points do not lie on unique circle")
    R = [[y2 - y3, x3 - x2] [y2 - y1, x1 - x2]] * [(x1^2 + y1^2)-(x2^2 + y2^2), (x2^2 + y2^2)-(x3^2 + y3^2)] / (2 * D)
    ρ = euclidean([x1, y1], R)
    (first(ind[2]) <= R[1] <= last(ind[2]) && first(ind[1]) <= R[2] <= last(ind[1])) || error("Center of circle is out of the bounds of image")
    ρ - 1 <= minimum(abs, [R[1] - first(ind[2]), R[1] - last(ind[2]), R[2] - first(ind[1]), R[2] - last(ind[1])]) || error("Circle is out of the bounds of image : Radius is too large")
    draw!(img, CirclePointRadius(Point(round(Int,R[1]), round(Int,R[2])), ρ), color)
end


text:
using DataFrames
using Gadfly
svg = true

tags = [
"CPAR100000lin49x50T49.0TC",
"CPAR100000lin49x20T49.0TC",
"CPAR100000lin49x10T49.0TC",
]

function rmean(v)
    cumsum(v)./[1:length(v)]
end

theme = Theme(panel_fill=color("light grey"), grid_color=color("white"))
geo = Geom.line

for tag in tags
    theme = Theme(panel_fill=color("light grey"), grid_color=color("white"), default_color=color("black"),
    default_point_size=0.01cm)
    geo = Geom.line

    thetas = readcsv("$tag/thetas$tag.csv")'
    N = min(size(thetas,2), 100000)
    thetas = thetas[:, 1:10:N]
    N = size(thetas,2)

    dfth = DataFrame(n=1:N, c1=thetas[1,:][:], c2=thetas[2,:][:], c3=thetas[3,:][:], c4=thetas[4,:][:],
                            c5=thetas[5,:][:], c6=thetas[6,:][:], c7=thetas[7,:][:], c8=thetas[8,:][:],
    tr="iterates")
    dftrue = DataFrame(n=1:N, c1=0.1, c2=0.7, c3=0.35, c4=0.2, c5=0.1, c6=0.9, c7=0.3, c8=0.1, tr="Truth")
    dfrmean = DataFrame(n=1:N, c1=rmean(thetas[1,:][:]), c2=rmean(thetas[2,:][:]), c3=rmean(thetas[3,:][:]), c4=rmean(thetas[4,:][:]),
                               c5=rmean(thetas[5,:][:]), c6=rmean(thetas[6,:][:]), c7=rmean(thetas[7,:][:]), c8=rmean(thetas[8,:][:]),
    tr="Running mean")

    dftrue2 = stack(dftrue, [2,3,4,5,6,7,8,9])
    dfth2 = stack(dfth, [2,3,4,5,6,7,8,9])
    dfrm2 = stack(dfrmean, [2,3,4,5,6,7,8,9])
    df = vcat(dftrue2, dfth2, dfrm2)

    p = plot(df, x="n", y="value", color="tr", ygroup="variable", Geom.subplot_grid(Stat.yticks,  geo, free_y_axis=true),
    Guide.xlabel("Iteration"), Guide.ylabel("Parameters"), Guide.title("Trace"),  Guide.colorkey(""),
    Scale.discrete_color_manual("green","red","grey12"), Scale.x_continuous(format=:plain), Scale.y_continuous(format=:plain),
    theme
    )
    draw(PNG("$tag/thetas$tag.png",18cm, 22cm), p)
    draw(PDF("$tag/thetas$tag.pdf",18cm, 22cm), p)

end


text:
module StrategyExample

abstract type Algo end
struct Memoized <: Algo end
struct Iterative <: Algo end

# Memoization startegy
using Memoize
@memoize function _fib(n)
    n <= 2 ? 1 : _fib(n-1) + _fib(n-2)
end

function fib(::Memoized, n)
    println("Using memoization algorithm")
    _fib(n)
end

# Iterative strategy
function fib(algo::Iterative, n)
    println("Using iterative algorithm")
    n <= 2 && return 1
    prev1, prev2 = 1, 1
    local curr
    for i in 3:n
        curr = prev1 + prev2
        prev1, prev2 = curr, prev1
    end
    return curr
end

# From client's perpsective, it just needs to pass an algorithm
# that is implemented for the `fib` function.
# function find_golden_ratio(n; algo = Memoized())
#     return fib(algo, n) / fib(algo, n-1)
# end

# auto-selection
function fib(n)
    algo = n > 50 ? Memoized() : Iterative()
    return fib(algo, n)
end

function test()
    @show fib(30)
    @show fib(60)
end

end #module

using .StrategyExample
StrategyExample.test()


text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta

mutable struct CustomHttpPattern <: ProtoType
    kind::AbstractString
    path::AbstractString
    CustomHttpPattern(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct CustomHttpPattern

mutable struct HttpRule <: ProtoType
    selector::AbstractString
    get::AbstractString
    put::AbstractString
    post::AbstractString
    delete::AbstractString
    patch::AbstractString
    custom::CustomHttpPattern
    body::AbstractString
    response_body::AbstractString
    additional_bindings::Base.Vector{HttpRule}
    HttpRule(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct HttpRule
const __fnum_HttpRule = Int[1,2,3,4,5,6,8,7,12,11]
const __oneofs_HttpRule = Int[0,1,1,1,1,1,1,0,0,0]
const __oneof_names_HttpRule = [Symbol("pattern")]
meta(t::Type{HttpRule}) = meta(t, ProtoBuf.DEF_REQ, __fnum_HttpRule, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, __oneofs_HttpRule, __oneof_names_HttpRule, ProtoBuf.DEF_FIELD_TYPES)

mutable struct Http <: ProtoType
    rules::Base.Vector{HttpRule}
    fully_decode_reserved_expansion::Bool
    Http(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct Http

export Http, HttpRule, CustomHttpPattern

max_stars_repo_name: sivakon/grpc-julia-gcp-example

max_issues_repo_name: sivakon/grpc-julia-gcp-example

max_forks_repo_name: sivakon/grpc-julia-gcp-example


text:
using JLD2, FileIO
using RCall
using Distributions

function loadSingleObj(objPath)
  data = load(objPath)
  return data[collect(keys(data))[1]]
end

include("PreProcess.jl")
include("../sim_study/util.jl")
cbDataPath = "data/cytof_cb_with_nan.jld2"
cbReducedDataPath = "data/reduced_cb.jld2"
table = R"table"

cbData = loadSingleObj(cbDataPath)
# Reduce Data by removing highly non-expressive / expressive columns
goodColumns, J = PreProcess.preprocess!(cbData, maxNanOrNegProp=.9, maxPosProp=.9)

@save cbReducedDataPath deepcopy(cbData)

# Histograms
legend = R"legend"
R"require(rcommon)"

I = length(cbData)
J_new = size(cbData[1], 2)

mar_tmp = [2, 5.1, 2, 2.1]
imgDir = "results/misc/img/"
mkpath(imgDir)

println("Create index of missing obs")
for i in 1:I
  idx_missing = findall(isnan.(cbData[i]))
  open("$(imgDir)/index_of_missing_obs_sample$(i).txt", "w") do f
    for idx in idx_missing
      n, j = (idx[1], idx[2])
      write(f, "i: $i, n: $n, j: $j\n")
    end
  end
end

util.plotPdf("$(imgDir)/reducedDataHist.pdf");
R"par(mfrow=c(4, 2), mar=$mar_tmp, cex.main=.8)";
for i in 1:I
  for j in 1:J_new
    yij = cbData[i][:, j]
    yij_obs = filter(x -> !isnan(x), yij)
    n = length(yij)
    util.hist(yij, ylab="density", main="sample $i, marker $(findall(goodColumns)[j])",
              xlab="", border="transparent", col="steelblue",
              nclass=80, prob=true, xlim=[-10, 4])
    prop_missing = mean(isnan.(yij))
    ymean = mean(yij_obs)
    ymean_pos = mean(filter(x -> x > 0, yij_obs))
    ymean_neg = mean(filter(x -> x < 0, yij_obs))
    ysd = std(yij_obs)
    ysd_pos = std(filter(x -> x > 0, yij_obs))
    ysd_neg = std(filter(x -> x < 0, yij_obs))
    legend("topleft", bty="n", cex=1,
           legend=["prop. missing: $(round(prop_missing, digits=5))",
                   "num. obs: $(n)",
                   "mean: $(round(mean(yij_obs), digits=2))",
                   "mean (+): $(round(ymean_pos, digits=2))",
                   "mean (-): $(round(ymean_neg, digits=2))",
                   "sd: $(round(ysd, digits=2))",
                   "sd (+): $(round(ysd_pos, digits=2))",
                   "sd (-): $(round(ysd_neg, digits=2))"])
    util.abline(v=[ymean, ymean_pos, ymean_neg], col="red", lty=[1, 2, 2])
  end
end
R"par(mfrow=c(1, 1), oma=rcommon::oma.default(), mar=rcommon::mar.default(), cex.main=1)"
util.devOff()

# Count number of missing values per cell per sample
N = size.(cbData, 1)
num_missing = [[sum(isnan.(cbData[i][n, :])) for n in 1:N[i]] for i in 1:I]
for i in 1:I
  open("$(imgDir)/missing_values_per_cell_sample$(i).txt", "w") do f
    for n in 1:N[i]
      write(f, "$(num_missing[i][n])\n")
    end
  end
end

for i in 1:I
  util.plotPdf("$(imgDir)/missing_values_per_cell_sample$(i).pdf")
  util.plot(R"table($(num_missing[i])) / $(N[i])",
            xlim=[0, J_new],
            ylab="Proportion of cells missing number of markers",
            xlab="number of missing marker expression levels")
  util.devOff()
end


text:
#=
The Fibonacci sequence is defined by the recurrence relation:

Fn = Fn−1 + Fn−2, where F1 = 1 and F2 = 1.
Hence the first 12 terms will be:

F1 = 1
F2 = 1
F3 = 2
F4 = 3
F5 = 5
F6 = 8
F7 = 13
F8 = 21
F9 = 34
F10 = 55
F11 = 89
F12 = 144
The 12th term, F12, is the first term to contain three digits.

What is the first term in the Fibonacci sequence to contain 1000 digits?
=#
function calc()
  a = BigInt[1, 1]
  while length("$(last(a))") < 1_000
    a1 = a[length(a)] + a[length(a)-1]
    push!(a, a1)
  end
  length(a)
end
@time println(calc())

max_stars_repo_name: daniel-beard/JuliaProjectEuler

max_issues_repo_name: daniel-beard/JuliaProjectEuler

max_forks_repo_name: daniel-beard/JuliaProjectEuler


text:
@inline function ExaTronKernels.dmid(n::Int, x::ROCDeviceArray{Float64,1},
                      xl::ROCDeviceArray{Float64,1}, xu::ROCDeviceArray{Float64,1})
    tx = workitemIdx().x

    if tx <= n
        @inbounds x[tx] = max(xl[tx], min(x[tx], xu[tx]))
    end
    AMDGPU.sync_workgroup()

    return
end


text:
# TODO: Finish this test file
# TODO: Pull in existing tests into this file
# TODO: Rename to TestDataArray
module TestDataArrays
	using DataArrays

	# DataArray{T, N}(d::Array{T, N}, m::BitArray{N} = falses(size(d)))
	DataArray([1, 2], falses(2))
	DataArray(repeat([1, 2], outer = [1, 2]), falses(2, 2))
	DataArray(repeat([1, 2], outer = [1, 2, 2]), falses(2, 2, 2))

	DataArray([1, 2], trues(2))
	DataArray(repeat([1, 2], outer = [1, 2]), trues(2, 2))
	DataArray(repeat([1, 2], outer = [1, 2, 2]), trues(2, 2, 2))

	DataArray(['a', 'b'], falses(2))
	DataArray(repeat(['a', 'b'], outer = [1, 2]), falses(2, 2))
	DataArray(repeat(['a', 'b'], outer = [1, 2, 2]), falses(2, 2, 2))

	DataArray(['a', 'b'], trues(2))
	DataArray(repeat(['a', 'b'], outer = [1, 2]), trues(2, 2))
	DataArray(repeat(['a', 'b'], outer = [1, 2, 2]), trues(2, 2, 2))

	# DataArray{T, N}(d::Array{T, N}, m::BitArray{N} = falses(size(d)))
	DataArray([1, 2])
	DataArray(repeat([1, 2], outer = [1, 2]))
	DataArray(repeat([1, 2], outer = [1, 2, 2]))

	DataArray(['a', 'b'])
	DataArray(repeat(['a', 'b'], outer = [1, 2]))
	DataArray(repeat(['a', 'b'], outer = [1, 2, 2]))

	# DataArray(d::Array, m::Array{Bool})
	DataArray([1, 2], [false, false])
	DataArray(['a', 'b'], [false, false])

	DataArray([1, 2], [true, true])
	DataArray(['a', 'b'], [true, true])

	# DataArray(d::BitArray, m::BitArray = falses(size(d)))
	convert(DataArray, falses(2))
	convert(DataArray, trues(2))

	# DataArray(d::BitArray, m::BitArray = falses(size(d)))
	convert(DataArray, falses(2))
	convert(DataArray, trues(2))

	# DataArray(d::Ranges, m::BitArray = falses(length(d)))
	convert(DataArray, 1:2)
	convert(DataArray, 1:2)

	# DataArray(t::Type, dims::Integer...)
	DataArray(Float64, 2)
	DataArray(Float64, 2, 2)
	DataArray(Float64, 2, 2, 2)

	DataArray(Char, 2)
	DataArray(Char, 2, 2)
	DataArray(Char, 2, 2, 2)

	# DataArray{N}(t::Type, dims::NTuple{N,Int})
	DataArray(Float64, (2, ))
	DataArray(Float64, (2, 2))
	DataArray(Float64, (2, 2, 2))

	DataArray(Char, (2, ))
	DataArray(Char, (2, 2))
	DataArray(Char, (2, 2, 2))

	# Base.copy(d::DataArray)
	copy(DataArray([1, 2], falses(2)))

	# Base.deepcopy(d::DataArray)
	deepcopy(DataArray([1, 2], falses(2)))

	# Base.copy!(dest::DataArray, src::Any)
	da = DataArray([1, 2], falses(2))
	copy!(da, [3, 4])
	da

	# function Base.similar(d::DataArray, T::Type, dims::Dims)
	similar(DataArray([1, 2], falses(2)), Float64, 2)
	similar(DataArray([1, 2], falses(2)), Float64, 2, 2)
	similar(DataArray([1, 2], falses(2)), Float64, 2, 2, 2)

	# Base.size(d::DataArray) = size(d.data)
	size(DataArray([1, 2], falses(2)))
	size(DataArray(repeat([1, 2], outer = [1, 2]), falses(2, 2)))
	size(DataArray(repeat([1, 2], outer = [1, 2, 2]), falses(2, 2, 2)))

	# Base.ndims(d::DataArray) = ndims(d.data)
	ndims(DataArray([1, 2], falses(2)))
	ndims(DataArray(repeat([1, 2], outer = [1, 2]), falses(2, 2)))
	ndims(DataArray(repeat([1, 2], outer = [1, 2, 2]), falses(2, 2, 2)))

	# Base.length(d::DataArray) = length(d.data)
	length(DataArray([1, 2], falses(2)))
	length(DataArray(repeat([1, 2], outer = [1, 2]), falses(2, 2)))
	length(DataArray(repeat([1, 2], outer = [1, 2, 2]), falses(2, 2, 2)))

	# Base.endof(d::DataArray) = endof(d.data)
	endof(DataArray([1, 2], falses(2)))
	endof(DataArray(repeat([1, 2], outer = [1, 2]), falses(2, 2)))
	endof(DataArray(repeat([1, 2], outer = [1, 2, 2]), falses(2, 2, 2)))

	# Base.eltype{T, N}(d::DataArray{T, N}) = T
	eltype(DataArray([1, 2], falses(2)))
	eltype(DataArray(repeat([1, 2], outer = [1, 2]), falses(2, 2)))
	eltype(DataArray(repeat([1, 2], outer = [1, 2, 2]), falses(2, 2, 2)))

	# Base.find(da::AbstractDataArray{Bool}) = find(array(da, false))
	find(DataArray([false, true, true, false]))

	# function array{T}(da::DataArray{T})
	# array(DataArray([1, 0, 3], [false, true, false]))
	# array(DataArray([1, 2, 3], [false, false, false]))

	# function array{T}(da::DataArray{T}, replacement::T)
	array(DataArray([1, 0, 3], [false, true, false]), -1)
	array(DataArray([1, 2, 3], [false, false, false]), -1)

	# dropna(da::DataArray)
	dropna(DataArray([1, 0, 3], [false, true, false]))
	dropna(DataArray([1, 2, 3], [false, false, false]))
	# dropna{T}(da::AbstractDataVector{T})
	# dropna(@data([1, NA, 3]))

	# Iterators

	# GetIndex
	# Base.getindex(x::Vector, inds::AbstractDataVector{Bool})
	dinds = @data([true, false, false])
	[1, 2, 3][dinds]

	# Base.getindex(x::Vector, inds::AbstractDataArray{Bool})
	dinds = @data([true, false, false])
	[1, 2, 3][dinds]

	# Base.getindex(x::Array, inds::AbstractDataVector{Bool})
	dinds = @data([true, false, false])
	[1 2; 3 4][dinds]

	# Base.getindex(x::Array, inds::AbstractDataArray{Bool})
	dinds = @data([true, false, false])
	[1 2; 3 4][dinds]

	# Base.getindex{S, T}(x::Vector{S}, inds::AbstractDataArray{T})
	dinds = @data([1, 2, NA])
	[1.0, 2.0, 3.0, 4.0][dinds]

	# Base.getindex{S, T}(x::Array{S}, inds::AbstractDataArray{T})
	dinds = @data([1, 2, NA])
	[1.0 2.0; 3.0 4.0][dinds]

	# Base.getindex(x::Array, inds::AbstractDataVector{Bool})
	dinds = @data([true, false, false])
	[1 2; 3 4][dinds]

	# Base.getindex(x::Array, inds::AbstractDataArray{Bool})
	dinds = @data([true, false, false])
	[1 2; 3 4][dinds]

	# Base.getindex(d::DataArray, i::SingleIndex)
	da = @data([1, 2, NA, 4])
	da[1]
	da[3]
	da[1.0]
	da[3.0]

	# Base.getindex(d::DataArray, inds::AbstractDataVector{Bool})
	da = @data([1, 2, NA, 4])
	dinds = @data([true, false, false, NA])
	da[dinds]

	# Base.getindex(d::DataArray, inds::AbstractDataVector)
	da = @data([1, 2, NA, 4])
	dinds = @data([1, 2, NA, 2])
	da[dinds]

	# Base.getindex{T <: Number, N}(d::DataArray{T,N}, inds::BooleanIndex)
	# da = @data([1, 2, NA, 4])
	# inds = [1, 2, NA, 2]
	# da[inds]

	# Base.getindex(d::DataArray, inds::BooleanIndex)
	# da = @data([1.0, 2.0, NA, 4.0])
	# inds = [1, 2, NA, 2]
	# da[inds]

	# Base.getindex{T <: Number, N}(d::DataArray{T, N}, inds::MultiIndex)
	da = @data([1.0, 2.0, NA, 4.0])
	inds = [1, 2, 2]
	da[inds]

	# Base.getindex(d::DataArray, inds::MultiIndex)
	da = @data([1.0, 2.0, NA, 4.0])
	inds = [1, 2, 2]
	da[inds]

	# Base.getindex{T <: Number, N}(d::DataArray{T, N}, inds::BooleanIndex)
	da = @data([1.0, 2.0, NA, 4.0])
	inds = [true, true, false, false]
	da[inds]

	# Base.getindex{T <: Number, N}(d::DataArray{T, N}, inds::MultiIndex)
	da = @data([1.0, 2.0, NA, 4.0])
	inds = [1, 2, 2]
	da[inds]

	# Base.setindex!(da::DataArray, val::NAtype, i::SingleIndex)
	da = @data([1.0, 2.0, NA, 4.0])
	da[1] = NA

	# Base.setindex!(da::DataArray, val::Any, i::SingleIndex)
	da = @data([1.0, 2.0, NA, 4.0])
	da[1] = 3.0

	# Base.setindex!(da::DataArray{NAtype}, val::NAtype, inds::AbstractVector{Bool})
	# da = DataArray([NA, NA], falses(2))
	# da[[true, false]] = NA

	# Base.setindex!(da::DataArray{NAtype}, val::NAtype, inds::AbstractVector)
	# da = DataArray([NA, NA], falses(2))
	# da[[1, 2]] = NA

	# Base.setindex!(da::DataArray, val::NAtype, inds::AbstractVector{Bool})
	da = @data([1, 2])
	da[[true, false]] = NA

	# Base.setindex!(da::DataArray, val::NAtype, inds::AbstractVector)
	da = @data([1, 2])
	da[[1, 2]] = NA

	# Base.setindex!(da::AbstractDataArray, vals::AbstractVector, inds::AbstractVector{Bool})
	da = @data([1, 2])
	da[[true, false]] = [3]

	# Base.setindex!(da::AbstractDataArray, vals::AbstractVector, inds::AbstractVector)
	da = @data([1, 2])
	da[[1, 2]] = [3, 4]

	# Base.setindex!{T}(da::AbstractDataArray{T}, val::Union(Number, String, T), inds::AbstractVector{Bool})
	da = @data([1, 2])
	da[[true, false]] = 5

	# Base.setindex!{T}(da::AbstractDataArray{T}, val::Union(Number, String, T), inds::AbstractVector)
	da = @data([1, 2])
	da[[1, 2]] = 5

	# Base.setindex!(da::AbstractDataArray, val::Any, inds::AbstractVector{Bool})
	da = @data([1, 2])
	da[[true, false]] = 5.0

	# Base.setindex!{T}(da::AbstractDataArray{T}, val::Any, inds::AbstractVector)
	da = @data([1, 2])
	da[[1, 2]] = 5

	# isna(a::AbstractArray)
	isna([1, 2])
	isna(repeat([1, 2], outer = [1, 2]))
	isna(repeat([1, 2], outer = [1, 2, 2]))

	# isna(da::DataArray)
	isna(DataArray([1, 2], falses(2)))
	isna(DataArray(repeat([1, 2], outer = [1, 2]), falses(2, 2)))
	isna(DataArray(repeat([1, 2], outer = [1, 2, 2]), falses(2, 2, 2)))

	# Base.isnan(da::DataArray)
	isnan(DataArray([1, 2], falses(2)))
	isnan(DataArray(repeat([1, 2], outer = [1, 2]), falses(2, 2)))
	isnan(DataArray(repeat([1, 2], outer = [1, 2, 2]), falses(2, 2, 2)))

	# Base.isfinite(da::DataArray)
	isfinite(DataArray([1, 2], falses(2)))
	isfinite(DataArray(repeat([1, 2], outer = [1, 2]), falses(2, 2)))
	isfinite(DataArray(repeat([1, 2], outer = [1, 2, 2]), falses(2, 2, 2)))

	# anyna(a::AbstractArray)
	anyna([1, 2])
	anyna(repeat([1, 2], outer = [1, 2]))
	anyna(repeat([1, 2], outer = [1, 2, 2]))

	# anyna(d::AbstractDataArray)
	anyna(DataArray([1, 2], falses(2)))
	anyna(DataArray(repeat([1, 2], outer = [1, 2]), falses(2, 2)))
	anyna(DataArray(repeat([1, 2], outer = [1, 2, 2]), falses(2, 2, 2)))

	# allna(a::AbstractArray)
	allna([1, 2])
	allna(repeat([1, 2], outer = [1, 2]))
	allna(repeat([1, 2], outer = [1, 2, 2]))

	# allna(d::AbstractDataArray)
	allna(DataArray([1, 2], falses(2)))
	allna(DataArray(repeat([1, 2], outer = [1, 2]), falses(2, 2)))
	allna(DataArray(repeat([1, 2], outer = [1, 2, 2]), falses(2, 2, 2)))

	# Base.start(x::AbstractDataArray)
	start(DataArray([1, 2], falses(2)))
	start(DataArray(repeat([1, 2], outer = [1, 2]), falses(2, 2)))
	start(DataArray(repeat([1, 2], outer = [1, 2, 2]), falses(2, 2, 2)))

	# Base.next(x::AbstractDataArray, state::Integer)
	next(DataArray([1, 2], falses(2)), 1)
	next(DataArray(repeat([1, 2], outer = [1, 2]), falses(2, 2)), 1)
	next(DataArray(repeat([1, 2], outer = [1, 2, 2]), falses(2, 2, 2)), 1)

	# Base.done(x::AbstractDataArray, state::Integer)
	done(DataArray([1, 2], falses(2)), 1)
	done(DataArray(repeat([1, 2], outer = [1, 2]), falses(2, 2)), 1)
	done(DataArray(repeat([1, 2], outer = [1, 2, 2]), falses(2, 2, 2)), 1)

	# Base.convert{N}(::Type{BitArray{N}}, d::DataArray{BitArray{N}, N})

	# Base.convert{T, N}(::Type{BitArray{N}}, d::DataArray{T, N})
	# convert(BitArray, DataArray([1, 2], falses(2)))
	# convert(BitArray, DataArray(repeat([1, 2], outer = [1, 2]), falses(2, 2)))
	# convert(BitArray, DataArray(repeat([1, 2], outer = [1, 2, 2]), falses(2, 2, 2)))

	# Base.convert{S, T, N}(::Type{Array{S, N}}, x::DataArray{T, N})
	convert(Array{Float64}, DataArray([1, 2], falses(2)))
	convert(Array{Float64}, DataArray(repeat([1, 2], outer = [1, 2]), falses(2, 2)))
	convert(Array{Float64}, DataArray(repeat([1, 2], outer = [1, 2, 2]), falses(2, 2, 2)))

	# Base.convert{T, N}(::Type{Array{T, N}}, x::DataArray{T, N})
	convert(Array, DataArray([1, 2], falses(2)))
	convert(Array, DataArray(repeat([1, 2], outer = [1, 2]), falses(2, 2)))
	convert(Array, DataArray(repeat([1, 2], outer = [1, 2, 2]), falses(2, 2, 2)))

	# Base.convert{S, T, N}(::Type{DataArray{S, N}}, x::Array{T, N})
	convert(DataArray{Float64}, [1, 2])
	convert(DataArray{Float64}, repeat([1, 2], outer = [1, 2]))
	convert(DataArray{Float64}, repeat([1, 2], outer = [1, 2, 2]))

	# Base.convert{T, N}(::Type{DataArray}, x::Array{T, N})
	convert(DataArray, [1, 2])
	convert(DataArray, repeat([1, 2], outer = [1, 2]))
	convert(DataArray, repeat([1, 2], outer = [1, 2, 2]))

	# Base.convert{S, T, N}(::Type{DataArray{S, N}}, x::DataArray{T, N})
	convert(DataArray{Float64}, DataArray([1, 2], falses(2)))
	convert(DataArray{Float64}, DataArray(repeat([1, 2], outer = [1, 2]), falses(2, 2)))
	convert(DataArray{Float64}, DataArray(repeat([1, 2], outer = [1, 2, 2]), falses(2, 2, 2)))

	# Base.convert{T, N}(::Type{DataArray}, x::DataArray{T, N})
	convert(DataArray, DataArray([1, 2], falses(2)))
	convert(DataArray, DataArray(repeat([1, 2], outer = [1, 2]), falses(2, 2)))
	convert(DataArray, DataArray(repeat([1, 2], outer = [1, 2, 2]), falses(2, 2, 2)))

	# int(da::DataArray)
	int(DataArray([1, 2], falses(2)))
	int(DataArray(repeat([1, 2], outer = [1, 2]), falses(2, 2)))
	int(DataArray(repeat([1, 2], outer = [1, 2, 2]), falses(2, 2, 2)))

	# float(da::DataArray)
	float(DataArray([1, 2], falses(2)))
	float(DataArray(repeat([1, 2], outer = [1, 2]), falses(2, 2)))
	float(DataArray(repeat([1, 2], outer = [1, 2, 2]), falses(2, 2, 2)))

	# bool(da::DataArray)
	bool(DataArray([1, 0], falses(2)))
	bool(DataArray(repeat([1, 0], outer = [1, 2]), falses(2, 2)))
	bool(DataArray(repeat([1, 0], outer = [1, 2, 2]), falses(2, 2, 2)))

	# Base.hash(a::AbstractDataArray)
	hash(DataArray([1, 2], falses(2)))
	hash(DataArray(repeat([1, 2], outer = [1, 2]), falses(2, 2)))
	hash(DataArray(repeat([1, 2], outer = [1, 2, 2]), falses(2, 2, 2)))
end


text:
@testset "Escaping" begin
    import InfluxDB: escape_field, escape_tag, escape_field_key, escape_tag_key, escape_measurement
    dontescape = join(['a':'z'; 'A':'Z'; '0':'9'; "!#\$%&'()*+-./:;<>?@[]^_`{|}~"; "ᕙ༼*◕_◕*༽ᕤ😆😎😵😗😈🐻💓🌼🎁🎂📱💻"])

    server = InfluxServer("http://localhost:8086")
    create_database(server, "escape_test")

    @testset "Escaping measurement" begin
        @test_throws ErrorException escape_measurement("\n")
        @test_throws ErrorException escape_measurement("a\\")

        tests = [
            dontescape => dontescape
            "a\\a"     => "a\\a"
            "b "       => "b\\ "
            "c,"       => "c\\,"
            "d="       => "d\\="
            "e\""      => "e\\\""
            "f\\\\f"   => "f\\\\f"
            "g\\ "     => "g\\\\\\ "
            "h\\,"     => "h\\\\\\,"
            "i\\="     => "i\\\\\\="
            "j\\\""    => "j\\\\\\\""
        ]
        write(server, "escape_test", Measurement.(first.(tests), Ref(Dict("a"=>"1")), Ref(Dict("b"=>"2"))))
        meas, = InfluxDB.list_measurements(server, "escape_test")

        for (str, exp) in tests
            @test escape_measurement(str) == exp
            @test str ∈ meas.name
        end
    end

    @testset "Escaping field" begin

        tests = [
            dontescape*'\n' => '"' * dontescape*'\n' * '"'
            "a\\"      => "\"a\\\\\""
            "b "       => "\"b \""
            "c,"       => "\"c,\""
            "d="       => "\"d=\""
            "e\""      => "\"e\\\"\""
            "f\\\\"    => "\"f\\\\\\\\\""
            "g\\ "     => "\"g\\\\ \""
            "h\\,"     => "\"h\\\\,\""
            "i\\="     => "\"i\\\\=\""
            "j\\\""    => "\"j\\\\\\\"\""
        ]
        for (i, (str, exp)) in enumerate(tests)
            @test escape_field(str) == exp
            write(server, "escape_test", [Measurement("test_field", Dict("a" => str, "i" => i), Dict("b"=>"1"), i)])
            df, = query(server, "escape_test", SELECT(;measurements=["test_field"], condition="\"i\"=$i"))
            @test df.a[1] == str
        end

        @test escape_field(123) == "123"
        @test escape_field(1.0) == "1.0"
        @test escape_field(true) == "true"

        m = Measurement(
            "test_field",
            Dict("c" => 123, "d" => 1.4, "e" => false, "i" => 1000),
            Dict("b" => "1"),
            1000
        )

        write(server, "escape_test", [m])
        df, = query(server, "escape_test", SELECT(;measurements=["test_field"], condition="\"i\"=1000"))
        @test df.c[1] == 123
        @test df.d[1] == 1.4
        @test df.e[1] == false
    end

    @testset "Escaping tag" begin
        @test_throws ErrorException escape_tag("\n")
        @test_throws ErrorException escape_tag("a\\")
        tests = [
            dontescape => dontescape
            "a\\a"     => "a\\a"
            "b "       => "b\\ "
            "c,"       => "c\\,"
            "d="       => "d\\="
            "e\""      => "e\""
            "f\\\\f"   => "f\\\\f"
            "g\\ "     => "g\\\\ "
            "h\\,"     => "h\\\\,"
            "i\\="     => "i\\\\="
            "j\\\""    => "j\\\""
        ]
        for (i, (str, exp)) in enumerate(tests)
            @test escape_tag(str) == exp
            write(server, "escape_test", [Measurement("test_tag", Dict("i" => i), Dict("a"=>str), i)])
            df, = query(server, "escape_test", SELECT(;measurements=["test_tag"], condition="\"i\"=$i"))
            @test df.a[1] == str
        end
    end

    @testset "Escaping field key" begin
        @test_throws ErrorException escape_field_key("\n")
        @test_throws ErrorException escape_field_key("a\\")

        tests = [
            dontescape => dontescape
            "a\\a"     => "a\\a"
            "b "       => "b\\ "
            "c,"       => "c\\,"
            "d="       => "d\\="
            "e\""      => "e\\\""
            "f\\\\f"   => "f\\\\f"
            # "g\\ "     => "g\\\\ " # ???
            # "h\\,"     => "h\\\\," # ???
            # "i\\="     => "i\\\\=" # ???
            # "j\\\""    => "j\\\"" # ???
        ]

        for (i, (str, exp)) in enumerate(tests)
            @test escape_field_key(str) == exp
            write(server, "escape_test", [Measurement("test_field_key", Dict("i" => i, str => true), Dict("a"=>"1"), i)])
            df, = query(server, "escape_test", SELECT(;measurements=["test_field_key"], condition="\"i\"=$i"))
            @test str ∈ names(df) && df[1,str] == true
        end
    end

    @testset "Escaping tag key" begin
        @test_throws ErrorException escape_tag_key("\n")
        @test_throws ErrorException escape_tag_key("a\\")

        tests = [
            dontescape => dontescape
            "a\\a"     => "a\\a"
            "b "       => "b\\ "
            "c,"       => "c\\,"
            "d="       => "d\\="
            "e\""      => "e\""
            "f\\\\f"   => "f\\\\f"
            # "g\\ "     => "g\\\\ " # ???
            # "h\\,"     => "h\\\\," # ???
            # "i\\="     => "i\\\\=" # ???
            # "j\\\""    => "j\\\"" # ???
        ]

        for (i, (str, exp)) in enumerate(tests)
            @test escape_tag_key(str) == exp
            write(server, "escape_test", [Measurement("test_tag_key", Dict("i" => i), Dict("a"=>"1", str => "1"), i)])
            df, = query(server, "escape_test", SELECT(;measurements=["test_tag_key"], condition="\"i\"=$i"))
            @test str ∈ names(df) && df[1,str] == "1"
        end
    end
end


text:
function beam01(E,rho,A,I,L,ne,k1,k2)
	#
	# mass and stiffness matrices for a beam using Hermite elements
	# degres of freedom of the element: [y1 rot1 y2 rot2]
	#
	# E:   elasticity modulus
	# rho: density
	# A:   element section
	# I:   inercia
	# L:   length of the beam
	# ne:  number of elements
	# k1:  added stiffness in vertical direction of first node
	# k2:  added stiffness in vertical direction of last node
	#
	# javier.cara@upm.es,
	#

	# local matrices
	# *******************************

	# element length
	le = L/ne

	# matriz de rigidez del elemento
	ke = [12 6*le -12 6*le;
	6*le 4*le^2 -6*le 2*le^2;
	-12 -6*le 12 -6*le;
	6*le 2*le^2 -6*le 4*le^2]
	ke = E*I/(le^3)*ke

	# matriz de masa del elemento (consistente)
	me = [156 22*le 54 -13*le;
	22*le 4*le^2 13*le -3*le^2;
	54 13*le 156 -22*le;
	-13*le -3*le^2 -22*le 4*le^2]
	me = rho*A*le/420*me

	# matrices globales
	# *******************************

	# numero de grados de libertad
	ndof = (ne+1)*2

	K = zeros(ndof,ndof)
	M = zeros(ndof,ndof)
	for e=1:ne
		 dof1 = (e-1)*2 + 1
		 dof2 = (e-1)*2 + 4
		 K[dof1:dof2,dof1:dof2] = K[dof1:dof2,dof1:dof2] + ke
		 M[dof1:dof2,dof1:dof2] = M[dof1:dof2,dof1:dof2] + me
	end

	# muelles en los extremos en direccion y
	K[1,1] = K[1,1] + k1
	K[2*ne+1,2*ne+1] = K[2*ne+1,2*ne+1] + k2

	return K,M

end

# ================================================================
function beam01_test()

	E = 2.1e11
	rho = 7850.
	b = 2.
	h = 0.20
	L = 20.
	ne = 18

	I = 1/12*b*h^3
	A = b*h

	k1 = 100*E*I
	k2 = 100*E*I

	(K,M) = beam01(E,rho,A,I,L,ne,k1,k2)

	# frecuencias naturales del model FEM
	(W,V) = eig(K,M)
	wo = sqrt(W) # (rad/s)

	# frecuencias naturales teoricas
	# wn=Cn*sqrt(E*I/(rho**L^4), donde Cn=(n*pi)^2
	wn = zeros(5)
	for m=1:5
		 wn[m] = (m*pi)^2*sqrt(E*I/(rho*A*L^4))
	end

	return wn,wo[1:5]

end


text:
struct Figure5
    colors::Vector{Symbol}
    styles::Vector{Symbol}
    names::Vector{String}
    θ_default::Vector{Float64}
    α::Vector{Float64}
    ωₚ::Vector{Float64}
    γ::Vector{Float64}
    r::Vector{Float64}
    ϕₘ::Vector{Float64}
    β::Vector{Float64}
    ν::Vector{Float64}
    σₗ::Vector{Float64}
    σᵣ::Vector{Float64}
    function Figure5(;
            colors = parametercolors(),
            styles = [:dot, :dash, :dashdot, :solid],
            names = ["\\alpha", "\\omega_p", "\\gamma", "r", "\\omega_p", "\\beta", "\\nu", "\\sigma_l", "\\sigma_r"],
            θ_default = [0.7,0.8,3.3,5,pi,4,2.7,0.55,0.26],
            α = collect(0.4:0.2:1.0),
            ωₚ = collect(0.8:0.1:1.1),
            γ = [1,2.5,4,5.5],
            r = collect(4:0.5:5.5),
            ϕₘ = collect(range(π/2,3π/2,length=4)),
            β = collect(3:6),
            ν = collect(1:4),
            σₗ = collect(0.45:0.05:0.6),
            σᵣ = collect(0.15:0.05:0.3)
        )
        @assert all(length(styles) == length(a) for a in (α, ωₚ, γ, r, ϕₘ, β, ν, σₗ, σᵣ))
        new(colors,styles,names,θ_default,α, ωₚ, γ, r, ϕₘ, β, ν, σₗ, σᵣ)
    end
end
ϕₛ(ω,θ) = θ[6]*exp(-θ[7] * (abs(ω) > θ[2] ? θ[2]/abs(ω) : 1))
σ(ω,θ) = θ[8] - θ[9]/3 * (4*(θ[2]/abs(ω))^2 - (θ[2]/abs(ω))^8)
function makefulltheta(θ_default,val,ind)
    θ = copy(θ_default)
    θ[ind] = val
    return θ
end

@recipe function f(f::Figure5)
    layout := (3,3)
    Ω = range(0,2,length=100)
    for (i,parameter) ∈ enumerate([:α,:ωₚ,:γ,:r]), (j,par) ∈ enumerate(getfield(f,parameter))
        @series begin
            subplot := i
            ylims := (0,1.2)
            seriescolor := f.colors[i]
            linestyle := f.styles[j]
            label := L"%$(f.names[i]) = %$par"
            if i < 4
                yguide := L"f(\omega)"
            elseif i == 4
                yguide := L"10\,\log_{10}\;f(\omega)"
                yaxis := :log10
                ylims := (0.001,1.2^10)
                return Ω[2:end], ω -> (real(sdf(JONSWAP{0}(makefulltheta(f.θ_default,par,i)[1:4]),ω))) ^ 10 # for 10log10
            end
            return Ω, ω -> real(sdf(JONSWAP{0}(makefulltheta(f.θ_default,par,i)[1:4]),ω))
        end
    end
    for (i,parameter) ∈ enumerate([:ϕₘ,:β,:ν]), (j,par) ∈ enumerate(getfield(f,parameter))
        @series begin
            subplot := i+4
            seriescolor := f.colors[i+4]
            linestyle := f.styles[j]
            label := L"%$(f.names[i+4]) = %$par"
            Ω, ω -> makefulltheta(f.θ_default, par, i+4)[5] + ϕₛ(ω,makefulltheta(f.θ_default, par, i+4))/2
        end
        @series begin
            subplot := i+4
            ylims := (pi-2.1,pi+2.1)
            seriescolor := f.colors[i+4]
            linestyle := f.styles[j]
            label := false
            yguide := L"$\phi_{m2}(\omega)$ and $\phi_{m1}(\omega)$"
            if i == 3
                xguide := "angular frequency"
            end
            Ω, ω -> makefulltheta(f.θ_default, par, i+4)[5] - ϕₛ(ω,makefulltheta(f.θ_default, par, i+4))/2
        end
    end
    for (i,parameter) ∈ enumerate([:σₗ,:σᵣ]), (j,par) ∈ enumerate(getfield(f,parameter))
        @series begin
            subplot := i+7
            ylims := (0.08,1)
            seriescolor := f.colors[i+7]
            linestyle := f.styles[j]
            label := L"%$(f.names[i+7]) = %$par"
            yguide := L"\sigma(\omega)"
            xguide := "angular frequency"
            Ω, ω -> σ(ω,makefulltheta(f.θ_default,par,i+7))
        end
    end
end


text:
inner(x) = dot(x,x)
outer(x) = x*x'
sym(x) = Symmetric(x)
sym(x::Number) = x
_logdet(Σ, d::Integer) = LinearAlgebra.logdet(Σ)
lchol(Σ) = cholesky(sym(Σ)).U'
lchol(Σ::Number) = sqrt(Σ)
cholesky_(x) = cholesky(x)
cholesky_(x::Number) = x

struct Zero <: AbstractVector{Bool}
end
Base.:+(x::AbstractVector, z::Zero) = x
Base.:+(z::Zero, x::AbstractVector) = x

Base.:*(x::AbstractMatrix, z::Zero) = z
Base.:-(x::AbstractVector, z::Zero) = x
flat(x) = collect(Iterators.flatten(x))


text:
abstract type distribution_hyper_params end
#Suff statistics must contain N which is the number of points associated with the cluster
abstract type sufficient_statistics end
abstract type distibution_sample end
import Base.copy

struct model_hyper_params
    distribution_hyper_params::distribution_hyper_params
    α::Float32
    total_dim::Int64
end

mutable struct cluster_parameters
    hyperparams::distribution_hyper_params
    distribution::distibution_sample
    suff_statistics::sufficient_statistics
    posterior_hyperparams::distribution_hyper_params
end

mutable struct splittable_cluster_params
    cluster_params::cluster_parameters
    cluster_params_l::cluster_parameters
    cluster_params_r::cluster_parameters
    lr_weights::AbstractArray{Float32, 1}
    splittable::Bool
    logsublikelihood_hist::AbstractArray{Float32,1}
end

mutable struct thin_cluster_params{T <: distibution_sample}
    cluster_dist::T
    l_dist::T
    r_dist::T
    lr_weights::AbstractArray{Float32, 1}
end

mutable struct thin_suff_stats
    cluster_suff::sufficient_statistics
    l_suff::sufficient_statistics
    r_suff::sufficient_statistics
end

mutable struct local_cluster
    cluster_params::splittable_cluster_params
    total_dim::Int64
    points_count::Int64
    l_count::Int64
    r_count::Int64
end

mutable struct local_group
    model_hyperparams::model_hyper_params
    points::AbstractArray{Float32,2}
    labels::AbstractArray{Int64,1}
    labels_subcluster::AbstractArray{Int64,1}
    local_clusters::Vector{local_cluster}
    weights::Vector{Float32}
end

mutable struct pts_less_group
    model_hyperparams::model_hyper_params
    labels::AbstractArray{Int64,1}
    labels_subcluster::AbstractArray{Int64,1}
    local_clusters::Vector{local_cluster}
    weights::Vector{Float32}
end

mutable struct local_group_stats
    labels::AbstractArray{Int64,1}
    labels_subcluster::AbstractArray{Int64,1}
    local_clusters::Vector{local_cluster}
end

mutable struct dp_parallel_sampling
    model_hyperparams::model_hyper_params
    group::local_group
end

function copy_local_cluster(c::local_cluster)
    return deepcopy(c)
end

function create_pts_less_group(group::local_group)
    return pts_less_group(group.model_hyperparams , Array(group.labels), Array(group.labels_subcluster), group.local_clusters, group.weights)
end

function create_model_from_saved_data(group::pts_less_group, points::AbstractArray{Float32,2},model_hyperparams::model_hyper_params)
    group = local_group(group.model_hyperparams, points , group.labels, group.labels_subcluster, group.local_clusters, group.weights)
    return dp_parallel_sampling(model_hyperparams, group)
end


text:
################################################################################
#
#    This file contains types for a variety of empirical atomic potentials
#
################################################################################

abstract type ArbitraryPotential end
abstract type EmpiricalPotential <: ArbitraryPotential end
abstract type FittedPotential <: ArbitraryPotential end
abstract type MixedPotential <: ArbitraryPotential end

############################### Empirical Potentials ################################################
include("EmpiricalPotentials/lj.jl")
# include("EmpiricalPotentials/bm.jl")
# include("EmpiricalPotentials/coulomb.jl")
include("EmpiricalPotentials/virials.jl")
include("EmpiricalPotentials/potential_energy.jl")
include("EmpiricalPotentials/force.jl")

################################ SNAP ##############################################################
include("SNAP/snap.jl")

################################ GaN ###############################################################
# include("GaN/gan.jl")


text:
using Distributions, Plots
gr(fontfamily = "Palatino", legend = false)
scatter([Tuple(c) for c in eachcol(rand(MvNormal([30.0,100], [15.0 12; 12 13]), 100))],xlabel = "ad spending", ylabel = "revenue")

savefig("ad-spending.jl")

max_stars_repo_name: seanrattana/courses

max_issues_repo_name: seanrattana/courses

max_forks_repo_name: seanrattana/courses


text:
# test global localization by mcl, particle filter
# calculate accuracy by 1000 times test

module TestGlobalMcl
  using Distributions

  include(joinpath(split(@__FILE__, "src")[1], "src/model/world/world.jl"))
  include(joinpath(split(@__FILE__, "src")[1], "src/model/robot/differential_wheeled_robot/differential_wheeled_robot.jl"))
  include(joinpath(split(@__FILE__, "src")[1], "src/localization/global_localization/global_mcl.jl"))

  function single_test(time_interval=0.1, test_time=30; particle_num=1000)
    w = World(-5.0, 5.0, -5.0, 5.0)

    m = Map()
    add_object(m, Object(-4.0, 2.0, id=1))
    add_object(m, Object(2.0, -3.0, id=2))
    add_object(m, Object(3.0, 3.0, id=3))

    s = Sensor(m, dist_noise_rate=0.1, dir_noise=pi/90,
               dist_bias_rate_stddev=0.1, dir_bias_stddev=pi/90)

    init_pose = [rand(Uniform(-5.0, 5.0)), rand(Uniform(-5.0, 5.0)), rand(Uniform(-pi, pi))]
    e = GlobalMcl(particle_num, env_map=m)
    circling = Agent(0.2, 10.0/180*pi, estimator=e)
    r = DifferentialWheeledRobot(init_pose, 0.2, "black",
                                 circling, time_interval,
                                 noise_per_meter=5,
                                 noise_std=pi/30,
                                 bias_rate_stds=[0.1, 0.1],
                                 sensor=s)

    for t in 0:time_interval:test_time
      draw(w)
      draw!(m)
      draw!(r)
    end

    actual = r.pose
    est = r.agent.estimator.estimated_pose
    diff = sqrt((actual[1] - est[1])^2 + (actual[2] - est[2])^2)
    println("Actual:$(actual), Est:$(est), Diff:$(diff)")
    if diff <= 1.0
      return true
    else
      return false
    end
  end

  function main(test_num=1000)
    convergence_count = 0
    result = false

    for test in 1:test_num
      println("test: $(test)")
      result = single_test(particle_num=1000)
      if result == true
        convergence_count += 1
      end
    end

    println("Global Mcl Accuracy:$(convergence_count)/$(test_num)")
  end
end

max_stars_repo_name: ShisatoYano/JuliaAutonomy

max_issues_repo_name: ShisatoYano/JuliaAutonomy

max_forks_repo_name: ShisatoYano/JuliaAutonomy


text:
using GraphDatasets
using Test

include("utils.jl")

if parse(Bool, get(ENV, "CI", "false"))
    @info "CI detected: skipping some test"
else
    include("TUDatasets.jl")
end


text:
using PooksoftAlphaVantageDataStore
using JSON

# ok, so lets build by user object -
my_current_dir = pwd()   # where am I?
path_to_data_file = my_current_dir*"/test/data/IBM-Intraday-04-14-2021.json"

# load -
json_dictionary = JSON.parsefile(path_to_data_file)


text:
abstract type AbstractModel end
mutable struct Model <: AbstractModel
    laser::AbstractLaser
    waveguide::Waveguide
end
mutable struct ToyModel{TF, TG, TL, TPow, TLen, TT} <: AbstractModel
    ω0::TF
    FSR::TF
    nonlinearcoeff::TG
    linearloss::TL
    coupling::Float64
    power_in::TPow
    length::TLen
    detuning::Float64
    betacoeff::Vector
    pulsetime::TT
    has_shock::Bool
    has_raman::Bool
end
ToyModel(;ω0=200.0THz, FSR=0.1THz, nonlinearcoeff=1.0/W/m, linearloss=0.009/m,
          coupling=0.009, power_in = 0.755W, length=(628e-6)m, detuning=0.0534,
          pulsetime=0.1ps, betacoeff=[0/m, 0ps/m, -0.05ps^2/m],
          has_shock=false, has_raman=false) = ToyModel(ω0, FSR, nonlinearcoeff, linearloss,
                                              coupling, power_in, length,
                                              detuning, betacoeff, pulsetime,
                                              has_shock, has_raman)

@inline currentstructure(model::Model, z) = currentstructure(model.waveguide, z)
@inline currentstructure(model::ToyModel, z) = (1, nothing)
@inline pathlength(m::Model) = pathlength(m.waveguide)
@inline pathlength(m::ToyModel) = m.length

function derive_pulse(m::Model, t)
    #TODO: Generic injection of pulse to arbitrary mode
    pulse = derive_pulse(m.laser, t)
    full_pulse = zeros(eltype(pulse), length(t), num_modes(m))
    full_pulse[:, 1] = pulse
    return full_pulse
end
function derive_pulse(m::ToyModel, t)
    sqrtpower = sqrt(m.power_in)
    pulse_0 = m.pulsetime/(2*log(1+sqrt(2)))
    return (1+0im)*sqrtpower*sech.(t/pulse_0)
end
function derive_pulse(laser::PulsedLaser, t)
    sqrtpower = sqrt(laser.power_in)
    pulse_0 = laser.pulsetime/(2*log(1+sqrt(2)))
    return (1+0im)*sqrtpower*sech.(t/pulse_0)
end

linearloss(model::ToyModel, ω) = model.linearloss
function linearloss(model::Model, ω)
    firstmode = first(first(model.waveguide.structures).medium.modes)
    TLoss = eltype(firstmode.linearloss.(ω))
    loss = zeros(TLoss, length(ω), num_modes(model), num_structures(model))
    for (i, structure) in enumerate(model.waveguide.structures)
        for (j, mode) in enumerate(structure.medium.modes)
            loss[:, j, i] = mode.linearloss.(ω)
        end
    end
    return loss
end

nonlinearcoeff(model::ToyModel, ω) = model.nonlinearcoeff
function nonlinearcoeff(model::Model, ω)
    firstmode = first(first(model.waveguide.structures).medium.modes)
    nl_n = first(model.waveguide.structures).medium.material.electronic.nl_index
    TNL = eltype(ω .* nl_n .* firstmode.corefraction(first(ω)) / firstmode.effectivearea(first(ω)) / c0)
    # pull n2 from structure
    nlcoeff = zeros(TNL, length(ω), num_modes(model), num_structures(model))
    for (i, structure) in enumerate(model.waveguide.structures)
        nl_n = structure.medium.material.electronic.nl_index
        for (j, mode) in enumerate(structure.medium.modes)
            nlcoeff[:, j, i] = ω .* nl_n .* mode.corefraction.(ω) ./ mode.effectivearea.(ω) / c0
        end
    end
    return nlcoeff
end

function stationarydispersion(model::ToyModel, ω)
    beta = model.betacoeff
    betas = beta .* [1/factorial(n) for n=0:(length(beta)-1)]
    dispersion = Poly(betas).(ω - model.ω0)
    return dispersion
end
function stationarydispersion(model::Model, ω)
    dispersion = zeros(eltype(ω / c0), length(ω), num_modes(model), num_structures(model))
    for (i, structure) in enumerate(model.waveguide.structures)
        for (j, mode) in enumerate(structure.medium.modes)
            dispersion[:, j, i] = get_stationarydispersion(mode, model.laser).(ω)
        end
    end
    return dispersion
end

has_raman(m::ToyModel) = m.has_raman
function has_raman(model::Model)
    has_raman = fill(false, num_modes(model))
    for structure in model.waveguide.structures
        for (i, mode) in enumerate(structure.medium.modes)
            has_raman[i] = mode.has_raman
        end
    end
    return has_raman
end

@inline num_structures(m::Model) = num_structures(m.waveguide)
@inline num_modes(m::Model) = num_modes(m.waveguide)

# need to wait for Tensors.jl to merge with rotation
function raman_tensor(structure, z)
    axis_vec, rot_angle = orientation(structure, z)
    #tensor = material.raman_tensor
    tensor = structure.medium.material.raman.tensor
    #tensor = rotate(tensor, axis_vec, rot_angle)
    return tensor
end
function electronic_tensor(structure, z)
    axis_vec, rot_angle = orientation(structure, z)
    #tensor = material.electronic_tensor
    tensor = structure.medium.material.electronic.tensor
    #tensor = rotate(tensor, axis_vec, rot_angle)
    return tensor
end

@inline function raman_fraction(structure)
    return structure.medium.material.raman.raman_fraction
end

@inline interacts_with_other_mode(mode) = mode.has_interaction
@inline get_paired_mode_idx(med::Medium, mode_idx) = med.interactions[mode_idx]
@inline get_polarization(mode) = mode.polarization == :TM ? 1 : 2


text:
include("../src/Tau.jl")
using Main.Tau
VERSION < v"0.7.0-beta2.199" ? using Base.Test : using Test

@testset "self-identity" begin
    @test isa(tau, Irrational)
    @test τ == τ
    @test τ == tau
end

@testset "tau vs. 2pi" begin

    @testset "symbols" begin
        @test τ ≠ 2π # tau is Irrational, can't be equal to an AbstractFloat
        @test 2π ≠ τ
        @test π ≠ τ/2
        @test τ/2 ≠ π
    end

    @testset "ascii" begin
        @test tau ≠ 2*pi # tau is Irrational, can't be equal to an AbstractFloat
        @test 2*pi ≠ tau
        @test pi ≠ tau/2 # pi is Irrational, can't be equal to an AbstractFloat
        @test tau/2 ≠ pi
    end

    @testset "explicit type conversions" begin
        @test Float32(tau) == 2 * Float32(pi)
        @test Float64(Float32(tau)) == Float64(2 * Float32(pi))
        @test BigFloat(tau) == 2 * BigFloat(pi)
    end

end

@testset "sintau/costau" begin

    @testset "approximate results for fractional inputs" begin

        @testset "real values" begin
            for T = (Float32, Float64), x = -3:0.1:3.0
                @test @inferred(sintau(T(x))) ≈ T(sinpi(2 * parse(BigFloat, string(x))))
                @test @inferred(costau(T(x))) ≈ T(cospi(2 * parse(BigFloat, string(x))))
            end
        end

        @testset "complex values" begin
            for x in -1.0:0.1:1.0, y in -1.0:0.1:1.0
                z = complex(x, y)
                @test @inferred(sintau(z)) ≈ sinpi(2 * z)
                @test @inferred(costau(z)) ≈ cospi(2 * z)
            end
        end

    end

    @testset "exact results for integer inputs" begin

        @testset "real and complex values passed as integer types" begin
            for T = (Int, Complex), x = -3:3
                @test @inferred(sintau(T(x))) == zero(T)
                @test @inferred(costau(T(x))) == one(T)
            end
        end

        @testset "real and complex values passed as floating point types" begin
            for T = (Float32, Float64, BigFloat, Complex), x = -3.0:3.0
                @test @inferred(sintau(T(x))) == sign(x) * zero(T)
                @test @inferred(costau(T(x))) == one(T)
            end
        end

    end

    @testset "corner cases for abnormal inputs" begin

        @testset "real values" begin
            for x in (Inf, -Inf)
                @test_throws DomainError sintau(x)
                @test_throws DomainError costau(x)
            end
            @test isequal(@inferred(sintau(NaN)), sinpi(NaN))
            @test isequal(@inferred(costau(NaN)), cospi(NaN))
        end

        @testset "complex values" begin
            for x in (Inf, NaN, 0), y in (Inf, NaN, 0)
                z = complex(x, y)
                @test isequal(@inferred(sintau(z)), sinpi(2 * z))
                @test isequal(@inferred(costau(z)), cospi(2 * z))
            end
        end

    end

    # Adapted from julia/test/math.jl
    @testset "type stability" begin
        for T = (Int, Float32, Float64, BigFloat), f = (sintau, costau)
            @test Base.return_types(f, Tuple{T}) == [T]
            @test Base.return_types(f, Tuple{Complex{T}}) == [Complex{float(T)}]
        end
    end
end

# Adapted from julia/test/mod2pi.jl
@testset "modtau" begin
    @test_throws ArgumentError modtau(Int64(2)^60 - 1)
    @test modtau(10) ≈ mod(10, tau)
    @test modtau(-10) ≈ mod(-10, tau)
    @test modtau(355) ≈ 3.1416227979431572
    @test modtau(Int32(355)) ≈ 3.1416227979431572
    @test modtau(355.0) ≈ 3.1416227979431572
    @test modtau(355.0f0) ≈ 3.1416228f0
    @test modtau(Int64(2)^60) == modtau(2.0^60)
end


text:
using Core.Intrinsics: llvmcall
using Base.Threads

import Base.Threads: atomictypes, llvmtypes, inttype, ArithmeticTypes, FloatTypes,
       atomic_cas!,
       atomic_xchg!,
       atomic_add!, atomic_sub!,
       atomic_and!, atomic_nand!, atomic_or!, atomic_xor!,
       atomic_max!, atomic_min!

import Base.Sys: ARCH, WORD_SIZE

for typ ∈ atomictypes
    lt = llvmtypes[typ]
    ilt = llvmtypes[inttype(typ)]
    # Note: atomic_cas! succeeded (i.e. it stored "new") if and only if the result is "cmp"
    if typ <: Integer
        @eval atomic_cas!(x::Ptr{$typ}, cmp::$typ, new::$typ) =
            llvmcall($"""
                     %ptr = inttoptr i$WORD_SIZE %0 to $lt*
                     %rs = cmpxchg $lt* %ptr, $lt %1, $lt %2 acq_rel acquire
                     %rv = extractvalue { $lt, i1 } %rs, 0
                     ret $lt %rv
                     """, $typ, Tuple{Ptr{$typ},$typ,$typ},
                     x, cmp, new)
    else
        @eval atomic_cas!(x::Ptr{$typ}, cmp::$typ, new::$typ) =
            llvmcall($"""
                     %iptr = inttoptr i$WORD_SIZE %0 to $ilt*
                     %icmp = bitcast $lt %1 to $ilt
                     %inew = bitcast $lt %2 to $ilt
                     %irs = cmpxchg $ilt* %iptr, $ilt %icmp, $ilt %inew acq_rel acquire
                     %irv = extractvalue { $ilt, i1 } %irs, 0
                     %rv = bitcast $ilt %irv to $lt
                     ret $lt %rv
                     """, $typ, Tuple{Ptr{$typ},$typ,$typ},
                     x, cmp, new)
    end

    arithmetic_ops = [:add, :sub]
    for rmwop in [arithmetic_ops..., :xchg, :and, :nand, :or, :xor, :max, :min]
        rmw = string(rmwop)
        fn = Symbol("atomic_", rmw, "!")
        if (rmw == "max" || rmw == "min") && typ <: Unsigned
            # LLVM distinguishes signedness in the operation, not the integer type.
            rmw = "u" * rmw
        end
        if rmwop in arithmetic_ops && !(typ <: ArithmeticTypes) continue end
        if typ <: Integer
            @eval $fn(x::Ptr{$typ}, v::$typ) =
                llvmcall($"""
                         %ptr = inttoptr i$WORD_SIZE %0 to $lt*
                         %rv = atomicrmw $rmw $lt* %ptr, $lt %1 acq_rel
                         ret $lt %rv
                         """, $typ, Tuple{Ptr{$typ}, $typ}, x, v)
        else
            rmwop == :xchg || continue
            @eval $fn(x::Ptr{$typ}, v::$typ) =
                llvmcall($"""
                         %iptr = inttoptr i$WORD_SIZE %0 to $ilt*
                         %ival = bitcast $lt %1 to $ilt
                         %irv = atomicrmw $rmw $ilt* %iptr, $ilt %ival acq_rel
                         %rv = bitcast $ilt %irv to $lt
                         ret $lt %rv
                         """, $typ, Tuple{Ptr{$typ}, $typ}, x, v)
        end
    end
end

const opnames = Dict{Symbol, Symbol}(:+ => :add, :- => :sub)
for op in [:+, :-, :max, :min]
    opname = get(opnames, op, op)
    @eval function $(Symbol("atomic_", opname, "!"))(var::Ptr{T}, val::T) where T<:FloatTypes
        IT = inttype(T)
        old = unsafe_load(var)
        while true
            new = $op(old, val)
            cmp = old
            old = atomic_cas!(var, cmp, new)
            reinterpret(IT, old) == reinterpret(IT, cmp) && return new
            # Temporary solution before we have gc transition support in codegen.
            ccall(:jl_gc_safepoint, Cvoid, ())
        end
    end
end


text:
using Documenter
using RobotDynamics

makedocs(
    sitename = "RobotDynamics",
    format = Documenter.HTML(prettyurls = false),
    pages = [
        "Introduction" => "index.md",
        "Getting Started" => [
            "models.md",
            "dynamics_evaluation.md"
        ],
        "Documentation" => [
            "discretization.md",
            "knotpoints.md",
            "liemodel.md",
            "rigidbody.md",
            "linearmodel.md",
            "finite_diff.md",
            "dynamics_api.md"
        ]
    ]
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(
    repo = "github.com/RoboticExplorationLab/RobotDynamics.jl.git",
)