text:
mutable struct JacobianCache{CacheType1,CacheType2,CacheType3,ColorType,SparsityType,fdtype,returntype}
    x1  :: CacheType1
    fx  :: CacheType2
    fx1 :: CacheType3
    colorvec :: ColorType
    sparsity :: SparsityType
end

function JacobianCache(
    x,
    fdtype     :: Union{Val{T1},Type{T1}} = Val(:forward),
    returntype :: Type{T2} = eltype(x);
    inplace    :: Union{Val{T3},Type{T3}} = Val(true),
    colorvec = 1:length(x),
    sparsity = nothing) where {T1,T2,T3}

    fdtype isa Type && (fdtype = fdtype())
    inplace isa Type && (inplace = inplace())
    if eltype(x) <: Real && fdtype==Val(:complex)
        x1  = false .* im .* x
        _fx = false .* im .* x
    else
        x1 = copy(x)
        _fx = copy(x)
    end

    if fdtype==Val(:complex)
        _fx1  = nothing
    else
        _fx1 = copy(x)
    end

    JacobianCache(x1,_fx,_fx1,fdtype,returntype;colorvec=colorvec,sparsity=sparsity)
end

function JacobianCache(
    x ,
    fx,
    fdtype     :: Union{Val{T1},Type{T1}} = Val(:forward),
    returntype :: Type{T2} = eltype(x);
    inplace    :: Union{Val{T3},Type{T3}} = Val(true),
    colorvec = 1:length(x),
    sparsity = nothing) where {T1,T2,T3}

    fdtype isa Type && (fdtype = fdtype())
    inplace isa Type && (inplace = inplace())
    if eltype(x) <: Real && fdtype==Val(:complex)
        x1  = false .* im .* x
    else
        x1 = copy(x)
    end

    if eltype(fx) <: Real && fdtype==Val(:complex)
        _fx = false .* im .* fx
    else
        _fx = copy(fx)
    end

    if fdtype==Val(:complex)
        _fx1  = nothing
    else
        _fx1 = copy(fx)
    end

    JacobianCache(x1,_fx,_fx1,fdtype,returntype;colorvec=colorvec,sparsity=sparsity)
end

function JacobianCache(
    x1 ,
    fx ,
    fx1,
    fdtype     :: Union{Val{T1},Type{T1}} = Val(:forward),
    returntype :: Type{T2} = eltype(fx);
    inplace    :: Union{Val{T3},Type{T3}} = Val(true),
    colorvec = 1:length(x1),
    sparsity = nothing) where {T1,T2,T3}

    fdtype isa Type && (fdtype = fdtype())
    inplace isa Type && (inplace = inplace())
    if fdtype==Val(:complex)
        !(returntype<:Real) && fdtype_error(returntype)

        if eltype(fx) <: Real
            _fx  = false .* im .* fx
        else
            _fx = fx
        end
        if eltype(x1) <: Real
            _x1  = false .* im .* x1
        else
            _x1 = x1
        end
    else
        _x1 = x1
        @assert eltype(fx) == T2
        @assert eltype(fx1) == T2
        _fx = fx
    end
    JacobianCache{typeof(_x1),typeof(_fx),typeof(fx1),typeof(colorvec),typeof(sparsity),fdtype,returntype}(_x1,_fx,fx1,colorvec,sparsity)
end

function _make_Ji(::SparseMatrixCSC, rows_index,cols_index,dx,colorvec,color_i,nrows,ncols)
    pick_inds = [i for i in 1:length(rows_index) if colorvec[cols_index[i]] == color_i]
    rows_index_c = rows_index[pick_inds]
    cols_index_c = cols_index[pick_inds]
    Ji = sparse(rows_index_c, cols_index_c, dx[rows_index_c],nrows,ncols)
    Ji
end

function _make_Ji(::AbstractArray, rows_index,cols_index,dx,colorvec,color_i,nrows,ncols)
    pick_inds = [i for i in 1:length(rows_index) if colorvec[cols_index[i]] == color_i]
    rows_index_c = rows_index[pick_inds]
    cols_index_c = cols_index[pick_inds]
    len_rows = length(pick_inds)
    unused_rows = setdiff(1:nrows,rows_index_c)
    perm_rows = sortperm(vcat(rows_index_c,unused_rows))
    cols_index_c = vcat(cols_index_c,zeros(Int,nrows-len_rows))[perm_rows]
    Ji = [j==cols_index_c[i] ? dx[i] : false for i in 1:nrows, j in 1:ncols]
    Ji
end

function _make_Ji(::SparseMatrixCSC, xtype, dx, color_i, nrows, ncols)
    Ji = sparse(1:nrows,fill(color_i,nrows),dx,nrows,ncols)
    Ji
end

function _make_Ji(::AbstractArray, xtype, dx, color_i, nrows, ncols)
    Ji = mapreduce(i -> i==color_i ? dx : zero(dx), hcat, 1:ncols)
    size(Ji) != (nrows, ncols) ? reshape(Ji, (nrows, ncols)) : Ji #branch when size(dx) == (1,) => size(Ji) == (1,) while size(J) == (1,1)
end

function finite_difference_jacobian(f, x,
    fdtype     = Val(:forward),
    returntype = eltype(x),
    f_in       = nothing;
    relstep=default_relstep(fdtype, eltype(x)),
    absstep=relstep,
    colorvec = 1:length(x),
    sparsity = nothing,
    jac_prototype = nothing,
    dir=true)

    if f_in isa Nothing
        fx = f(x)
    else
        fx = f_in
    end
    cache = JacobianCache(x, fx, fdtype, returntype)
    finite_difference_jacobian(f, x, cache, fx; relstep=relstep, absstep=absstep, colorvec=colorvec, sparsity=sparsity, jac_prototype=jac_prototype, dir=dir)
end

void_setindex!(args...) = (setindex!(args...); return)

function finite_difference_jacobian(
    f,
    x,
    cache::JacobianCache{T1,T2,T3,cType,sType,fdtype,returntype},
    f_in=nothing;
    relstep=default_relstep(fdtype, eltype(x)),
    absstep=relstep,
    colorvec = cache.colorvec,
    sparsity = cache.sparsity,
    jac_prototype = nothing,
    dir=true) where {T1,T2,T3,cType,sType,fdtype,returntype}

    x1, fx, fx1 = cache.x1, cache.fx, cache.fx1

    if !(f_in isa Nothing)
        vecfx = _vec(f_in)
    elseif fdtype == Val(:forward)
        vecfx = _vec(f(x))
    elseif fdtype == Val(:complex) && returntype <: Real
        vecfx = real(fx)
    else
        vecfx = _vec(fx)
    end
    vecx = _vec(x)
    vecx1 = _vec(x1)
    J = jac_prototype isa Nothing ? (sparsity isa Nothing ? Array{eltype(x),2}(undef, length(vecfx), 0) : zeros(eltype(x),size(sparsity))) : zero(jac_prototype)
    nrows, ncols = size(J)

    if !(sparsity isa Nothing)
        rows_index, cols_index = ArrayInterface.findstructralnz(sparsity)
        rows_index = [rows_index[i] for i in 1:length(rows_index)]
        cols_index = [cols_index[i] for i in 1:length(cols_index)]
    end

    if fdtype == Val(:forward)

        function calculate_Ji_forward(i)
            x_save = ArrayInterface.allowed_getindex(vecx, i)
            epsilon = compute_epsilon(Val(:forward), x_save, relstep, absstep, dir)
            _vecx1 = Base.setindex(vecx, x_save+epsilon, i)
            _x1 = reshape(_vecx1, axes(x))
            vecfx1 = _vec(f(_x1))
            dx = (vecfx1-vecfx) / epsilon
            return dx
        end

        if jac_prototype isa Nothing && sparsity isa Nothing
            J = mapreduce(calculate_Ji_forward, hcat, 1:maximum(colorvec))
            J = _mat(J)
        else
            @inbounds for color_i ∈ 1:maximum(colorvec)
                if sparsity isa Nothing
                    dx = calculate_Ji_forward(color_i)
                    J = J + _make_Ji(J, eltype(x), dx, color_i, nrows, ncols)
                else
                    tmp = norm(vecx .* (colorvec .== color_i))
                    epsilon = compute_epsilon(Val(:forward), sqrt(tmp), relstep, absstep, dir)
                    _vecx = @. vecx + epsilon * (colorvec == color_i)
                    _x = reshape(_vecx, axes(x))
                    vecfx1 = _vec(f(_x))
                    dx = (vecfx1-vecfx)/epsilon
                    Ji = _make_Ji(J,rows_index,cols_index,dx,colorvec,color_i,nrows,ncols)
                    J = J + Ji
                end
            end
        end
    elseif fdtype == Val(:central)

        function calculate_Ji_central(i)
            x1_save = ArrayInterface.allowed_getindex(vecx1,i)
            x_save = ArrayInterface.allowed_getindex(vecx,i)
            epsilon = compute_epsilon(Val(:forward), x1_save, relstep, absstep, dir)
            _vecx1 = Base.setindex(vecx1,x1_save+epsilon,i)
            _vecx = Base.setindex(vecx,x_save-epsilon,i)
            _x1 = reshape(_vecx1, axes(x))
            _x = reshape(_vecx, axes(x))
            vecfx1 = _vec(f(_x1))
            vecfx = _vec(f(_x))
            dx = (vecfx1-vecfx)/(2epsilon)
            return dx
        end

        if jac_prototype isa Nothing && sparsity isa Nothing
            J = mapreduce(calculate_Ji_central, hcat, 1:maximum(colorvec))
            J = _mat(J)
        else
            @inbounds for color_i ∈ 1:maximum(colorvec)
                if sparsity isa Nothing
                    dx = calculate_Ji_central(color_i)
                    J = J + _make_Ji(J, eltype(x), dx, color_i, nrows, ncols)
                else
                    tmp = norm(vecx1 .* (colorvec .== color_i))
                    epsilon = compute_epsilon(Val(:forward), sqrt(tmp), relstep, absstep, dir)
                    _vecx1 = @. vecx1 + epsilon * (colorvec == color_i)
                    _vecx = @. vecx - epsilon * (colorvec == color_i)
                    _x1 = reshape(_vecx1, axes(x))
                    _x = reshape(_vecx, axes(x))
                    vecfx1 = _vec(f(_x1))
                    vecfx = _vec(f(_x))
                    dx = (vecfx1-vecfx)/(2epsilon)
                    Ji = _make_Ji(J,rows_index,cols_index,dx,colorvec,color_i,nrows,ncols)
                    J = J + Ji
                end
            end
        end
    elseif fdtype == Val(:complex) && returntype <: Real
        epsilon = eps(eltype(x))

        function calculate_Ji_complex(i)
            x_save = ArrayInterface.allowed_getindex(vecx,i)
            _vecx = Base.setindex(complex.(vecx),x_save+im*epsilon,i)
            _x = reshape(_vecx, axes(x))
            vecfx = _vec(f(_x))
            dx = imag(vecfx)/epsilon
            return dx
        end

        if jac_prototype isa Nothing && sparsity isa Nothing
            J = mapreduce(calculate_Ji_complex, hcat, 1:maximum(colorvec))
            J = _mat(J)
        else
            @inbounds for color_i ∈ 1:maximum(colorvec)
                if sparsity isa Nothing
                    dx = calculate_Ji_complex(color_i)
                    J = J + _make_Ji(J, eltype(x), dx, color_i, nrows, ncols)
                else
                    _vecx = @. vecx + im * epsilon * (colorvec == color_i)
                    _x = reshape(_vecx, axes(x))
                    vecfx = _vec(f(_x))
                    dx = imag(vecfx)/epsilon
                    Ji = _make_Ji(J,rows_index,cols_index,dx,colorvec,color_i,nrows,ncols)
                    J = J + Ji
                end
            end
        end
    else
        fdtype_error(returntype)
    end
    J
end

function finite_difference_jacobian!(J,
    f,
    x,
    fdtype     = Val(:forward),
    returntype = eltype(x),
    f_in       = nothing;
    relstep=default_relstep(fdtype, eltype(x)),
    absstep=relstep,
    colorvec = 1:length(x),
    sparsity = ArrayInterface.has_sparsestruct(J) ? J : nothing)
    if f_in isa Nothing && fdtype == Val(:forward)
        if size(J,1) == length(x)
            fx = zero(x)
        else
            fx = zeros(returntype,size(J,1))
        end
        f(fx,x)
        cache = JacobianCache(x, fx, fdtype, returntype)
    elseif f_in isa Nothing
        cache = JacobianCache(x, fdtype, returntype)
    else
        cache = JacobianCache(x, f_in, fdtype, returntype)
    end
    finite_difference_jacobian!(J, f, x, cache, cache.fx; relstep=relstep, absstep=absstep, colorvec=colorvec, sparsity=sparsity)
end

function finite_difference_jacobian!(
    J,
    f,
    x,
    cache::JacobianCache{T1,T2,T3,cType,sType,fdtype,returntype},
    f_in = nothing;
    relstep = default_relstep(fdtype, eltype(x)),
    absstep = relstep,
    colorvec = cache.colorvec,
    sparsity = cache.sparsity,
    dir = true) where {T1,T2,T3,cType,sType,fdtype,returntype}

    m, n = size(J)
    _color = reshape(colorvec, axes(x)...)

    x1, fx, fx1 = cache.x1, cache.fx, cache.fx1
    copyto!(x1, x)
    vfx = _vec(fx)

    rows_index = nothing
    cols_index = nothing
    if _use_findstructralnz(sparsity)
        rows_index, cols_index = ArrayInterface.findstructralnz(sparsity)
    end

    if sparsity !== nothing
        fill!(J,false)
    end

    if fdtype == Val(:forward)
        vfx1 = _vec(fx1)

        if f_in isa Nothing
            f(fx, x)
            vfx = _vec(fx)
        else
            vfx = _vec(f_in)
        end

        @inbounds for color_i ∈ 1:maximum(colorvec)
            if sparsity isa Nothing
                x1_save = ArrayInterface.allowed_getindex(x1,color_i)
                epsilon = compute_epsilon(Val(:forward), x1_save, relstep, absstep, dir)
                ArrayInterface.allowed_setindex!(x1, x1_save + epsilon, color_i)
                f(fx1, x1)
                # J is dense, so either it is truly dense or this is the
                # compressed form of the coloring, so write into it.
                @. J[:,color_i] = (vfx1 - vfx) / epsilon
                # Now return x1 back to its original value
                ArrayInterface.allowed_setindex!(x1, x1_save, color_i)
            else # Perturb along the colorvec vector
                @. fx1 = x1 * (_color == color_i)
                tmp = norm(fx1)
                epsilon = compute_epsilon(Val(:forward), sqrt(tmp), relstep, absstep, dir)
                @. x1 = x1 + epsilon * (_color == color_i)
                f(fx1, x1)
                # J is a sparse matrix, so decompress on the fly
                @. vfx1 = (vfx1 - vfx) / epsilon
                if ArrayInterface.fast_scalar_indexing(x1)
                    _colorediteration!(J,sparsity,rows_index,cols_index,vfx1,colorvec,color_i,n)
                else
                    #=
                    J.nzval[rows_index] .+= (colorvec[cols_index] .== color_i) .* vfx1[rows_index]
                    or
                    J[rows_index, cols_index] .+= (colorvec[cols_index] .== color_i) .* vfx1[rows_index]
                    += means requires a zero'd out start
                    =#
                    if J isa SparseMatrixCSC
                        @. void_setindex!((J.nzval,), getindex((J.nzval,), rows_index) + (getindex((_color,), cols_index) == color_i) * getindex((vfx1,), rows_index), rows_index)
                    else
                        @. void_setindex!((J,), getindex((J,), rows_index, cols_index) + (getindex((_color,), cols_index) == color_i) * getindex((vfx1,), rows_index), rows_index, cols_index)
                    end
                end
                # Now return x1 back to its original value
                @. x1 = x1 - epsilon * (_color == color_i)
            end
        end #for ends here
    elseif fdtype == Val(:central)
        vfx1 = _vec(fx1)
        @inbounds for color_i ∈ 1:maximum(colorvec)
            if sparsity isa Nothing
                x_save = ArrayInterface.allowed_getindex(x, color_i)
                epsilon = compute_epsilon(Val(:central), x_save, relstep, absstep, dir)
                ArrayInterface.allowed_setindex!(x1, x_save + epsilon, color_i)
                f(fx1, x1)
                ArrayInterface.allowed_setindex!(x1, x_save - epsilon, color_i)
                f(fx, x1)
                @. J[:,color_i] = (vfx1 - vfx) / 2epsilon
                ArrayInterface.allowed_setindex!(x1, x_save, color_i)
            else # Perturb along the colorvec vector
                @. fx1 = x1 * (_color == color_i)
                tmp = norm(fx1)
                epsilon = compute_epsilon(Val(:central), sqrt(tmp), relstep, absstep, dir)
                @. x1 = x1 + epsilon * (_color == color_i)
                @. x  = x  - epsilon * (_color == color_i)
                f(fx1, x1)
                f(fx, x)
                @. vfx1 = (vfx1 - vfx) / 2epsilon
                if ArrayInterface.fast_scalar_indexing(x1)
                    _colorediteration!(J,sparsity,rows_index,cols_index,vfx1,colorvec,color_i,n)
                else
                    if J isa SparseMatrixCSC
                        @. void_setindex!((J.nzval,), getindex((J.nzval,), rows_index) + (getindex((_color,), cols_index) == color_i) * getindex((vfx1,), rows_index), rows_index)
                    else
                        @. void_setindex!((J,), getindex((J,), rows_index, cols_index) + (getindex((_color,), cols_index) == color_i) * getindex((vfx1,), rows_index), rows_index, cols_index)
                    end
                end
                @. x1 = x1 - epsilon * (_color == color_i)
                @. x  = x  + epsilon * (_color == color_i)
            end
        end
    elseif fdtype==Val(:complex) && returntype<:Real
        epsilon = eps(eltype(x))
        @inbounds for color_i ∈ 1:maximum(colorvec)
            if sparsity isa Nothing
                x1_save = ArrayInterface.allowed_getindex(x1, color_i)
                ArrayInterface.allowed_setindex!(x1, x1_save + im*epsilon, color_i)
                f(fx,x1)
                @. J[:,color_i] = imag(vfx) / epsilon
                ArrayInterface.allowed_setindex!(x1, x1_save,color_i)
            else # Perturb along the colorvec vector
                @. x1 = x1 + im * epsilon * (_color == color_i)
                f(fx,x1)
                @. vfx = imag(vfx) / epsilon
                if ArrayInterface.fast_scalar_indexing(x1)
                    _colorediteration!(J,sparsity,rows_index,cols_index,vfx,colorvec,color_i,n)
                else
                   if J isa SparseMatrixCSC
                        @. void_setindex!((J.nzval,), getindex((J.nzval,), rows_index) + (getindex((_color,), cols_index) == color_i) * getindex((vfx,),rows_index), rows_index)
                    else
                        @. void_setindex!((J,), getindex((J,), rows_index, cols_index) + (getindex((_color,), cols_index) == color_i) * getindex((vfx,), rows_index), rows_index, cols_index)
                    end
                end
                @. x1 = x1 - im * epsilon * (_color == color_i)
            end
        end
    else
        fdtype_error(returntype)
    end
    nothing
end

function resize!(cache::JacobianCache, i::Int)
    resize!(cache.x1,  i)
    resize!(cache.fx,  i)
    cache.fx1 != nothing && resize!(cache.fx1, i)
    cache.colorvec = 1:i
    nothing
end


text:
function Queue(conn::Connection, name::AbstractString)
    queue_handle_ref = Ref{Ptr{Cvoid}}()
    result = dpiConn_newQueue(conn.handle, name, C_NULL, queue_handle_ref)
    error_check(context(conn), result)
    return Queue(conn, queue_handle_ref[], name)
end

function enqueue(queue::Queue, msg::Message)
    result = dpiQueue_enqOne(queue.handle, msg.handle)
    error_check(context(queue), result)
    nothing
end

function enqueue(queue::Queue, msgs::Vector{T}) where {T<:Message}
    msg_handles = [ msg.handle for msg in msgs ]
    result = dpiQueue_enqMany(queue.handle, length(msg_handles), pointer(msg_handles))
    error_check(context(queue), result)
    nothing
end

function dequeue(queue::Queue) :: Message
    msg_props_handle_ref = Ref{Ptr{Cvoid}}()
    result = dpiQueue_deqOne(queue.handle, msg_props_handle_ref)
    error_check(context(queue), result)
    return Message(queue, msg_props_handle_ref[])
end

function dequeue(queue::Queue, max_msgs::Integer) :: Vector{Message}
    num_props = Ref{UInt32}(max_msgs)
    msg_handle_array_ref = Ref{Ptr{Cvoid}}()

    # dpiMsgProps **props
    # an array of references to message properties which will be populated upon successful completion of this function.
    # Each of these references should be released when they are no longer needed by calling the function dpiMsgProps_release().
    buffer = Vector{Ptr{Cvoid}}(undef, max_msgs)
    err = dpiQueue_deqMany(queue.handle, num_props, pointer(buffer))
    error_check(context(queue), err)

    if iszero(num_props[])
        return Vector{Message}()
    end

    result = Vector{Message}(undef, num_props[])
    for i in 1:num_props[]
        result[i] = Message(queue, buffer[i])
    end

    return result
end

ext: jl


text:
#Kozachenko-Leonenko estimator
import NearestNeighbors
import Distances
const NN = NearestNeighbors

function nn_entropy{T<:Real}(X::Array{T,2})
	dims = size(X,1)
	N = size(X,2)
    S = spherical_volume(dims)
	t = NN.NaiveNeighborTree(X,NN.Euclidean())
	ee = 0.0
	for i=1:N
		didx,dd = NN.nearest(t,X[:,i],1,i)
		ee += log2(dd[end])
	end
    H = (dims/N)*ee + log2(S*(N-1)/dims) + eulergamma/log(2)
end

function ksg_information{T<:Real}(X1::Array{T,2}, X2::Array{T,2};k::Integer=1)

	dims = size(X1,1)
	N1 = size(X1,2)
	N2 = size(X2,2)
	N = N1 + N2
	#contruct the full matrix
	X = cat(2,X1,X2)
    S = spherical_volume(dims)

	#compute the joint entropy as well as the conditional entropy using the Kozachenko-Leonenko estimator
	#get the distances
	t = NN.NaiveNeighborTree(X,NN.Euclidean())
	eex = 0.0
	eey = 0.0
	ee = 0.0
	for i=1:N
		didx,dd = NN.nearest(t,X[:,i],k,i)
		ee += log2(dd[end])
		#get the number of points within this distance for x and y separately
		qidx,bdd = NN.inball(t,X[:,i],dd[end],1)
		if i <= N1
			n = sum(qidx .<= N1) #x
			eex += digamma(n+1)
		else
			n = sum(qidx .> N1) #y
			eey += digamma(n+1)
		end
	end
	ee = digamma(k) - (eex/N1 + eey/N2) + digamma(N) #Eq. 8
    H = (dims/N)*ee + log2(S*(N-1)/dims) + eulergamma/log(2)
	return H
end

function information(::Type{NSBEntropy},X::Array{Int64,2}, S::Array{Int64,1})
	@assert size(X,ndims(X)) == length(S)
	K = (maximum(X)+1)^size(X,1)
	labels = sort(unique(S))
	Hs = zeros(length(labels))
	y = hash(X)
	n = counts(y)
	H = nsb_entropy(n,K)
	for (i,l) in enumerate(labels)
		n = counts(hash(y[S.==l]))
		Hs[i] = nsb_entropy(n,K)
	end
	return H,Hs
end

function GroupedTemporalEntropy(N::Array{Int64,2}, bins::Array{Float64,1}, trial_labels::Array{Int64,1},word_size::Integer)
	ulabels = sort(unique(trial_labels))
	nlabels = length(ulabels)
	H = NaN*zeros(size(N,1)-word_size)
	dH = NaN*zeros(size(N,1)-word_size)
	Hc = NaN*zeros(size(N,1)-word_size,nlabels)
	dHc = NaN*zeros(size(N,1)-word_size,nlabels)
	trials_per_label = zeros(Int64,nlabels)
	for (k,l) in enumerate(ulabels)
		trials_per_label[k] = sum(trial_labels.==l)
	end
	for i=1:size(N,1)-word_size
		y = Entropy.hash(N[i:i+word_size-1,:])
		n = Entropy.counts(y)
		try
			H[i],dH[i] = Entropy.nsb_entropy(n,maximum(y))
		catch ee
			println("Problem computing entropy for bin $i")
		end

		H[i] /= log(2)
		dH[i] /= log(2)
		for (k,l) in enumerate(ulabels)
			n = Entropy.counts(y[trial_labels.==l])
			try
				Hc[i,k],dHc[i,k] = Entropy.nsb_entropy(n,maximum(y))
				Hc[i,k] /= log(2)
				dHc[i,k] /= log(2)
			catch ee
				println("Problem computing entropy for bin $i location $l")
			end
		end
	end
	return GroupedTemporalEntropy(Hc,dHc,H,dH,bins,word_size, ulabels,trials_per_label)
end

function GroupedCountEntropy(N::Array{Int64,2}, bins::Array{Float64,1}, trial_labels::Array{Int64,1},word_size::Integer)
	ntrials = size(N,2)
	ulabels = sort(unique(trial_labels))
	nlabels = length(ulabels)
	H = NaN*zeros(size(N,1)-word_size)
	dH = NaN*zeros(size(N,1)-word_size)
	Hc = NaN*zeros(size(N,1)-word_size,nlabels)
	dHc = NaN*zeros(size(N,1)-word_size,nlabels)
	trials_per_label = zeros(Int64,nlabels)
	for (k,l) in enumerate(ulabels)
		trials_per_label[k] = sum(trial_labels.==l)
	end
	y = zeros(Int64,ntrials)
	for i=1:size(N,1)-word_size
		fill!(y,0)
		for j=1:ntrials
			for k=1:word_size
				y[j] += N[i+k-1,j]
			end
		end
		n = Entropy.counts(y)
		try
			H[i],dH[i] = Entropy.nsb_entropy(n,maximum(y))
		catch ee
			println("Problem computing entropy for bin $i")
		end
		H[i] /= log(2)
		dH[i] /= log(2)
		for (k,l) in enumerate(ulabels)
			n = Entropy.counts(y[trial_labels.==l])
			try
				Hc[i,k],dHc[i,k] = Entropy.nsb_entropy(n,maximum(y))
				Hc[i,k] /= log(2)
				dHc[i,k] /= log(2)
			catch ee
				println("Problem computing entropy for bin $i location $l")
			end
		end
	end
	return GroupedCountEntropy(Hc,dHc,H,dH,bins,word_size, ulabels,trials_per_label)
end

function GroupedPopulationEntropy(N::Array{Int64,3}, bins::Array{Float64,1}, trial_labels::Array{Int64,1})
	nbins,ntrials,ncells = size(N)
	ulabels = sort(unique(trial_labels))
	nlabels = length(ulabels)
	H = NaN*zeros(nbins)
	dH = NaN*zeros(H)
	Hc = NaN*zeros(nbins,nlabels)
	dHc = NaN*zeros(Hc)
	trials_per_label = zeros(Int64,nlabels)
	for (k,l) in enumerate(ulabels)
		trials_per_label[k] = sum(trial_labels.==l)
	end
	for i=1:size(N,1)
		y = Entropy.hash(squeeze(N[i,:,:],1)')
		n = Entropy.counts(y)
		try
			H[i],dH[i] = Entropy.nsb_entropy(n,maximum(y))
		catch ee
			println("Problem computing entropy for bin $i")
		end

		H[i] /= log(2)
		dH[i] /= log(2)
		for (k,l) in enumerate(ulabels)
			n = Entropy.counts(y[trial_labels.==l])
			try
				Hc[i,k],dHc[i,k] = Entropy.nsb_entropy(n,maximum(y))
				Hc[i,k] /= log(2)
				dHc[i,k] /= log(2)
			catch ee
				println("Problem computing entropy for bin $i location $l")
			end
		end
	end
	I,dI = compute_information(Hc,dHc, H, dH, trials_per_label)
	return GroupedPopulationEntropy(Hc,dHc,H,dH,I,dI,bins,ncells, ulabels,trials_per_label)
end

function compute_information(Hc::Array{Float64,2}, dHc::Array{Float64,2},H::Array{Float64,1},dH::Array{Float64,1}, trials_per_group::Array{Int64,1})
	ps = trials_per_group/sum(trials_per_group)
	_hc = zeros(size(Hc,1))
	_dhc = zeros(size(Hc,1))
	for i=1:size(Hc,2)
		for j=1:size(Hc,1)
			_hc[j] += ps[i]*Hc[j,i]
			_dhc[j] += ps[i]*ps[i]*dHc[j,i]*dHc[j,i] #variance
		end
	end
	I = H - _hc
	dI = sqrt(dH.*dH + _dhc)
	return I, dI
end


text:
# PROJECT_SOC    Project a point onto the second-order cone.
#
#   Suppose v is an (n+1)-dimensional vector, so v = [t;v0], where t is a
#   scalar and v0 is n-dimensional. Then project_soc(v) is the projection of
#   (t,v0) onto the second-order cone; the result is also (n+1)-dimensional.

function project_soc(v)
    # assume v = (t, v0)
    nv = norm(v[2:end])
    x  = nan(size(v))

    if nv <= -v[1]
        x = zeros(size(v))
    elseif nv <= v[1]
        x = v
    else
        r = 0.5 * (1 + v[1] / nv)
        x[1] = r * nv
        x[2:end] = r * v[2:end]
    end

    return x
end


text:
mutable struct Trace{T}
    save_every::Int
    iteration::Int
    way_x::Vector{T}
    way_l::Vector{Float64}
    isshow::Bool
end

Trace(u0::T, l::Float64; save_every::Int=1, isshow::Bool=false) where T = Trace(save_every, 0, [copy(u0)], [l], isshow)
Trace(prob; save_every::Int=1, isshow::Bool=false) = Trace(prob.u0, prob.f(prob.u0, prob.p); save_every, isshow)

Base.show(io::IO, tr::Trace) = show(io, "i = $(tr.iteration), loss = $(tr.way_l[end])")

function update_trace!(tr::Trace{T}, x::T, l) where T
    tr.iteration += 1
    if tr.iteration % tr.save_every == 0
        push!(tr.way_x, copy(x))
        push!(tr.way_l, l)
    end
end

function (trace::Trace)(x, l)
    update_trace!(trace, x, l)
    trace.isshow && println(trace)
    false
end

Base.iterate(t::Trace, state=1) = state > t.iteration + 1 ? nothing : (t.way_x[state], state + 1)
Base.eltype(::Type{Trace{T}}) where T = T
Base.length(t::Trace) = t.iteration + 1
Base.getindex(t::Trace, i) = t.way_x[i]
Base.firstindex(t::Trace) = 1
Base.lastindex(t::Trace) = t.iteration + 1


text:
# Necessary to have initialize set u_modified to false if all don't do anything
# otherwise unnecessary save
INITIALIZE_DEFAULT(cb,u,t,integrator) = u_modified!(integrator, false)

struct ContinuousCallback{F1,F2,F3,F4,T,T2,I} <: AbstractContinuousCallback
  condition::F1
  affect!::F2
  affect_neg!::F3
  initialize::F4
  idxs::I
  rootfind::Bool
  interp_points::Int
  save_positions::BitArray{1}
  abstol::T
  reltol::T2
  ContinuousCallback(condition::F1,affect!::F2,affect_neg!::F3,
                     initialize::F4,idxs::I,rootfind,
                     interp_points,save_positions,abstol::T,reltol::T2) where {F1,F2,F3,F4,T,T2,I} =
                       new{F1,F2,F3,F4,T,T2,I}(condition,
                                               affect!,affect_neg!,
                                               initialize,idxs,rootfind,interp_points,
                                               BitArray(collect(save_positions)),
                                               abstol,reltol)
end

ContinuousCallback(condition,affect!,affect_neg!;
                   initialize = INITIALIZE_DEFAULT,
                   idxs = nothing,
                   rootfind=true,
                   save_positions=(true,true),
                   interp_points=10,
                   abstol=10eps(),reltol=0) = ContinuousCallback(
                              condition,affect!,affect_neg!,initialize,
                              idxs,
                              rootfind,interp_points,
                              save_positions,abstol,reltol)

function ContinuousCallback(condition,affect!;
                   initialize = INITIALIZE_DEFAULT,
                   idxs = nothing,
                   rootfind=true,
                   save_positions=(true,true),
                   affect_neg! = affect!,
                   interp_points=10,
                   abstol=10eps(),reltol=0)

 ContinuousCallback(
            condition,affect!,affect_neg!,initialize,idxs,
            rootfind,interp_points,
            collect(save_positions),abstol,reltol)

end

struct DiscreteCallback{F1,F2,F3} <: AbstractDiscreteCallback
  condition::F1
  affect!::F2
  initialize::F3
  save_positions::BitArray{1}
  DiscreteCallback(condition::F1,affect!::F2,
                   initialize::F3,save_positions) where {F1,F2,F3} = new{F1,F2,F3}(condition,
                                                                                affect!,initialize,
                                                                                BitArray(collect(save_positions)))
end
DiscreteCallback(condition,affect!;
        initialize = INITIALIZE_DEFAULT,save_positions=(true,true)) = DiscreteCallback(condition,affect!,initialize,save_positions)

# DiscreteCallback(condition,affect!,save_positions) = DiscreteCallback(condition,affect!,save_positions)

struct CallbackSet{T1<:Tuple,T2<:Tuple} <: DECallback
  continuous_callbacks::T1
  discrete_callbacks::T2
end

CallbackSet(callback::AbstractDiscreteCallback) = CallbackSet((),(callback,))
CallbackSet(callback::AbstractContinuousCallback) = CallbackSet((callback,),())
CallbackSet() = CallbackSet((),())
CallbackSet(cb::Nothing) = CallbackSet()

# For Varargs, use recursion to make it type-stable
CallbackSet(callbacks::Union{DECallback,Nothing}...) = CallbackSet(split_callbacks((), (), callbacks...)...)

@inline split_callbacks(cs, ds) = cs, ds
@inline split_callbacks(cs, ds, c::Nothing, args...) = split_callbacks(cs, ds, args...)
@inline split_callbacks(cs, ds, c::AbstractContinuousCallback, args...) = split_callbacks((cs..., c), ds, args...)
@inline split_callbacks(cs, ds, d::AbstractDiscreteCallback, args...) = split_callbacks(cs, (ds..., d), args...)
@inline function split_callbacks(cs, ds, d::CallbackSet, args...)
  split_callbacks((cs...,d.continuous_callbacks...), (ds..., d.discrete_callbacks...), args...)
end

# Recursively apply initialize! and return whether any modified u
function initialize!(cb::CallbackSet,u,t,integrator::DEIntegrator)
  initialize!(u,t,integrator,false,cb.continuous_callbacks...,cb.discrete_callbacks...)
end
initialize!(cb::CallbackSet{Tuple{},Tuple{}},u,t,integrator::DEIntegrator) = false
function initialize!(u,t,integrator::DEIntegrator,any_modified::Bool,
                     c::DECallback,cs::DECallback...)
  c.initialize(c,u,t,integrator)
  initialize!(u,t,integrator,any_modified || integrator.u_modified,cs...)
end
function initialize!(u,t,integrator::DEIntegrator,any_modified::Bool,
                     c::DECallback)
  c.initialize(c,u,t,integrator)
  any_modified || integrator.u_modified
end

# Helpers
Base.isempty(cb::CallbackSet) = isempty(cb.continuous_callbacks) && isempty(cb.discrete_callbacks)
Base.isempty(cb::AbstractContinuousCallback) = false
Base.isempty(cb::AbstractDiscreteCallback) = false

#======================================================#
# Callback handling
#======================================================#

function get_tmp(integrator::DEIntegrator, callback)
  _tmp = get_tmp_cache(integrator)
  _tmp === nothing && return nothing
  _cache = first(_tmp)
  if callback.idxs isa Nothing
    tmp = _cache
  elseif !(callback.idxs isa Number)
    tmp = @view _cache[callback.idxs]
  else
    tmp = nothing
  end
  return tmp
end

function get_condition(integrator::DEIntegrator, callback, abst)
  tmp = get_tmp(integrator, callback)
  ismutable = !(tmp === nothing)
  ismutable && !(typeof(callback.idxs) isa Number) ? integrator(tmp,abst,Val{0},idxs=callback.idxs) :
                                                     tmp = integrator(abst,Val{0},idxs=callback.idxs)
  return callback.condition(tmp,abst,integrator)
end

# Use Recursion to find the first callback for type-stability

# Base Case: Only one callback
function find_first_continuous_callback(integrator, callback::AbstractContinuousCallback)
  (find_callback_time(integrator,callback,1)...,1,1)
end

# Starting Case: Compute on the first callback
function find_first_continuous_callback(integrator, callback::AbstractContinuousCallback, args...)
  find_first_continuous_callback(integrator,find_callback_time(integrator,callback,1)...,1,1,args...)
end

function find_first_continuous_callback(integrator,tmin::Number,upcrossing::Number,
                                        event_occured::Bool,idx::Int,counter::Int,
                                        callback2)
  counter += 1 # counter is idx for callback2.
  tmin2,upcrossing2,event_occurred2 = find_callback_time(integrator,callback2,counter)

  if event_occurred2 && (tmin2 < tmin || !event_occured)
    return tmin2,upcrossing2,true,counter,counter
  else
    return tmin,upcrossing,event_occured,idx,counter
  end
end

function find_first_continuous_callback(integrator,tmin::Number,upcrossing::Number,event_occured::Bool,idx::Int,counter::Int,callback2,args...)
  find_first_continuous_callback(integrator,find_first_continuous_callback(integrator,tmin,upcrossing,event_occured,idx,counter,callback2)...,args...)
end

@inline function determine_event_occurance(integrator,callback,counter)
  event_occurred = false
  if callback.interp_points!=0
    addsteps!(integrator)
  end
  Θs = range(typeof(integrator.t)(0), stop=typeof(integrator.t)(1), length=callback.interp_points)
  interp_index = 0
  # Check if the event occured
  if typeof(callback.idxs) <: Nothing
    previous_condition = callback.condition(integrator.uprev,integrator.tprev,integrator)
  else
    @views previous_condition = callback.condition(integrator.uprev[callback.idxs],integrator.tprev,integrator)
  end

  if integrator.event_last_time == counter && ODE_DEFAULT_NORM(previous_condition,integrator.t) < 100ODE_DEFAULT_NORM(integrator.last_event_error,integrator.t)

    # If there was a previous event, utilize the derivative at the start to
    # chose the previous sign. If the derivative is positive at tprev, then
    # we treat the value as positive, and derivative is negative then we
    # treat the value as negative, reguardless of the postiivity/negativity
    # of the true value due to it being =0 sans floating point issues.

    # Only due this if the discontinuity did not move it far away from an event
    # Since near even we use direction instead of location to reset

    if callback.interp_points==0
      addsteps!(integrator)
    end

    abst = integrator.tprev+sign(integrator.dt)*100*eps(integrator.tprev) # Evaluate condition slightly in future
    tmp_condition = get_condition(integrator, callback, abst)
    prev_sign = tmp_condition > previous_condition ? 1.0 : -1.0
  else
    prev_sign = sign(previous_condition)
  end

  prev_sign_index = 1
  if typeof(callback.idxs) <: Nothing
    next_sign = sign(callback.condition(integrator.u,integrator.t,integrator))
  else
    @views next_sign = sign(callback.condition(integrator.u[callback.idxs],integrator.t,integrator))
  end

  if ((prev_sign<0 && !(typeof(callback.affect!)<:Nothing)) || (prev_sign>0 && !(typeof(callback.affect_neg!)<:Nothing))) && prev_sign*next_sign<=0
    event_occurred = true
    interp_index = callback.interp_points
  elseif callback.interp_points!=0 && !isdiscrete(integrator.alg) # Use the interpolants for safety checking
    for i in 2:length(Θs)
      abst = integrator.tprev+integrator.dt*Θs[i]
      new_sign = get_condition(integrator, callback, abst)
      if ((prev_sign<0 && !(typeof(callback.affect!)<:Nothing)) || (prev_sign>0 && !(typeof(callback.affect_neg!)<:Nothing))) && prev_sign*new_sign<0
        event_occurred = true
        interp_index = i
        break
      else
        prev_sign_index = i
      end
    end
  end

  event_occurred,interp_index,Θs,prev_sign,prev_sign_index
end

function find_callback_time(integrator,callback,counter)
  event_occurred,interp_index,Θs,prev_sign,prev_sign_index = determine_event_occurance(integrator,callback,counter)
  if event_occurred
    if typeof(callback.condition) <: Nothing
      new_t = zero(typeof(integrator.t))
    else
      if callback.interp_points!=0
        top_Θ = Θs[interp_index] # Top at the smallest
        bottom_θ = Θs[prev_sign_index]
      else
        top_Θ = typeof(integrator.t)(1)
        bottom_θ = typeof(integrator.t)(0)
      end
      if callback.rootfind && !isdiscrete(integrator.alg)
        zero_func = (Θ) -> begin
          abst = integrator.tprev+integrator.dt*Θ
          return get_condition(integrator, callback, abst)
        end
        if zero_func(top_Θ) == 0
          Θ = top_Θ
        else
          if integrator.event_last_time == counter &&
            abs(zero_func(bottom_θ)) < 100abs(integrator.last_event_error) &&
            prev_sign_index == 1

            # Determined that there is an event by derivative
            # But floating point error may make the end point negative

            sign_top = sign(zero_func(top_Θ))
            bottom_θ += 2eps(typeof(bottom_θ))
            iter = 1
            while sign(zero_func(bottom_θ)) == sign_top && iter < 12
              bottom_θ *= 5
              iter += 1
            end
            iter == 12 && error("Double callback crossing floating pointer reducer errored. Report this issue.")
          end
          Θ = prevfloat(find_zero(zero_func, (bottom_θ,top_Θ), Roots.AlefeldPotraShi(), atol = callback.abstol/100))
          integrator.last_event_error = ODE_DEFAULT_NORM(zero_func(Θ),integrator.t+integrator.dt*Θ)
        end
        #Θ = prevfloat(...)
        # prevfloat guerentees that the new time is either 1 floating point
        # numbers just before the event or directly at zero, but not after.
        # If there's a barrier which is never supposed to be crossed,
        # then this will ensure that
        # The item never leaves the domain. Otherwise Roots.jl can return
        # a float which is slightly after, making it out of the domain, causing
        # havoc.
        new_t = integrator.dt*Θ
      elseif interp_index != callback.interp_points && !isdiscrete(integrator.alg)
        new_t = integrator.dt*Θs[interp_index]
      else
        # If no solve and no interpolants, just use endpoint
        new_t = integrator.dt
      end
    end
  else
    new_t = zero(typeof(integrator.t))
  end

  new_t,prev_sign,event_occurred
end

function apply_callback!(integrator,callback::ContinuousCallback,cb_time,prev_sign)
  if cb_time == zero(typeof(integrator.t))
    error("Event repeated at the same time. Please report this error")
  end
  change_t_via_interpolation!(integrator,integrator.tprev+cb_time)

  # handle saveat
  _, savedexactly = savevalues!(integrator)
  saved_in_cb = true

  @inbounds if callback.save_positions[1]
    # if already saved then skip saving
    savedexactly || savevalues!(integrator,true)
  end

  integrator.u_modified = true

  if prev_sign < 0
    if typeof(callback.affect!) <: Nothing
      integrator.u_modified = false
    else
      callback.affect!(integrator)
    end
  elseif prev_sign > 0
    if typeof(callback.affect_neg!) <: Nothing
      integrator.u_modified = false
    else
      callback.affect_neg!(integrator)
    end
  end

  if integrator.u_modified
    reeval_internals_due_to_modification!(integrator)
    @inbounds if callback.save_positions[2]
      savevalues!(integrator,true)
      saved_in_cb = true
    end
    return true,saved_in_cb
  end
  false,saved_in_cb
end

#Base Case: Just one
@inline function apply_discrete_callback!(integrator,callback::DiscreteCallback)
  saved_in_cb = false
  if callback.condition(integrator.u,integrator.t,integrator)
    # handle saveat
    _, savedexactly = savevalues!(integrator)
    saved_in_cb = true
    @inbounds if callback.save_positions[1]
      # if already saved then skip saving
      savedexactly || savevalues!(integrator,true)
    end
    integrator.u_modified = true
    callback.affect!(integrator)
    @inbounds if callback.save_positions[2]
      savevalues!(integrator,true)
      saved_in_cb = true
    end
  end
  integrator.u_modified,saved_in_cb
end

#Starting: Get bool from first and do next
@inline function apply_discrete_callback!(integrator,callback::DiscreteCallback,args...)
  apply_discrete_callback!(integrator,apply_discrete_callback!(integrator,callback)...,args...)
end

@inline function apply_discrete_callback!(integrator,discrete_modified::Bool,saved_in_cb::Bool,callback::DiscreteCallback,args...)
  bool,saved_in_cb2 = apply_discrete_callback!(integrator,apply_discrete_callback!(integrator,callback)...,args...)
  discrete_modified || bool, saved_in_cb || saved_in_cb2
end

@inline function apply_discrete_callback!(integrator,discrete_modified::Bool,saved_in_cb::Bool,callback::DiscreteCallback)
  bool,saved_in_cb2 = apply_discrete_callback!(integrator,callback)
  discrete_modified || bool, saved_in_cb || saved_in_cb2
end


text:
### A Pluto.jl notebook ###
# v0.15.1

using Markdown
using InteractiveUtils

# ╔═╡ 8f9a3074-f92f-11eb-33dd-41eeab6f1ece
using Pkg

# ╔═╡ 15d15589-759d-4490-91dd-033f0eff550e
Pkg.add("LightGraphs")

# ╔═╡ 50dbc3e2-cbdf-470a-b71b-488a8423877b
begin
	Pkg.add("PyPlot")
end

# ╔═╡ c3cd9a02-a051-49be-a558-b57ecedbad5d
using GraphRecipes

# ╔═╡ 33e4b886-01f1-4438-889b-f4ebdefed9ce
using LightGraphs

# ╔═╡ 79a55dae-cd3e-4f90-9973-b710c744b6a6
using GraphPlot

# ╔═╡ 78cc9dcf-50d8-4284-9669-7323a05cd8f1
using Plots

# ╔═╡ 2ac3ef7b-7a9c-42c9-ae49-694094cbba10
#pgfplotsx()

# ╔═╡ 29f9ee4a-40a1-406a-bcdd-d763165a4e48
pyplot()

# ╔═╡ 2b807a12-be01-4dee-82f0-55e8e538406a
positions = [
	0 0;
	-1 1;
	1 1
]

# ╔═╡ 0bca08d0-20d5-4b2f-b682-7e4b3f7c9c78
n=3

# ╔═╡ 64cd57ea-3146-4d05-9e40-d9111971c502
pp = Dict(
	:a => Dict(
		:b => 1,
		:c => 2
	)
)

# ╔═╡ 7832eb0c-35fa-49c7-8ec4-e75059314d0c
const A = Float64[ rand() < 0.5 ? 0 : rand() for i=1:n, j=1:n]

# ╔═╡ abb9e1b1-20b2-4f6b-a374-a1d9bd67ae93
names = ["A", "B", "C"]

# ╔═╡ a046ced9-90b7-47c0-bf42-f5b1c50d2ca1
begin
	nvertices = 15
	G1 = Graph(nvertices)
	labels = [
		"¹H+¹H → ²H + e",
		"¹H+¹H+e → ²H",
		"²H+H → ³H",
		"³He + ³He → ⁴He + 2 ¹H",
		"³He + ⁴He → ⁷Be",
		"³He + ⁴He → ⁷Be",
		"⁷Be + e → ⁷Li + ν_e",
		"⁷Be + ¹H → ⁸B + γ",
		"⁷Li + ¹H -> ⁴He + ⁴He",
		"⁸B -> ⁷Be +e + ν",
		"⁸Be -> ⁴He + ⁴He",
	]
	#labels = ["1a", "b", "c²", "d", "e", "f", "g", "h", "i", "j", "k"]
	edgelabels = Dict()
	x = [
		-1,
		1,
		0,
		-2,
		0,
		2,
		0,
		2,
		0,
		2,
		2
	]
	y = [
		8,
		8,
		6,
		3,
		2,
		3,
		1,
		2,
		0,
		1,
		0
	]

	add_edge!(G1, 1, 3)
	edgelabels[(1,3)] = "99,75%         "
	add_edge!(G1, 2, 3)
	edgelabels[(2,3)] = "    0,25%"

	add_edge!(G1, 3, 4)
	edgelabels[(3,4)] = " I\n91%"
	add_edge!(G1, 3, 5)
	edgelabels[(3,5)] = "II\n9%"
	add_edge!(G1, 3, 6)
	edgelabels[(3,6)] = "III\n0.1% "

	add_edge!(G1, 5, 7)
	add_edge!(G1, 6, 8)

	add_edge!(G1, 7, 9)
	add_edge!(G1, 8, 10)

	add_edge!(G1, 10, 11)
end

# ╔═╡ 6d50a9a9-d650-4997-9257-0f4da6136914
plot = graphplot(
	G1,
	names=labels,
	fontsize=10,
	nodeshape=:rect,
	nodesize=0.0,
	x=x,#.*400,
	y=y,#.*400,
	size=(800,1000),
	dpi=100,
	curvature_scalar=0,
	nodecolor=:white,
	shorten=0.1,
	edgelabel=edgelabels,
	edge_label_box=true,
	edgelabel_offset=0.2,
)

# ╔═╡ b7c805b7-cd32-4248-bb8a-5ee067b1183e
typeof(G1)

# ╔═╡ f431b7f9-4c36-4435-85fc-3102e1f2dea7
labels

# ╔═╡ ee189abc-2707-4281-b5cb-73b2890ae023
savefig(plot, "pp.png")

# ╔═╡ ea90e75c-ed5b-48b2-8e7e-f4af4d0d666c
savefig(plot, "pp.svg")

# ╔═╡ 153589fc-5060-4fc2-b6b5-1834450b6847
nv(G1)

# ╔═╡ 7bc0c196-f54b-4c11-8f6e-54ad4636df3d
# set positions manually to avoid randomness!

# ╔═╡ 25cebcd0-a75a-4373-96a6-a0b778fb4f90
savefig(plot, "test.svg")

# ╔═╡ 1d325978-a695-4db9-a55a-353e8b216430
length(labels)

# ╔═╡ cf1815ee-4219-4cac-8225-200b37d13a08
ne(G1)

# ╔═╡ 18247747-32a6-4eb0-87cc-1595c2582e0e
begin
	G2 = Graph(10)
	add_edge!(G2, 1, 3)
	add_edge!(G2, 1, 2)

	# first fusion
	add_edge!(G2, 2, 4)
	add_edge!(G2, 3, 4)

	# left branch
	add_edge!(G2, 4, 5)

	# middle branch
	#add_edge!(G2, 4, 6)

	# right branch
	#add_edge!(G2, 4, 7)
	#x = [0, -1, 1, 0, 0]
	#y = [3, 2, 2, 1, 0]
	#z = [0, 0, 0, 0, 0]
end

# ╔═╡ 054e3631-29b9-45b3-8d62-0ecc5646fe87
graphplot(G2, method=:tree)#, x=x, y=y, z=z)

# ╔═╡ e60e8737-1494-4b34-bc48-ca843c0f17f9
gplot(G2)

# ╔═╡ edf2e096-01e2-4697-8a87-05c188e5b03b

# ╔═╡ 6f7611e6-210a-4b64-b44a-e649b5fb0fa4
a  = Dict()

# ╔═╡ 0e915d78-676f-4e62-b524-72817c48c710
a["1"] = 2

# ╔═╡ 1b44df22-ed90-43ef-a6d5-92f9609b7187
a

# ╔═╡ Cell order:
# ╠═8f9a3074-f92f-11eb-33dd-41eeab6f1ece
# ╠═15d15589-759d-4490-91dd-033f0eff550e
# ╠═c3cd9a02-a051-49be-a558-b57ecedbad5d
# ╠═33e4b886-01f1-4438-889b-f4ebdefed9ce
# ╠═79a55dae-cd3e-4f90-9973-b710c744b6a6
# ╠═78cc9dcf-50d8-4284-9669-7323a05cd8f1
# ╠═2ac3ef7b-7a9c-42c9-ae49-694094cbba10
# ╠═50dbc3e2-cbdf-470a-b71b-488a8423877b
# ╠═29f9ee4a-40a1-406a-bcdd-d763165a4e48
# ╠═2b807a12-be01-4dee-82f0-55e8e538406a
# ╠═0bca08d0-20d5-4b2f-b682-7e4b3f7c9c78
# ╠═64cd57ea-3146-4d05-9e40-d9111971c502
# ╠═7832eb0c-35fa-49c7-8ec4-e75059314d0c
# ╠═abb9e1b1-20b2-4f6b-a374-a1d9bd67ae93
# ╠═a046ced9-90b7-47c0-bf42-f5b1c50d2ca1
# ╠═6d50a9a9-d650-4997-9257-0f4da6136914
# ╠═b7c805b7-cd32-4248-bb8a-5ee067b1183e
# ╠═f431b7f9-4c36-4435-85fc-3102e1f2dea7
# ╠═ee189abc-2707-4281-b5cb-73b2890ae023
# ╠═ea90e75c-ed5b-48b2-8e7e-f4af4d0d666c
# ╠═153589fc-5060-4fc2-b6b5-1834450b6847
# ╠═7bc0c196-f54b-4c11-8f6e-54ad4636df3d
# ╠═25cebcd0-a75a-4373-96a6-a0b778fb4f90
# ╠═1d325978-a695-4db9-a55a-353e8b216430
# ╠═cf1815ee-4219-4cac-8225-200b37d13a08
# ╠═18247747-32a6-4eb0-87cc-1595c2582e0e
# ╠═054e3631-29b9-45b3-8d62-0ecc5646fe87
# ╠═e60e8737-1494-4b34-bc48-ca843c0f17f9
# ╠═edf2e096-01e2-4697-8a87-05c188e5b03b
# ╠═6f7611e6-210a-4b64-b44a-e649b5fb0fa4
# ╠═0e915d78-676f-4e62-b524-72817c48c710
# ╠═1b44df22-ed90-43ef-a6d5-92f9609b7187


text:
push!(LOAD_PATH, "../src/")
using Documenter , UHDBindings

makedocs(sitename="UHDBindings.jl",
		 format = Documenter.HTML(),
		 pages    = Any[
						"Introduction to UHDBindings"   => "index.md",
						"Function list"         => "base.md",
						"Examples"              => Any[
														 "Examples/example_setup.md"
														 "Examples/example_parameters.md"
														 "Examples/example_benchmark.md"
														 ],
						],
		 );

#makedocs(sitename="My Documentation", format = Documenter.HTML(prettyurls = false))

deploydocs(
    repo = "github.com/JuliaTelecom/UHDBindings.jl",
)


text:
"""
    radial_dimension(setup)

Get dimension of radial direction.
"""
function radial_dimension end

radial_dimension(::PlateSetup) = 1
radial_dimension(::CylinderSetup) = 2
radial_dimension(::SphereSetup) = 3


text:
module TestMulitDim

using Test
using ThreadPools

include("util.jl")

@testset "generators" begin

    @testset "tmap" begin
        N = 2 * Threads.nthreads()
        primary = Threads.nthreads() == 1
        fn! = (x) -> begin
            Threads.threadid() == 1 && (primary = true)
            x.data
        end
        @test tmap(fn!, (TestObj(x) for x in 1:N)) == collect(1:N)
        @test primary
    end

    @testset "bmap" begin
        N = 2 * Threads.nthreads()
        fn! = (x) -> begin
            Threads.nthreads() == 1 || Threads.threadid() == 1 && error("Task on primary")
            x.data
        end
        @test bmap(fn!, (TestObj(x) for x in 1:N)) == collect(1:N)
        @inferred bmap(fn!, (TestObj(x) for x in 1:N))
    end

    @testset "qmap" begin
        N = 2 * Threads.nthreads()
        primary = Threads.nthreads() == 1
        fn! = (x) -> begin
            Threads.threadid() == 1 && (primary = true)
            x.data
        end
        @test qmap(fn!, (TestObj(x) for x in 1:N)) == collect(1:N)
        @test primary
        @inferred qmap(fn!, (TestObj(x) for x in 1:N))
    end

    @testset "qbmap" begin
        N = 2 * Threads.nthreads()
        fn! = (x) -> begin
            Threads.nthreads() == 1 || Threads.threadid() == 1 && error("Task on primary")
            x.data
        end
        @test qbmap(fn!, (TestObj(x) for x in 1:N)) == collect(1:N)
        @inferred qbmap(fn!, (TestObj(x) for x in 1:N))
    end

    @testset "logtmap" begin
        N = 2 * Threads.nthreads()
        primary = Threads.nthreads() == 1
        fn! = (x) -> begin
            Threads.threadid() == 1 && (primary = true)
            x.data
        end
        p, r = logtmap(fn!, (TestObj(x) for x in 1:N))
        @test r == collect(1:N)
        @test length(p.recs) == N*2
        @test primary
        @inferred logtmap(fn!, (TestObj(x) for x in 1:N))
    end

    @testset "logbmap" begin
        N = 2 * Threads.nthreads()
        fn! = (x) -> begin
            Threads.nthreads() == 1 || Threads.threadid() == 1 && error("Task on primary")
            x.data
        end
        p, r = logbmap(fn!, (TestObj(x) for x in 1:N))
        @test r == collect(1:N)
        @test length(p.recs) == N*2
        @inferred logbmap(fn!, (TestObj(x) for x in 1:N))
    end

    @testset "logqmap" begin
        N = 2 * Threads.nthreads()
        primary = Threads.nthreads() == 1
        fn! = (x) -> begin
            Threads.threadid() == 1 && (primary = true)
            x.data
        end
        p, r = logqmap(fn!, (TestObj(x) for x in 1:N))
        @test r == collect(1:N)
        @test length(p.recs) == N*2
        @test primary
        @inferred logqmap(fn!, (TestObj(x) for x in 1:N))
    end

    @testset "logqbmap" begin
        N = 2 * Threads.nthreads()
        fn! = (x) -> begin
            Threads.nthreads() == 1 || Threads.threadid() == 1 && error("Task on primary")
            x.data
        end
        p, r = logqbmap(fn!, (TestObj(x) for x in 1:N))
        @test r == collect(1:N)
        @test length(p.recs) == N*2
        @inferred logqbmap(fn!, (TestObj(x) for x in 1:N))
    end

    @testset "bthreads" begin
        function test()
            N = 2 * Threads.nthreads()
            @bthreads for x in (TestObj(x) for x in 1:N)
                Threads.nthreads() == 1 || Threads.threadid() == 1 && error("Task on primary")
                x.data += 1
            end
            @test [x.data for x in (TestObj(x) for x in 1:N)] == collect(1:N)
        end
        test()
        test()
    end

    @testset "qthreads" begin
        function test()
            N = 2 * Threads.nthreads()
            primary = Threads.nthreads() == 1
            @qthreads for x in (TestObj(x) for x in 1:N)
                Threads.threadid() == 1 && (primary = true)
                x.data += 1
                sleep(0.001)
            end
            @test [x.data for x in (TestObj(x) for x in 1:N)] == collect(1:N)
            @test primary
        end
        test()
        test()
    end

end

@testset "2d array" begin

    @testset "tmap" begin
        N = 2 * Threads.nthreads()
        objs = [collect(TestObj(x) for x in 1:N/2) collect(TestObj(x) for x in N/2+1:N)]
        primary = Threads.nthreads() == 1
        fn! = (x) -> begin
            Threads.threadid() == 1 && (primary = true)
            x.data
        end
        @test tmap(fn!, objs) == [collect(1:N/2) collect(N/2+1:N)]
        @test primary
        @inferred tmap(fn!, objs)
    end

    @testset "bmap" begin
        N = 2 * Threads.nthreads()
        objs = [collect(TestObj(x) for x in 1:N/2) collect(TestObj(x) for x in N/2+1:N)]
        fn! = (x) -> begin
            Threads.nthreads() == 1 || Threads.threadid() == 1 && error("Task on primary")
            x.data
        end
        @test bmap(fn!, objs) == [collect(1:N/2) collect(N/2+1:N)]
        @inferred bmap(fn!, objs)
    end

    @testset "qmap" begin
        N = 2 * Threads.nthreads()
        objs = [collect(TestObj(x) for x in 1:N/2) collect(TestObj(x) for x in N/2+1:N)]
        primary = Threads.nthreads() == 1
        fn! = (x) -> begin
            Threads.threadid() == 1 && (primary = true)
            x.data
        end
        @test qmap(fn!, objs) == [collect(1:N/2) collect(N/2+1:N)]
        @test primary
        @inferred qmap(fn!, objs)
    end

    @testset "qbmap" begin
        N = 2 * Threads.nthreads()
        objs = [collect(TestObj(x) for x in 1:N/2) collect(TestObj(x) for x in N/2+1:N)]
        fn! = (x) -> begin
            Threads.nthreads() == 1 || Threads.threadid() == 1 && error("Task on primary")
            x.data
        end
        @test qbmap(fn!, objs) == [collect(1:N/2) collect(N/2+1:N)]
        @inferred qbmap(fn!, objs)
    end

    @testset "logtmap" begin
        N = 2 * Threads.nthreads()
        objs = [collect(TestObj(x) for x in 1:N/2) collect(TestObj(x) for x in N/2+1:N)]
        primary = Threads.nthreads() == 1
        fn! = (x) -> begin
            Threads.threadid() == 1 && (primary = true)
            x.data
        end
        p, r = logtmap(fn!, objs)
        @test r == [collect(1:N/2) collect(N/2+1:N)]
        @test length(p.recs) == N*2
        @test primary
        @inferred logtmap(fn!, objs)
    end

    @testset "logbmap" begin
        N = 2 * Threads.nthreads()
        objs = [collect(TestObj(x) for x in 1:N/2) collect(TestObj(x) for x in N/2+1:N)]
        fn! = (x) -> begin
            Threads.nthreads() == 1 || Threads.threadid() == 1 && error("Task on primary")
            x.data
        end
        p, r = logbmap(fn!, objs)
        @test r == [collect(1:N/2) collect(N/2+1:N)]
        @test length(p.recs) == N*2
        @inferred logbmap(fn!, objs)
    end

    @testset "logqmap" begin
        N = 2 * Threads.nthreads()
        objs = [collect(TestObj(x) for x in 1:N/2) collect(TestObj(x) for x in N/2+1:N)]
        primary = Threads.nthreads() == 1
        fn! = (x) -> begin
            Threads.threadid() == 1 && (primary = true)
            x.data
        end
        p, r = logqmap(fn!, objs)
        @test r == [collect(1:N/2) collect(N/2+1:N)]
        @test length(p.recs) == N*2
        @test primary
        @inferred logqmap(fn!, objs)
    end

    @testset "logqbmap" begin
        N = 2 * Threads.nthreads()
        objs = [collect(TestObj(x) for x in 1:N/2) collect(TestObj(x) for x in N/2+1:N)]
        fn! = (x) -> begin
            Threads.nthreads() == 1 || Threads.threadid() == 1 && error("Task on primary")
            x.data
        end
        p, r = logqbmap(fn!, objs)
        @test r == [collect(1:N/2) collect(N/2+1:N)]
        @test length(p.recs) == N*2
        @inferred logqbmap(fn!, objs)
    end

    @testset "bthreads" begin
        function test()
            N = 2 * Threads.nthreads()
            objs = [collect(TestObj(x) for x in 1:N/2) collect(TestObj(x) for x in N/2+1:N)]
            @bthreads for x in objs
                Threads.nthreads() == 1 || Threads.threadid() == 1 && error("Task on primary")
                x.data += 1
            end
            @test [x.data for x in objs] == [collect(2:N/2+1) collect(N/2+2:N+1)]
        end
        test()
        test()
    end

    @testset "qthreads" begin
        function test()
            N = 2 * Threads.nthreads()
            objs = [collect(TestObj(x) for x in 1:N/2) collect(TestObj(x) for x in N/2+1:N)]
            primary = Threads.nthreads() == 1
            @qthreads for x in objs
                Threads.threadid() == 1 && (primary = true)
                x.data += 1
                sleep(0.001)
            end
            @test [x.data for x in objs] == [collect(2:N/2+1) collect(N/2+2:N+1)]
            @test primary
        end
        test()
        test()
    end

end

end # module


text:
using Clustering
using Test
using Random
using LinearAlgebra
using SparseArrays
using Statistics

tests = ["seeding",
         "kmeans",
         "kmedoids",
         "affprop",
         "dbscan",
         "fuzzycmeans",
         "counts",
         "silhouette",
         "varinfo",
         "randindex",
         "hclust",
         "mcl",
         "vmeasure",
         "mutualinfo"]

println("Runing tests:")
for t in tests
    fp = "$(t).jl"
    println("* $fp ...")
    include(fp)
end


text:
module EcoBase

import Base: show, view
import RecipesBase

include("DataTypes.jl")
include("Interface.jl")
include("PlotRecipes.jl")

export nthings, nplaces, occupancy, richness, records, placenames, thingnames
export occurring, noccurring, occupied, noccupied, occurrences
export placeoccurrences, thingoccurrences, cooccurring, places, things
export indices, coordinates, xcells, ycells, cells, xmin, xmax, ymin, ymax
export xrange, yrange, xcellsize, ycellsize, cellsize, getcoords

end # module


text:
module CoinMPInterface

import Cbc

const CONTINUOUS = 1
const INTEGER = 2 # including binary

export CONTINUOUS, INTEGER,
    optionList,
    setOption,
    getOptionValue,
    getOptionDefault,
    getOptionMin,
    getOptionMax,
    CoinProblem,
    GetSolverName,
    GetSolverVersionStr,
    GetSolverVersion,
    GetInfinity,
    LoadMatrix,
    LoadInteger,
    LoadPriority,
    LoadSemiCont,
    CheckProblem,
    GetColCount,
    GetRowCount,
    RegisterMsgLogCallback,
    RegisterLPIterCallback,
    RegisterMipNodeCallback,
    OptimizeProblem,
    GetSolutionStatus,
    GetSolutionText,
    GetObjectValue,
    GetMipBestBound,
    GetIterCount,
    GetMipNodeCount,
    GetSolutionValues,
    ReadFile,
    WriteFile,
    OpenLogFile,
    CloseLogFile

# helper macros/functions

macro coin_ccall(func, args...)
    f = "Coin$(func)"
    quote
        @unix_only ret = ccall(($f,Cbc.libcoinmp), $(args...))
        @windows_only ret = ccall(($f,Cbc.libcoinmp), stdcall, $(args...))
        ret
    end
end

macro coin_checkedccall(func, args...)
    f = "Coin$(func)"
    quote
        @unix_only ret = ccall(($f,Cbc.libcoinmp), Int32, $(args...))
        @windows_only ret = ccall(($f,Cbc.libcoinmp), stdcall, Int32, $(args...))
        if ret != 0
            error("Internal error in ", $f)
        end
    end
end

type CoinProblem
    p::Ptr{Void}
    function CoinProblem()
        p = @coin_ccall CreateProblem Ptr{Void} (Ptr{Uint8},) bytestring("")
        prob = new(p)
        finalizer(prob, delete_problem)
        return prob
    end
end

function delete_problem(prob::CoinProblem)
    if prob.p == C_NULL
        return
    end
    @coin_checkedccall UnloadProblem (Ptr{Void},) prob.p
    prob.p = C_NULL
    return
end

function check_problem(prob::CoinProblem)
   if prob.p == C_NULL
       error("Invalid CoinProblem")
   end
   return true
end

# copied from Clp.jl
typealias VecOrNothing Union(Vector,Nothing)
function vec_or_null{T}(::Type{T}, a::VecOrNothing, len::Integer)
    if isequal(a, nothing) || isa(a, Array{None}) || length(a) == 0
        return C_NULL
    else # todo: helpful message if convert fails
        if length(a) != len
            error("Expected vector to have length $len")
        end
        return convert(Vector{T},a)
    end
end

## Convenience methods for options
# There's no documentation for the meanings of various options and their values
# The only way is to dig into the CoinMP source or find examples...

# returns the list of options
function optionList()
    dummy = CoinProblem()
    l = ASCIIString[]
    for i in 1:200
        s = GetOptionName(dummy,i)
        if s != ""
            push!(l,s)
        end
    end
    return l
end

function setOption(prob::CoinProblem, name::ASCIIString, value)
    id = LocateOptionName(prob,name)
    t = GetOptionType(prob,id)
    if t == 2 || t == 1
        SetIntOption(prob, id, value)
    elseif t == 4
        SetRealOption(prob, id, value)
    else
        error("Unkown option type, this is a bug")
    end
end

for (n1, n2) in ( ("Value",""), ("Default","DefaultValue"),
    ("Min","MinValue"), ("Max","MaxValue") )
    fname = symbol(string("getOption",n1))
    coinint = symbol(string("GetIntOption",n2))
    coinreal = symbol(string("GetRealOption",n2))
    @eval begin function $(fname)(prob::CoinProblem, name::ASCIIString)
        id = LocateOptionName(prob,name)
        t = GetOptionType(prob,id)
        if t == 2 || t == 1
            $coinint(prob, id)
        elseif t == 4
            $coinreal(prob, id)
        else
            error("Unkown option type, this is a bug")
        end
    end end
end

## actual interface

# CoinInitSolver and CoinFreeSolver don't do anything

function GetSolverName()
    s = @coin_ccall GetSolverName Ptr{Uint8} ()
    bytestring(s)
end

function GetSolverVersionStr()
    s = @coin_ccall GetSolverVersionStr Ptr{Uint8} ()
    bytestring(s)
end

function GetSolverVersion()
    @coin_ccall GetSolverVersion Float64 ()
end

# not sure the point of CoinGetFeatures and CoinGetMethods

function GetInfinity()
    @coin_ccall GetInfinity Float64 ()
end

function LoadMatrix(prob::CoinProblem, objective_sense::Integer,
    objective_offset::Float64, objective_coeffs::VecOrNothing,
    col_lb::VecOrNothing, col_ub::VecOrNothing,
    row_lb::VecOrNothing, row_ub::VecOrNothing,
    constraint_matrix::AbstractMatrix)
    check_problem(prob)

    mat = convert(SparseMatrixCSC{Float64,Int32},constraint_matrix)
    nrow,ncol = mat.m, mat.n

    @coin_checkedccall LoadMatrix (Ptr{Void}, Int32, Int32, Int32, Int32, Int32,
        Float64, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Uint8}, Ptr{Float64},
        Ptr{Float64}, Ptr{Int32}, Ptr{Int32}, Ptr{Int32},
        Ptr{Float64}) prob.p ncol nrow nnz(mat) 0 objective_sense objective_offset vec_or_null(Float64, objective_coeffs, ncol) vec_or_null(Float64, col_lb, ncol) vec_or_null(Float64, col_ub, ncol) C_NULL vec_or_null(Float64, row_lb, nrow) vec_or_null(Float64, row_ub, nrow) mat.colptr.-int32(1) C_NULL vec_or_null(Int32, mat.rowval.-int32(1),nnz(mat)) vec_or_null(Float64, mat.nzval, nnz(mat))
end

# TODO: CoinLoadNames

# CoinLoadInitValues doesn't do anything

# 'B' - binary
# 'I' - integer
# 'C' - continuous
function LoadInteger(prob::CoinProblem, column_type::Vector{Uint8})
    check_problem(prob)
    if length(column_type) != GetColCount(prob)
        error("Length of vector must equal number of variables in the problem")
    end
    @coin_checkedccall LoadInteger (Ptr{Void}, Ptr{Uint8}) prob.p column_type
end

function LoadPriority(prob::CoinProblem,
    PriorIndex::Vector{Int32}, PriorValues::Vector{Int32})
    # PriorBranch argument is ignored?
    check_problem(prob)
    len = length(PriorIndex)
    PriorIndex = PriorIndex .- int32(1) # input is 1-based
    if len != length(PriorValues)
        error("All input vectors must be same length")
    end
    if len > GetColCount(prob)
        error("More priorities than variables")
    end
    @coin_checkedccall LoadPriority (Ptr{Void}, Ptr{Int32}, Ptr{Int32}, Ptr{Int32}) prob.p PriorIndex PriorValues C_NULL
end

# TODO: CoinLoadSos

function LoadSemiCont(prob::CoinProblem, SemiIndex::Vector{Int32})
    check_problem(prob)
    SemiIndex = SemiIndex - int32(1)
    @coin_checkedccall LoadSemiCont (Ptr{Void}, Int32, Ptr{Int32}) prob.p length(SemiIndex) SemiIndex
end

# CoinLoadQuadratic and CoinLoadNonlinear are not implemented by CoinMP

# This returns 0 if the problem data is ok, otherwise it returns
# one of the SOLV_CHECK_* error codes defined in CoinMP.h
function CheckProblem(prob::CoinProblem)
    check_problem(prob)
    @coin_ccall CheckProblem Int32 (Ptr{Void},) prob.p
end

# TODO: CoinGetProblemName -- is this useful?

function GetColCount(prob::CoinProblem)
    check_problem(prob)
    @coin_ccall GetColCount Int32 (Ptr{Void},) prob.p
end

function GetRowCount(prob::CoinProblem)
    check_problem(prob)
    @coin_ccall GetRowCount Int32 (Ptr{Void},) prob.p
end

# TODO: CoinGetColName, CoinGetRowName

# Pass in c-pointer to callback function. (from cfunction())
# See CoinMP.h for signature
# TODO: test these
function RegisterMsgLogCallback(prob::CoinProblem, func::Ptr{Void})
    check_problem(prob)
    @coin_checkedccall RegisterMsgLogCallback (Ptr{Void},Ptr{Void},Ptr{Void}) prob.p func C_NULL
end

function RegisterLPIterCallback(prob::CoinProblem, func::Ptr{Void})
    check_problem(prob)
    @coin_checkedccall RegisterLPIterCallback (Ptr{Void},Ptr{Void},Ptr{Void}) prob.p func C_NULL
end

function RegisterMipNodeCallback(prob::CoinProblem, func::Ptr{Void})
    check_problem(prob)
    @coin_checkedccall RegisterMipNodeCallback (Ptr{Void},Ptr{Void},Ptr{Void}) prob.p func C_NULL
end

function OptimizeProblem(prob::CoinProblem)
    check_problem(prob)
    # Method argument is ignored?
    @coin_checkedccall OptimizeProblem (Ptr{Void},Int32) prob.p 0
end

# no documentation for what these values mean...
function GetSolutionStatus(prob::CoinProblem)
    check_problem(prob)
    @coin_ccall GetSolutionStatus Int32 (Ptr{Void},) prob.p
end

function GetSolutionText(prob::CoinProblem)
    check_problem(prob)
    s = @coin_ccall GetSolutionText Ptr{Uint8} (Ptr{Void},) prob.p
    return bytestring(s)
end

function GetObjectValue(prob::CoinProblem)
    check_problem(prob)
    @coin_ccall GetObjectValue Float64 (Ptr{Void},) prob.p
end

function GetMipBestBound(prob::CoinProblem)
    check_problem(prob)
    @coin_ccall GetMipBestBound Float64 (Ptr{Void},) prob.p
end

function GetIterCount(prob::CoinProblem)
    check_problem(prob)
    @coin_ccall GetIterCount Int32 (Ptr{Void},) prob.p
end

function GetMipNodeCount(prob::CoinProblem)
    check_problem(prob)
    @coin_ccall GetMipNpdeCount Int32 (Ptr{Void},) prob.p
end

# return vector of primal solutions
# we could return more in the future
function GetSolutionValues(prob::CoinProblem)
    check_problem(prob)
    ncol = GetColCount(prob)
    x = Array(Float64,ncol)
    @coin_checkedccall GetSolutionValues (Ptr{Void},
        Ptr{Float64},Ptr{Float64},Ptr{Float64},Ptr{Float64}) prob.p x C_NULL C_NULL C_NULL
    return x
end

# TODO: GetSolutionRanges -- what does this mean? sensitivity analysis?

# TODO: GetSolutionBasis -- pretty useless since you can't input a starting basis

# the only file type currently supported is MPS (3)
function ReadFile(prob::CoinProblem, FileType::Integer, FileName::ASCIIString)
    check_problem(prob)
    @coin_checkedccall ReadFile (Ptr{Void},Int32,Ptr{Uint8}) prob.p FileType FileName
end

function WriteFile(prob::CoinProblem, FileType::Integer, FileName::ASCIIString)
    check_problem(prob)
    @coin_checkedccall WriteFile (Ptr{Void},Int32,Ptr{Uint8}) prob.p FileType FileName
end

function OpenLogFile(prob::CoinProblem, FileName::ASCIIString)
    check_problem(prob)
    @coin_checkedccall OpenLogFile (Ptr{Void},Ptr{Uint8}) prob.p FileName
end

function CloseLogFile(prob::CoinProblem)
    check_problem(prob)
    @coin_checkedccall CloseLogFile (Ptr{Void},) prob.p
end

function GetOptionCount(prob::CoinProblem)
    check_problem(prob)
    @coin_ccall GetOptionCount Int32 (Ptr{Void},) prob.p
end

function LocateOptionID(prob::CoinProblem, OptionID::Integer)
    check_problem(prob)
    @coin_ccall LocateOptionID Int32 (Ptr{Void},Int32) prob.p OptionID
end

function LocateOptionName(prob::CoinProblem, OptionName::ASCIIString)
    check_problem(prob)
    @coin_ccall LocateOptionName Int32 (Ptr{Void},Ptr{Uint8}) prob.p OptionName
end

# Don't know what this does
function GetOptionID(prob::CoinProblem, OptionNr::Integer)
    check_problem(prob)
    @coin_ccall GetOptionID Int32 (Ptr{Void},Int32) prob.p OptionNr
end

# TODO: CoinGetOptionInfo, CoinGetIntOptionMinMax, CoinGetRealOptionMinMax

# TODO: CoinGetOptionGroup

# 1 -- 0/1, 2 -- INT, 4 -- REAL,
function GetOptionType(prob::CoinProblem, OptionID::Integer)
    check_problem(prob)
    @coin_ccall GetOptionType Int32 (Ptr{Void},Int32) prob.p OptionID
end

for s in (:GetIntOptionDefaultValue, :GetIntOptionMinValue, :GetIntOptionMaxValue)
    @eval begin
        function $(s)(prob::CoinProblem, OptionID::Integer)
            check_problem(prob)
            @coin_ccall $s Int32 (Ptr{Void},Int32) prob.p OptionID
        end
    end
end

for s in (:GetRealOptionDefaultValue, :GetRealOptionMinValue, :GetRealOptionMaxValue)
    @eval begin
        function $(s)(prob::CoinProblem, OptionID::Real)
            check_problem(prob)
            @coin_ccall $s Float64 (Ptr{Void},Int32) prob.p OptionID
        end
    end
end

for s in (:GetOptionName, :GetOptionShortName, :GetStringOption)
    @eval begin
        function $(s)(prob::CoinProblem, OptionID::Integer)
            check_problem(prob)
            bytestring(@coin_ccall GetOptionName Ptr{Uint8} (Ptr{Void},Int32) prob.p OptionID)
        end
    end
end

# TODO: CoinGetOptionChanged

function GetIntOption(prob::CoinProblem, OptionID::Integer)
    check_problem(prob)
    @coin_ccall GetIntOption Int32 (Ptr{Void},Int32) prob.p OptionID
end

function SetIntOption(prob::CoinProblem, OptionID::Integer, Value::Integer)
    check_problem(prob)
    @coin_checkedccall SetIntOption (Ptr{Void},Int32,Int32) prob.p OptionID Value
end

function GetRealOption(prob::CoinProblem, OptionID::Integer)
    check_problem(prob)
    @coin_ccall GetRealOption Float64 (Ptr{Void},Int32) prob.p OptionID
end

function SetRealOption(prob::CoinProblem, OptionID::Integer, Value::Real)
    check_problem(prob)
    @coin_checkedccall SetRealOption (Ptr{Void},Int32,Float64) prob.p OptionID Value
end

function SetStringOption(prob::CoinProblem, OptionID::Integer, Value::ASCIIString)
    check_problem(prob)
    @coin_checkedccall SetStringOption (Ptr{Void},Int32,Ptr{Uint8}) prob.p OptionID Value
end

end


text:
# Different shader string literals- usage: e.g. frag" my shader code"
macro frag_str(source::AbstractString)
    quote
        ($source, GL_FRAGMENT_SHADER)
    end
end
macro vert_str(source::AbstractString)
    quote
        ($source, GL_VERTEX_SHADER)
    end
end
macro geom_str(source::AbstractString)
    quote
        ($source, GL_GEOMETRY_SHADER)
    end
end
macro comp_str(source::AbstractString)
    quote
        ($source, GL_COMPUTE_SHADER)
    end
end

function getinfolog(obj::GLuint)
    # Return the info log for obj, whether it be a shader or a program.
    isShader    = glIsShader(obj)
    getiv       = isShader == GL_TRUE ? glGetShaderiv : glGetProgramiv
    get_log     = isShader == GL_TRUE ? glGetShaderInfoLog : glGetProgramInfoLog

    # Get the maximum possible length for the descriptive error message
    maxlength = GLint[0]
    getiv(obj, GL_INFO_LOG_LENGTH, maxlength)
    maxlength = first(maxlength)
    # Return the text of the message if there is any
    if maxlength > 0
        buffer = zeros(GLchar, maxlength)
        sizei = GLsizei[0]
        get_log(obj, maxlength, sizei, buffer)
        length = first(sizei)
        return unsafe_string(pointer(buffer), length)
    else
        return "success"
    end
end

function iscompiled(shader::GLuint)
    success = GLint[0]
    glGetShaderiv(shader, GL_COMPILE_STATUS, success)
    return first(success) == GL_TRUE
end
islinked(program::GLuint) = glGetProgramiv(program, GL_LINK_STATUS) == GL_TRUE

function createshader(shadertype::GLenum)
    shaderid = glCreateShader(shadertype)
    @assert shaderid > 0 "opengl context is not active or shader type not accepted. Shadertype: $(GLENUM(shadertype).name)"
    shaderid::GLuint
end
function createprogram()
    program = glCreateProgram()
    @assert program > 0 "couldn't create program. Most likely, opengl context is not active"
    program::GLuint
end

shadertype(s::Shader) = s.typ
function shadertype(ext::AbstractString)
    ext == ".comp" && return GL_COMPUTE_SHADER
    ext == ".vert" && return GL_VERTEX_SHADER
    ext == ".frag" && return GL_FRAGMENT_SHADER
    ext == ".geom" && return GL_GEOMETRY_SHADER
    error("$ext not a valid shader extension")
end

function uniformlocations(nametypedict::Dict{Symbol, GLenum}, program)
    result = Dict{Symbol, Tuple}()
    texturetarget = -1 # start -1, as texture samplers start at 0
    for (name, typ) in nametypedict
        loc = get_uniform_location(program, name)
        str_name = string(name)
        if istexturesampler(typ)
            texturetarget += 1
            result[name] = (loc, texturetarget)
        else
            result[name] = (loc,)
        end
    end
    return result
end

abstract type AbstractLazyShader end
struct LazyShader <: AbstractLazyShader
    paths::Tuple
    kw_args::Dict{Symbol, Any}
    function LazyShader(paths...; kw_args...)
        args = Dict{Symbol, Any}(kw_args)
        get!(args, :view, Dict{String, String}())
        new(paths, args)
    end
end

gl_convert(shader::GLProgram, data) = shader

# caching templated shaders is a pain -.-

# cache for template keys per file
# path --> template keys
const _template_cache = Dict{String, Vector{String}}()
# path --> Dict{template_replacements --> Shader)
const _shader_cache = Dict{String, Dict{Any, Shader}}()
const _program_cache = Dict{Any, GLProgram}()

function empty_shader_cache!()
    empty!(_template_cache)
    empty!(_shader_cache)
    empty!(_program_cache)
end

# TODO remove this silly constructor
function compile_shader(source::Vector{UInt8}, typ, name)
    shaderid = createshader(typ)
    glShaderSource(shaderid, source)
    glCompileShader(shaderid)
    if !GLAbstraction.iscompiled(shaderid)
        GLAbstraction.print_with_lines(String(source))
        @warn("shader $(name) didn't compile. \n$(GLAbstraction.getinfolog(shaderid))")
    end
    Shader(name, source, typ, shaderid)
end

function compile_shader(path, source_str::AbstractString)
    typ = shadertype(splitext(path)[2])
    source = Vector{UInt8}(source_str)
    name = Symbol(path)
    compile_shader(source, typ, name)
end

function get_shader!(path, template_replacement, view, attributes)
    # this should always be in here, since we already have the template keys
    shader_dict = _shader_cache[path]
    get!(shader_dict, template_replacement) do
        template_source = read(path, String)
        source = mustache_replace(template_replacement, template_source)
        compile_shader(path, source)
    end::Shader
end

function get_template!(path, view, attributes)
    get!(_template_cache, path) do
        _, ext = splitext(path)

        typ = shadertype(ext)
        template_source = read(path, String)
        source, replacements = template2source(
            template_source, view, attributes
        )
        s = compile_shader(path, source)
        template_keys = collect(keys(replacements))
        template_replacements = collect(values(replacements))
        # can't yet be in here, since we didn't even have template keys
        _shader_cache[path] = Dict(template_replacements => s)

        template_keys
    end
end

function compile_program(shaders, fragdatalocation)
    # Remove old shaders
    program = createprogram()
    #attach new ones
    foreach(shaders) do shader
        glAttachShader(program, shader.id)
    end

    #Bind frag data
    for (location, name) in fragdatalocation
        glBindFragDataLocation(program, location, ascii(name))
    end

    #link program
    glLinkProgram(program)
    if !GLAbstraction.islinked(program)
        error(
            "program $program not linked. Error in: \n",
            join(map(x-> string(x.name), shaders), " or "), "\n", getinfolog(program)
        )
    end
    # Can be deleted, as they will still be linked to Program and released after program gets released
    #foreach(glDeleteShader, shader_ids)
    # generate the link locations
    nametypedict = uniform_name_type(program)
    uniformlocationdict = uniformlocations(nametypedict, program)
    GLProgram(program, shaders, nametypedict, uniformlocationdict)
end

function get_view(kw_dict)
    _view = kw_dict[:view]
    extension = Sys.isapple() ? "" : "#extension GL_ARB_draw_instanced : enable\n"
    _view["GLSL_EXTENSION"] = extension*get(_view, "GLSL_EXTENSIONS", "")
    _view["GLSL_VERSION"] = glsl_version_string()
    _view
end

function gl_convert(lazyshader::AbstractLazyShader, data)
    kw_dict = lazyshader.kw_args
    paths = lazyshader.paths
    if all(x-> isa(x, Shader), paths)
        fragdatalocation = get(kw_dict, :fragdatalocation, Tuple{Int, String}[])
        return compile_program([paths...], fragdatalocation)
    end
    v = get_view(kw_dict)
    fragdatalocation = get(kw_dict, :fragdatalocation, Tuple{Int, String}[])

    # Tuple(Source, ShaderType)
    if all(paths) do x
            isa(x, Tuple) && length(x) == 2 &&
            isa(first(x), String) &&
            isa(last(x), GLenum)
        end
        # we don't cache view & templates for shader strings!
        shaders = map(paths) do source_typ
            source, typ = source_typ
            src, _ = template2source(source, v, data)
            compile_shader(Vector{UInt8}(src), typ, :from_string)
        end
        return compile_program([shaders...], fragdatalocation)
    end
    if !all(x-> isa(x, String), paths)
        error("Please supply only paths or tuples of (source, typ) for Lazy Shader
            Found: $paths"
        )
    end
    template_keys = Vector{Vector{String}}(undef, length(paths))
    replacements = Vector{Vector{String}}(undef, length(paths))
    for (i, path) in enumerate(paths)
        template = get_template!(path, v, data)
        template_keys[i] = template
        replacements[i] = String[mustache2replacement(t, v, data) for t in template]
    end
    program = get!(_program_cache, (paths, replacements)) do
        # when we're here, this means there were uncached shaders, meaning we definitely have
        # to compile a new program
        shaders = Vector{Shader}(undef, length(paths))
        for (i, path) in enumerate(paths)
            tr = Dict(zip(template_keys[i], replacements[i]))
            shaders[i] = get_shader!(path, tr, v, data)
        end
        compile_program(shaders, fragdatalocation)
    end
end

function insert_from_view(io, replace_view::Function, keyword::AbstractString)
    print(io, replace_view(keyword))
    nothing
end

function insert_from_view(io, replace_view::Dict, keyword::AbstractString)
    if haskey(replace_view, keyword)
        print(io, replace_view[keyword])
    end
    nothing
end
"""
Replaces
{{keyword}} with the key in `replace_view`, or replace_view(key)
in a string
"""
function mustache_replace(replace_view::Union{Dict, Function}, string)
    io = IOBuffer()
    replace_started = false
    open_mustaches = 0
    closed_mustaches = 0
    i = 0
    replace_begin = i
    last_char = SubString(string, 1, 1)
    len = lastindex(string)
    while i <= len
        i = nextind(string, i)
        i > len && break
        char = string[i]
        if replace_started
            # ignore, or wait for }
            if char == '}'
                closed_mustaches += 1
                if closed_mustaches == 2 # we found a complete mustache!
                    insert_from_view(io, replace_view, SubString(string, replace_begin+1, i-2))
                    open_mustaches = 0
                    closed_mustaches = 0
                    replace_started = false
                end
            else
                closed_mustaches = 0
                continue
            end
        elseif char == '{'
            open_mustaches += 1
            if open_mustaches == 2
                replace_begin = i
                replace_started = true
            end
        else
            if open_mustaches == 1
                print(io, last_char)
            end
            print(io, char) # just copy all the rest
            open_mustaches = 0
            closed_mustaches = 0
        end
        last_char = char
    end
    String(take!(io))
end

function mustache2replacement(mustache_key, view, attributes)
    haskey(view, mustache_key) && return view[mustache_key]
    for postfix in ("_type", "_calculation")
        keystring = replace(mustache_key, postfix => "")
        keysym = Symbol(keystring)
        if haskey(attributes, keysym)
            val = attributes[keysym]
            if !isa(val, AbstractString)
                if postfix == "_type"
                    return toglsltype_string(val)::String
                else  postfix == "_calculation"
                    return glsl_variable_access(keystring, val)
                end
            end
        end
    end
    return ""
    # error("No match found: $(mustache_key)")
end

# Takes a shader template and renders the template and returns shader source
template2source(source::Vector{UInt8}, view, attributes::Dict{Symbol, Any}) = template2source(String(source), attributes, view)
function template2source(source::AbstractString, view, attributes::Dict{Symbol, Any})
    replacements = Dict{String, String}()
    source = mustache_replace(source) do mustache_key
        r = mustache2replacement(mustache_key, view, attributes)
        replacements[mustache_key] = r
        return r
    end
    return source, replacements
end


text:
type SwapForwardBasisSystem <: MarketModelBasisSystem
  rateTimes::Vector{Float64}
  exerciseTimes::Vector{Float64}
  currentIndex::Int
  rateIndex::Vector{Int}
  evolution::EvolutionDescription
end

function SwapForwardBasisSystem(rateTimes::Vector{Float64}, exerciseTimes::Vector{Float64})
  evolution = EvolutionDescription(rateTimes, exerciseTimes)
  rateIndex = Vector{Int}(length(exerciseTimes))
  j = 1
  for i in eachindex(exerciseTimes)
    while j <= length(rateTimes) && rateTimes[j] < exerciseTimes[i]
      j += 1
    end

    rateIndex[i] = j
  end

  return SwapForwardBasisSystem(rateTimes, exerciseTimes, -1, rateIndex, evolution)
end

reset!(bs::SwapForwardBasisSystem) = bs.currentIndex = 1

next_step!(bs::SwapForwardBasisSystem, ::CurveState) = bs.currentIndex += 1

function number_of_functions(bs::SwapForwardBasisSystem)
  n = length(bs.exerciseTimes)
  sizes = fill(10, n)
  if bs.rateIndex[n] == length(bs.rateIndex) - 3
    sizes[end] = 6
  end
  if bs.rateIndex[n] == length(bs.rateIndex) - 2
    sizes[end] = 3
  end

  return sizes
end

function set_values!(bs::SwapForwardBasisSystem, currentState::CurveState, results::Vector{Float64})
  rateIndex = bs.rateIndex[bs.currentIndex - 1]

  if rateIndex < length(bs.rateTimes) - 2
    resize!(results, 10)

    x = forward_rate(currentState, rateIndex)
    y = coterminal_swap_rate(currentState, rateIndex + 1)
    z = discount_ratio(currentState, rateIndex, length(bs.rateTimes))

    results[1] = 1.0
    results[2] = x
    results[3] = y
    results[4] = z
    results[5] = x*y
    results[6] = y*z
    results[7] = z*x
    results[8] = x*x
    results[9] = y*y
    results[10] = z*z
  else
    if rateIndex == length(bs.rateTimes) - 2
      x = forward_rate(currentState, rateIndex)
      y = forward_rate(currentState, rateIndex + 1)

      resize!(results, 6)

      results[1] = 1.0
      results[2] = x
      results[3] = y
      results[4] = x*x
      results[5] = x*y
      results[6] = y*y
    else
      x = forward_rate(currentState, rateIndex)

      resize!(results, 3)
      results[1] = 1.0
      results[2] = x
      results[3] = x*x
    end
  end

  return results
end

clone(bs::SwapForwardBasisSystem) = SwapForwardBasisSystem(copy(bs.rateTimes), copy(bs.exerciseTimes), bs.currentIndex, copy(bs.rateIndex), clone(bs.evolution))


text:
"""
    struct ChebyshevNodes{T} <: AbstractIntervalGrid{T}

A grid with chebyshev nodes on [-1,1].

# Example
```jldocs
julia> ChebyshevExtremae(4)
4-element ChebyshevExtremae{Float64}:
  1.0
  0.5000000000000001
 -0.4999999999999998
 -1.0
```
"""
struct ChebyshevTNodes{T} <: AbstractIntervalGrid{T}
    n   ::  Int
end

const ChebyshevNodes = ChebyshevTNodes
const ChebyshevPoints = ChebyshevTNodes
unsafe_getindex(g::ChebyshevTNodes{T}, i::Int) where {T} = T(-1)*cos((i-T(1)/2) * T(pi) / (g.n) )

"""
    struct ChebyshevExtremae{T} <: AbstractIntervalGrid{T}

A grid with chebyshev extrema on [-1,1].

# Example
```jldocs
julia> ChebyshevExtremae(4)
4-element ChebyshevExtremae{Float64}:
  1.0
  0.5000000000000001
 -0.4999999999999998
 -1.0
```
"""
struct ChebyshevExtremae{T} <: AbstractIntervalGrid{T}
    n   ::  Int
end

const ChebyshevPointsOfTheSecondKind = ChebyshevExtremae

# TODO: flip the values so that they are sorted
unsafe_getindex(g::ChebyshevExtremae{T}, i::Int) where {T} = i == 0 ? T(0) : cos((i-1)*T(pi) / (g.n-1) )

struct ChebyshevUNodes{T} <: AbstractIntervalGrid{T}
    n   :: Int
end

unsafe_getindex(nodes::ChebyshevUNodes{T}, i::Int) where T = cos((nodes.n + 1 - i) * convert(T,π) / (nodes.n + 1))

struct LegendreNodes{T} <: AbstractIntervalGrid{T}
    n   :: Int
    nodes   :: Vector{T}
    LegendreNodes(vector::Vector{T}) where {T} = new{T}(length(vector), vector)
end

struct LaguerreNodes{T} <: AbstractIntervalGrid{T}
    n   :: Int
    α   :: T
    nodes   :: Vector{T}
    LaguerreNodes(α::T, vector::Vector{T}) where {T} = new{T}(length(vector), α, vector)
end

struct HermiteNodes{T} <: AbstractIntervalGrid{T}
    n   :: Int
    nodes   :: Vector{T}
    HermiteNodes(vector::Vector{T}) where {T} = new{T}(length(vector), vector)
end

struct JacobiNodes{T} <: AbstractIntervalGrid{T}
    n   ::  Int
    α   ::  T
    β   ::  T
    nodes   ::Vector
    JacobiNodes(α::T, β::T, vector::Vector{T}) where {T} = new{T}(length(vector), α, β, vector)
end

for GRID in (:LegendreNodes,:HermiteNodes,:LaguerreNodes,:JacobiNodes)
    @eval unsafe_getindex(grid::$GRID, i::Int) = unsafe_getindex(grid.nodes, i)
end

for (Tp,fun) in zip((:ChebyshevTWeights,), (:chebyshevtweights_fun, ))
    @eval begin
        struct $Tp{T,N,Axes} <:FillArrays.AbstractFill{T,N,Axes}
            axes::Axes
            @inline $Tp{T, N}(sz::Axes) where Axes<:Tuple{Vararg{AbstractUnitRange,N}} where {T, N} =
                new{T,N,Axes}(sz)
            @inline $Tp{T,0,Tuple{}}(sz::Tuple{}) where T = new{T,0,Tuple{}}(sz)
        end

        @inline $Tp{T, 0}(sz::Tuple{}) where {T} = $Tp{T,0,Tuple{}}(sz)
        @inline $Tp{T, N}(sz::Tuple{Vararg{<:Integer, N}}) where {T, N} = $Tp{T,N}(Base.OneTo.(sz))
        @inline $Tp{T, N}(sz::Vararg{<:Integer, N}) where {T, N} = $Tp{T,N}(sz)
        @inline $Tp{T}(sz::Vararg{Integer,N}) where {T, N} = $Tp{T, N}(sz)
        @inline $Tp{T}(sz::SZ) where SZ<:Tuple{Vararg{Any,N}} where {T, N} = $Tp{T, N}(sz)
        @inline $Tp(sz::Vararg{Any,N}) where N = $Tp{Float64,N}(sz)
        @inline $Tp(sz::SZ) where SZ<:Tuple{Vararg{Any,N}} where N = $Tp{Float64,N}(sz)

        @inline $Tp{T,N}(A::AbstractArray{V,N}) where{T,V,N} = $Tp{T,N}(size(A))
        @inline $Tp{T}(A::AbstractArray) where{T} = $Tp{T}(size(A))
        @inline $Tp(A::AbstractArray) = $Tp(size(A))

        @inline axes(Z::$Tp) = Z.axes
        @inline size(Z::$Tp) = length.(Z.axes)
        @inline FillArrays.getindex_value(Z::$Tp) = $fun(Z)

        AbstractArray{T}(F::$Tp{T}) where T = F
        AbstractArray{T,N}(F::$Tp{T,N}) where {T,N} = F
        AbstractArray{T}(F::$Tp) where T = $Tp{T}(F.axes)
        AbstractArray{T,N}(F::$Tp{V,N}) where {T,V,N} = $Tp{T}(F.axes)
        convert(::Type{AbstractArray{T}}, F::$Tp{T}) where T = AbstractArray{T}(F)
        convert(::Type{AbstractArray{T,N}}, F::$Tp{T,N}) where {T,N} = AbstractArray{T,N}(F)
        convert(::Type{AbstractArray{T}}, F::$Tp) where T = AbstractArray{T}(F)
        convert(::Type{AbstractArray{T,N}}, F::$Tp) where {T,N} = AbstractArray{T,N}(F)

        getindex(F::$Tp{T,0}) where T = getindex_value(F)
        function getindex(F::$Tp{T}, kj::Vararg{AbstractVector{II},N}) where {T,II<:Integer,N}
            checkbounds(F, kj...)
            Fill{T}(FillArrays.getindex_value(F),length.(kj))
        end

        function getindex(A::$Tp{T}, kr::AbstractVector{Bool}) where T
            length(A) == length(kr) || throw(DimensionMismatch())
            Fill{T}(FillArrays.getindex_value(F),count(kr))
        end
        function getindex(A::$Tp{T}, kr::AbstractArray{Bool}) where T
            size(A) == size(kr) || throw(DimensionMismatch())
            Fill{T}(FillArrays.getindex_value(F),count(kr))
        end
    end
end

chebyshevtweights_fun(Z::ChebyshevTWeights{T}) where T =  convert(T,π) / convert(T,length(Z))

abstract type VectorAlias{T} <: AbstractVector{T} end
getindex(vector::VectorAlias, i::Int) = getindex(vector.vector, i)
unsafe_getindex(vector::VectorAlias, i::Int) = unsafe_getindex(vector.vector, i)
size(vector::VectorAlias) = size(vector.vector)

abstract type FunctionVector{T} <: AbstractVector{T} end
length(vector::FunctionVector) = vector.n
size(vector::FunctionVector) = (length(vector),)

function getindex(vector::FunctionVector, i::Int)
    @boundscheck (1 <= i <= length(vector)) || throw(BoundsError())
    @inbounds unsafe_getindex(vector, i)
end

struct ChebyshevUWeights{T} <: FunctionVector{T}
    n   :: Int
end

unsafe_getindex(weights::ChebyshevUWeights{T}, i::Int) where {T} =
    convert(T,π)/(weights.n + 1) * sin(convert(T,weights.n + 1 -i) / (weights.n + 1) * convert(T,π))^2

gausschebyshev(::Type{T}, n::Int) where T =
    ChebyshevTNodes{T}(n), ChebyshevTWeights{T}(n)
gausschebyshev(n::Int) = gausschebyshev(Float64, n)

gausschebyshevu(::Type{T}, n::Int) where T =
    ChebyshevUNodes{T}(n), ChebyshevUWeights{T}(n)
gausschebyshevu(n::Int) = gausschebyshevu(Float64,n)

struct LegendreWeights{T} <: VectorAlias{T}
    vector  ::  Vector{T}
end

function gausslegendre(::Type{Float64}, n::Int)
    x,w = FastGaussQuadrature.gausslegendre(n)
    LegendreNodes(x), LegendreWeights(w)
end
function gausslegendre(::Type{T}, n::Int) where T
    x,w = GaussQuadrature.legendre(T, n)
    LegendreNodes(x), LegendreWeights(w)
end
gausslegendre(n::Int) = gausslegendre(Float64, n)

struct LaguerreWeights{T} <: VectorAlias{T}
    α       ::  T
    vector  ::  Vector{T}
end

function gausslaguerre(::Type{Float64}, n::Int, α::Float64)
    x,w = FastGaussQuadrature.gausslaguerre(n, α)
    LaguerreNodes(α, x), LaguerreWeights(α, w)
end
function gausslaguerre(::Type{T}, n::Int, α::T) where T
    x,w = GaussQuadrature.laguerre(n, α)
    LaguerreNodes(α, x), LaguerreWeights(α, w)
end
gausslaguerre(n::Int, α::T) where T = gausslaguerre(T, n, α)

struct HermiteWeights{T} <: VectorAlias{T}
    vector  ::  Vector{T}
end
function gausshermite(::Type{Float64}, n::Int)
    x,w = FastGaussQuadrature.gausshermite(n)
    HermiteNodes(x), HermiteWeights(w)
end
function gausshermite(::Type{T}, n::Int) where T
    x,w = GaussQuadrature.hermite(T, n)
    HermiteNodes(x), HermiteWeights(w)
end
gausshermite(n::Int) = gausshermite(Float64, n)

struct JacobiWeights{T} <: VectorAlias{T}
    α       ::  T
    β       ::  T
    vector  ::  Vector{T}
end
function gaussjacobi(::Type{Float64}, n::Int, α::Float64, β::Float64)
    x,w = FastGaussQuadrature.gaussjacobi(n, α, β)
    JacobiNodes(α, β, x), JacobiWeights(α, β, w)
end
function gaussjacobi(::Type{T}, n::Int, α::T, β::T) where T
    x,w = GaussQuadrature.jacobi(n, α, β)
    JacobiNodes(α, β, x), JacobiWeights(α, β, w)
end
gaussjacobi(n::Int, α::T, β::T) where T = gaussjacobi(T, n, α, β)


text:
# -- PRIVATE FUNCTIONS ------------------------------------------- #
function _iv_objective_function(x, assetSet::Set{PSAbstractAsset}, parameters::PSBinaryLatticeModel, baseAssetPriceValue::Float64, optionPriceValue::Float64)

    volatility::Float64
    timeToExercise::Float64
    riskFreeRate::Float64
    dividendRate::Float64
    numberOfLevels::Int64

    # create a new paramter w/new volatility -
    perturbedVolatility = x[1]
    perturbedLatticeModel = PSBinaryLatticeModel(perturbedVolatility, parameters.timeToExercise, parameters.riskFreeRate,
        parameters.dividendRate; numberOfLevels=parameters.numberOfLevels)

    #= We are calling:
    option_contract_price(contractSet::Set{PSAbstractAsset}, latticeModel::PSBinaryLatticeModel, baseUnderlyingPrice::Float64;
        earlyExercise::Bool = true)::PooksoftBase.PSResult
    =#

    # re-compute the option price -
    result = option_contract_price(assetSet, perturbedLatticeModel, baseAssetPriceValue)
    results_tuple = result.value
    estimatedPriceValue = results_tuple.cost_calculation_result.option_contract_price_array[1];

    # compute the error -
    error_term = (optionPriceValue - estimatedPriceValue)

    # return -
    return error_term*error_term
end
# ---------------------------------------------------------------- #

# -- PUBLIC FUNCTIONS -------------------------------------------- #
function estimate_implied_volatility(contract::PSAbstractAsset, parameters::PSBinaryLatticeModel, baseAssetPriceValue::Float64;
    earlyExercise::Bool = false)

    # setup asset set -
    assetSet = Set{PSAbstractAsset}()
    push!(assetSet, contract)

    # what is the premium?
    optionPriceValue = contract.premimumValue

    # what is the initial volatility -
    initialVolatility = [parameters.volatility]

    # setup the objective function -
    OF(p) = _iv_objective_function(p, assetSet, parameters, baseAssetPriceValue, optionPriceValue)

    # call the optimizer -
    opt_result = Optim.optimize(OF, initialVolatility, BFGS())

    # return the result -
    return Optim.minimizer(opt_result)[1]
end

function compute_weighted_volatility(data::DataFrame, weightKey::Symbol, volatilityKey::Symbol)::PSResult

    # initialize -
    tmp_array = Array{Float64,1}()
    volatility_value = nothing

    # sum the wght key -
    iv_array = data[!,volatilityKey]
    wght_array = data[!,weightKey]

    # main loop -
    for (index, iv_value) in enumerate(iv_array)

        # get weight value -
        wght_value = wght_array[index]

        # compute weighted iv value numerator -
        value = wght_value*iv_array[index]
        push!(tmp_array,value)
    end

    # compute value -
    wght_sum = sum(wght_array)
    volatility_value = sum(tmp_array)*(1/wght_sum)

    # return -
    return PSResult(volatility_value)
end
# ---------------------------------------------------------------- #


text:
using SnakePlot, GLMakie

path = "SnakePlot/example/piezo2.txt"
tms = [13:24, 31:43, 51:76, 123:141, 222:237, 241:258, 265:287, 336:355,
		493:514, 520:531, 536:562, 584:614, 690:703, 710:728, 738:757, 790:811,
		958:973, 980:989, 998:1018, 1075:1099, 1141:1155, 1158:1171, 1183:1202,
		1240:1260, 1315:1327, 1334:1346, 1356:1381, 1431:1447, 1992:2006, 2014:2025,
		2032:2053, 2087:2105, 2260:2279, 2302:2322, 2327:2350, 2360:2382, 2468:2491,
		2740:2760];

seq = readFasta(path)
seqdf = fasta2df(seq)
structdf = assignStructure(seqdf, tms)
positionTMs!(structdf, -1)

loops = @view structdf[structdf.element .== "loop", :]
loopNums = unique(loops.structMap)

yOffset = 0.2
for i in loopNums
    loop = @view loops[loops.structMap .== i, :]
	startRes = loop[1,:resNum]
	n = size(loop, 1)
    prevX = loop.structMap[1] - 0.25
    nxtX = prevX + 1

    ## Check if cytoplasmic or extracellular loop
    if structdf[startRes-1, :y][1] >= 0.5
		side = 1
		prevY = 1
	else
		side = -1
		prevY = 0
	end

    ## Just a line if it is small
	if n <= 5
		dX = (nxtX - prevX)/(n-1)
		for i in 1:size(loop, 1)
			structdf[loop[i, :resNum], :x] = prevX + dX[1] * (i-1)
			structdf[loop[i, :resNum], :y] = prevY + side * yOffset
		end

    ## Square parabola if it is fairly small
    elseif n > 5 && n <= 20
        ## Check if the loop has an odd or even length
        if n%2 != 0
            nBot = 5
        else
            nBot = 4
        end

        ## Divide loop among structure
        leg = div(n - nBot, 2)
        lft = 1:leg
        bot = leg+1:leg+1+nBot
        rgt = n-(leg-1):n

        ## Assign positions
        dY = 0.15
        dX = (nxtX-prevX)/(nBot-1)
        for i in lft
            structdf[loop[i, :resNum], :x] = prevX
            structdf[loop[i, :resNum], :y] = prevY + side*(yOffset + (i-1) * dY)
        end
        for i in bot
            structdf[loop[i, :resNum], :x] = prevX + dX * (i-bot[1])
            structdf[loop[i, :resNum], :y] = prevY + side*(yOffset + (bot[1]-1) * dY)
        end
        structdf[loop[rgt, :resNum], :y] = reverse(loop[lft, :y])
        structdf[loop[rgt, :resNum], :x] .= nxtX
    elseif n > 40 && n <= 120
        loopLayout!(loop, 4, n, prevX, nxtX, prevY, side)
    elseif n > 120 && n <= 240
        loopLayout!(loop, 6, n, prevX, nxtX, prevY, side, expand=2)
    else
        loopLayout!(loop, 12, n, prevX, nxtX, prevY, side, expand=6)
    end
end

f = Figure(resolution=(2000,800))
ax = Axis(f[1, 1], xlabel = "", ylabel = "",
    title = "", grid=:none)
hlines!(ax, 1.1, color=:black, linewidth=2)
hlines!(ax, -0.1, color=:black, linewidth=2)
lines!(structdf.x, structdf.y, strokewidth=1, color = :black)
scatter!(structdf.x, structdf.y, color=:orange, strokewidth=1, strokecolor = :black, markersize=8)
ylims!(ax, -5, 5)
display(f)