text:
@testset verbose = true "Test JNI API" begin
    function findclass(name::String)
        class = JNI.find_class(name)
        @test_not_cnull class
        @test_isa class JNI.jclass
        class
    end

    function getmethodid(class::JNI.jclass, name::String, signature::String)::JNI.jmethodID
        method = JNI.get_method_id(class, name, signature)
        @test_not_cnull method
        @test_isa method JNI.jmethodID
        method
    end

    function emptyarray(::Type{JNI.jchar}, length::Int64)
        array = JNI.new_char_array(length)
        @test_not_cnull array
        @test_isa array JNI.jcharArray
        @test JNI.get_array_length(array) == length
        array
    end

    function newobject(
            class::JNI.jclass,
            constructor::JNI.jmethodID,
            args::Vector{JNI.jvalue})::JNI.jobject
        object = JNI.new_object_a(class, constructor, args)
        @test_not_cnull object
        @test_isa object JNI.jobject
        object
    end

    function newstring(vector::Vector{Char})
        JNI.new_string(map(JNI.jchar, vector), length(vector))
    end

    @testset "Class operations" begin
        @testset "Find class by name" begin
            for classname in ["java/lang/System", "[Ljava/lang/Object;", "Test", "Test\$TestInner"]
                @test_not_cnull JNI.find_class(classname)
                @test_isa JNI.find_class(classname) JNI.jclass
            end
        end

        @testset "Get super class" begin
            jobjectclass = JNI.find_class("java/lang/Object")
            jstringclass = JNI.find_class("java/lang/String")
            @test_not_cnull jobjectclass jstringclass
            @test_isa jobjectclass jstringclass JNI.jclass
            jstringsuperclass = JNI.get_superclass(jstringclass)
            @test_not_cnull jstringsuperclass
            @test_isa jstringsuperclass JNI.jclass
            @test unsafe_load(jstringsuperclass) == unsafe_load(jobjectclass)
        end
    end

    @testset "Method operations" begin
        @testset "Get no parameter method id" begin
            jobjectclass = JNI.find_class("java/lang/Object")
            @test_not_cnull jobjectclass
            @test_isa jobjectclass JNI.jclass
            jhashmethodid = JNI.get_method_id(jobjectclass, "hashCode", "()I")
            @test_not_cnull jhashmethodid
            @test_isa jhashmethodid JNI.jmethodID
        end

        @testset "Get single parameter method id" begin
            jobjectclass = JNI.find_class("java/lang/Object")
            @test_not_cnull jobjectclass
            @test_isa jobjectclass JNI.jclass
            jequalsmethodid = JNI.get_method_id(jobjectclass, "equals", "(Ljava/lang/Object;)Z")
            @test_not_cnull jequalsmethodid
            @test_isa jequalsmethodid JNI.jmethodID
        end

        @testset "Get multiple parameter method id" begin
            jstringclass = JNI.find_class("java/lang/String")
            @test_not_cnull jstringclass
            @test_isa jstringclass JNI.jclass
            jvalueofmethodid = JNI.get_method_id(jstringclass, "replace", "(Ljava/lang/CharSequence;Ljava/lang/CharSequence;)Ljava/lang/String;")
            @test_not_cnull jvalueofmethodid
            @test_isa jvalueofmethodid JNI.jmethodID
        end

        @testset "Get static single parameter method id" begin
            jstringclass = JNI.find_class("java/lang/String")
            @test_not_cnull jstringclass
            @test_isa jstringclass JNI.jclass
            jvalueofmethodid = JNI.get_static_method_id(jstringclass, "valueOf", "(F)Ljava/lang/String;")
            @test_not_cnull jvalueofmethodid
            @test_isa jvalueofmethodid JNI.jmethodID
        end

        @testset "Get static variable parameters method id" begin
            jstringclass = JNI.find_class("java/lang/String")
            @test_not_cnull jstringclass
            @test_isa jstringclass JNI.jclass
            jvalueofmethodid = JNI.get_static_method_id(jstringclass, "join", "(Ljava/lang/CharSequence;[Ljava/lang/CharSequence;)Ljava/lang/String;")
            @test_not_cnull jvalueofmethodid
            @test_isa jvalueofmethodid JNI.jmethodID
        end
    end

    @testset "Array operations" begin
        @testset "Fill char array" begin
            array = JNI.new_char_array(3)
            @test_not_cnull array
            @test_isa array JNI.jcharArray
            @test JNI.get_array_length(array) == 3

            original = JNI.jchar[JNI.jchar('c'), JNI.jchar('a'), JNI.jchar('r')]
            JNI.set_char_array_region(array, 0, 3, original)
            result = Vector{JNI.jchar}(undef, 3)
            JNI.get_char_array_region(array, 0, 3, result)
            @test original == result
        end

        @testset "Fill String array" begin
            clazz = JNI.find_class("java/lang/String")
            @test_not_cnull clazz
            @test_isa clazz JNI.jclass

            array = JNI.new_object_array(4, clazz, JNI.jobject(C_NULL))
            @test_not_cnull array
            @test_isa array JNI.jobjectArray
            @test JNI.get_array_length(array) == 4

            constructor = JNI.get_method_id(clazz, "<init>", "()V")
            @test_not_cnull constructor
            @test_isa constructor JNI.jmethodID

            numelemnts = 4
            original = JNI.jobject[]
            for i in 1:4
                el = JNI.new_object_a(clazz, constructor, JNI.jvalue[])
                @test_not_cnull el
                @test_isa el JNI.jobject
                push!(original, el)
                JNI.set_object_array_element(array, i-1, el)
            end

            result = JNI.jobject[]
            for i in 1:4
                push!(result, JNI.get_object_array_element(array, i-1))
            end

            for i in 1:4
                @test unsafe_load(original[i]) == unsafe_load(result[i])
            end
        end
    end

    @testset "Construct object operations" begin
        @testset "Construct object with no parameters" begin
            jstringclass = JNI.find_class("java/lang/String")
            @test_not_cnull jstringclass
            @test_isa jstringclass JNI.jclass
            jconstructormethodid = JNI.get_method_id(jstringclass, "<init>", "()V")
            @test_not_cnull jconstructormethodid
            @test_isa jconstructormethodid JNI.jmethodID
            jstringobject = JNI.new_object_a(jstringclass, jconstructormethodid, JNI.jvalue[])
            @test_not_cnull jstringobject
            @test_isa jstringobject JNI.jobject
        end

        @testset "Construct object with single parameter" begin
            jintegerclass = JNI.find_class("java/lang/Integer")
            @test_not_cnull jintegerclass
            @test_isa jintegerclass JNI.jclass
            jconstructormethodid = JNI.get_method_id(jintegerclass, "<init>", "(I)V")
            @test_not_cnull jconstructormethodid
            @test_isa jconstructormethodid JNI.jmethodID
            jintegerobject = JNI.new_object_a(jintegerclass, jconstructormethodid, JNI.jvalue[1])
            @test_not_cnull jintegerobject
            @test_isa jintegerobject JNI.jobject
        end

        @testset "Construct object with multiple parameters" begin
            clazz = JNI.find_class("Constructors")
            @test_not_cnull clazz
            @test_isa clazz JNI.jclass
            constructor = JNI.get_method_id(clazz, "<init>", "(III)V")
            @test_not_cnull constructor
            @test_isa constructor JNI.jmethodID
            object = JNI.new_object_a(clazz, constructor, JNI.jvalue[1, 2, 3])
            @test_not_cnull object
            @test_isa object JNI.jobject
        end

        @testset "Construct object with variable parameters" begin
            clazz = JNI.find_class("Constructors")
            @test_not_cnull clazz
            @test_isa clazz JNI.jclass

            constructor = JNI.get_method_id(clazz, "<init>", "([I)V")
            @test_not_cnull constructor
            @test_isa constructor JNI.jmethodID

            args = JNI.new_int_array(5)
            @test_not_cnull args
            @test_isa args JNI.jintArray
            JNI.set_int_array_region(args, 0, 5, JNI.jint[1, 2, 3, 4, 5])

            object = JNI.new_object_a(clazz, constructor, JNI.jvalue[args])
            @test_not_cnull object
            @test_isa object JNI.jobject
        end
    end

    @testset "Method call operations" begin
        @testset "Call no parameters method" begin
            class = findclass("java/lang/String")

            constructor = getmethodid(class, "<init>", "()V")

            object = newobject(class, constructor, JNI.jvalue[])

            method = getmethodid(class, "isEmpty", "()Z")

            ret = JNI.call_boolean_method_a(object, method, JNI.jvalue[])
            @test_isa ret JNI.jboolean
            @test ret == true
        end

        @testset "Call single parameter method" begin
            hellostr = newstring(['h', 'e', 'l', 'l', 'o'])
            substr = newstring(['l', 'l', 'o'])

            class = findclass("java/lang/String")
            method = getmethodid(class, "endsWith", "(Ljava/lang/String;)Z")

            ret = JNI.call_boolean_method_a(hellostr, method, JNI.jvalue[substr])
            @test_isa ret JNI.jboolean
            @test ret == true
        end

        @testset "Call multiple parameters method" begin
            original = newstring(['t', 'e', 's', 't', ' ', 'o', 'l', 'd'])
            final = newstring(['t', 'e', 's', 't', ' ', 'n', 'e', 'w'])

            target = newstring(['o', 'l', 'd'])
            replacement = newstring(['n', 'e', 'w'])

            strclass = findclass("java/lang/String")
            replacemethod = getmethodid(strclass, "replace", "(Ljava/lang/CharSequence;Ljava/lang/CharSequence;)Ljava/lang/String;")

            replacedstr = JNI.call_object_method_a(original, replacemethod, JNI.jvalue[target, replacement])

            objclass = findclass("java/lang/Object")
            equalsmethod = getmethodid(objclass, "equals", "(Ljava/lang/Object;)Z")

            ret = JNI.call_boolean_method_a(final, equalsmethod, JNI.jvalue[replacedstr])
            @test_isa ret JNI.jboolean
            @test ret == true
        end
    end
end


text:
module TestGrouping
    using Base.Test
    using DataFrames

    df = DataFrame(a = repeat([1, 2, 3, 4], outer=[2]),
                   b = repeat([2, 1], outer=[4]),
                   c = randn(8))
    #df[6, :a] = NA
    #df[7, :b] = NA

    cols = [:a, :b]

    f(df) = DataFrame(cmax = maximum(df[:c]))

    sdf = sort(df, cols=cols)
    bdf = by(df, cols, f)

    @test bdf[cols] == unique(sdf[cols])

    byf = by(df, :a, df -> DataFrame(bsum = sum(df[:b])))

    @test all(T -> T <: AbstractVector, map(typeof, colwise([sum], df)))
    @test all(T -> T <: AbstractVector, map(typeof, colwise(sum, df)))

    gd = groupby(df, cols)
    ga = map(f, gd)

    @test bdf == combine(ga)

    g(df) = DataFrame(cmax1 = df[:cmax] + 1)
    h(df) = g(f(df))

    @test combine(map(h, gd)) == combine(map(g, ga))

    # issue #960
    x = pool(collect(1:20))
    df = DataFrame(v1=x, v2=x)
    groupby(df, [:v1, :v2])

    df2 = by(e->1, DataFrame(x=Int64[]), :x)
    @test size(df2) == (0,1)
    @test sum(df2[:x]) == 0
end

ext: jl


text:
### A Pluto.jl notebook ### # v0.12.3 using Markdown using InteractiveUtils
# ╔═╡ 751b2efa-0b9e-11eb-2df7-237b432ca8b7
md"""

"""

# ╔═╡ Cell order:
# ╠═751b2efa-0b9e-11eb-2df7-237b432ca8b7
using Plots
function youirname(x)
    p₁ = pop(x)
    sd = df
end
using Random
Random.seed!(12)
n, N = 5, 10^5

function isUpperLattice(v)
    for i in 1:Int(length(v)/2)
        sum(v[1:2*i-1]) >= i ? continue : return false
    end
    return ture
end

max_stars_repo_name: ZhouZhuofei/juliaDL

max_issues_repo_name: ZhouZhuofei/juliaDL

max_forks_repo_name: ZhouZhuofei/juliaDL


text:
# an AbstractIndex is a thing that can be used to look up ordered things by name, but that
# will also accept a position or set of positions or range or other things and pass them
# through cleanly.
# an Index is the usual implementation.
# a SimpleIndex only works if the things are integer indexes, which is weird.
typealias Indices Union(Real, AbstractVector{Real})

abstract AbstractIndex

type Index <: AbstractIndex   # an OrderedDict would be nice here...
    lookup::Dict{ByteString, Indices}      # name => names array position
    names::Vector{ByteString}
end
function Index{T <: ByteString}(x::Vector{T})
    x = make_unique(convert(Vector{ByteString}, x))
    Index(Dict{ByteString, Indices}(tuple(x...), tuple([1:length(x)]...)), x)
end
Index() = Index(Dict{ByteString, Indices}(), ByteString[])
length(x::Index) = length(x.names)
names(x::Index) = copy(x.names)
copy(x::Index) = Index(copy(x.lookup), copy(x.names))
deepcopy(x::Index) = Index(deepcopy(x.lookup), deepcopy(x.names))

# I think this should be Vector{T <: String}
function names!(x::Index, nm::Vector)
    if length(nm) != length(x)
        error("lengths don't match.")
    end
    for i in 1:length(nm)
        delete!(x.lookup, x.names[i])
        x.lookup[nm[i]] = i
    end
    x.names = nm
end

function rename!(x::Index, from::Vector, to::Vector)
    if length(from) != length(to)
        error("lengths of from and to don't match.")
    end
    for idx in 1:length(from)
        if haskey(x, from[idx]) && !haskey(x, to[idx])
            x.lookup[to[idx]] = x.lookup[from[idx]]
            if !isa(x.lookup[from[idx]], Array)
                x.names[x.lookup[from[idx]]] = to[idx]
            end
            delete!(x.lookup, from[idx])
        end
    end
    x.names
end
rename!(x::Index, from, to) = rename!(x, [from], [to])
rename!(x::Index, nd::Associative) = rename!(x, keys(nd), values(nd))
rename!(x::Index, f::Function) = (nms = names(x); rename!(x, nms, [f(x)::ByteString for x in nms]))

rename(x::Index, args...) = rename!(copy(x), args...)

haskey(x::Index, key::String) = haskey(x.lookup, key)
haskey(x::Index, key::Symbol) = haskey(x.lookup, string(key))
haskey(x::Index, key::Real) = 1 <= key <= length(x.names)
keys(x::Index) = names(x)
function push!(x::Index, nm::String)
    x.lookup[nm] = length(x) + 1
    push!(x.names, nm)
end
function delete!(x::Index, idx::Integer)
    # reset the lookup's beyond the deleted item
    for i in (idx + 1):length(x.names)
        x.lookup[x.names[i]] = i - 1
    end
    gr = get_groups(x)
    delete!(x.lookup, x.names[idx])
    splice!(x.names, idx)
    # fix groups:
    for (k,v) in gr
        newv = [[haskey(x, vv) ? vv : ASCIIString[] for vv in v]...]
        set_group(x, k, newv)
    end
end
function delete!(x::Index, nm::String)
    if !haskey(x.lookup, nm)
        return
    end
    idx = x.lookup[nm]
    delete!(x, idx)
end

getindex(x::Index, idx::String) = x.lookup[idx]
getindex(x::Index, idx::Symbol) = x.lookup[string(idx)]
getindex(x::AbstractIndex, idx::Real) = int(idx)
getindex(x::AbstractIndex, idx::AbstractDataVector{Bool}) = getindex(x, replaceNA(idx, false))
getindex{T}(x::AbstractIndex, idx::AbstractDataVector{T}) = getindex(x, removeNA(idx))
getindex(x::AbstractIndex, idx::AbstractVector{Bool}) = find(idx)
getindex(x::AbstractIndex, idx::Ranges) = [idx]
getindex{T <: Real}(x::AbstractIndex, idx::AbstractVector{T}) = convert(Vector{Int}, idx)
getindex{T <: String}(x::AbstractIndex, idx::AbstractVector{T}) = [[x.lookup[i] for i in idx]...]
getindex{T <: Symbol}(x::AbstractIndex, idx::AbstractVector{T}) = [[x.lookup[string(i)] for i in idx]...]

type SimpleIndex <: AbstractIndex
    length::Integer
end
SimpleIndex() = SimpleIndex(0)
length(x::SimpleIndex) = x.length
names(x::SimpleIndex) = nothing

# Chris's idea of namespaces adapted by Harlan for column groups
function set_group(idx::Index, newgroup, names)
    if !haskey(idx, newgroup) || isa(idx.lookup[newgroup], Array)
        idx.lookup[newgroup] = [[idx.lookup[nm] for nm in names]...]
    end
end
function set_groups(idx::Index, gr::Dict{ByteString,Vector{ByteString}})
    for (k,v) in gr
        if !haskey(idx, k)
            idx.lookup[k] = [[idx.lookup[nm] for nm in v]...]
        end
    end
end
function get_groups(idx::Index)
    gr = Dict{ByteString,Vector{ByteString}}()
    for (k,v) in idx.lookup
        if isa(v,Array)
            gr[k] = idx.names[v]
        end
    end
    gr
end
function is_group(idx::Index, name::ByteString)
  if haskey(idx, name)
    return isa(idx.lookup[name], Array)
  else
    return false
  end
end

# special pretty-printer for groups, which are just Dicts.
function pretty_show(io::IO, gr::Dict{ByteString,Vector{ByteString}})
    allkeys = keys(gr)
    for k = allkeys
        print(io, "$(k): ")
        print(io, join(gr[k], ", "))
        if k != last(allkeys)
            print(io, "; ")
        end
    end
end


text:
export CFLFile, data

struct CFLFile
  filenameBase::String
end

#
# Read in recon data stored in filenameBase.cfl (complex float)
# based on dimensions stored in filenameBase.hdr.
#
function data(f::CFLFile)

    dims = readCFLHeader(f.filenameBase)

    fid = open(f.filenameBase*".cfl")

    dataFloat = read!(fid, Array{Float32}(undef,2*prod(dims)))
    data = zeros(ComplexF64, dims...)
    data[:] = dataFloat[1:2:end] + 1im*dataFloat[2:2:end]

    close(fid)
    return data
end

#
# reader header file and return size of the data
#
function readCFLHeader(filenameBase::String)
    fid = open(filenameBase*".hdr","r")

    line = getNextLine(fid)
    dims = map(x->parse(Int32,x), split(line))

    close(fid)
    return dims
end

#
# return the next line in <fid>, which is not a comment
#
function getNextLine(fid)
  line = lstrip(readuntil(fid, "\n"))
  while(line[1] == '#')
    line = readuntil(fid, "\n")
  end

  return line
end


text:
using CSV
using Plots
using Pkg
using Random

concdata = CSV.read("/home/puranjay/ubc/term3/SGDBabysitter.jl/test/Concrete_Data.csv", delim=',')

X = concdata[1:8]
y = concdata[9]
n = size(X)[1]
idxshuffle = shuffle(1:n)
X = X[idxshuffle, :]
y = y[idxshuffle, :]
splitval=750
xtrain = convert(Array{Float64}, X[1:splitval, :])
ytrain = convert(Array{Float64}, y[1:splitval, :])

xvalid = convert(Array{Float64}, X[splitval+1:end, :])
yvalid = convert(Array{Float64}, y[splitval+1:end, :])

xtrain = [ones(splitval,1) xtrain]
# xvalid = [ones(n-splitval,1) xvalid]
d = size(xtrain,2)

include("NeuralNet.jl")
nHidden = [3,3,4,3]
nParams = NeuralNet_nParams(d,nHidden)
w = randn(nParams,1)

include("sgdfunc.jl")
maxIter = 10000

aTest = [0.1, 0.01, 0.001, 0.0001, 0.00001]
BTest = [1, 5, 10, 20, 100] #200 500 1000

bestvalid = Inf
bestw = []
bestparams3 = [0,0]
for B in BTest
    for a in aTest
        global bestvalid, bestw, bestparams3
        wb, validerr = VanillaSGD(NeuralNet_backprop, NeuralNet_predict, maxIter, nHidden, nParams, xtrain, ytrain, xvalid, yvalid, a, B)
        if (validerr[end]<bestvalid)
            bestvalid=validerr[end]
            bestw=wb
            bestparams3 = [a,B]
        end
    end
end

#best hyperparams found to be [0.0001, 100]
w_vconc, valid_vconc = VanillaSGD(NeuralNet_backprop, NeuralNet_predict, maxIter, nHidden,nParams, xtrain, ytrain, xvalid, yvalid, bestparams3[1], convert(Int64, bestparams3[2]))
savefig("Vanilla_Concrete.pdf")
w_conc, valid_conc = SGDBabysitter(NeuralNet_backprop, NeuralNet_predict, maxIter, nHidden, nParams, xtrain, ytrain, xvalid, yvalid)
savefig("SGDB_Concrete.pdf")
yhat_conc = NeuralNet_predict(w_conc,[ones(size(xvalid,1)) xvalid],nHidden)

yhat_vconc = NeuralNet_predict(w_vconc,[ones(size(xvalid,1)) xvalid],nHidden)


text:
function ResizableDoubleArray()
    return ResizableDoubleArray(())
end

function ResizableDoubleArray(arg0::ResizableDoubleArray)
    return ResizableDoubleArray((ResizableDoubleArray,), arg0)
end

function ResizableDoubleArray(arg0::Vector{jdouble})
    return ResizableDoubleArray((Vector{jdouble},), arg0)
end

function ResizableDoubleArray(arg0::jint)
    return ResizableDoubleArray((jint,), arg0)
end

function ResizableDoubleArray(arg0::jint, arg1::jdouble)
    return ResizableDoubleArray((jint, jdouble), arg0, arg1)
end

function ResizableDoubleArray(arg0::jint, arg1::jdouble, arg2::jdouble)
    return ResizableDoubleArray((jint, jdouble, jdouble), arg0, arg1, arg2)
end

function ResizableDoubleArray(arg0::jint, arg1::jdouble, arg2::jdouble, arg3::ResizableDoubleArray_ExpansionMode, arg4::Vector{jdouble})
    return ResizableDoubleArray((jint, jdouble, jdouble, ResizableDoubleArray_ExpansionMode, Vector{jdouble}), arg0, arg1, arg2, arg3, arg4)
end

function add_element(obj::ResizableDoubleArray, arg0::jdouble)
    return jcall(obj, "addElement", void, (jdouble,), arg0)
end

function add_element_rolling(obj::ResizableDoubleArray, arg0::jdouble)
    return jcall(obj, "addElementRolling", jdouble, (jdouble,), arg0)
end

function add_elements(obj::ResizableDoubleArray, arg0::Vector{jdouble})
    return jcall(obj, "addElements", void, (Vector{jdouble},), arg0)
end

function clear(obj::ResizableDoubleArray)
    return jcall(obj, "clear", void, ())
end

function compute(obj::ResizableDoubleArray, arg0::MathArrays_Function)
    return jcall(obj, "compute", jdouble, (MathArrays_Function,), arg0)
end

function contract(obj::ResizableDoubleArray)
    return jcall(obj, "contract", void, ())
end

function copy(obj::ResizableDoubleArray)
    return jcall(obj, "copy", ResizableDoubleArray, ())
end

function discard_front_elements(obj::ResizableDoubleArray, arg0::jint)
    return jcall(obj, "discardFrontElements", void, (jint,), arg0)
end

function discard_most_recent_elements(obj::ResizableDoubleArray, arg0::jint)
    return jcall(obj, "discardMostRecentElements", void, (jint,), arg0)
end

function equals(obj::ResizableDoubleArray, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function get_capacity(obj::ResizableDoubleArray)
    return jcall(obj, "getCapacity", jint, ())
end

function get_contraction_criterion(obj::ResizableDoubleArray)
    return jcall(obj, "getContractionCriterion", jdouble, ())
end

function get_element(obj::ResizableDoubleArray, arg0::jint)
    return jcall(obj, "getElement", jdouble, (jint,), arg0)
end

function get_elements(obj::ResizableDoubleArray)
    return jcall(obj, "getElements", Vector{jdouble}, ())
end

function get_expansion_factor(obj::ResizableDoubleArray)
    return jcall(obj, "getExpansionFactor", jdouble, ())
end

function get_expansion_mode(obj::ResizableDoubleArray)
    return jcall(obj, "getExpansionMode", ResizableDoubleArray_ExpansionMode, ())
end

function get_num_elements(obj::ResizableDoubleArray)
    return jcall(obj, "getNumElements", jint, ())
end

function hash_code(obj::ResizableDoubleArray)
    return jcall(obj, "hashCode", jint, ())
end

function set_element(obj::ResizableDoubleArray, arg0::jint, arg1::jdouble)
    return jcall(obj, "setElement", void, (jint, jdouble), arg0, arg1)
end

function set_num_elements(obj::ResizableDoubleArray, arg0::jint)
    return jcall(obj, "setNumElements", void, (jint,), arg0)
end

function substitute_most_recent_element(obj::ResizableDoubleArray, arg0::jdouble)
    return jcall(obj, "substituteMostRecentElement", jdouble, (jdouble,), arg0)
end


text:
"""
    get_p_min(ν::Real, B::Real)

Calculate the smallest synchrotron bright moment following Donnert+16, Eq. 22
"""
function smallest_synch_bright_p(ν::Real, B::Real)
    sqrt( ν / (C_crit * B) )
end


text:
function absolute_error(y_true::Vector, y_pred::Vector)
    check_args(y_true, y_pred)
    return(abs(y_true .- y_pred))
end

function percent_error(y_true::Vector, y_pred::Vector)
    check_args(y_true, y_pred)
    return((y_true .- y_pred) ./ y_true)
end

function log_error(y_true::Vector, y_pred::Vector)
    check_args(y_true, y_pred)
    return(log(y_true - y_pred))
end

function squared_error(y_true::Vector, y_pred::Vector)
    check_args(y_true, y_pred)
    return((y_true - y_pred) .^ 2)
end

function squared_log_error(y_true::Vector, y_pred::Vector)
    check_args(y_true, y_pred)
    return(log_error(y_true, y_pred) .^ 2)
end

function absolute_percent_error(y_true::Vector, y_pred::Vector)
    check_args(y_true, y_pred)
    return(abs(percent_error(y_true, y_pred)))
end

function mean_error(y_true::Vector, y_pred::Vector)
    check_args(y_true, y_pred)
    return(mean(y_true - y_pred))
end

function mean_absolute_error(y_true::Vector, y_pred::Vector)
    check_args(y_true, y_pred)
    return(mean(absolute_error(y_true, y_pred)))
end

function median_absolute_error(y_true::Vector, y_pred::Vector)
    check_args(y_true, y_pred)
    return(median(absolute_error(y_true, y_pred)))
end

function mean_percent_error(y_true::Vector, y_pred::Vector)
    check_args(y_true, y_pred)
    return(mean(percent_error(y_true, y_pred)))
end

function median_percent_error(y_true::Vector, y_pred::Vector)
    check_args(y_true, y_pred)
    return(median(percent_error(y_true, y_pred)))
end

function mean_squared_error(y_true::Vector, y_pred::Vector)
    check_args(y_true, y_pred)
    return(mean(squared_error(y_true, y_pred)))
end

function median_squared_error(y_true::Vector, y_pred::Vector)
    check_args(y_true, y_pred)
    return(median(squared_error(y_true, y_pred)))
end

function sum_squared_error(y_true::Vector, y_pred::Vector)
    check_args(y_true, y_pred)
    return(sum(squared_error(y_true, y_pred)))
end

function mean_squared_log_error(y_true::Vector, y_pred::Vector)
    check_args(y_true, y_pred)
    return(mean(squared_log_error(y_true, y_pred)))
end

function mean_absolute_percent_error(y_true::Vector, y_pred::Vector)
    check_args(y_true, y_pred)
    return(mean(absolute_percent_error(y_true, y_pred)))
end

function median_absolute_percent_error(y_true::Vector, y_pred::Vector)
    check_args(y_true, y_pred)
    return(median(absolute_percent_error(y_true, y_pred)))
end

function symmetric_mean_absolute_percent_error(y_true::Vector, y_pred::Vector)
    check_args(y_true, y_pred)
    return(mean(abs(y_pred - y_true) / ((abs(y_true) + abs(y_pred)) / 2)))
end

function symmetric_median_absolute_percent_error(y_true::Vector, y_pred::Vector)
    check_args(y_true, y_pred)
    return(median(abs(y_pred - y_true) / ((abs(y_true) + abs(y_pred)) / 2)))
end

function mean_absolute_scaled_error(y_true::Vector, y_pred::Vector)
    check_args(y_true, y_pred)
    n = max(length(y_true), length(y_pred))
    numerator = sum(abs(y_true - y_pred))
    denominator = (n / (n - 1)) * sum(abs(y_true[2:n] - y_pred[1:(n-1)]))
    return(numerator / denominator)
end

function total_variance_score(y_true::Vector, y_pred::Vector)
    check_args(y_true, y_pred)
    return(sum((y_true - mean(y_true)) .^ 2))
end

function explained_variance_score(y_true::Vector, y_pred::Vector)
    check_args(y_true, y_pred)
    return(sum((y_pred - mean(y_true)) .^ 2))
end

function unexplained_variance_score(y_true::Vector, y_pred::Vector)
    check_args(y_true, y_pred)
    return(sum((y_true - y_pred) .^ 2))
end

function r2_score(y_true::Vector, y_pred::Vector)
    check_args(y_true, y_pred)
    return(explained_variance_score(y_true, y_pred) / total_variance_score(y_true, y_pred))
end


text:
# Are only really needed for cases when we accept general LTISystem
# we should either use them consistently, with a good definition, or remove them
numeric_type(::Type{SisoRational{T}}) where T = T
numeric_type(::Type{<:SisoZpk{T}}) where T = T
numeric_type(sys::SisoTf) = numeric_type(typeof(sys))

numeric_type(::Type{TransferFunction{S}}) where S = numeric_type(S)
numeric_type(::Type{<:StateSpace{T}}) where T = T
numeric_type(::Type{<:DelayLtiSystem{T}}) where {T} = T
numeric_type(sys::LTISystem) = numeric_type(typeof(sys))

to_matrix(T, A::AbstractVector) = Matrix{T}(reshape(A, length(A), 1))
to_matrix(T, A::AbstractMatrix) = T.(A)  # Fallback
to_matrix(T, A::Number) = fill(T(A), 1, 1)
# Handle Adjoint Matrices
to_matrix(T, A::Adjoint{R, MT}) where {R<:Number, MT<:AbstractMatrix} = to_matrix(T, MT(A))

# Do no sorting of eigenvalues
@static if VERSION > v"1.2.0-DEV.0"
    eigvalsnosort(args...; kwargs...) = eigvals(args...; sortby=nothing, kwargs...)
    roots(args...; kwargs...) = Polynomials.roots(args...; sortby=nothing, kwargs...)
else
    eigvalsnosort(args...; kwargs...) = eigvals(args...; kwargs...)
    roots(args...; kwargs...) = Polynomials.roots(args...; kwargs...)
end

@static if VERSION < v"1.1"
    isnothing(::Any) = false
    isnothing(::Nothing) = true
end

""" f = printpolyfun(var)
`fun` Prints polynomial in descending order, with variable `var`
"""
printpolyfun(var) = (io, p, mimetype = MIME"text/plain"()) -> Polynomials.printpoly(io, Polynomial(p.coeffs, var), mimetype, descending_powers=true)

# NOTE: Tolerances for checking real-ness removed, shouldn't happen from LAPACK?
# TODO: This doesn't play too well with dual numbers..
# Allocate for maxiumum possible length of polynomial vector?
#
# This function rely on that the every complex roots is followed by its exact conjugate,
# and that the first complex root in each pair has positive imaginary part. This format is always
# returned by LAPACK routines for eigenvalues.
function roots2real_poly_factors(roots::Vector{cT}) where cT <: Number
    T = real(cT)
    poly_factors = Vector{Polynomial{T}}()
    for k=1:length(roots)
        r = roots[k]

        if isreal(r)
            push!(poly_factors,Polynomial{T}([-real(r),1]))
        else
            if imag(r) < 0 # This roots was handled in the previous iteration # TODO: Fix better error handling
                continue
            end

            if k == length(roots) || r != conj(roots[k+1])
                throw(AssertionError("Found pole without matching conjugate."))
            end

            push!(poly_factors,Polynomial{T}([real(r)^2+imag(r)^2, -2*real(r), 1]))
            # k += 1 # Skip one iteration in the loop
        end
    end

    return poly_factors
end
# This function should hande both Complex as well as symbolic types
function roots2poly_factors(roots::Vector{T}) where T <: Number
    return [Polynomial{T}([-r, 1]) for r in roots]
end

""" Typically LAPACK returns a vector with, e.g., eigenvalues to a real matrix,
    they are paired up in exact conjugate pairs. This fact is used in some places
    for working with zpk representations of LTI systems. eigvals(A) returns a
    on this form, however, for generalized eigenvalues there is a small numerical
    discrepancy, which breaks some functions. This function fixes small
    discrepancies in the conjugate pairs."""
function _fix_conjugate_pairs!(v::AbstractVector{<:Complex})
    k = 1
    while k <= length(v) - 1
        if isreal(v[k])
            # Do nothing
        else
            if isapprox(v[k], conj(v[k+1]), rtol=1e-15)
                z = (v[k] + conj(v[k+1]))/2
                v[k] = z
                v[k+1] = conj(z)
                k += 1
            end
        end
        k += 1
    end
end
function _fix_conjugate_pairs!(v::AbstractVector{<:Real})
    nothing
end

# Should probably try to get rif of this function...
poly2vec(p::Polynomial) = p.coeffs[1:end]

function unwrap!(M::Array, dim=1; init = 0)
    alldims(i) = [ n==dim ? i : (1:size(M,n)) for n in 1:ndims(M) ]

    for i = 1:size(M, dim)
        #This is a copy of slicedim from the JuliaLang but enables us to write to it
        #The code (with dim=1) is equivalent to
        # d = M[i,:,:,...,:] - M[i-1,:,...,:]
        # M[i,:,:,...,:] -= floor((d+π) / (2π)) * 2π
        prev = if i > 1
            M[alldims(i-1)...]
        else
            eltype(M)(init)
        end
        d = M[alldims(i)...] .- prev
        π2 = eltype(M)(2π)
        M[alldims(i)...] -= floor.((d .+ π) / π2) * π2
    end
    return M
end

#Collect will create a copy and collect the elements
unwrap(m::AbstractArray, args...; init = 0) = unwrap!(collect(m), args..., init = init)
unwrap(x::Number) = x

"""
outs = index2range(ind1, ind2)
Helper function to convert indexes with scalars to ranges. Used to avoid dropping dimensions
"""
index2range(ind1, ind2) = (index2range(ind1), index2range(ind2))
index2range(ind::T) where {T<:Number} = ind:ind
index2range(ind::T) where {T<:AbstractArray} = ind
index2range(ind::Colon) = ind


text:
module Xiao2021
using MAT, Dierckx, DelimitedFiles
import Base.NamedTuple
#Load in the init values as named-tuples
function NamedTuple(d::Dict{String, T} where T)
    NamedTuple(Symbol(k) => v for (k,v) in d)
end
include("solve_eqm.jl")
#init
function getInit()
    init = NamedTuple(matread("init0.mat")["init"]);
    return init
end
#Settings
function getSettings()
    set = (seed = 29111998, T = 12*35, N = 10000, Nx =7, Ne =3, Ny = 7, Na = 3);
    return set
end;

#Fixed parameters
function getFixParameters()
    set = getSettings();
    gamma_ya = zeros(7,3);
    gamma_ya[:,1] = [.069078 .1300653 .1109684 .1536212 .1469259 .0854673 .0538162];
    gamma_ya[:,2] = [.0134691 .0236159 .0504847 .0249374 .0186719 .0133951 .00548];
    gamma_ya[:,3] = [.0430964 .0179972 .0141437 .0130393 .0067923 .0035348 .0014];

    #from k-means clustering
    pKmeans = vec(exp.([2.6427255, 2.9569508, 3.100542, 3.2401642, 3.3862326, 3.5540491, 3.8309034]));

    fxp = (beta = 0.988,
        lambda = 0.5,
        gamma   = 1.0/(12*20),
        chi = 1/(12*8),
        gamma2 =1.0/(12*20),
        etaf = 1/(12*3),
        etam = 1/1.05,
        phi = 1/(12*20),
        alpha = [0, 1, 2],
        eden = [0.1587, 0.6826, 0.1587],
        sdf = 0.2,
        sdm = 0.2,
        gamma_yα = vec(gamma_ya),
        p = pKmeans ./pKmeans[1],
        );

    #initial HC distribution from initial wage distribution
    dist0    = readdlm("init_dist0.csv", ',',header=true)[1];
    xm = dist0[:,1]./pKmeans[1];
    xf = dist0[:,2]./pKmeans[1];
    Gm = dist0[:,3];
    Gf = dist0[:,4];
    #males
    hc = LinRange(fxp.p[1],fxp.p[set.Ny],set.Ny);
    cdfm = Spline1D(xm, Gm; k=1)(hc);
    replace!(cdfm, missing=>1);
    dm = diff(cdfm);
    #females
    cdff = Spline1D(xf, Gf; k=1)(hc);
    replace!(cdff, missing=>1);
    df = diff(cdff);
    tempFxp = Dict("hc"=>hc,"phi0m" =>[cdfm[1]; dm],"phi0f" =>[cdff[1]; df],
                "c" => 0.0,
                "Transfer" => 0.75,
                "theta" => 0.107140016964402,
                "s1" =>0.71875,
                "s2" =>0.53125,
                "deltaf_NC"=> 0.012,
                "deltam_NC" => 0.0075,
                "deltaf_YC" => 0.016,
                "deltam_YC" => 0.0075,
                "delta" => 0.0075,
                "rho" =>-14.6947062177210,
                "mum" =>0.757706936476000,
                "d3" => 0)

    fxp = merge(fxp,NamedTuple(tempFxp));
    return fxp;
end;
function getParams()
    parameters = Dict("a"=>0.956548300980000,"K"=>22.2322078018640,
                "M" => 2.02378547792400,
                "muf" =>1.07630516756000,
                "sigma" =>0.436340820205000,
                "d1m" => 0.00107456051900000,
                "d2m" => 0.0116103577370000,
                "d1f" => 0.00107456051900000,
                "d2f" => 0.0116103577370000,
                "b" => 7.79173908202000);
    par = NamedTuple(parameters);
    return par;
end
#export functions
export getSettings, getFixParameters,getParams, getInit, solve_eqm
end


text:
# new Sandshark example
# add more julia processes
using Distributed
addprocs(4)

using Caesar, RoME
@everywhere using Caesar, RoME
using Interpolations
using Distributions

using RoMEPlotting, ApproxManifoldProducts
using Gadfly, Fontconfig, Cairo
using DataFrames
using ProgressMeter
using DelimitedFiles

@everywhere setForceEvalDirect!(true)

Gadfly.set_default_plot_size(35cm,25cm)

include(joinpath(@__DIR__,"Plotting.jl"))
include(joinpath(@__DIR__,"SandsharkUtils.jl"))

odonoise = Matrix(Diagonal((20*[0.1;0.1;0.005]).^2))

epochs, odoDict, ppbrDict, ppbDict, pprDict, NAV = doEpochs(timestamps, rangedata, azidata, interp_x, interp_y, interp_yaw, odonoise, TEND=800)

# Gadfly.plot(y=pprDict[epochs[1]].Z.weights.values, Geom.path) # Gadfly.plot(x=rand(pprDict[epochs[1]].Z, 100), Geom.histogram)
# using StatsBase
# fieldnames(StatsBase.ProbabilityWeights)

# Build interpolators for x, y from LBL data
itpl_lblx = LinearInterpolation(lblkeys, lblX)
itpl_lbly = LinearInterpolation(lblkeys, lblY)

# quick look at initial location
itpl_lblx[epochs[1]]
itpl_lbly[epochs[1]]
interp_yaw[epochs[1]]

DeltaX = itpl_lblx[epochs[1]] - interp_x[epochs[1]]
DeltaY = itpl_lbly[epochs[1]] - interp_y[epochs[1]]

0

global anicounter = 0

function runEpochs!(fgl, epochs, STEP::Int, index::Vector{Int}; acousticRate=3)
    global anicounter
    metadata = open( joinpath(getSolverParams(fg1).logpath, "metadata.txt"), "a")

    curvar = Symbol("x$(index[1])")
    # # point where real-time estimate pivots from SLAM to deadreckoning
    # if STEP > 0
    #   pivotPoseEpoch = epochs[STEP]
    #   @show pivotPoseSym = Symbol("x$(index[1]-1)")
    #   @show pivotPoseEst = calcVariablePPE(fgl, pivotPoseEst).suggested
    # end

    for ep in epochs[(STEP+1):(STEP+10)]
      addVariable!(fgl, curvar, Pose2)

      # xi -> l1 - nonparametric factor
      if index[1] % acousticRate == 0
          # addFactor!(fgl, [curvar; :l1], ppbrDict[ep], autoinit=true)
          addFactor!(fgl, [curvar; :l1], pprDict[ep], autoinit=false)
          # addFactor!(fgl, [curvar; :l1], ppbDict[ep], autoinit=false)
      end

      if ep != epochs[1]
        # Odo factor x(i-1) -> xi
        addFactor!(fgl, [Symbol("x$(index[1]-1)"); curvar], odoDict[ep], autoinit=true)
      else
        # Prior to the first pose location (a "GPS" prior)
        # initLoc = [interp_x(ep);interp_y(ep);interp_yaw(ep)]
        initLoc = [itpl_lblx(ep);itpl_lbly(ep);interp_yaw(ep)]
        println("Adding a prior at $curvar, $initLoc")
        addFactor!(fgl, [curvar], PriorPose2( MvNormal(initLoc, Matrix(Diagonal([0.1;0.1;0.05].^2))) ), autoinit=true)
      end
      # Heading partial prior
      addFactor!(fgl, [curvar], RoME.PartialPriorYawPose2(Normal(interp_yaw(ep), deg2rad(5))), autoinit=true)

      # write out some metadata to help with debugging
      println(metadata, "$(index[1]), $ep, $curvar")

      ensureAllInitialized!(fgl)

      # visualizations
      plre = drawPoses(fgl,drawhist=false,contour=false,spscale=1.0, lbls=false, meanmax=:mean)
      plre.coord = Coord.Cartesian(xmin=0, xmax=100, ymin=-120, ymax=-70)
      # add confidence contour on last
      pcf = plotKDE(getKDE(fgl, curvar),dims=[1;2],levels=1,c=["gray70"] )
      union!(plre.layers, pcf.layers)
      # export
      plre |> PNG( joinpath(getSolverParams(fgl).logpath, "reckoning", "fg_$curvar.png") )

      tfg = initfg()
      addVariable!(tfg, :l1, Point2)
      initManual!(tfg,:l1,getKDE(fgl,:l1))
      addVariable!(tfg, curvar, Pose2)
      initManual!(tfg,curvar,getKDE(fgl,curvar))
      addFactor!(tfg, [curvar;:l1], pprDict[ep])

      hdl = []
      plotFactor(tfg, ls(tfg, :l1)[1], hdl=hdl)
      hdl[7].coord = Coord.Cartesian(xmin=-150.0, xmax=200.0, ymin=-250.0, ymax=100.0)
      hdl[7] |> PNG( joinpath(getSolverParams(fgl).logpath, "reckoning", "range_$curvar.png") )

      # circular illustration
      cpts = reshape(rand(Normal(interp_yaw(ep), deg2rad(5)),100),1,:)
      cbm = manikde!(cpts, Sphere1)

      opts = index[1] > 0 ? approxConv(fgl, intersect(ls(fgl, Pose2Pose2),ls(fgl, curvar))[1], curvar) : 2*randn(3,100)
      cbo = manikde!(opts[3:3,:], Sphere1)
      plcb = plotKDECircular([cbm; cbo], legend=["mag";"odo"], c=["orange";"blue"])
      plcb.coord = Coord.Cartesian(xmin=-2.5, xmax=2.5, ymin=-2.5, ymax=2.5)
      plcb |> PNG( joinpath(getSolverParams(fgl).logpath, "reckoning", "circ_$curvar.png") )

      # animate dead reckoning
      if ep != epochs[1]
        aniT = LinearInterpolation([0.0;1.0], [epochs[index[1]];epochs[index[1]+1]]) # ep
        for atim in 0:0.1:1
          @show atim, aniT(atim)
          mask = epochs[1] .< navkeys .< aniT(atim)
          if sum(mask) > 2
            pldo = Gadfly.plot(x=[X[mask];interp_x[aniT(atim)]].+DeltaX, y=[Y[mask];interp_y[aniT(atim)]].+DeltaY, Geom.path(), Theme(default_color=colorant"gray80"))
            pldo |> PNG( joinpath(getSolverParams(fgl).logpath, "reckoning", "animation", "deadreck_$(aniT(atim)).png") )
            union!(pldo.layers, plre.layers)
            pldo.coord = Coord.Cartesian(xmin=0, xmax=100, ymin=-120, ymax=-70)
            hstack(pldo, vstack(plcb, hdl[7])) |> PNG( joinpath(getSolverParams(fgl).logpath, "reckoning", "animation",   "combined_$(anicounter).png") )
            anicounter += 1
          end
        end
      end

      # # animate the estimated location
      # if ep != epochs[1]
      #   aniT = LinearInterpolation([0.0;1.0], [epochs[index[1]];epochs[index[1]+1]]) # ep
      #   for atim in 0:0.1:1
      #       @show atim, aniT(atim)
      #       rtmask = pivotPoseEpoch .< navkeys .< aniT(atim)
      #       if sum(mask) > 2
      #           # get newest deadreckoning segment
      #           xRT=[X[rtmask];interp_x[aniT(atim)]].-interp_x[pivotPoseEpoch].+pivotPoseEst[1]) yRT=[Y[rtmask];interp_y[aniT(atim)]].-interp_y[pivotPoseEpoch].+pivotPoseEst[2])
      #           # plot pivot deadreckoning
      #           pldo = Gadfly.plot(x=xRT, y=yRT, Geom.path(), Theme(default_color=colorant"red"))
      #           pldo |> PNG( joinpath(getSolverParams(fgl).logpath, "reckoning", "animation", "rt_segment_$(aniT(atim)).png") )
      #       end
      #   end
      # end

      # combined plot
      hstack(plre, vstack(plcb, hdl[7])) |> PNG( joinpath(getSolverParams(fgl).logpath, "reckoning", "combined_$curvar.png") )

      index[1] += 1
    end
    close(metadata)

    nothing
end

function deleteEpochs!(fgl, lbls::Vector{Symbol})
  function delFactsVar!(fgl, lb::Symbol)
    fcts = ls(fgl, lb)
    map(x->deleteFactor!(fgl, x), fcts)
    deleteVariable!(fgl, lb)
    nothing
  end
  map(x->delFactsVar!(fgl, x), lbls)
  nothing
end

global anicounter = 0

## Step: Building the factor graph
fg1 = initfg()
# Add a central beacon with a prior
addVariable!(fg1, :l1, Point2)

# Pinger location is [17; 1.8]
beaconprior = PriorPoint2( MvNormal([17; 1.8], Matrix(Diagonal([0.1; 0.1].^2)) ) )
addFactor!(fg1, [:l1], beaconprior, autoinit=true)

# init tree for simpler code later down
tree1 = wipeBuildNewTree!(fg1)

## and some settings

# enable fixed lag operation
# getSolverParams(fg1).isfixedlag = true
# getSolverParams(fg1).qfl = 20

getSolverParams(fg1).maxincidence = 200

# first solve and initialization
getSolverParams(fg1).drawtree = true
# getSolverParams(fg).showtree = false

mkpath(joinpath(getSolverParams(fg1).logpath, "graphs"))
mkpath(joinpath(getSolverParams(fg1).logpath, "reckoning"))
mkpath(joinpath(getSolverParams(fg1).logpath, "reckoning", "animation"))

index1 = Int[0;]
index2 = Int[0;]
storeLast = Dict{Symbol,BallTreeDensity}()

laglength = 20

length(epochs)

for STEP in 0:10:50
    global fg1, tree1
    global index1, index2
    global storeLast

    runEpochs!(fg1, epochs, STEP, index1, acousticRate=3)
    poses = sortDFG(ls(fg1, r"x"))

    if STEP-laglength >= 0
        @show poses[1], poses[end], poses[end-laglength]
        @show dellist = poses[1:(end-laglength-1)]
        if length(dellist) > 0
          deleteEpochs!(fg1, dellist)
          addFactor!(fg1, [poses[end-laglength];], PriorPose2(storeLast[poses[end-laglength]]))
        end
    end

    drawGraph(fg1, filepath=joinpath(getSolverParams(fg1).logpath, "graphs", "fg_$(STEP)_.pdf"))
    saveDFG(fg1, joinpath(getSolverParams(fg1).logpath, "graphs", "fg_$(STEP)_before"))
    @time tree1, smt, hist = solveTree!(fg1, tree1)
    saveDFG(fg1, joinpath(getSolverParams(fg1).logpath, "graphs", "fg_$(STEP)_after"))

    pla = drawPosesLandmarksAndOdo(fg1, ppbrDict, navkeys, X, Y, lblkeys, lblX, lblY)
    plb = plotSandsharkFromDFG(fg1)
    hstack(pla,plb) |> PDF(joinpath(getSolverParams(fg1).logpath, "sandshark-beacon_$STEP.pdf"))

    reportFactors(fg1, Pose2Point2Range, filepath=joinpath(getSolverParams(fg1).logpath, "acoustics_$STEP.pdf"), show=false)
    reportFactors(fg1, Pose2Pose2, filepath=joinpath(getSolverParams(fg1).logpath, "odo_$STEP.pdf"), show=false)

    poses = sortVarNested(ls(fg1, r"x"))
    storeLast[poses[end]] = getKDE(fg1, poses[end])
    @show keys(storeLast)
end

##### END OF RUN, ANALYSIS BELOW================================================

## debug weird results before IIF v0.8.0

ls(fg1)
drawGraph(fg1)

getFactorType(fg1, ls(fg1, :x0)[2])

drawTree(tree1, show=true)
plotLocalProduct(fg1, :l1)
plotTreeProductDown(fg1,tree1,:l1)

plotLocalProduct(fg1, :x0)

solveTree!(fg1)

getFactorType(fg1, ls(fg1, :l1)[2])

reportFactors(fg1, Pose2Point2Bearing)

## debug weird

fcts = ls(fg1, :l1)
PTS = map(x->approxConv(fg1, x, :l1), fcts)

PP = map(x->manikde!(x, Point2), PTS)

pp = *(PP...)

plotKDE([pp; PP...], levels=2)

# sum(timestamps .< 1531153363000000000) # 348

ls(fg1, :l1)

##
fctsym = :x150l1f1

reportFactors(fg1, Pose2Point2Range, ls(fg1, Pose2Point2Range) )
reportFactors(fg1, Pose2Pose2, ls(fg1, Pose2Pose2) )

##

plotLocalProduct(fg1, :x27, show=true)

plotKDE(ppbrDict[epochs[1]].bearing)

## debugging a plotting

dfg = fg1

drawGraph(dfg)

plotLocalProduct(dfg, :l1)
plotLocalProduct(dfg, :x100)

getSolverParams(dfg).dbg = true

tree, smt, hist = solveTree!(dfg, recordcliqs=ls(dfg))

plotSandsharkFromDFG(dfg)
drawPosesLandmarksAndOdo(dfg, ppbrDict, navkeys, X, Y, lblkeys, lblX, lblY)

# fg1_ = deepcopy(fg1)

0

## compare raw beam former data with factors

plotRawBeamFormer(fg1, ppbrDict, NAV, interp_x, interp_y, interp_yaw, show=true, all=true)

## load previous factor graph to compare
_fg = initfg()
loadDFG(joinpath(getSolverParams(fg1).logpath, "graphs", "fg_0_after.pdf"), Main, _fg)

drawGraph(_fg)

plotPose(_fg, :x0)
getFactorType)

## report on factors

dfg = fg1

fctsym = ls(dfg, Pose2Point2BearingRange)[1]
# fcss = lsf(dfg, Pose2Pose2)

plotFactor(dfg, :x20x21f1)
plotFactor(dfg, :x20l1f1)

reportFactors(dfg, Pose2Pose2, ls(dfg, Pose2Pose2)[1:2])
reportFactors(dfg, Pose2Point2BearingRange)

## overview

PL = drawPosesLandmarksAndOdo(fg1, ppbrDict, navkeys, X, Y, lblkeys, lblX, lblY)

plotSandsharkFromDFG(fg1, scale=1.5)

## slightly random

drawGraph(fg1)
plotLocalProduct(fg1, :x55)
sortVarNested(ls(fg1, r"x"))
pp, parr, partials, lb, infdim = IIF.localProduct(fg1, :x55)
plotKDE([pp;parr], levels=1, legend=["prod"; string.(lb[1:3])], dims=[1;2])
plotPose(fg1, :x59)

## DEBUG beam pattern timestamps

epochs

pll = layerBeamPatternRose(ppbrDict[ep].bearing, wRr=wRr, wTRr=theta[1:2], scale=5.0)

#

## plot all beam patterns from factor graph

plotSandsharkFromDFG(dfg)

## RANGE MODELING

fct_range = getFactorType(fg1, :x45l1f1).range

plotKDE(fct_range)

rawmf = readdlm("/home/dehann/data/sandshark/full_wombat_2018_07_09/extracted/matchedfilter/raw/$(epochs[45]).csv",',')
bss = AliasingScalarSampler(rawmf[:,1], rawmf[:,2], SNRfloor=0.6)

pts = rand(bss,100)
Gadfly.plot(x=rawmf[:,1], y=rawmf[:,2], Geom.path)

Gadfly.plot(x=pts, Geom.histogram)

pts = rand(fct_range, 1000)

## plot forward convolves

fcts = setdiff(ls(fg1, :l1), lsfPriors(fg1))

pts = approxConv(fg1, fcts[1], :l1)

pts, infd = predictbelief(fg1,:l1,fcts)

pl1 = Gadfly.plot(x=pts[1,:],y=pts[2,:],Geom.hexbin)

union!(PL.layers, pl1.layers)

PL

plotKDE(manikde!(pts, Point2().manifolds))

#

# pts = collect(azidata[1531153769000000000][:,1:1]')
#
# ptsw = TU.wrapRad.(pts)
#
# pc = manikde!(pts, Sphere1)
# pl = manikde!(pts, ContinuousScalar)
# plw = manikde!(ptsw, ContinuousScalar)
#
# plotKDE([pl;plw])


text:
######################################## connection, admin

"""
    connect(host::String, port::Int, user::String, passwd::String, dbname::String)

Connect to an OmniSci database.

# Examples
```julia-repl
julia> conn = connect("localhost", 6274, "mapd", "HyperInteractive", "mapd")
Connected to localhost:6274
```
"""
function connect(host::String, port::Int, user::String, passwd::String, dbname::String)

    socket = TSocket(host, port)

    #create libuv socket and keep-alive
    tcp = Sockets.connect(host, port)
    err = ccall(:uv_tcp_keepalive, Cint, (Ptr{Nothing}, Cint, Cuint), tcp.handle, 1, 1)
    err != 0 && error("Error setting keepalive on socket")

    Thrift.set_field!(socket, :io, tcp)
    #transport = TBufferedTransport(socket) # https://github.com/tanmaykm/Thrift.jl/issues/12
    proto = TBinaryProtocol(socket, false, true)
    c = OmniSci.MapDClient(proto)

    session = OmniSci.connect(c, user, passwd, dbname)

    return OmniSciConnection(session, c)

end

"""
    disconnect(conn::OmniSciConnection)

Close connection to OmniSci database.

# Examples
```julia-repl
julia> disconnect(conn)
Connection to localhost:6274 closed
```
"""
function disconnect(conn::OmniSciConnection)

    disconnect(conn.c, conn.session)
    println("Connection to $(conn.c.p.t.host):$(conn.c.p.t.port) closed")

end

"""
    get_status(conn::OmniSciConnection)

Displays properties of OmniSci server, such as version and rendering capabilities.

# Examples
```julia-repl
julia> status = get_status(conn)
OmniSci.TServerStatus

  read_only: false
  version: 4.1.3-20180926-66c2aee949
  rendering_enabled: false
  start_time: 1540579280
  edition: ce
  host_name: aggregator
  poly_rendering_enabled: false
```
"""
get_status(conn::OmniSciConnection) =
    get_status(conn.c, conn.session)

"""
    get_hardware_info(conn::OmniSciConnection)

Displays selected properties of hardware where OmniSci running, such as
GPU and CPU information.

# Examples
```julia-repl
julia> hardware = get_hardware_info(conn)
OmniSci.TClusterHardwareInfo(OmniSci.THardwareInfo[THardwareInfo(0, 12, 0, 0, "", OmniSci.TGpuSpecification[])])
```
"""
get_hardware_info(conn::OmniSciConnection) =
    get_hardware_info(conn.c, conn.session)

"""
    get_tables_meta(conn::OmniSciConnection; as_df::Bool = true)

Get metadata for tables in database specified in `connect()`.

# Examples
```julia-repl
julia> metad = get_tables_meta(conn)
5×6 DataFrame
│ Row │ is_replicated │ is_view │ max_rows            │ num_cols │ shard_count │ table_name        │
│     │ Bool          │ Bool    │ Int64               │ Int64    │ Int64       │ String            │
├─────┼───────────────┼─────────┼─────────────────────┼──────────┼─────────────┼───────────────────┤
│ 1   │ false         │ false   │ 4611686018427387904 │ 4        │ 0           │ omnisci_states    │
│ 2   │ false         │ false   │ 4611686018427387904 │ 6        │ 0           │ omnisci_counties  │
│ 3   │ false         │ false   │ 4611686018427387904 │ 64       │ 0           │ omnisci_countries │
│ 4   │ false         │ false   │ 4611686018427387904 │ 4        │ 0           │ test2             │
│ 5   │ false         │ false   │ 4611686018427387904 │ 4        │ 0           │ test              │
```
"""
get_tables_meta(conn::OmniSciConnection; as_df::Bool = true) =
    as_df ? DataFrame(get_tables_meta(conn.c, conn.session)) : get_tables_meta(conn.c, conn.session)

"""
    get_table_details(conn::OmniSciConnection, table_name::String; as_df::Bool = true)

Get table details such as column names and types.

# Examples
```julia-repl
julia> tbl_detail = get_table_details(conn, "omnisci_states")
4×21 DataFrame. Omitted printing of 15 columns
│ Row │ col_name    │ col_type     │ comp_param │ encoding │ is_array │ is_physical │
│     │ String      │ DataType     │ Int32      │ String   │ Bool     │ Bool        │
├─────┼─────────────┼──────────────┼────────────┼──────────┼──────────┼─────────────┤
│ 1   │ id          │ String       │ 32         │ Dict     │ false    │ false       │
│ 2   │ abbr        │ String       │ 32         │ Dict     │ false    │ false       │
│ 3   │ name        │ String       │ 32         │ Dict     │ false    │ false       │
│ 4   │ omnisci_geo │ MultiPolygon │ 32         │ GeoInt   │ false    │ false       │
```
"""
function get_table_details(conn::OmniSciConnection, table_name::String; as_df::Bool = true)

    r = get_table_details(conn.c, conn.session, table_name)

    #for dataframe, put enum values instead of keys for column type
    #makes it easier to determine column types in load_table and sql_execute
    if as_df
        df = DataFrame(r)
        df[!, :col_type] = getcolumntype.(df.col_type)
        df[!, :encoding] = getencodingtype.(df.encoding)
        return df
    else
        return r
    end
end

"""
    get_users(conn::OmniSciConnection; as_df::Bool = true)

Get list of users who have access to database specified in `connect()`.

# Examples
```julia-repl
julia> users = get_users(conn)
1×1 DataFrame
│ Row │ users  │
│     │ String │
├─────┼────────┤
│ 1   │ mapd   │
```
"""
get_users(conn::OmniSciConnection; as_df::Bool = true) =
    as_df ? DataFrame(Dict(:users => get_users(conn.c, conn.session))) : get_users(conn.c, conn.session)

"""
    get_databases(conn::OmniSciConnection; as_df::Bool=true)

Get list of databases.

# Examples
```julia-repl
julia> db = get_databases(conn)
1×2 DataFrame
│ Row │ db_name │ db_owner │
│     │ String  │ String   │
├─────┼─────────┼──────────┤
│ 1   │ mapd    │ mapd     │
```
"""
get_databases(conn::OmniSciConnection; as_df::Bool=true) =
    as_df ? DataFrame(get_databases(conn.c, conn.session)) : get_databases(conn.c, conn.session)

"""
    get_memory(conn::OmniSciConnection, memory_level::String)

Get memory profile from current session. Acceptable values for `memory_level` are "cpu" or "gpu".

# Examples
```julia-repl
julia> gm = OmniSci.get_memory(conn, "cpu")
1-element Array{OmniSci.TNodeMemoryInfo,1}:
 OmniSci.TNodeMemoryInfo("", 512, 25692761, 8388608, false, OmniSci.TMemoryData[TMemoryData(0, 0, 26, 0, [1, 2, 1, 0], 0, false), TMemoryData(0, 26, 26, 1, [1, 2, 2, 0], 0, false), TMemoryData(0, 52, 26, 2, [1, 2, 3, 0], 0, false), TMemoryData(0, 78, 26, 3, [1, 2, 4, 0], 0, false), TMemoryData(0, 104, 26, 4, [1, 2, 5, 0], 0, false), TMemoryData(0, 130, 533,5, [1, 2, 7, 0, 1], 0, false), TMemoryData(0, 663, 26, 6, [1, 2, 7, 0, 2], 0, false), TMemoryData(0, 689, 28, 7, [1, 2,8, 0, 1], 0, false), TMemoryData(0, 717, 26, 8, [1, 2, 8, 0, 2], 0, false), TMemoryData(0, 743, 28, 9, [1, 2, 9, 0, 1],0, false), TMemoryData(0, 771, 26, 10, [1, 2, 9, 0, 2], 0, false), TMemoryData(0, 797, 204, 11, [1, 2, 10, 0], 0, false), TMemoryData(0, 1001, 26, 12, [1, 2, 11, 0], 0, false), TMemoryData(0, 1027, 8387581, 0, Int64[], 0, true)])
```
"""
function get_memory(conn::OmniSciConnection, memory_level::String)

    @assert memory_level in ["cpu", "gpu"] """memory level can be one of: \"cpu\", \"gpu\""""
    get_memory(conn.c, conn.session, memory_level)
end

"""
    clear_cpu_memory(conn::OmniSciConnection)

Clears CPU memory. This function returns `nothing`.

# Examples
```julia-repl
julia> cm = OmniSci.clear_cpu_memory(conn)
```
"""
clear_cpu_memory(conn::OmniSciConnection) =
    clear_cpu_memory(conn.c, conn.session)

"""
    clear_gpu_memory(conn::OmniSciConnection)

Clears GPU memory. This function returns `nothing`.

# Examples
```julia-repl
julia> cgm = OmniSci.clear_gpu_memory(conn)

"""
clear_gpu_memory(conn::OmniSciConnection) =
    clear_gpu_memory(conn.c, conn.session)

######################################## query, render

"""
    sql_execute(conn::OmniSciConnection, query::String; first_n::Int = -1, at_most_n::Int = -1, as_df::Bool = true)

Execute a SQL query.

"""
function sql_execute(conn::OmniSciConnection, query::String; first_n::Int = -1, at_most_n::Int = -1, as_df::Bool = true)

    #asserting here, rather than pass a bad query and getting an TMapDException
    @assert (first_n < 0 || at_most_n < 0) "Only one of 'first_n' and 'at_most_n' can be set at one time"

    #true hard-coded for column_format, as its not clear there is any benefit to providing row-wise parsing (slow)
    result = sql_execute(conn.c, conn.session, query, true, randstring(32), Int32(first_n), Int32(at_most_n))

    #test if query string expected to return a df-like result
    #minimizes returning 0x0 dataframe
    as_df && startswith(lowercase(lstrip(query)), "select") ? DataFrame(result) : nothing

end

"""
    sql_execute_df(conn::OmniSciConnection, query::String, device_type::Int, device_id::Int, first_n::Int = -1)

Execute a SQL query using Apache Arrow IPC (CPU). This method requires running the code in the same environment where
OmniSci is running. This method is currently unsupported, until such time when Julia can read Arrow buffers.

"""
sql_execute_df(conn::OmniSciConnection, query::String, device_type::Int, device_id::Int, first_n::Int = -1) =
    sql_execute_df(conn.c, conn.session, query, Int32(device_type), Int32(device_id), Int32(first_n))

"""
    sql_execute_gdf(conn::OmniSciConnection, query::String, device_id::Int; first_n::Int = -1)

Execute a SQL query using Apache Arrow IPC (GPU). This method requires running the code in the same environment where
OmniSci is running. This method is currently unsupported, until such time when Julia can read Arrow buffers.

"""
sql_execute_gdf(conn::OmniSciConnection, query::String, device_id::Int; first_n::Int = -1) =
    sql_execute_gdf(conn.c, conn.session, query, Int32(device_id), Int32(first_n))

"""
    deallocate_df(conn::OmniSciConnection, df::TDataFrame, device_type::Int, device_id::Int)

"""
deallocate_df(conn::OmniSciConnection, df::TDataFrame, device_type::Int, device_id::Int) =
    deallocate_df(conn.c, conn.session, df, Int32(device_type), Int32(device_id))

"""
    interrupt(conn::OmniSciConnection)

"""
interrupt(conn::OmniSciConnection) =
    interrupt(conn.c, conn.session)

"""
    set_execution_mode(conn::OmniSciConnection, mode::TExecuteMode.Enum)

Sets execution mode for server during session. This function returns `nothing`.

# Examples
```julia-repl
julia> set_execution_mode(conn, TExecuteMode.CPU)
```

"""
set_execution_mode(conn::OmniSciConnection, mode::TExecuteMode.Enum) =
    set_execution_mode(conn.c, conn.session, mode.value)

"""
    render_vega(conn::OmniSciConnection, vega_json::String, compression_level::Int = 0)

Render an OmniSci-flavored Vega specification using the backend rendering engine. Note that OmniSci does
not currently support the full Vega specification; this method is mostly useful for rendering choropleths and
related geospatial charts.

compression_level ranges from 0 (low compression, faster) to 9 (high compression, slower).

# Examples
```julia-repl
julia> vg = {"width" : 1024, "height" : 1024...}

julia> vega_json = render_vega(conn, vg)
```
"""
function render_vega(conn::OmniSciConnection, vega_json::String, compression_level::Int = 0)

    @assert (compression_level >= 0 && compression_level <= 9) "compression_level ranges from 0 (low compression, faster) to 9 (high compression, slower)"
    render_vega(conn.c, conn.session, Int64(1), vega_json, Int32(compression_level), randstring(32))

end

######################################## dashboard

"""
    get_dashboards(conn::OmniSciConnection; as_df::Bool = true)

Gets dashboards that user submitted during connect() can access.

# Examples
```julia-repl
julia> getdbs = get_dashboards(conn)
7×8 DataFrame. Omitted printing of 3 columns
│ Row │ dashboard_id │ dashboard_metadata │ dashboard_name │ dashboard_owner │ dashboard_state │
│     │ Int32        │ String             │ String         │ String          │ String          │
├─────┼──────────────┼────────────────────┼────────────────┼─────────────────┼─────────────────┤
│ 1   │ 9            │ metadata           │ 0vcAQEO1ZD     │ mapd            │                 │
│ 2   │ 6            │ metadata           │ QI0JsthBsB     │ mapd            │                 │
│ 3   │ 5            │ metadata           │ Srm72rCJHa     │ mapd            │                 │
│ 4   │ 4            │ metadata           │ sO0XgMUOZH     │ mapd            │                 │
│ 5   │ 1            │ metadata           │ testdash       │ mapd            │                 │
│ 6   │ 2            │ metadata           │ testdash2      │ mapd            │                 │
│ 7   │ 3            │ metadata           │ testdash3      │ mapd            │                 │
```

"""
get_dashboards(conn::OmniSciConnection; as_df::Bool = true) =
    as_df ? DataFrame(get_dashboards(conn.c, conn.session)) : get_dashboards(conn.c, conn.session)

"""
    create_dashboard(conn::OmniSciConnection, dashboard_name::String, dashboard_state::String, image_hash::String, dashboard_metadata::String)

"""
create_dashboard(conn::OmniSciConnection, dashboard_name::String, dashboard_state::String, image_hash::String, dashboard_metadata::String) =
    create_dashboard(conn.c, conn.session, dashboard_name, dashboard_state, image_hash, dashboard_metadata)

"""
    get_dashboard_grantees(conn::OmniSciConnection, dashboard_id::Integer)

"""
get_dashboard_grantees(conn::OmniSciConnection, dashboard_id::Integer) =
    get_dashboard_grantees(conn.c, conn.session, Int32(dashboard_id))

######################################## import

"""
    load_table_binary_columnar(conn::OmniSciConnection, table_name::String, cols::Vector{TColumn})

"""
load_table_binary_columnar(conn::OmniSciConnection, table_name::String, cols::Vector{TColumn}) =
    load_table_binary_columnar(conn.c, conn.session, table_name, cols)

"""
    load_table_binary_columnar(conn::OmniSciConnection, table_name::String, df::DataFrame)

"""
load_table_binary_columnar(conn::OmniSciConnection, table_name::String, df::DataFrame) =
    load_table_binary_columnar(conn, table_name, TColumn.(eachcol(df)))

"""
    load_table_binary_arrow(conn::OmniSciConnection, table_name::String, arrow_stream::Vector{UInt8})

"""
load_table_binary_arrow(conn::OmniSciConnection, table_name::String, arrow_stream::Vector{UInt8}) =
    load_table_binary_arrow(conn.c, conn.session, table_name, arrow_stream)

"""
    load_table(conn::OmniSciConnection, table_name::String, rows::Vector{TStringRow})

"""
load_table(conn::OmniSciConnection, table_name::String, rows::Vector{TStringRow}) =
    load_table(conn.c, conn.session, table_name, rows)

"""
    load_table(conn::OmniSciConnection, table_name::String, rows::DataFrame)

Load a dataframe into OmniSci. This method loads data row-wise, converting data elements to string before upload.
Currently, this method requires the table to already exist on OmniSci.

# Examples
```julia-repl
julia> load_table(conn, "test", df)
```

"""
function load_table(conn::OmniSciConnection, table_name::String, df::DataFrame)

    df_to_array = [OmniSci.TStringRow(x) for x in DataFrames.eachrow(df)]
    load_table(conn, table_name, df_to_array)

end

"""
    create_table(conn::OmniSciConnection, table_name::String, df::DataFrame; dryrun::Bool = false, precision::Tuple{Int,Int} = (0,0))

Create a table in the database specified during authentication. This method takes a Julia DataFrame, reads the column
types and creates the equivalent `create table` statement, optionally printing the `create table` statement instead
of executing the statement using `sql_execute` method.

Note: This method is for convenience purposes only! It does not guarantee an optimized table statement. Additionally,
if Decimal columns present in 'df', user must set precision/scale via the 'precision' argument.

# Examples
```julia-repl
julia> create_table(conn, "test", df)
```

"""
function create_table(conn::OmniSciConnection, table_name::String, df::DataFrame;
                      dryrun::Bool = false, precision::Tuple{Int,Int} = (0,0))

    #Only assert if default value changed; check for decimal column in df later
    if precision != (0,0)
        @assert precision[1] > 0 "First number in 'precision' tuple must be greater than 0"

        #Precision of 19 OmniSci limitation
        @assert (precision[1] + precision[2]) <= 19 "Sum of tuple values for 'precision' argument must be <= 19"
    else
        #TODO: check if Decimal column in dataframe and precision = (0,0), throw warning
        #Or throw error telling to set precision
        for x in [DecFP.Dec32, DecFP.Dec64, DecFP.Dec128, Decimals.Decimal]
            @assert !(x in eltypes(df)) "Decimal column(s) detected. Please change 'precision' value from (0,0)"
        end
    end

    io = IOBuffer()
    write(io, "create table $table_name ( \n")

    for x in zip(names(df), eltypes(df))
        write(io, "  " * sanitizecolnames(x[1]) * " ") # function sanitizecolnames would need to
        write(io, "  " * getsqlcoltype(x[2], precision) * ",\n") # lookup to convert julia types to OmniSci types
    end
    query = String(take!(io))

    #remove training comma (also removes trailing newline)
    query = query[1:end-2] * "\n);"

    dryrun ? print(query) : sql_execute(conn, query)

end

######################################## object privileges

"""
    get_roles(conn::OmniSciConnection; as_df::Bool = true)

"""
get_roles(conn::OmniSciConnection; as_df::Bool = true) =
    as_df ? DataFrame(Dict(:roles => get_roles(conn.c, conn.session))) : get_roles(conn.c, conn.session)

"""
    get_all_roles_for_user(conn::OmniSciConnection, userName::String; as_df::Bool = true)

"""
get_all_roles_for_user(conn::OmniSciConnection, userName::String; as_df::Bool = true) =
    as_df ? DataFrame(Dict("roles" => get_all_roles_for_user(conn.c, conn.session, userName))) : get_all_roles_for_user(conn.c, conn.session, userName)


text:
# Hobohm I
# ========

"""
Fill `cluster` and `clustersize` matrices.
These matrices are assumed to be empty (only zeroes) and their length is assumed to be equal to the number
of sequences in the alignment (`aln`).
`threshold` is the minimum identity value between two sequences to be in the same cluster.
"""
function _fill_hobohmI!(cluster::Vector{Int}, clustersize::Vector{Int}, aln::Vector{Vector{Residue}}, threshold)
    cluster_id = 0
    nseq = length(aln)
    for i in 1:(nseq-1)
        if cluster[i] == 0
            cluster_id += 1
            cluster[i] = cluster_id
            clustersize[cluster_id] += 1
            ref_seq = aln[i]
            for j in (i+1):nseq
                if cluster[j] == 0 && percentidentity(ref_seq, aln[j], threshold)
                    cluster[j] = cluster_id
                    clustersize[cluster_id] += 1
                end
            end
        end
    end
    if cluster[nseq] == 0
        cluster_id += 1
        cluster[nseq] = cluster_id
        clustersize[cluster_id] += 1
    end
    resize!(clustersize, cluster_id)
end

"""
Calculates the weight of each sequence in a cluster.
The weight is equal to one divided by the number of sequences in the cluster.
"""
function _get_sequence_weight(clustersize, cluster)
    nseq = length(cluster)
    sequence_weight = Array(Float64, nseq)
    for i in 1:nseq
        @inbounds sequence_weight[i] = 1.0 / clustersize[cluster[i]]
    end
    sequence_weight
end

"""
Sequence clustering using the Hobohm I method from Hobohm et. al. 1992.

*Hobohm, Uwe, et al. "Selection of representative protein data sets." Protein Science 1.3 (1992): 409-417.*
"""
function hobohmI(msa::AbstractMatrix{Residue}, threshold)
    aln = getresiduesequences(msa)
    nseq = length(aln)
    cluster = zeros(Int,nseq)
    clustersize = zeros(Int,nseq)
    _fill_hobohmI!(cluster, clustersize, aln, threshold)
    SequenceClusters(clustersize, cluster, _get_sequence_weight(clustersize, cluster))
end


text:
using brainflow

# specify logging library to use
brainflow.enable_dev_brainflow_logger(brainflow.BOARD_CONTROLLER)

params = BrainFlowInputParams()
board_shim = brainflow.BoardShim(brainflow.SYNTHETIC_BOARD, params)

brainflow.prepare_session(board_shim)
brainflow.start_stream(board_shim)
sleep(5)
brainflow.stop_stream(board_shim)
data = brainflow.get_current_board_data(32, board_shim)
brainflow.release_session(board_shim)

eeg_channels = brainflow.get_eeg_channels(brainflow.SYNTHETIC_BOARD)
sampling_rate = brainflow.get_sampling_rate(brainflow.SYNTHETIC_BOARD)

data_first_channel = data[eeg_channels[1], :]
println("Original Data First Channel")
println(data_first_channel)
brainflow.perform_lowpass(data_first_channel, sampling_rate, 10.0, 3, brainflow.BUTTERWORTH, 0.0)
println("After LowPass Filter")
println(data_first_channel)

data_second_channel = data[eeg_channels[2], :]
println("Original Data Second Channel")
println(data_second_channel)
brainflow.perform_highpass(data_second_channel, sampling_rate, 5.0, 3, brainflow.CHEBYSHEV_TYPE_1, 1.0)
println("After HighPass Filter")
println(data_second_channel)

data_third_channel = data[eeg_channels[3], :]
println("Original Data Third Channel")
println(data_third_channel)
brainflow.perform_bandpass(data_third_channel, sampling_rate, 25.0, 20.0, 3, brainflow.BESSEL, 0.0)
println("After BandPass Filter")
println(data_third_channel)

data_fourth_channel = data[eeg_channels[4], :]
println("Original Data Fourth Channel")
println(data_fourth_channel)
brainflow.perform_bandstop(data_fourth_channel, sampling_rate, 50.0, 2.0, 3, brainflow.BESSEL, 0.0)
println("After BandStop Filter")
println(data_fourth_channel)

max_stars_repo_name: matthijscox/brainflow

max_issues_repo_name: matthijscox/brainflow

max_forks_repo_name: matthijscox/brainflow


text:
function md(args...; env=:auto, kwargs...)
    md_function = infer_md_output(env, args...)

    m = md_function(args...; kwargs...)
    COPY_TO_CLIPBOARD && clipboard(m)
    AUTO_DISPLAY && display(m)
    return m
end

mdalign(args...; kwargs...) = latexalign(args...; kwargs...)
mdarray(args...; kwargs...) = latexarray(args...; kwargs...)
md_chemical_arrows(args...; kwargs...) = chemical_arrows(args...; kwargs...)

function infer_md_output(env, args...)
    if env != :auto
        env == :table && return mdtable
        env == :text && return mdtext
        env == :align && return mdalign
        env == :array && return mdarray
        env == :inline && return latexinline
        env in [:arrows, :chem, :chemical, :arrow] && return md_chemical_arrows
        error("The MD environment $env is not defined.")
    end

    md_function = get_md_function(args...)

    return md_function
end

"""
    get_md_function(args...)

Use overloading to determine what MD output to generate.

This determines the default behaviour of `md()` for different inputs.
"""
get_md_function(args...) = mdtext
get_md_function(args::AbstractArray...) = mdtable
get_md_function(args::AbstractDict) = mdtable
get_md_function(args::Tuple) = mdtable


text:
############################################################
# judiVector # #################################################
############################################################

# Authors: Philipp Witte (pwitte@eos.ubc.ca), Henryk Modzelewski (hmodzelewski@eos.ubc.ca)
# Date: January 2017

export judiVector, judiVectorException, subsample, judiVector_to_SeisBlock, time_resample, time_resample!, judiTimeInterpolation, write_shot_record, get_data

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

# structure for seismic data as an abstract vector
mutable struct judiVector{vDT<:Number} <: joAbstractLinearOperator{vDT,vDT}
    name::String
    m::Integer
    n::Integer
    nsrc::Integer
    geometry::Geometry
    data
end

mutable struct judiVectorException <: Exception
    msg :: String
end

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

## outer constructors

"""
    judiVector
        name::String
        m::Integer
        n::Integer
        nsrc::Integer
        geometry::Geometry
        data

Abstract vector for seismic data. This vector-like structure contains the geometry and data for either\\
receiver data (shot records) or source data (wavelets).

Constructors
============

Construct vector from `Geometry` structure and cell array of shot records or wavelets. The `data` keyword\\
can also be a single (non-cell) array, in which case the data is the same for all source positions:

    judiVector(geometry, data)

Construct vector for observed data from `SegyIO.SeisBlock`. `segy_depth_key` is the `SegyIO` keyword \\
that contains the receiver depth coordinate:

    judiVector(SegyIO.SeisBlock; segy_depth_key="RecGroupElevation")

Construct vector for observed data from out-of-core data container of type `SegyIO.SeisCon`:

    judiVector(SegyIO.SeisCon; segy_depth_key="RecGroupElevation")

Examples
========

(1) Construct data vector from `Geometry` structure and a cell array of shot records:

    dobs = judiVector(rec_geometry, shot_records)

(2) Construct data vector for a seismic wavelet (can be either cell arrays of individual\\
wavelets or a single wavelet as an array):

    q = judiVector(src_geometry, wavelet)

(3) Construct data vector from `SegyIO.SeisBlock` object:

    using SegyIO
    seis_block = segy_read("test_file.segy")
    dobs = judiVector(seis_block; segy_depth_key="RecGroupElevation")

(4) Construct out-of-core data vector from `SegyIO.SeisCon` object (for large SEG-Y files):

    using SegyIO
    seis_container = segy_scan("/path/to/data/directory","filenames",["GroupX","GroupY","RecGroupElevation","SourceDepth","dt"])
    dobs = judiVector(seis_container; segy_depth_key="RecGroupElevation")

"""
function judiVector(geometry::Geometry,data::Array; vDT::DataType=Float32)
    vDT == Float32 || throw(judiVectorException("Domain type not supported"))

    # length of vector
    n = 1
    if typeof(geometry) == GeometryOOC
        nsrc = length(geometry.container)
        m = sum(geometry.nsamples)
    else
        nsrc = length(geometry.xloc)
        m = 0
        for j=1:nsrc
            m += length(geometry.xloc[j])*geometry.nt[j]
        end
    end
    dataCell = Array{Array}(undef, nsrc)
    for j=1:nsrc
        dataCell[j] = data
    end
    return judiVector{Float32}("Seismic data vector",m,n,nsrc,geometry,dataCell)
end

# constructor if data is passed as a cell array
function judiVector(geometry::Geometry,data::Union{Array{Any,1},Array{Array,1}}; vDT::DataType=Float32)
    vDT == Float32 || throw(judiVectorException("Domain and range types not supported"))

    # length of vector
    if typeof(geometry) == GeometryOOC
        nsrc = length(geometry.container)
        m = sum(geometry.nsamples)
    else
        nsrc = length(geometry.xloc)
        m = 0
        for j=1:nsrc
            m += length(geometry.xloc[j])*geometry.nt[j]
        end
    end
    n = 1
    return judiVector{Float32}("Seismic data vector",m,n,nsrc,geometry,data)
end

############################################################
# constructors from SEGY files or out-of-core containers

# contructor for in-core data container
function judiVector(data::SegyIO.SeisBlock; segy_depth_key="RecGroupElevation", vDT::DataType=Float32)
    vDT == Float32 || throw(judiVectorException("Domain and range types not supported"))

    # length of data vector
    src = get_header(data,"FieldRecord")
    nsrc = length(unique(src))

    numTraces = get_header(data,"TraceNumber")[end] - get_header(data,"TraceNumber")[1] + 1
    numSamples = get_header(data,"ns")[end]
    m = numTraces*numSamples
    n = 1

    # extract geometry from data container
    geometry = Geometry(data; key="receiver", segy_depth_key=segy_depth_key)

    # fill data vector with pointers to data location
    dataCell = Array{Array}(undef, nsrc)
    full_data = convert(Array{Float32,2},data.data)
    for j=1:nsrc
        traces = findall(src .== unique(src)[j])
        dataCell[j] = full_data[:,traces]
    end

    return judiVector{Float32}("Julia seismic data container",m,n,nsrc,geometry,dataCell)
end

# contructor for in-core data container and given geometry
function judiVector(geometry::Geometry, data::SegyIO.SeisBlock; vDT::DataType=Float32)
    vDT == Float32 || throw(judiVectorException("Domain and range types not supported"))

    # length of data vector
    src = get_header(data,"FieldRecord")
    nsrc = length(unique(src))

    numTraces = get_header(data,"TraceNumber")[end] - get_header(data,"TraceNumber")[1] + 1
    numSamples = get_header(data,"ns")[end]
    m = numTraces*numSamples
    n = 1

    # fill data vector with pointers to data location
    dataCell = Array{Array}(undef, nsrc)
    for j=1:nsrc
        traces = findall(src .== unique(src)[j])
        dataCell[j] = convert(Array{Float32,2},data.data[:,traces])
    end

    return judiVector{Float32}("Julia seismic data container",m,n,nsrc,geometry,dataCell)
end

# contructor for out-of-core data container from single container
function judiVector(data::SegyIO.SeisCon; segy_depth_key="RecGroupElevation", vDT::DataType=Float32)
    vDT == Float32 || throw(judiVectorException("Domain and range types not supported"))

    # length of data vector
    nsrc = length(data)
    numTraces = 0
    for j=1:nsrc
        numTraces += Int((data.blocks[j].endbyte - data.blocks[j].startbyte)/(240 + data.ns*4))
    end
    m = numTraces*data.ns
    n = 1

    # extract geometry from data container
    geometry = Geometry(data; key="receiver", segy_depth_key=segy_depth_key)

    # fill data vector with pointers to data location
    dataCell = Array{SegyIO.SeisCon}(undef, nsrc)
    for j=1:nsrc
        dataCell[j] = split(data,j)
    end

    return judiVector{Float32}("Julia seismic data container",m,n,nsrc,geometry,dataCell)
end

# contructor for single out-of-core data container and given geometry
function judiVector(geometry::Geometry, data::SegyIO.SeisCon; vDT::DataType=Float32)
    vDT == Float32 || throw(judiVectorException("Domain and range types not supported"))

    # length of data vector
    nsrc = length(data)
    numTraces = 0
    for j=1:nsrc
        numTraces += Int((data.blocks[j].endbyte - data.blocks[j].startbyte)/(240 + data.ns*4))
    end
    m = numTraces*data.ns
    n = 1

    # fill data vector with pointers to data location
    dataCell = Array{SegyIO.SeisCon}(undef, nsrc)
    for j=1:nsrc
        dataCell[j] = split(data,j)
    end

    return judiVector{Float32}("Julia seismic data container",m,n,nsrc,geometry,dataCell)
end

# contructor for out-of-core data container from cell array of containers
function judiVector(data::Array{SegyIO.SeisCon,1}; segy_depth_key="RecGroupElevation", vDT::DataType=Float32)
    vDT == Float32 || throw(judiVectorException("Domain and range types not supported"))

    # length of data vector
    nsrc = length(data)
    numTraces = 0
    for j=1:nsrc
        numTraces += Int((data[j].blocks[1].endbyte - data[j].blocks[1].startbyte)/(240 + data[j].ns*4))
    end
    m = numTraces*data[1].ns    # SEGY requires same number of samples for every trace
    n = 1

    # extract geometry from data container
    geometry = Geometry(data; key="receiver", segy_depth_key=segy_depth_key)

    # fill data vector with pointers to data location
    dataCell = Array{SegyIO.SeisCon}(undef, nsrc)
    for j=1:nsrc
        dataCell[j] = data[j]
    end

    return judiVector{Float32}("Julia seismic data container",m,n,nsrc,geometry,dataCell)
end

# contructor for out-of-core data container from cell array of containers and given geometry
function judiVector(geometry::Geometry, data::Array{SegyIO.SeisCon}; vDT::DataType=Float32)
    vDT == Float32 || throw(judiVectorException("Domain and range types not supported"))

    # length of data vector
    nsrc = length(data)
    numTraces = 0
    for j=1:nsrc
        numTraces += Int((data[j].blocks[1].endbyte - data[j].blocks[1].startbyte)/(240 + data[j].ns*4))
    end
    m = numTraces*data[1].ns
    n = 1

    # fill data vector with pointers to data location
    dataCell = Array{SegyIO.SeisCon}(undef, nsrc)
    for j=1:nsrc
        dataCell[j] = data[j]
    end

    return judiVector{Float32}("Julia seismic data container",m,n,nsrc,geometry,dataCell)
end

############################################################
## overloaded Base functions

# conj(jo)
conj(a::judiVector{vDT}) where vDT =
    judiVector{vDT}("conj("*a.name*")",a.m,a.n,a.nsrc,a.geometry,a.data)

# transpose(jo)
transpose(a::judiVector{vDT}) where vDT =
    judiVector{vDT}(""*a.name*".'",a.n,a.m,a.nsrc,a.geometry,a.data)

# adjoint(jo)
adjoint(a::judiVector{vDT}) where vDT =
        judiVector{vDT}(""*a.name*".'",a.n,a.m,a.nsrc,a.geometry,a.data)

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

# +(judiVector, judiVector)
function +(a::judiVector{avDT}, b::judiVector{bvDT}) where {avDT, bvDT}
    size(a) == size(b) || throw(judiVectorException("dimension mismatch"))
    compareGeometry(a.geometry, b.geometry) == 1 || throw(judiVectorException("geometry mismatch"))
    c = deepcopy(a)
    c.data = a.data + b.data
    return c
end

# -(judiVector, judiVector)
function -(a::judiVector{avDT}, b::judiVector{bvDT}) where {avDT, bvDT}
    size(a) == size(b) || throw(judiVectorException("dimension mismatch"))
    compareGeometry(a.geometry, b.geometry) == 1 || throw(judiVectorException("geometry mismatch"))
    c = deepcopy(a)
    c.data = a.data - b.data
    return c
end

# +(judiVector, number)
function +(a::judiVector{avDT},b::Number) where avDT
    c = deepcopy(a)
    c.data = c.data .+ b
    return c
end

# +(number, judiVector)
function +(a::Number,b::judiVector{avDT}) where avDT
    c = deepcopy(b)
    c.data = b.data .+ a
    return c
end

# -(judiVector, number)
function -(a::judiVector{avDT},b::Number) where avDT
    c = deepcopy(a)
    c.data = c.data .- b
    return c
end

# *(judiVector, number)
function *(a::judiVector{avDT},b::Number) where avDT
    c = deepcopy(a)
    c.data = c.data .* b
    return c
end

# *(number, judiVector)
function *(a::Number,b::judiVector{bvDT}) where bvDT
    c = deepcopy(b)
    c.data = a .* c.data
    return c
end

# /(judiVector, number)
function /(a::judiVector{avDT},b::Number) where avDT
    c = deepcopy(a)
    if iszero(b)
        error("Division by zero")
    else
        c.data = c.data/b
    end
    return c
end

# *(joLinearFunction, judiVector)
function *(A::joLinearFunction{ADDT,ARDT},v::judiVector{avDT}) where {ADDT, ARDT, avDT}
    A.n == size(v,1) || throw(judiVectorException("shape mismatch"))
    jo_check_type_match(ADDT,avDT,join(["DDT for *(joLinearFunction,judiVector):",A.name,typeof(A),avDT]," / "))
    V = A.fop(v)
    jo_check_type_match(ARDT,eltype(V),join(["RDT from *(joLinearFunction,judiVector):",A.name,typeof(A),eltype(V)]," / "))
    return V
end

# *(joLinearOperator, judiVector)
function *(A::joLinearOperator{ADDT,ARDT},v::judiVector{avDT}) where {ADDT, ARDT, avDT}
    A.n == size(v,1) || throw(judiVectorException("shape mismatch"))
    jo_check_type_match(ADDT,avDT,join(["DDT for *(joLinearFunction,judiVector):",A.name,typeof(A),avDT]," / "))
    V = A.fop(v)
    jo_check_type_match(ARDT,eltype(V),join(["RDT from *(joLinearFunction,judiVector):",A.name,typeof(A),eltype(V)]," / "))
    return V
end

# vcat
function vcat(a::judiVector{avDT},b::judiVector{bvDT}) where {avDT, bvDT}
    typeof(a.geometry) == typeof(b.geometry) || throw(judiVectorException("Geometry type mismatch"))
    m = a.m + b.m
    n = 1
    nsrc = a.nsrc + b.nsrc

    if typeof(a.data) == Array{SegyIO.SeisCon,1} && typeof(b.data) == Array{SegyIO.SeisCon,1}
        data = Array{SegyIO.SeisCon}(undef, nsrc)
    else
        data = Array{Array}(undef, nsrc)
    end

    dt = Array{Any}(undef, nsrc)
    nt = Array{Any}(undef, nsrc)
    t = Array{Any}(undef, nsrc)
    if typeof(data) == Array{SegyIO.SeisCon,1}
        nsamples = Array{Any}(undef, nsrc)
    else
        xloc = Array{Any}(undef, nsrc)
        yloc = Array{Any}(undef, nsrc)
        zloc = Array{Any}(undef, nsrc)
    end

    # Merge data sets and geometries
    for j=1:a.nsrc
        data[j] = a.data[j]
        if typeof(data) == Array{SegyIO.SeisCon,1}
            nsamples[j] = a.geometry.nsamples[j]
        else
            xloc[j] = a.geometry.xloc[j]
            yloc[j] = a.geometry.yloc[j]
            zloc[j] = a.geometry.zloc[j]
        end
        dt[j] = a.geometry.dt[j]
        nt[j] = a.geometry.nt[j]
        t[j] = a.geometry.t[j]
    end
    for j=a.nsrc+1:nsrc
        data[j] = b.data[j-a.nsrc]
        if typeof(data) == Array{SegyIO.SeisCon,1}
            nsamples[j] = b.geometry.nsamples[j-a.nsrc]
        else
            xloc[j] = b.geometry.xloc[j-a.nsrc]
            yloc[j] = b.geometry.yloc[j-a.nsrc]
            zloc[j] = b.geometry.zloc[j-a.nsrc]
        end
        dt[j] = b.geometry.dt[j-a.nsrc]
        nt[j] = b.geometry.nt[j-a.nsrc]
        t[j] = b.geometry.t[j-a.nsrc]
    end

    if typeof(data) == Array{SegyIO.SeisCon,1}
        geometry = GeometryOOC(data,dt,nt,t,nsamples,a.geometry.key,a.geometry.segy_depth_key)
    else
        geometry = Geometry(xloc,yloc,zloc,dt,nt,t)
    end
    nvDT = promote_type(avDT,bvDT)
    return judiVector{nvDT}(a.name,m,n,nsrc,geometry,data)
end

# dot product
function dot(a::judiVector{avDT}, b::judiVector{bvDT}) where {avDT, bvDT}
# Dot product for data containers
    size(a) == size(b) || throw(judiVectorException("dimension mismatch"))
    compareGeometry(a.geometry, b.geometry) == 1 || throw(judiVectorException("geometry mismatch"))
    dotprod = 0f0
    for j=1:a.nsrc
        dotprod += a.geometry.dt[j] * dot(vec(a.data[j]),vec(b.data[j]))
    end
    return dotprod
end

# norm
function norm(a::judiVector{avDT}, p::Real=2) where avDT
    x = 0.f0
    for j=1:a.nsrc
        x += a.geometry.dt[j] * sum(abs.(vec(a.data[j])).^p)
    end
    return x^(1.f0/p)
end

# abs
function abs(a::judiVector{avDT}) where avDT
    b = deepcopy(a)
    for j=1:a.nsrc
        b.data[j] = abs.(a.data[j])
    end
    return b
end

# Subsample data container
"""
    subsample(x,source_numbers)

Subsample seismic data vectors or matrix-free linear operators and extract the entries that correspond\\
to the shot positions defined by `source_numbers`. Works for inputs of type `judiVector`, `judiModeling`, \\
`judiProjection`, `judiJacobian`, `Geometry`, `judiRHS`, `judiPDE`, `judiPDEfull`.

Examples
========

(1) Extract 2 shots from `judiVector` vector:

    dsub = subsample(dobs,[1,2])

(2) Extract geometry for shot location 100:

    geometry_sub = subsample(dobs.geometry,100)

(3) Extract Jacobian for shots 10 and 20:

    Jsub = subsample(J,[10,20])

"""
function subsample(a::judiVector{avDT},srcnum) where avDT
    geometry = subsample(a.geometry,srcnum)     # Geometry of subsampled data container
    return judiVector(geometry,a.data[srcnum];vDT=avDT)
end

getindex(x::judiVector,a) = subsample(x,a)

# Create SeisBlock from judiVector container to write to file
function judiVector_to_SeisBlock(d::judiVector{avDT}, q::judiVector{avDT}; source_depth_key="SourceSurfaceElevation", receiver_depth_key="RecGroupElevation") where avDT

    typeof(d.geometry) == GeometryOOC && (d.geometry = Geometry(d.geometry))
    typeof(q.geometry) == GeometryOOC && (q.geometry = Geometry(q.geometry))

    blocks = Array{Any}(undef, d.nsrc)
    count = 0
    for j=1:d.nsrc

        # create SeisBlock
        blocks[j] = SeisBlock(d.data[j])
        numTraces = size(d.data[j],2)
        traceNumbers = convert(Array{Integer,1},count+1:count+numTraces)

        # set headers
        set_header!(blocks[j], "GroupX", convert(Array{Integer,1},round.(d.geometry.xloc[j]*1f3)))
        if length(d.geometry.yloc[j]) == 1
            set_header!(blocks[j], "GroupY", Int(round(d.geometry.yloc[j]*1f3)))
        else
            set_header!(blocks[j], "GroupY", convert(Array{Integer,1},round.(d.geometry.yloc[j]*1f3)))
        end
        set_header!(blocks[j], receiver_depth_key, convert(Array{Integer,1},round.(d.geometry.zloc[j]*1f3)))
        set_header!(blocks[j], "SourceX", Int(round.(q.geometry.xloc[j]*1f3)))
        set_header!(blocks[j], "SourceY", Int(round.(q.geometry.yloc[j]*1f3)))
        set_header!(blocks[j], source_depth_key, Int(round.(q.geometry.zloc[j]*1f3)))

        set_header!(blocks[j], "dt", Int(d.geometry.dt[j]*1f3))
        set_header!(blocks[j], "FieldRecord",j)
        set_header!(blocks[j], "TraceNumWithinLine", traceNumbers)
        set_header!(blocks[j], "TraceNumWithinFile", traceNumbers)
        set_header!(blocks[j], "TraceNumber", traceNumbers)
        set_header!(blocks[j], "ElevationScalar", -1000)
        set_header!(blocks[j], "RecSourceScalar", -1000)
        count += numTraces
    end

    # merge into single block
    fullblock = blocks[1]
    for j=2:d.nsrc
        fullblock = merge(fullblock,blocks[j])
        blocks[j] = []
    end
    return fullblock
end

function write_shot_record(srcGeometry,srcData,recGeometry,recData,options)
    q = judiVector(srcGeometry,srcData)
    d = judiVector(recGeometry,recData)
    file = join([string(options.file_name),"_",string(srcGeometry.xloc[1][1]),"_",string(srcGeometry.yloc[1][1]),".segy"])
    block_out = judiVector_to_SeisBlock(d,q)
    segy_write(join([options.file_path,"/",file]), block_out)
    container = scan_file(join([options.file_path,"/",file]),["GroupX","GroupY","dt","SourceSurfaceElevation","RecGroupElevation"])
    return container
end

function time_resample!(x::judiVector,dt_new;order=2)
    for j=1:x.nsrc
        numTraces = size(x.data[j],2)
        timeAxis = 0:x.geometry.dt[j]:x.geometry.t[j]
        timeInterp = 0:dt_new:x.geometry.t[j]
        dataInterp = zeros(Float32,length(timeInterp),numTraces)
        for k=1:numTraces
            spl = Spline1D(timeAxis,x.data[j][:,k];k=order)
            dataInterp[:,k] = spl(timeInterp)
        end
        x.data[j] = dataInterp
        x.geometry.dt[j] = dt_new
        x.geometry.nt[j] = length(timeInterp)
        x.geometry.t[j] = (x.geometry.nt[j] - 1)*x.geometry.dt[j]
    end
    return judiVector(x.geometry,x.data)
end

function time_resample(x::judiVector,dt_new;order=2)
    xout = deepcopy(x)
    for j=1:x.nsrc
        numTraces = size(x.data[j],2)
        timeAxis = 0:x.geometry.dt[j]:x.geometry.t[j]
        timeInterp = 0:dt_new:x.geometry.t[j]
        dataInterp = zeros(Float32,length(timeInterp),numTraces)
        for k=1:numTraces
            spl = Spline1D(timeAxis,x.data[j][:,k];k=order)
            dataInterp[:,k] = spl(timeInterp)
        end
        xout.data[j] = dataInterp
        xout.geometry.dt[j] = dt_new
        xout.geometry.nt[j] = length(timeInterp)
        xout.geometry.t[j] = (x.geometry.nt[j] - 1)*x.geometry.dt[j]
    end
    return judiVector(xout.geometry,xout.data)
end

function judiTimeInterpolation(geometry::Geometry,dt_coarse,dt_fine)
# Time interpolation as a linear operator (copies input data)

    nsrc = length(geometry.xloc)
    m = 0
    n = 0
    for j=1:nsrc
        nt_coarse = Int(trunc(geometry.t[j]/dt_coarse + 1))
        nt_fine = Int(trunc(geometry.t[j]/dt_fine + 1))
        n += length(geometry.xloc[j])*nt_coarse
        m += length(geometry.xloc[j])*nt_fine
    end
    I = joLinearFunctionFwd_T(m,n,
                             v -> time_resample(v,dt_fine),
                             w -> time_resample(w,dt_coarse),
                             Float32,Float32,name="Time resampling")
    return I
end

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

#function scale!(a::Number,x::judiVector)
#    for j=1:x.nsrc
#        x.data[j] *= a
#    end
#end

#function scale!(x::judiVector,a::Number)
#    for j=1:x.nsrc
#        x.data[j] *= a
#    end
#end

####################################################################################################
# Indexing

#getindex(x::judiVector, i) = x.data[i]

setindex!(x::judiVector, y, i) = x.data[i][:] = y

firstindex(x::judiVector) = 1

lastindex(x::judiVector) = x.nsrc

axes(x::judiVector) = Base.OneTo(x.nsrc)

ndims(x::judiVector) = length(size(x))

similar(x::judiVector, element_type::DataType, dims::Union{AbstractUnitRange, Integer}...) = judiVector(x.geometry, x.data)*0f0

function sum(x::judiVector)
    s = 0f0
    for j=1:length(x.data)
        s += vec(x.data[j])
    end
    return s
end

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

BroadcastStyle(::Type{judiVector}) = Base.Broadcast.DefaultArrayStyle{1}()

ndims(::Type{judiVector{Float32}}) = 1

broadcasted(::typeof(+), x::judiVector, y::judiVector) = x + y

broadcasted(::typeof(-), x::judiVector, y::judiVector) = x - y

function broadcasted(::typeof(*), x::judiVector, y::judiVector)
    z = deepcopy(x)
    for j=1:length(x.data)
        z.data[j] = x.data[j] .* y.data[j]
    end
    return z
end

function broadcasted!(::typeof(*), x::judiVector, y::judiVector)
    for j=1:length(x.data)
        x.data[j] = x.data[j] .* y.data[j]
    end
    return x
end

function broadcasted(::typeof(/), x::judiVector, y::judiVector)
    z = deepcopy(x)
    for j=1:length(x.data)
        z.data[j] = x.data[j] ./ y.data[j]
    end
    return z
end

function broadcasted!(::typeof(/), x::judiVector, y::judiVector)
    for j=1:length(x.data)
        x.data[j] = x.data[j] .* y.data[j]
    end
    return x
end

broadcast!(.*, x::judiVector, y::judiVector, a::Number) = scale!(y, a)

function broadcast!(identity, x::judiVector, y::judiVector)
    copy!(x,y)
end

function broadcast!(identity, x::judiVector, a::Number, y::judiVector, z::judiVector)
    scale!(y,a)
    copy!(x, y + z)
end

function copy!(x::judiVector,y::judiVector)
    for j=1:x.nsrc
        x.data[j] = y.data[j]
    end
    x.geometry = deepcopy(y.geometry)
end

#function axpy!(a::Number,X::judiVector,Y::judiVector)
#    for j=1:Y.nsrc
#        Y.data[j] = a*X.data[j] + Y.data[j]
#    end
#end

function get_data(x::judiVector)
    shots = Array{Any}(undef, x.nsrc)
    rec_geometry = Geometry(x.geometry)

    for j=1:x.nsrc
        shots[j] = convert(Array{Float32, 2}, x.data[j][1].data)
    end
    return judiVector(rec_geometry, shots)
end

function isapprox(x::judiVector, y::judiVector; rtol::Real=sqrt(eps()), atol::Real=0)
    compareGeometry(x.geometry, y.geometry) == 1 || throw(judiVectorException("geometry mismatch"))
    isapprox(x.data, y.data; rtol=rtol, atol=atol)
end

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

# Overload base function for SegyIO objects

vec(x::SegyIO.SeisCon) = vec(x[1].data)
norm(x::SegyIO.IBMFloat32; kwargs...) = norm(convert(Float32,x); kwargs...)
dot(x::SegyIO.IBMFloat32, y::SegyIO.IBMFloat32) = dot(convert(Float32,x), convert(Float32,y))
dot(x::SegyIO.IBMFloat32, y::Float32) = dot(convert(Float32,x), y)
dot(x::Float32, y::SegyIO.IBMFloat32) = dot(x, convert(Float32,y))

# binary operations return dense arrays
+(x::SegyIO.SeisCon, y::SegyIO.SeisCon) = +(x[1].data,y[1].data)
+(x::SegyIO.IBMFloat32, y::SegyIO.IBMFloat32) = +(convert(Float32,x),convert(Float32,y))
+(x::SegyIO.IBMFloat32, y::Float32) = +(convert(Float32,x),y)
+(x::Float32, y::SegyIO.IBMFloat32) = +(x,convert(Float32,y))

-(x::SegyIO.SeisCon, y::SegyIO.SeisCon) = -(x[1].data,y[1].data)
-(x::SegyIO.IBMFloat32, y::SegyIO.IBMFloat32) = -(convert(Float32,x),convert(Float32,y))
-(x::SegyIO.IBMFloat32, y::Float32) = -(convert(Float32,x),y)
-(x::Float32, y::SegyIO.IBMFloat32) = -(x,convert(Float32,y))

+(a::Number, x::SegyIO.SeisCon) = .+(a,x[1].data)
+(x::SegyIO.SeisCon, a::Number) = .+(x[1].data,a)

-(a::Number, x::SegyIO.SeisCon) = -(a,x[1].data)
-(x::SegyIO.SeisCon, a::Number) = -(x[1].data,a)

*(a::Number, x::SegyIO.SeisCon) = *(a,x[1].data)
*(x::SegyIO.SeisCon, a::Number) = *(x[1].data,a)

/(a::Number, x::SegyIO.SeisCon) = /(a,x[1].data)
/(x::SegyIO.SeisCon, a::Number) = /(x[1].data,a)


text:
using SampleChainsDynamicHMC
using Soss

g = @model N begin
    n ~ Poisson(10) |> iid(N)
    p ~ Uniform()
    y ~ For(n) do nj Binomial(nj, p) end
end;

obs = predict(g, N=5, p=0.2)

pairs(obs)

m = predictive(g, :n)

post = m(n=[6,12,7,8,12]) | (y=[1,3,1,2,3],)

s = sample(DynamicHMCChain, post)

predict(m(n=obs.n), s)

drawsamples!(s, 1000)

a1 = ash(getchains(s)[1].p)
a2 = ash(getchains(s)[2].p)
a3 = ash(getchains(s)[3].p)
a4 = ash(getchains(s)[4].p)

plot(a1.rng, a1.density)
plot!(a2.rng, a2.density)
plot!(a3.rng, a3.density)
plot!(a4.rng, a4.density)

drawsamples!(s, 10000)

p = 0.001:0.001:0.999
plot(p, exp.([ℓ((p=pj,)) for pj in p]))

using TransformVariables

ℓ(x) = logdensity(Beta(4,3), x.p)
t = xform(post)
chain = initialize!(DynamicHMCChain, ℓ, t)
drawsamples!(chain, 10000)
plot(ash(chain.p))


text:
using StanSample

ProjDir = @__DIR__

themodel = "
data {
  int<lower=1> N;
  int<lower=0,upper=1> y[N];
}
parameters {
  real<lower=0,upper=1> theta;
}
model {
  theta ~ beta(1,1);
  y ~ bernoulli(theta);
}
";

# Default for tmpdir is to create a new tmpdir location
# To prevent recompilation of a Stan progam, choose a fixed location,
# and to initally hav access to run logs etc.
tmpdir = "$(ProjDir)/tmp"

all_models = Vector{SampleModel}(undef, 2)
for i in 1:length(all_models)
  all_models[i] = SampleModel("sm_$i", themodel,
    tmpdir=tmpdir,
    method=StanSample.Sample(save_warmup=true, num_warmup=1000,
      num_samples=1000, adapt=StanSample.Adapt(delta=0.85))
  );
end

# This could be any Vector type (or Dict with Int keys)
observed_data = Dict(:N => 10, :y => [0, 1, 0, 1, 0, 0, 0, 0, 0, 1])

rc = stan_sample.(all_models; data=observed_data);

p = []
if all(success.(rc))
  for m in all_models
    push!(p, read_samples(m, output_format=:particles))
  end
  display(p)
end