text:
push!(LOAD_PATH, "/Users/czhang/Desktop/Projects/dw_/julialib/")
import DimmWitted
DimmWitted.set_libpath("/Users/czhang/Desktop/Projects/dw_/libdw_julia")

immutable DoublePair
	d1::Cdouble
	d2::Cdouble
end

######################################
# The following piece of code creates a
# synthetic data set:
#    - Data type is Cdouble
#    - Modle type is Array{Cdouble}
#
nexp = 100000
nfeat = 1024
examples = Array(DoublePair, nexp, nfeat+1)
for row = 1:nexp
	for col = 1:nfeat
		examples[row, col] = DoublePair(1.0,1.0)
	end
	if rand() > 0.8
		examples[row, nfeat+1] = DoublePair(1.0,0.0)
	else
		examples[row, nfeat+1] = DoublePair(0.0,1.0)
	end
end
model = DoublePair[DoublePair(0.0,0.0) for i = 1:nfeat]

######################################
# Define the loss function and gradient
# function for logistic regression
#
function loss(row::Array{DoublePair,1}, model::Array{DoublePair,1})
	const label1 = row[length(row)].d1
	const label2 = row[length(row)].d2
	const nfeat = length(model)
	d1 = 0.0
	d2 = 0.0
	for i = 1:nfeat
		d1 = d1 + row[i].d1*model[i].d1
		d2 = d2 + row[i].d2*model[i].d2
	end
	v1 = (-label1 * d1 + log(exp(d1) + 1.0))
	v2 = (-label2 * d2 + log(exp(d2) + 1.0))
	return v1 + v2
end

function grad(row::Array{DoublePair,1}, model::Array{DoublePair,1})
	const label1 = row[length(row)].d1
	const label2 = row[length(row)].d2
	const nfeat = length(model)
	d1 = 0.0
	d2 = 0.0
	for i = 1:nfeat
		d1 = d1 + row[i].d1*model[i].d1
		d2 = d2 + row[i].d2*model[i].d2
	end
	d1 = exp(-d1)
	d2 = exp(-d2)
	Z1 = 0.00001 * (-label1 + 1.0/(1.0+d1))
	Z2 = 0.00001 * (-label2 + 1.0/(1.0+d2))
  	for i = 1:nfeat
  		model[i] = DoublePair(model[i].d1 - row[i].d1 * Z1, model[i].d2 - row[i].d2 * Z2)
  	end
	return 1.0
end

######################################
# Create a DimmWitted object using data
# and model. You do not need to specify
# the type, they are infer'ed by the
# open() function, which is parametric.
#
dw = DimmWitted.open(examples, model,
                DimmWitted.MR_PERMACHINE,
                DimmWitted.DR_SHARDING,
                DimmWitted.AC_ROW)

######################################
# Register functions.
#
handle_loss = DimmWitted.register_row(dw, loss, true)
handle_grad = DimmWitted.register_row(dw, grad, true)
#DimmWitted.register_model_avg(dw, handle_loss, avg, true)
#DimmWitted.register_model_avg(dw, handle_grad, avg, true)

######################################
# Run 10 epoches.
#
for iepoch = 1:10
	rs = DimmWitted.exec(dw, handle_loss)
	println("LOSS: ", rs/nexp)
	rs = DimmWitted.exec(dw, handle_grad)
end

sum1 = 0.0
sum2 = 0.0
for i = 1:length(model)
	sum1 = sum1 + model[i].d1
	sum2 = sum2 + model[i].d2
end
println("SUM OF MODEL1: ", sum1)
println("SUM OF MODEL2: ", sum2)

max_stars_repo_name: GloriaGo/dw

max_issues_repo_name: GloriaGo/dw

max_forks_repo_name: GloriaGo/dw


text:
"""
    get_int_param(m::Model, param::Integer)

get the value of some integer valued parameter
"""
get_int_param(m::Model, param::Integer) = get_int_param(m::Model, Cint(param))
function get_int_param(m::Model, param::Cint)

    ipar = convert(Cint,param)
    igval = Array{Cint}(undef, 1)

    ret = @xprs_ccall(getintcontrol, Cint, (Ptr{Nothing},Cint,Ptr{Cint}),
        m.ptr_model, ipar, igval)
    if ret != 0
        throw(XpressError(m))
    end
    convert(Int, igval[1])
end

"""
    get_dbl_param(m::Model, param::Integer)

get the value of some double valued parameter
"""
get_dbl_param(m::Model, param::Integer) = get_dbl_param(m::Model, Cint(param))
function get_dbl_param(m::Model, param::Cint)

    ipar = convert(Cint,param)
    dgval = Array{Float64}(undef, 1)

    ret = @xprs_ccall(getdblcontrol, Cint, (Ptr{Nothing},Cint,Ptr{Cdouble}),
        m.ptr_model, ipar, dgval)
    if ret != 0
        throw(XpressError(m))
    end
    convert(Float64, dgval[1])
end

"""
    get_str_param(m::Model, param::Integer)

get the value of some string valued parameter
"""
get_str_param(m::Model, param::Integer) = get_str_param(m::Model, Cint(param))
function get_str_param(m::Model, param::Cint)

    ipar = convert(Cint,param)
    cgval = Array{Cchar}(undef,  256)

    ret = @xprs_ccall(getstrcontrol, Cint, (Ptr{Nothing},Cint,Ptr{Cchar}),
        m.ptr_model, ipar, cgval)
    if ret != 0
        throw(XpressError(m))
    end

    unsafe_string(pointer(out))
end

"""
    set_int_param(m::Model, ipar::Integer, isval::Integer)

set intger valued parameter
"""
set_int_param(m::Model, ipar::Integer, isval::Integer) = set_int_param(m::Model, Cint(ipar), Cint(isval))
function set_int_param(m::Model, ipar::Cint, isval::Cint)

    ipar = convert(Cint,ipar)
    isval = convert(Cint,isval)
    ret = @xprs_ccall(setintcontrol, Cint, (Ptr{Nothing},Cint,Cint),
        m.ptr_model, ipar, isval)
    if ret != 0
        throw(XpressError(m))
    end
    return nothing
end

"""
    set_dbl_param(m::Model, ipar::Int, dsval::Float64)

set double valued parameter
"""
set_dbl_param(m::Model, ipar::Int, dsval::Float64) = set_dbl_param(m::Model, Cint(ipa), Float64(dsval))
function set_dbl_param(m::Model, ipar::Cint, dsval::Float64)

    ipar = convert(Cint,ipar)
    dsval = convert(Cdouble,dsval)
    ret = @xprs_ccall(setdblcontrol, Cint, (Ptr{Nothing},Cint,Cdouble),
        m.ptr_model, ipar, dsval)
    if ret != 0
        throw(XpressError(m))
    end
    return nothing
end

"""
    set_str_param(m::Model, ipar::Integer, csval::String)

Set string valued parameter
"""
set_str_param(m::Model, ipar::Integer, csval::String) = set_str_param(m::Model, Cint(ipar), csval)
function set_str_param(m::Model, ipar::Cint, csval::String)

    ipar = convert(Cint,ipar)
    csval = convert(String,csval)
    ret = @xprs_ccall(setstrcontrol, Cint, (Ptr{Nothing},Cint,Ptr{Cchar}),
        m.ptr_model, ipar, csval)
    if ret != 0
        throw(XpressError(m))
    end
    return nothing
end

"""
    getparam(m::Model, param::Integer)

get parameter of any type
"""
function getparam(m::Model, param::Integer)
    if convert(Int,param) in XPRS_INT_CONTROLS
        return get_int_param(m, param)
    elseif convert(Int,param) in XPRS_DBL_CONTROLS
        return get_dbl_param(m, param)
    elseif convert(Int,param) in XPRS_STR_CONTROLS
        return get_str_param(m, param)
    else
        error("Unrecognized parameter number: $(param).")
    end
end

"""
    setparam!(m::Model, param::Symbol, val::Any)
    setparam!(m::Model, param::Integer, val::Any)

set parameter of any type
"""
setparam!(m::Model, param::Symbol, val::Any) = setparam!(m,Cint(XPRS_CONTROLS_DICT[param]),val)
setparam!(m::Model, param::Integer, val::Any) = setparam!(m::Model, Cint(param), val::Any)
function setparam!(m::Model, param::Cint, val::Any)
    if convert(Int,param) in XPRS_INT_CONTROLS
        set_int_param(m, param, val)
    elseif convert(Int,param) in XPRS_DBL_CONTROLS
        set_dbl_param(m, param, val)
    elseif convert(Int,param) in XPRS_STR_CONTROLS
        set_str_param(m, param, val)
    else
        error("Unrecognized parameter number: $(param).")
    end
end

"""
    setparams!(m::Model;args...)

set multiple parameters of any type
"""
function setparams!(m::Model;args...)
    for (param,val) in args
        setparam!(m,XPRS_CONTROLS_DICT[param],val)
    end
end

ext: jl


text:
# This file was formerly a part of Julia. License is MIT: https://julialang.org/license

import Base: put!, wait, isready, take!, fetch

mutable struct DictChannel <: AbstractChannel
    d::Dict
    cond_take::Condition    # waiting for data to become available
    DictChannel() = new(Dict(), Condition())
end

function put!(D::DictChannel, k, v)
    D.d[k] = v
    notify(D.cond_take)
    D
end

function take!(D::DictChannel, k)
    v=fetch(D,k)
    delete!(D.d, k)
    v
end

isready(D::DictChannel) = length(D.d) > 1
isready(D::DictChannel, k) = haskey(D.d,k)
function fetch(D::DictChannel, k)
    wait(D,k)
    D.d[k]
end

function wait(D::DictChannel, k)
    while !isready(D, k)
        wait(D.cond_take)
    end
end

max_stars_repo_name: vtjnash/Examples

max_issues_repo_name: JuliaArchive/Examples

max_forks_repo_name: JuliaArchive/Examples


text:
module StochMCMC

using Distributions

include("MH.jl")
include("HMC.jl")
include("SGHMC.jl")

export mcmc, MH, HMC, SGHMC

end


text:
module CoreDefs

export SubSystem, substitute, check_subset, empirical_frequency, dilate, id, draw, embed_aff, DGroupElem
export collar_in, is_interior, frequency, total_collaring, UnrecognizedCollar, transition_matrix
export embed_center
export autocorrelation

using Luxor
using LinearAlgebra

"""
    DGroupElem

A subtype of `DGroupElem` is the type of elements of some group
with a dilation.
It should implement the functions `id`, `*`, `inv` and `dilate`.
"""
abstract type DGroupElem end

"""
    dilate(λ, g)

`embed_aff` dilates an group element `g` by the coefficient `λ`.
"""
function dilate end

"""
    id(g)

Returns the identity element of the group `g` belongs to.
"""
function id end

function Base.:*(g :: G, x :: Pair) where {G<:DGroupElem}
    return g*x[1] => x[2]
end
function Base.:*(g :: G, X :: AbstractVector) where {G<:DGroupElem}
    return map(x -> g*x, X)
end
function Base.:*(g :: G, X :: AbstractSet) where {G<:DGroupElem}
    return map(x -> g*x, X)
end
function Base.:*(g :: G, X :: AbstractDict) where {G<:DGroupElem}
    return typeof(X)(g*x for x in X)
end

function dilate(λ, x :: Pair)
    return dilate(λ,x[1]) => x[2]
end
function dilate(λ, X :: AbstractVector)
    return map(x -> dilate(λ,x), X)
end
function dilate(λ, X :: AbstractSet)
    return map(x -> dilate(λ,x), X)
end

"""
    embed_aff(g)

`embed_aff` should return the embedding of `g` as an element of `Aff(2)`.
This should be a vector length `6` with the three columns of the augmented matrix of the transformation, ignoring its last row.
For example, if `g` should be embedded as a reflection along the x axis followed by a translation by `(1,0)`,
the output should be `[1., 0., 0., -1., 1., 0.]`.
This is the format Luxor uses to represent affine transformations.
"""
function embed_aff end

"""
    draw(t :: T, action :: Symbol)

`T` should be a type of prototiles and `s` should be an action Luxor recognizes when drawing shapes (such as `:fill` or `:stroke`).
`draw(t,s)` should then draw some geometric primitives in a Luxor context using the action `s`.
"""
function draw end

"""
    draw(t :: Pair{<:DGroupElem, T}, sc, hue, action)

If `t` is a tile, `draw(t, sc, hue, action)` draws it using
`draw(t[1], action)` after applying the transformation `embed_aff(t[2])`
with scale `sc`, hue `hue` and action p `action`.

"""
function draw(t :: Pair{<:DGroupElem, T}, sc, hue, action :: Symbol) where T
    scale(sc)
    sethue(hue)
    transform(embed_aff(t[1]))
    draw(t[2], action)
end

"""
    embed_center(t :: Pair{<:DGroupElem, T})

If `t` is a tile with coordinates that can be embedded in the euclidean plane,
outputs the center of the tile as a Luxor point
"""

function embed_center(g :: DGroupElem)
    return Point(embed_aff(g)[5:6]...)
end

"""
    SubSystem{G, D, L}

The description of a substitution system with coordinates in `G`,
labels of type `L`
and dilation coefficients of type `D`.
`G` should be a subset of `DGroupElem` implementing `dilate(λ :: D, g :: G)`.
"""
struct SubSystem{G, D, L}
    sub::Dict{L, Vector{Pair{G, L}}}
    λ :: D
end
VectorTiling{G,L} = AbstractVector{Pair{G, L}}
DictTiling{G,L} = AbstractDict{G, L}

"""
    substitute(S, tiling, n, [in_bounds, window])

Applies `n` inflation steps of the tiling:
each step dilates the tiling and then substitutes each tile according to the substitution rule of S.
This version computes the substitution depth first: this is more efficient unless the substitution rule has lots of overlaps
(see `substitute_bf`)

Optionally, one can filter which tiles are computed by providing `in_bounds` and `window`.
If provided, `in_bounds` should be a test function `in_bounds(tile, n, window)`
which recieves a tile, a depth and the supplied window,
and outputs whether the tile should be substituted further or discarded.
This can be used to greatly reduce compute times when rendering images of tilings.
"""
function substitute(S, tiling, n, in_bounds = nothing, window=nothing)
    result = typeof(tiling)()
    for tile in tiling
        substitute_df_inner!(S, n, tile, result, in_bounds, window)
    end
    return result
end
function substitute_df_inner!(S, n, tile, result, in_bounds = nothing, window=nothing)
    if isnothing(in_bounds) || in_bounds(tile, n, window)
        if n == 0
            push!(result, tile)
            return
        else
            for new_tile in S.sub[tile[2]]
                substitute_df_inner!(S, n-1, (dilate(S.λ, tile[1]) * new_tile), result, in_bounds, window)
            end
        end
    end
end

"""
    substitute_bf(S, tiling,n=1)

Set version of the substitute function: use when the substitution rule has overlaps
"""
function substitute_bf(S, tiling, n=1)
    if n == 0
        return tiling
    else
        step = mapreduce(tile -> (dilate(S.λ, tile[1]) * S.sub[tile[2]]), union, tiling)
        return substitute(S, Set(step), n-1)
    end
end

"""
    center_patch

Translates the patch so that the first tile is located at the identity.
"""
function center_patch(patch)
    centered = inv(patch[1][1])*patch
    @assert centered[1][1] == id(patch[1][1])
    return centered
end

"""
    empirical_frequency(patch, tiling)

Computes the empirical relative frequency of `patch` in the `tiling`:
that is, computes how many translates of `patch` are subsets of `tiling`
divided by the total amount of tiles in `tiling`.
"""
function empirical_frequency(patch :: Dict{G, L}, tiling :: Dict{G, L}) where {G<:DGroupElem, L}
    freq = 0//1
    n = 0

    origin_ptile = patch[id(G)]
    for g in keys(tiling)
        n += 1
        if origin_ptile == tiling[g]
            translated_patch = g*patch
            if translated_patch ⊆ tiling
                freq += 1
            end
        end
    end
    freq = freq//n
    return freq
end
function empirical_frequency(patch :: AbstractArray, tiling :: Dict)
    return empirical_frequency((Dict(patch)), tiling)
end
function empirical_frequency(patch :: AbstractArray, tiling :: AbstractArray)
    freq = 0//1
    n = 0

    origin_ptile = patch[1]
    for (g,l) in tiling
        n += 1
        if (g => l) in tiling
            translated_patch = g*patch
            if translated_patch ⊆ tiling
                freq += 1
            end
        end
    end
    freq = freq//n
    return freq
end

"""
    transition_matrix(S, order)

Given a substitution system `S` and a list of labels `order`,
computes the transition matrix of `S` in that order
"""
function transition_matrix(S, order)
    n = length(order)
    return [count(t -> t[2]==order[i],S.sub[order[j]]) for i=1:n, j=1:n ]
end

struct UnrecognizedCollar <: Exception end
function collar_in end
function is_interior end

function collar_class(tiling, g)
    return inv(g)*collar_in(tiling, g)
end

function center_label(patch)
    e = id(typeof(collect(keys(patch))[1]))
    return patch[e]
end
function center_tile(patch)
    e = id(typeof(collect(keys(patch))[1] ))
    return e => patch[e]
end

"""
    total_collaring(S, initial_collar)

Given a (primitive) substitution system `S` and some collar `initial_collar` (centered at the origin),
construct its total collaring.
The result is given as a tuple `(collars, Sc)`
where `collars` is the list of all collars that happen in `S`-legal tilings
and `Sc` is the corresponding substitution system,
where the labels are integers constituting indexes in `collars`.
"""
function total_collaring(S :: SubSystem{G,D,L}, initial_collar) where {G,D,L}
    collars = [initial_collar]
    visited = 0
    sub = Dict{Int,Vector{Pair{G, Int}}}([])
    while length(collars) > visited
        new_visited = length(collars)
        for i in (visited+1):length(collars)
            sub[i] = []
            patch = substitute(S, collars[i], 1)
            for t in substitute(S, [center_tile(collars[i])], 1)
                class = collar_class(patch, t[1])
                search = findfirst(c -> issetequal(class, c), collars)
                if !isnothing(search)
                    push!(sub[i], t[1] => search)
                else
                    push!(collars, class)
                    push!(sub[i], t[1] => length(collars))
                end
            end
        end
        visited = new_visited
    end
    return (collars, SubSystem(sub, S.λ))
end

"""
Given a (primitive) substitution system `S`,
some collar `initial_collar` and some patch `patch` (both centered at the origin) and some depth `depth`
calculates the frequency of `patch` in `S`.
`initial_collar` is used to construct a total collaring and doesn't matter as long as it's a legal collar.
`depth` determines how the level at which the frequency is calculated:
the result is only exact if `depth` is high enough, but lower values result in faster computation.

"""
function frequency(S :: SubSystem{G, D, L}, initial_collar, patch, depth) where {G, D, L}
    (collars, Sc) = total_collaring(S, initial_collar)
    n = length(collars)
    A_tr = transition_matrix(Sc, 1:n)
    (eigenvalues, eigenvectors) = eigen(A_tr)
    λ_PF = eigenvalues[n]
    v_PF = eigenvectors[:,n]/sum(eigenvectors[:,n])

    ptiles = unique([collar[id(G)] for collar in collars])
    collars_of_ptile = Dict([ptile => Int[] for ptile in ptiles])
    for i = 1:n
        push!(collars_of_ptile[collars[i][id(G)]], i)
    end

    patch_c = Dict([])
    for (g,l) in patch
        if is_interior(patch, g)
            patch_c[g] = (:interior, l)
        else
            patch_c[g] = (:exterior, collars_of_ptile[l])
        end
    end

    freq = 0.0
    for label in 1:n
        domain = substitute(Sc, [id(G) => label], depth-1)
        forced_uncollared_domain = substitute(S, collars[label], depth-1)

        for tile in domain
            translated_patch_c = tile[1] * patch_c
            is_subset = true
            for (g, (kind, detect)) in translated_patch_c
                if kind == :interior
                    if (g => detect) ∉ forced_uncollared_domain
                        is_subset = false
                        break
                    end
                else
                    if !any(c -> c ⊆ forced_uncollared_domain, Ref(g).*collars[detect])
                        is_subset = false
                        break
                    end
                end
            end
            if is_subset
                freq += v_PF[label]/λ_PF^(depth-1)
            end
        end
    end
    return freq
end

# Radius should be forces by the chosen prototile depth!
function autocorrelation(S :: SubSystem{G, D, L}, initial_collar, depth) where {G, D, L}
    (collars, Sc) = total_collaring(S, initial_collar)
    n = length(collars)
    A_tr = transition_matrix(Sc, 1:n)
    (eigenvalues, eigenvectors) = eigen(A_tr)
    λ_PF = eigenvalues[n]
    v_PF = eigenvectors[:,n]/sum(eigenvectors[:,n])

    ptiles = unique([collar[id(G)] for collar in collars])
    collars_of_ptile = Dict([ptile => Int[] for ptile in ptiles])
    for i = 1:n
        push!(collars_of_ptile[collars[i][id(G)]], i)
    end

    measure = Dict{G, Float64}()
    for label in 1:n
        domain = substitute(S, [id(G) => center_label(collars[label])], depth-1)
        forced_domain = substitute(S, collars[label], depth-1)

        for tile in domain
            translated_f_domain = inv(tile[1]) * forced_domain
            for t in translated_f_domain
                g = t[1]
                if g in keys(measure)
                    measure[g] += v_PF[label]/λ_PF^(depth-1)
                else
                    measure[g] = v_PF[label]/λ_PF^(depth-1)
                end
            end
        end
    end
    return measure
end

end


text:
using Documenter
using SustainableFinance

push!(LOAD_PATH,"../src/")
makedocs(
    sitename = "Julia for Sustainable Finance",
    format = Documenter.HTML(),
    modules = [SustainableFinance];
    pages = [
        "Home" => "index.md",
        "ESG Valuation" => Any[
            "Introduction"=>"ESG Valuation/intro.md",
            "GHG Emissions"=>"ESG Valuation/Environmentals/Climate/impacts.md",
            "Water Consumption"=>"ESG Valuation/Environmentals/Water Consumption/impacts.md",
        ],

    ]
)

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

max_stars_repo_name: tlorans/SustainableFinance

max_issues_repo_name: tlorans/SustainableFinance

max_forks_repo_name: tlorans/SustainableFinance


text:
symmetrize(P::Matrix) = (P+P')/2 #+ 1e-15eye(P) #TODO: Hacky. Square root filter?

function likelihood{T}(y::Vector{T}, C::Matrix{T}, σϵ²::T, Q::Matrix{T})

    n   = length(y)
    nₓ  = size(C, 1)

    x  = zeros(nₓ)
    P∞ = eye(nₓ)
    P = zeros(nₓ, nₓ)

    ll = -n*log(2π)

    d = 0
    diffuseP = true

    while diffuseP

        d += 1
        Cₜ = C[:,d]

        ν⁰ = y[d] - dot(Cₜ, x)

        F∞  = dot(Cₜ, P∞*Cₜ)
        F   = dot(Cₜ, P*Cₜ) + σϵ²
        F¹  = 1/F∞
        F²  = - F¹ * F * F¹
        posdefF∞ = !isapprox(0, F∞, atol=1e-10)

        M∞  = P∞ * Cₜ
        M   = P * Cₜ

        K⁰  = posdefF∞ ? M∞*F¹ : M/F
        K¹  = M * F¹ + M∞ * F²

        L⁰  = I - K⁰ * Cₜ'
        L¹  = -K¹ * Cₜ'

        # Calculate loglikelihood contribution
        ll-= posdefF∞ ? log(F∞) : log(F) + ν⁰*ν⁰/F

        # Calculate next values
        x   = x + K⁰*ν⁰
        P   = posdefF∞ ? symmetrize(P∞*L¹' + P*L⁰' + Q) : symmetrize(P*L⁰' + Q)
        P∞  = posdefF∞ ? symmetrize(P∞*L⁰') : P∞

        diffuseP = findfirst(round(P∞,10)) > 0

        d == n && error("Failed to substantiate diffuse initial conditions")

    end #while

    for t in d+1:n-1

        Cₜ = C[:,t]

        νₜ = y[t] - dot(Cₜ, x)
        Fₜ = dot(Cₜ, P*Cₜ) + σϵ²
        Kₜ = P * Cₜ / Fₜ

        # Calculate loglikelihood contribution
        ll -= Cₜ[1] == 0 ? 0 : log(Fₜ) + νₜ*νₜ/Fₜ

        # Predict next values
        x += Kₜ * νₜ
        P = P * (I - Cₜ * Kₜ') + Q

    end #for

    Cₜ = C[:,n]
    νₜ = y[n] - dot(Cₜ, x)
    Fₜ = dot(Cₜ, P*Cₜ) + σϵ²
    ll -= log(Fₜ) + νₜ*νₜ/Fₜ

    return -ll/2

end #likelihood

function smooth{T}(y::Vector{T}, C::Matrix{T}, σϵ²::T, σμ²::T)

    n   = length(y)
    nₓ  = size(C, 1)
    Q = zeros(nₓ, nₓ)
    Q[1] = σμ²

    K   = Array{T}(nₓ, n)
    F⁻¹ = Array{T}(n)
    ν   = Array{T}(n)

    P∞  = Matrix{T}[]
    L⁰  = Matrix{T}[]
    L¹  = Matrix{T}[]
    F∞  = T[]
    F²  = T[]

    x   = Array{T}(nₓ, n)
    P   = Array{T}(nₓ, nₓ, n)
    ϵ   = Array{T}(n)
    η   = Array{T}(nₓ, n)

    xₜ  = zeros(nₓ)
    P∞ₜ = eye(nₓ)
    Pₜ  = zeros(nₓ, nₓ)

    ll = -n*log(2π)

    d           = 0
    x[:,1]      = xₜ
    P[:,:,1]    = Pₜ
    push!(P∞, P∞ₜ)
    diffuseP    = true

    while diffuseP

        d += 1
        Cₜ = C[:,d]
        ν⁰ = y[d] - dot(Cₜ, xₜ)

        F∞ₜ = dot(Cₜ, P∞ₜ * Cₜ)
        Fₜ  = dot(Cₜ, Pₜ * Cₜ) + σϵ²
        F¹ₜ = 1/F∞ₜ
        F²ₜ = - F¹ₜ * Fₜ * F¹ₜ

        M∞  = P∞ₜ * Cₜ
        M   = Pₜ * Cₜ

        K⁰  = F∞ₜ == 0 ? M/Fₜ : M∞*F¹ₜ
        K¹ = M * F¹ₜ + M∞ * F²ₜ

        L⁰ₜ = I - K⁰ * Cₜ'
        L¹ₜ = -K¹ * Cₜ'

        ν[d]    = ν⁰
        F⁻¹[d]  = 1/Fₜ
        K[:, d] = K⁰
        push!(F∞, F∞ₜ)
        push!(F², F²ₜ)
        push!(L⁰, L⁰ₜ)
        push!(L¹, L¹ₜ)

        # Calculate next values
        xₜ  = xₜ + K⁰ * ν⁰
        Pₜ  = F∞ₜ == 0 ? Pₜ*L⁰ₜ' + Q : P∞ₜ*L¹ₜ' + Pₜ*L⁰ₜ' + Q
        P∞ₜ = F∞ₜ == 0 ? P∞ₜ         : P∞ₜ*L⁰ₜ'

        x[:,d+1]    = xₜ
        P[:,:,d+1]  = Pₜ
        push!(P∞, P∞ₜ)

        diffuseP = findfirst(round(P∞ₜ,12)) > 0

    end #while

    for t in d+1:n-1

        Cₜ = C[:,t]
        νₜ = y[t] - dot(Cₜ, xₜ)
        Fₜ = dot(Cₜ, Pₜ*Cₜ) + σϵ²
        Kₜ = Pₜ * Cₜ / Fₜ

        ν[t]    = νₜ
        F⁻¹[t]  = 1/Fₜ
        K[:,t]  = Kₜ

        # Predict next values
        xₜ += Kₜ * νₜ
        Pₜ = Pₜ * (I - Cₜ * Kₜ') + Q

        x[:,t+1]   = xₜ
        P[:,:,t+1] = Pₜ

    end #for

    Cₜ = C[:,n]
    νₜ = y[n] - dot(Cₜ, xₜ)
    Fₜ = dot(Cₜ, Pₜ*Cₜ) + σϵ²
    Kₜ = Pₜ * Cₜ / Fₜ

    ν[n]    = νₜ
    F⁻¹[n]  = 1/Fₜ
    K[:,n]  = Kₜ

    rₜ, Nₜ = zeros(nₓ), zeros(nₓ, nₓ)
    uₜ = νₜ/Fₜ - dot(Kₜ, rₜ)
    Dₜ = 1/Fₜ + dot(Kₜ, Nₜ * Kₜ)

    # Calculate smoothed disturbances
    ϵ[n] = σϵ² * uₜ
    η[:, n] = Q * rₜ

    # Calculate helper recursions
    rₜ = Cₜ * uₜ + rₜ
    Θ = Cₜ * Kₜ' * Nₜ
    Nₜ = Cₜ * Dₜ * Cₜ' + Nₜ - Θ - Θ'

    for t in n-1:-1:d+1

        Cₜ    = C[:,t]
        F⁻¹ₜ  = F⁻¹[t]
        νₜ    = ν[t]
        Kₜ    = K[:,t]
        Pₜ    = P[:,:,t]

        uₜ = F⁻¹ₜ * νₜ - dot(Kₜ, rₜ)
        Dₜ = F⁻¹ₜ + dot(Kₜ, Nₜ * Kₜ)

        # Calculate smoothed disturbances
        ϵ[t] = σϵ² * uₜ
        η[:, t] = Q * rₜ

        # Calculate helper recursions
        rₜ = Cₜ * uₜ + rₜ
        Θ = Cₜ * Kₜ' * Nₜ
        Nₜ = Cₜ * Dₜ * Cₜ' + Nₜ - Θ - Θ'

        # Calculate smoothed states
        x[:, t] += Pₜ * rₜ
        P[:,:,t] -= Pₜ * Nₜ * Pₜ

    end #for

    r⁰ₜ, r¹ₜ      = rₜ, zeros(rₜ)
    N⁰ₜ, N¹ₜ, N²ₜ = Nₜ, zeros(Nₜ), zeros(Nₜ)

    for t in d:-1:1

        Cₜ        = C[:,t]
        F∞ₜ       = F∞[t]
        #F¹ₜ       = F⁻¹[t]
        F²ₜ       = F²[t]
        F⁻¹ₜ      = F∞ₜ == 0 ? F⁻¹[t] : 0
        ν⁰ₜ       = ν[t]
        K⁰ₜ       = K[:,t]
        L⁰ₜ, L¹ₜ  = L⁰[t], L¹[t]
        Pₜ, P∞ₜ   = P[:,:,t], P∞[t]

        uₜ        = F⁻¹ₜ*ν⁰ₜ - dot(K⁰ₜ, r⁰ₜ)
        Dₜ        = F⁻¹ₜ + dot(K⁰ₜ, N⁰ₜ * K⁰ₜ)

        # Calculate smoothed disturbances
        ϵ[t] = σϵ² * uₜ
        η[:, t] = Q * rₜ

        # Calculate helper recursions
        r¹ₜ = Cₜ*ν⁰ₜ/F∞ₜ + L⁰ₜ'*r¹ₜ + L¹ₜ'*r⁰ₜ
        r⁰ₜ += Cₜ * uₜ
        Δ = L⁰ₜ' * N¹ₜ * L¹ₜ
        N²ₜ = Cₜ*F²ₜ*Cₜ' + L⁰ₜ'*N²ₜ*L⁰ₜ + Δ + Δ' + L¹ₜ'*N⁰ₜ*L¹ₜ
        N¹ₜ = Cₜ*Cₜ'/F∞ₜ + L⁰ₜ'*N¹ₜ*L⁰ₜ + L¹ₜ'*N⁰ₜ*L⁰ₜ
        Θ = Cₜ*K⁰ₜ'*N⁰ₜ
        N⁰ₜ += Cₜ*Dₜ*Cₜ' - Θ - Θ'

        # Calculate smoothed states
        Γ = P∞ₜ * N¹ₜ * Pₜ
        x[:,t] += Pₜ*r⁰ₜ + P∞ₜ*r¹ₜ
        P[:,:,t] -= Pₜ*N⁰ₜ*Pₜ + P∞ₜ*N²ₜ*P∞ₜ + Γ + Γ'

    end #for

    return ϵ, x, P, η

end #statesmooth


text:
function computetables(exponents::Matrix{T}) where {T<:Integer}
    m, n = size(exponents)
    lookuptable = permutedims(exponents, (2, 1))

    for j=1:m
        sort!(@view lookuptable[:, j])
    end

    j = 1
    maxk = 0
    for j=1:m
        k = 1
        for i=2:n
            if lookuptable[i, j] > lookuptable[k, j]
                k += 1
                lookuptable[k, j] = lookuptable[i, j]
            end
            if k < i
                lookuptable[i, j] = 0
            end
        end
        maxk = max(maxk, k)
    end
    # now we copy the relevant subsection
    # we do not compute the differences now since we need to compute the lookuptable first
    differences = Matrix{UInt8}(undef, m, maxk)
    for i=1:m, j=1:maxk
        differences[i, j] = lookuptable[j, i]
    end

    lookuptable = permutedims(lookuptable, (2, 1))
    for j=1:n, i=1:m
        for k = 1:maxk
            if differences[i, k] == exponents[i, j]
                lookuptable[i, j] = k
                break
            end
        end
    end

    # now we can compute the differences
    for i=1:m
        lastval = differences[i, 1]
        for j=2:maxk
            val = differences[i, j]
            if val > 0
                differences[i, j] = val - lastval
                lastval = val
            end
        end
    end
    differences, lookuptable
end

function computetables(exponents_vec::Vector{Matrix{T}}) where {T<:Integer}
    widths = size.(exponents_vec, 2)
    # we fake our exponents_vec as a single big matrix, compute the differences and a
    # big lookuptable and then split it again
    allexponents = hcat(exponents_vec...)
    differences, biglookuptable = computetables(allexponents)

    m = size(biglookuptable, 1)

    total = 1
    lookuptables = Vector{Matrix{T}}()
    for n in widths
        push!(lookuptables, biglookuptable[:, total : total + n - 1])
        total += n
    end
    differences, lookuptables
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 42549fee-2c45-11ec-3ddc-8db7e6f0c642
begin
	using PlutoUI
	using DataFrames
end

# ╔═╡ 5e04b370-80e4-4f41-8c60-27d4cc76488b
md"""
#### Solution: Modified VCRC Practice Problem
"""

# ╔═╡ b0a49c0e-476c-4aea-b36f-6a21a0ac5976
md"""
$(PlutoUI.LocalResource("./figs/Fig-ModVCRC-Work-Recycle-F20.png"))
"""

# ╔═╡ d863255d-217b-491a-8939-23c95a90de02
md"""
a) Compute the missing state values in the state table.

Look up values on the PH-diagram for the particular working fluid in this system (R134A).

"""

# ╔═╡ 36521d17-881c-4ce0-9483-cb4040ad4369
md"""$(PlutoUI.LocalResource("./figs/Fig-R134A-PH-Diagram-MVCRC.png"))"""

# ╔═╡ 9d05fdfc-bc82-43dd-838d-dae815b29926
md"""
b) Compute the missing heat and work rate values in the path table.

Use the open first and second laws around each of the process units to compute the missing values in the path table
"""

# ╔═╡ 4e55dacc-c131-4f29-ad74-42838539d332
md"""
c) Compute the reversible (ideal) coefficient of performance for the modified cycle using the values in the path table

The coefficient of performance (COP) for the VCRC is given by:

$$\omega = \dot{Q}_{c}/\dot{W}_{s,C}$$

We know all these values from the path table.
"""

# ╔═╡ f0b29ce1-28ec-46bf-a8e5-d77df196dbb5
md"""
d) Would you expect the ideal coefficient of performance to increase, decrease or stay the same if we did not recycle $\dot{W}_{s,T}$?

You should expect the COP to _decrease_ if we did not recycle the work from the turbine. Why? Because the recycled work _decreases_ the amount of additional work that we need to put to move a unit of heat. Thus, if we do _not_ have the recycled work, then you have to put more external work in (to move the same $\dot{Q}$), hence, $\omega$ will decrease.
"""

# ╔═╡ ddfdfb31-35e1-4ac0-8703-a4a75b67d40d
md"""
##### Julia functions for this problem below are defined below here
"""

# ╔═╡ ad1283c7-a2d9-4bc8-881c-f70a5a42f723
function _compute_COP(pathTable::DataFrame)

    QC = pathTable[1,:Qs_dot]
    WC = pathTable[2,:Wsc_dot]
    COP = (QC/WC)
    return COP
end

# ╔═╡ cd4e7dd8-a34b-4ec3-a4d4-2831f67797e2
function _build_default_state_table()

    # initialize with problem data -
    operating_points = [1,2,3,4]
    T = [
        -60.0       ;   # O1
        -60.0       ;   # O2
        95.0        ;   # O3
        70.0        ;   # O4
    ];

    P = [
        15.89       ;   # O1
        15.89       ;   # O2
        2117.34     ;   # O3
        2117.34     ;   # O4
    ];

    H = [
        0.0         ;   # O1
        361.7       ;   # O2
        464.8       ;   # O3
        304.8       ;   # O4
    ];

    S = [
        1.3347      ;   # O1
        1.8024      ;   # O2
        1.8024      ;   # O3
        1.3347      ;   # O4
    ];

    quality = [
        0.0         ;   # O1
        1.0         ;   # O2
        1.0         ;   # O3
        0.0         ;   # O4
    ];

    # build DataFrame -
    df = DataFrame(O=operating_points,T=T,P=P,H=H,S=S,quality=quality);

    # return -
    return df
end

# ╔═╡ b58cd3e8-31bb-40a3-a7ea-bba17f64f381
function _build_final_state_table(dataTable::DataFrame)

    # update -
    final_state_table = copy(dataTable)

    # we need to compute the quality
    # what are the endpoints at -60C
    Hl = 124.0
    Hv = 361.7
    Sl = 0.6873
    Sv = 1.8024

    # quality -
    S1 = dataTable[4,:S] # H1 = H4
    theta = (S1 - Sl)/(Sv - Sl)

    # update the quality -
    final_state_table[1,:quality] = theta

    # Compute S1 -
    H1 = (1.0 - theta)*Hl+theta*Hv
    final_state_table[1,:H] = H1

    # return -
    return final_state_table
end

# ╔═╡ 940da569-1f17-48e4-9e40-b0588ea47b46
# compute the remainder of the missing elements in the state table -
state_table = _build_final_state_table(_build_default_state_table())

# ╔═╡ f0328fc0-859f-463f-af1c-19c18dc93eef
function _build_default_path_table()

    # setup -
    operating_points = ["P12","P23","P34","P41","Cycle"]
    Qs_dot = [
        0.0     ;   # P1
        0.0     ;   # P2
        0.0     ;   # P3
        0.0     ;   # P4
        0.0     ;   # Cycle
    ];

    Wst_dot = [
        0.0     ;   # P1
        0.0     ;   # P2
        0.0     ;   # P3
        0.0     ;   # P4
        0.0     ;   # Cycle
    ];

    Wsc_dot = [
        0.0     ;   # P1
        0.0     ;   # P2
        0.0     ;   # P3
        0.0     ;   # P4
        0.0     ;   # Cycle
    ];

    Wst_star_dot = [
        0.0     ;   # P1
        0.0     ;   # P2
        0.0     ;   # P3
        0.0     ;   # P4
        0.0     ;   # Cycle
    ];

    Wsc_star_dot = [
        0.0     ;   # P1
        0.0     ;   # P2
        0.0     ;   # P3
        0.0     ;   # P4
        0.0     ;   # Cycle
    ];

    # build DataFrame -
    df = DataFrame(P=operating_points,Qs_dot=Qs_dot,Wst_dot=Wst_dot,Wsc_dot=Wsc_dot,
            Wst_star_dot=Wst_star_dot,Wsc_star_dot=Wsc_star_dot);

    # return -
    return df
end

# ╔═╡ 4e8a57d9-5c8e-421b-b823-4b1a4cff1eba
mutable struct VCRCProblem

    mdot::Float64
    eta_T::Float64
    eta_C::Float64

    VCRCProblem(mdot,eta_T,eta_C) = new(mdot,eta_T,eta_C)
end

# ╔═╡ 46103655-b977-4874-b861-7447095ac877
# before we do anything, lets build an instance of the problem -
problem = VCRCProblem(0.167,0.75,0.81)

# ╔═╡ 3585b45e-3f0a-4150-86d2-5ef204a816ff
function _build_final_path_table(pathTable::DataFrame, stateTable::DataFrame, problem::VCRCProblem)

    # initialize -
    final_path_table = copy(pathTable)

    # get stuff from problem -
    mdot = problem.mdot
    eta_T = problem.eta_T
    eta_C = problem.eta_C

    # fill in missing values -
    # P12 - qdot
    H2 = stateTable[2,:H]
    H1 = stateTable[1,:H]
    final_path_table[1,:Qs_dot] = mdot*(H2 - H1)

    # P34 - qdot
    H3 = stateTable[3,:H]
    H4 = stateTable[4,:H]
    final_path_table[3,:Qs_dot] = mdot*(H4 - H3)

    # P41 - Wst_dot
    H1 = stateTable[1,:H]
    H4 = stateTable[4,:H]
    final_path_table[4,:Wst_dot] = mdot*(H1 - H4)

    # P23 - Wst_dot
    H1 = stateTable[1,:H]
    H4 = stateTable[4,:H]
    final_path_table[2,:Wst_dot] = -1*mdot*(H1 - H4)

    # P23 - Wsc_dot
    H2 = stateTable[2,:H]
    H3 = stateTable[3,:H]
    final_path_table[2,:Wsc_dot] = mdot*(H3 - H2) - final_path_table[2,:Wst_dot]

    # real work -
    W_ideal = final_path_table[4,:Wst_dot]
    final_path_table[4,:Wst_star_dot] = eta_T*W_ideal
    final_path_table[2,:Wst_star_dot] = -1*eta_T*W_ideal

    W_ideal = final_path_table[2,:Wsc_dot]
    final_path_table[2,:Wsc_star_dot] = (1/eta_C)*W_ideal

    # cycle Q -
    cycle_Q = sum(final_path_table[1:4,:Qs_dot])
    final_path_table[5,:Qs_dot] = cycle_Q

    # cycle W -
    cycle_W = sum(final_path_table[1:4,:Wst_dot])
    final_path_table[5,:Wst_dot] = cycle_W

    cycle_W = sum(final_path_table[1:4,:Wsc_dot])
    final_path_table[5,:Wsc_dot] = cycle_W

    # return -
    return final_path_table
end

# ╔═╡ dccb6edb-3077-4d59-bbdb-a87a8528d4dc
begin

	# build the default path table -
    default_path_table = _build_default_path_table();

    # compute the remainder of the path table -
    path_table = _build_final_path_table(default_path_table, state_table, problem)
end

# ╔═╡ 1b5d5556-f391-4865-99c8-a169a24f5d37
COP = _compute_COP(path_table)

# ╔═╡ 182da5e8-6648-4058-92b4-93f4bb89541c
html"""
<style>
main {
    max-width: 1200px;
    width: 90%;
    margin: auto;
    font-family: "Roboto, monospace";
}

a {
    color: blue;
    text-decoration: none;
}

.H1 {
    padding: 0px 30px;
}
</style>
"""

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"

[compat]
DataFrames = "~1.2.2"
PlutoUI = "~0.7.16"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

[[ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[Compat]]
deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
git-tree-sha1 = "31d0151f5716b655421d9d75b7fa74cc4e744df2"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.39.0"

[[Crayons]]
git-tree-sha1 = "3f71217b538d7aaee0b69ab47d9b7724ca8afa0d"
uuid = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f"
version = "4.0.4"

[[DataAPI]]
git-tree-sha1 = "cc70b17275652eb47bc9e5f81635981f13cea5c8"
uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
version = "1.9.0"

[[DataFrames]]
deps = ["Compat", "DataAPI", "Future", "InvertedIndices", "IteratorInterfaceExtensions", "LinearAlgebra", "Markdown", "Missings", "PooledArrays", "PrettyTables", "Printf", "REPL", "Reexport", "SortingAlgorithms", "Statistics", "TableTraits", "Tables", "Unicode"]
git-tree-sha1 = "d785f42445b63fc86caa08bb9a9351008be9b765"
uuid = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
version = "1.2.2"

[[DataStructures]]
deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
git-tree-sha1 = "7d9d316f04214f7efdbb6398d545446e246eff02"
uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
version = "0.18.10"

[[DataValueInterfaces]]
git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6"
uuid = "e2d170a0-9d28-54be-80f0-106bbe20a464"
version = "1.0.0"

[[Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[DelimitedFiles]]
deps = ["Mmap"]
uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"

[[Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

[[Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[Formatting]]
deps = ["Printf"]
git-tree-sha1 = "8339d61043228fdd3eb658d86c926cb282ae72a8"
uuid = "59287772-0a20-5a39-b81b-1366585eb4c0"
version = "0.4.2"

[[Future]]
deps = ["Random"]
uuid = "9fa8497b-333b-5362-9e8d-4d0656e87820"

[[Hyperscript]]
deps = ["Test"]
git-tree-sha1 = "8d511d5b81240fc8e6802386302675bdf47737b9"
uuid = "47d2ed2b-36de-50cf-bf87-49c2cf4b8b91"
version = "0.0.4"

[[HypertextLiteral]]
git-tree-sha1 = "f6532909bf3d40b308a0f360b6a0e626c0e263a8"
uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
version = "0.9.1"

[[IOCapture]]
deps = ["Logging", "Random"]
git-tree-sha1 = "f7be53659ab06ddc986428d3a9dcc95f6fa6705a"
uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89"
version = "0.2.2"

[[InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[InvertedIndices]]
git-tree-sha1 = "bee5f1ef5bf65df56bdd2e40447590b272a5471f"
uuid = "41ab1584-1d38-5bbf-9106-f11c6c58b48f"
version = "1.1.0"

[[IteratorInterfaceExtensions]]
git-tree-sha1 = "a3f24677c21f5bbe9d2a714f95dcd58337fb2856"
uuid = "82899510-4779-5014-852e-03e436cf321d"
version = "1.0.0"

[[JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "8076680b162ada2a031f707ac7b4953e30667a37"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.2"

[[LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[LinearAlgebra]]
deps = ["Libdl"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[Missings]]
deps = ["DataAPI"]
git-tree-sha1 = "bf210ce90b6c9eed32d25dbcae1ebc565df2687f"
uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
version = "1.0.2"

[[Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[OrderedCollections]]
git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
version = "1.4.1"

[[Parsers]]
deps = ["Dates"]
git-tree-sha1 = "a8709b968a1ea6abc2dc1967cb1db6ac9a00dfb6"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.0.5"

[[Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[PlutoUI]]
deps = ["Base64", "Dates", "Hyperscript", "HypertextLiteral", "IOCapture", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "UUIDs"]
git-tree-sha1 = "4c8a7d080daca18545c56f1cac28710c362478f3"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.16"

[[PooledArrays]]
deps = ["DataAPI", "Future"]
git-tree-sha1 = "a193d6ad9c45ada72c14b731a318bedd3c2f00cf"
uuid = "2dfb63ee-cc39-5dd5-95bd-886bf059d720"
version = "1.3.0"

[[PrettyTables]]
deps = ["Crayons", "Formatting", "Markdown", "Reexport", "Tables"]
git-tree-sha1 = "69fd065725ee69950f3f58eceb6d144ce32d627d"
uuid = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
version = "1.2.2"

[[Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[Random]]
deps = ["Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[SharedArrays]]
deps = ["Distributed", "Mmap", "Random", "Serialization"]
uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"

[[Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[SortingAlgorithms]]
deps = ["DataStructures"]
git-tree-sha1 = "b3363d7460f7d098ca0912c69b082f75625d7508"
uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c"
version = "1.0.1"

[[SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[TableTraits]]
deps = ["IteratorInterfaceExtensions"]
git-tree-sha1 = "c06b2f539df1c6efa794486abfb6ed2022561a39"
uuid = "3783bdb8-4a98-5b6b-af9a-565f29a5fe9c"
version = "1.0.1"

[[Tables]]
deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "TableTraits", "Test"]
git-tree-sha1 = "fed34d0e71b91734bf0a7e10eb1bb05296ddbcd0"
uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
version = "1.6.0"

[[Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╟─5e04b370-80e4-4f41-8c60-27d4cc76488b
# ╟─b0a49c0e-476c-4aea-b36f-6a21a0ac5976
# ╠═46103655-b977-4874-b861-7447095ac877
# ╟─d863255d-217b-491a-8939-23c95a90de02
# ╠═940da569-1f17-48e4-9e40-b0588ea47b46
# ╟─36521d17-881c-4ce0-9483-cb4040ad4369
# ╟─9d05fdfc-bc82-43dd-838d-dae815b29926
# ╠═dccb6edb-3077-4d59-bbdb-a87a8528d4dc
# ╟─4e55dacc-c131-4f29-ad74-42838539d332
# ╠═1b5d5556-f391-4865-99c8-a169a24f5d37
# ╟─f0b29ce1-28ec-46bf-a8e5-d77df196dbb5
# ╟─ddfdfb31-35e1-4ac0-8703-a4a75b67d40d
# ╠═ad1283c7-a2d9-4bc8-881c-f70a5a42f723
# ╠═cd4e7dd8-a34b-4ec3-a4d4-2831f67797e2
# ╠═b58cd3e8-31bb-40a3-a7ea-bba17f64f381
# ╠═3585b45e-3f0a-4150-86d2-5ef204a816ff
# ╠═f0328fc0-859f-463f-af1c-19c18dc93eef
# ╠═4e8a57d9-5c8e-421b-b823-4b1a4cff1eba
# ╠═42549fee-2c45-11ec-3ddc-8db7e6f0c642
# ╟─182da5e8-6648-4058-92b4-93f4bb89541c
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
using jInv.InverseSolve

n = 100 * 3

dc   = rand(n).+2
dobs = rand(n).+3
wd   = rand(n).+1
dd   = rand(n)

for misfitFun in [ SSDFun HuberFun ]

	# Gradient test for ", misfitFun, " (real)")

	function testFun(x,v=[])
		mis,dmis,d2mis = misfitFun(x,dobs,wd)
		if isempty(v)
			return mis
		else
			return mis,dot(dmis,v)
		end
	end

	chkDer, = checkDerivative(testFun,dc,out=false)
	@test chkDer

end

for misfitFun in [ SSDFun ]

	# Hessian test for ", misfitFun, " (real)")

	function testFun(x,v=[])
		mis,dmis,d2mis = misfitFun(x,dobs,wd)
		if isempty(v)
			return dmis
		else
			return dmis,d2mis.*v
		end
	end

	chkDer, = checkDerivative(testFun,dc,out=false)
	@test chkDer
end

### Complex

dc   += 1im*rand(n)
dobs += 1im*rand(n)
wd   += 1im*rand(n)
dd   += 1im*rand(n)

for misfitFun in [ SSDFun ]

	# Gradient test for ", misfitFun, " (complex)")

	function testFun(x,v=[])
		mis,dmis,d2mis = misfitFun(x,dobs,wd)
		if isempty(v)
			return mis
		else
			return mis,real(dot(dmis,v))
		end
	end

	chkDer, = checkDerivative(testFun,dc,out=false)
	@test chkDer

end

# for misfitFun in [ SSDFun ]

	# Hessian test for ", misfitFun, " (complex)")
	# warn("skipped"); break
# 	function testFun(x,v=[])
# 		mis,dmis,d2mis = misfitFun(x,dobs,wd)
# 		if isempty(v)
# 			return real(dmis)
# 		else
# 			return real(dmis),real(d2mis.*v)
# 		end
# 	end
#
# 	chkDer, = checkDerivative(testFun,dc,out=true)
# 	@test chkDer
#
# end


text:
module Router

import OrderedCollections

import Dance.Configuration
import Dance.Logger
import Dance.URIUtils

include("mime_types.jl")

const GET= "GET"
const OPTIONS = "OPTIONS"
const POST= "POST"
const METHODS = [GET, OPTIONS, POST]

const EP_HTML = "HTML"
const EP_JSON = "JSON"
const EP_STATIC = "STATIC"
const ENDPOINTS = [EP_HTML, EP_JSON, EP_STATIC]

STATIC_DIR = ""
STATIC_ROUTE_PREFIX = ""

export GET, POST, EP_HTML, EP_JSON, EP_STATIC, output_file_as_string, route, route_group, static_dir

mutable struct Route
    endpoint::String
    method::String
    path::String
    has_regex::Bool
    action::Function
    html_file::String

    """
        Route(;endpoint::String, method::String, path::String, has_regex::Bool, action::Function, html_file::String)

    Create new Route object, by validating supplied Route parameters
    """
    function Route(;endpoint::String, method::String, path::String, has_regex::Bool, action::Function, html_file::String) :: Route
        error::String = ""

        if !(endpoint in ENDPOINTS)
            error = "Invalid endpoint. Must be \"EP_HTML\" or \"EP_JSON\" or \"EP_STATIC\""
        elseif !(method in METHODS)
            error = "Invalid method. Must be \"GET\" or \"POST\""
        elseif path_already_exists(path)
            error = "Path already exists"
        else
            try
                for item in split(path, "/")
                    isa("r\"" * path * "\"", Regex)
                end
            catch e
                error = "Invalid route regex format"
            end
            try
                methods(action)
            catch e
                error = "Action function not defined"
            end
        end

        if length(error)==0
            return new(endpoint, method, path, has_regex, action, html_file)
        else
            Logger.log("Route construction for `$path`: $error")
        end
    end
end

const ROUTES = OrderedCollections.OrderedDict{String, Route}()

"""
    delete_routes!()

Empty ROUTES ordered dict

Only called from tests as of now
"""
function delete_routes!() :: OrderedCollections.OrderedDict{String, Route}
    empty!(ROUTES)
end

"""
    get_route(route_path::String)

Search ROUTES ordered dict for route component from specified path
"""
function get_route(route_path::String) :: Union{Route, Nothing}
    found_route::Union{Route, Nothing} = nothing

    # Remove trailing slash (if not index url)
    route_path = remove_trailing_slash(route_path)

    for (path, route) in ROUTES
        if !route.has_regex
            if path==route_path
                found_route = route
            end
        else
            request_path_number_of_segments::Int8 = length(URIUtils.get_path_segments(route_path))[1]
            route_number_of_params::Int8 = length(collect(eachmatch(URIUtils.ROUTE_REGEX_BASIC_STRUCTURE, path)))
            if route.has_regex && request_path_number_of_segments==route_number_of_params
                try
                    m::Union{RegexMatch, Nothing} = match(Regex(path), route_path)
                    if isa(m, RegexMatch) && length(m.match)==length(route_path)
                        found_route = route
                    end
                catch e
                    nothing
                end
            end
        end
    end

    if isnothing(found_route)
        Logger.log("Getting Route: route with path `$route_path` not defined")
    end

    return found_route
end

"""
    output_file_as_string(file_path::String)

Read file_path as string and obtain corresponding mime-type based on file extension

First try searching in STATIC_DIR then try in project root, each time relative to route_path
`favicon.ico` is special exception case
"""
function output_file_as_string(file_path::String) :: Tuple{String, Dict}
    filename::String = ""
    file_output_as_string::String = ""
    headers::Dict{String, String} = Dict("Content-Type" => "")

    file_type::String = split(file_path, '.')[end]
    if file_type in keys(mime_types)
        headers["Content-Type"] = mime_types[file_type]
    else
        Logger.log("HTTP mime-type not yet supported for file: `$file_path`\nPlease file Github issue at: https://github.com/DanceJL/Dance.jl/issues")
    end

    if file_path=="/favicon.ico"
        filename = Configuration.Settings[:html_favicon_name] * ".ico"
    else
        try
            filename = STATIC_DIR * split(file_path, STATIC_ROUTE_PREFIX)[2]
        catch e
            filename = strip(file_path, '/')
        end
    end
    file_output_as_string = read(filename, String)
    headers["Content-Length"] = string(filesize(filename))

    return file_output_as_string, headers
end

"""
    path_already_exists(path::String)

Verify if ROUTES ordered dict already has existing entry for specified path
"""
function path_already_exists(path::String) :: Bool
    path_exists = false

    for route_path in keys(ROUTES)
        if route_path == path
            path_exists = true
            break
        end
    end

    return path_exists
end

"""
    populate(file_path)

Register Routes from `Configuration.Settings[:routes_filename]` file
"""
function populate(file_path::String) :: Bool
    is_success::Bool = false
    routes_filepath::String = joinpath(file_path, Configuration.Settings[:routes_filename]*".jl")

    if isfile(routes_filepath)
        is_success = true
    else
        @error "Populating Routes: file not found at $routes_filepath"
    end

    # Add favicon.ico
    route("/favicon.ico", output_file_as_string; method=GET, endpoint=EP_STATIC)

    return is_success
end

"""
    remove_trailing_slash(path::String)

Remove eventual trailing slash, if not index url
"""
function remove_trailing_slash(path::String) :: String
    if length(path)>1 && path[end]=='/'
        path = chop(path)
    end

    return path
end

"""
    route(path::Union{Regex, String}, action::Function; method::String=POST, endpoint=EP_JSON, html_file::String=Configuration.Settings[:html_base_filename]*".html")

Create new Route and add to ROUTES ordered dict
- If path is Regex, convert to String for storage (Route has `has_regex` field)
- Remove trailing slash
"""
function route(path::Union{Regex, String}, action::Function; method::String=POST, endpoint::String=EP_JSON, html_file::String=Configuration.Settings[:html_base_filename]*".html") :: OrderedCollections.OrderedDict{String, Route}
    has_regex::Bool = isa(path, Regex)

    if has_regex
        path::String = replace(rstrip(string(path), '"'), "r\""=>"")
    end
    path = remove_trailing_slash(path)

    # Prevent route paths contaning dots at start (static files relative path)
    if startswith(path, "/..")
        @error "Route paths cannot go higher than project route in directory structure"
    end

    # Route path must have leading slash
    if !startswith(path, "/")
        @error "Route must start with leading slash"
    end

    push!(ROUTES, path => Route(endpoint=endpoint, method=method, path=path, has_regex=has_regex, action=action, html_file=html_file))
end

"""
    route_group(routes::Array; route_prefix::String="", method::String="", endpoint::String="", html_file::String=Configuration.Settings[:html_base_filename]*".html")

Loop through array of named tuple routes, calling `route` function
- automatically prepend route_prefix to route path (if supplied)
- if other parameters supplied as function kwargs, set route parameters accordingly
"""
function route_group(routes::Array; route_prefix::String="", method::String="", endpoint::String="", html_file::String=Configuration.Settings[:html_base_filename]*".html") :: Union{Nothing, OrderedCollections.OrderedDict{String, Route}}
    for item in routes
        if isa(item, NamedTuple)
            path::Union{Regex, String} = ""

            for idx in [:path, :action]
                if !(idx in keys(item))
                    @error "Please make sure $item contains at least `path` and `action` keys"
                end
            end

            if length(route_prefix)>0 && route_prefix!="/"
                if isa(item.path, Regex)
                    path = Regex(route_prefix * string(item.path)[3:end-1])
                else
                    path = route_prefix * item.path
                end
            else
                path = item.path
                route_prefix = ""
            end

            if :method in keys(item)
                method = item.method
            end

            if :endpoint in keys(item)
                endpoint = item.endpoint
            end

            if :html_file in keys(item)
                html_file = item.html_file
            end

            route(path, item.action; method=method, endpoint=endpoint, html_file=html_file)
        else
            @error "Please supply route $item as NamedTuple"
        end
    end
end

"""
    static_dir(route_prefix::String, dir_path::String)

Parse supplied directory path and create routes for each item
"""
function static_dir(route_prefix::String, dir_path::String) :: Nothing
    global STATIC_DIR = dir_path
    global STATIC_ROUTE_PREFIX = route_prefix

    for (root, dirs, files) in walkdir(dir_path)
        for file in files
            if !occursin(".DS_Store", file)
                path::String = STATIC_ROUTE_PREFIX * split(joinpath(root, file), STATIC_DIR)[2]

                # Windows path issue (replace backslash by forward slash)
                if Sys.iswindows()
                    path = replace(path, "\\" => "/")
                end

                route(path, output_file_as_string; method=GET, endpoint=EP_STATIC)
            end
        end
    end
end

end


text:
using Adapt
using Test

# custom array type

struct CustomArray{T,N} <: AbstractArray{T,N}
    arr::Array
end

CustomArray(x::Array{T,N}) where {T,N} = CustomArray{T,N}(x)
Adapt.adapt_storage(::Type{<:CustomArray}, xs::Array) = CustomArray(xs)

Base.size(x::CustomArray, y...) = size(x.arr, y...)
Base.getindex(x::CustomArray, y...) = getindex(x.arr, y...)
Base.count(x::CustomArray) = count(x.arr)

const mat = CustomArray{Float64,2}(rand(2,2))
const vec = CustomArray{Float64,1}(rand(2))

const mat_bools = CustomArray{Bool,2}(rand(Bool,2,2))

macro test_adapt(to, src_expr, dst_expr, typ=nothing)
    quote
        src = $(esc(src_expr))
        dst = $(esc(dst_expr))

        @test adapt($(esc(to)), src) == dst
        @test typeof(adapt($(esc(to)), src)) == typeof(dst)
        if $(esc(typ)) !== nothing
            @test typeof(dst) <: $(esc(typ))
        end
    end
end

AnyCustomArray{T,N} = Union{CustomArray,WrappedArray{T,N,CustomArray,CustomArray{T,N}}}

# basic adaption
@test_adapt CustomArray mat.arr mat

# idempotency
@test_adapt CustomArray mat mat

# custom wrapper
struct Wrapper{T}
    arr::T
end
Wrapper(x::T) where T = Wrapper{T}(x)
Adapt.adapt_structure(to, xs::Wrapper) = Wrapper(adapt(to, xs.arr))
@test_adapt CustomArray Wrapper(mat.arr) Wrapper(mat)

@testset "base structures" begin

@test @inferred(adapt(nothing, NamedTuple())) == NamedTuple()
@test_adapt CustomArray (mat.arr,) (mat,)
@test @allocated(adapt(nothing, ())) == 0
@test @allocated(adapt(nothing, (1,))) == 0
@test @allocated(adapt(nothing, (1,2,3,4,5,6,7,8,9,10))) == 0

@test_adapt CustomArray (a=mat.arr,) (a=mat,)

end

# NOTE: if we put this in the preceding testset, unrelated tests start to allocate
@testset "closures" begin

# basic capture of a variable
function closure1(x)
    function foobar()
        x
    end
    foobar
end

f = closure1(mat)
@test f() == mat
f′ = adapt(CustomArray, f)
@test f′() == mat.arr
@test @inferred(adapt(nothing, f)()) == f()

# closure with sparams
function closure2(A::CustomArray{T}=zeros(1), b::Number=0) where {T}
    function f()
        convert(T, 0)
        A, A[1] == b
    end
    f
end

f = closure2(mat)
@test f() == (mat, false)
f′ = adapt(CustomArray, f)
@test f′() == (mat.arr, false)

end

@testset "array wrappers" begin

@test_adapt CustomArray view(mat.arr,:,:) view(mat,:,:) AnyCustomArray
inds = CustomArray{Int,1}([1,2])
@test_adapt CustomArray view(mat.arr,inds.arr,:) view(mat,inds,:) AnyCustomArray

# NOTE: manual creation of PermutedDimsArray because permutedims collects
@test_adapt CustomArray PermutedDimsArray(mat.arr,(2,1)) PermutedDimsArray(mat,(2,1)) AnyCustomArray

# NOTE: manual creation of ReshapedArray because Base.Array has an optimized `reshape`
@test_adapt CustomArray Base.ReshapedArray(mat.arr,(2,2),()) reshape(mat,(2,2)) AnyCustomArray

@test_adapt CustomArray Base.LogicalIndex(mat_bools.arr) Base.LogicalIndex(mat_bools) AnyCustomArray

@test_adapt CustomArray reinterpret(Int64,mat.arr) reinterpret(Int64,mat) AnyCustomArray

## doubly-wrapped

@test_adapt CustomArray reinterpret(Int64,view(mat.arr,:,:)) reinterpret(Int64,view(mat,:,:)) AnyCustomArray

@test_adapt CustomArray reshape(view(mat.arr,:,:), (2,2)) reshape(view(mat,:,:), (2,2)) AnyCustomArray
@test_adapt CustomArray reshape(reinterpret(Int64,mat.arr), (2,2)) reshape(reinterpret(Int64,mat), (2,2)) AnyCustomArray
@test_adapt CustomArray reshape(reinterpret(Int64,view(mat.arr,:,:)), (2,2)) reshape(reinterpret(Int64,view(mat,:,:)), (2,2)) AnyCustomArray

@test_adapt CustomArray view(reinterpret(Int64,mat.arr), :, :) view(reinterpret(Int64,mat), :, :) AnyCustomArray
@test_adapt CustomArray view(reinterpret(Int64,view(mat.arr,:,:)), :, :) view(reinterpret(Int64,view(mat,:,:)), :, :) AnyCustomArray
@test_adapt CustomArray view(Base.ReshapedArray(mat.arr,(2,2),()), :, :) view(reshape(mat, (2,2)), :, :) AnyCustomArray
@test_adapt CustomArray view(reshape(view(mat.arr,:,:), (2,2)), :, :) view(reshape(view(mat,:,:), (2,2)), :, :) AnyCustomArray
@test_adapt CustomArray view(reshape(reinterpret(Int64,mat.arr), (2,2)), :, :) view(reshape(reinterpret(Int64,mat), (2,2)), :, :) AnyCustomArray
@test_adapt CustomArray view(reshape(reinterpret(Int64,view(mat.arr,:,:)), (2,2)), :, :) view(reshape(reinterpret(Int64,view(mat,:,:)), (2,2)), :, :) AnyCustomArray

using LinearAlgebra

@test_adapt CustomArray mat.arr' mat' AnyCustomArray

@test_adapt CustomArray transpose(mat.arr) transpose(mat) AnyCustomArray

@test_adapt CustomArray LowerTriangular(mat.arr) LowerTriangular(mat) AnyCustomArray
@test_adapt CustomArray UnitLowerTriangular(mat.arr) UnitLowerTriangular(mat) AnyCustomArray
@test_adapt CustomArray UpperTriangular(mat.arr) UpperTriangular(mat) AnyCustomArray
@test_adapt CustomArray UnitUpperTriangular(mat.arr) UnitUpperTriangular(mat) AnyCustomArray

@test_adapt CustomArray Diagonal(vec.arr) Diagonal(vec) AnyCustomArray

dl = CustomArray{Float64,1}(rand(2))
du = CustomArray{Float64,1}(rand(2))
d = CustomArray{Float64,1}(rand(3))
@test_adapt CustomArray Tridiagonal(dl.arr, d.arr, du.arr) Tridiagonal(dl, d, du) AnyCustomArray

end

@testset "type information" begin
    @test Adapt.ndims(LinearAlgebra.Transpose{Float64,Array{Float64,1}}) == 2
    @test Adapt.ndims(Adapt.WrappedSubArray{Float64,3,Array{Float64,3}}) == 3

    @test Adapt.parent(LinearAlgebra.Transpose{Float64,Array{Float64,1}}) == Array
    @test Adapt.parent(Adapt.WrappedSubArray{Float64,3,Array{Float64,3}}) == Array
end


text:
import Revise
import GLRenderer as GL
import ThreeDP3 as T
import Images as I
import MiniGSG as S
import PoseComposition: Pose
using Gen
import MeshCatViz

import MeshCatViz
MeshCatViz.setup_visualizer()

import PyCall
PyCall.py"""

import os
import numpy as np
import h5py as h5
from imageio import imread
import IPython
import math

def im2col(im, psize):
    n_channels = 1 if len(im.shape) == 2 else im.shape[0]
    (n_channels, rows, cols) = (1,) * (3 - len(im.shape)) + im.shape

    im_pad = np.zeros((n_channels,
                       int(math.ceil(1.0 * rows / psize) * psize),
                       int(math.ceil(1.0 * cols / psize) * psize)))
    im_pad[:, 0:rows, 0:cols] = im

    final = np.zeros((im_pad.shape[1], im_pad.shape[2], n_channels,
                      psize, psize))
    for c in range(n_channels):
        for x in range(psize):
            for y in range(psize):
                im_shift = np.vstack(
                    (im_pad[c, x:], im_pad[c, :x]))
                im_shift = np.column_stack(
                    (im_shift[:, y:], im_shift[:, :y]))
                final[x::psize, y::psize, c] = np.swapaxes(
                    im_shift.reshape(im_pad.shape[1] / psize, psize,
                                     im_pad.shape[2] / psize, psize), 1, 2)

    return np.squeeze(final[0:rows - psize + 1, 0:cols - psize + 1])

def filterDiscontinuities(depthMap):
    filt_size = 7
    thresh = 1000

    # Ensure that filter sizes are okay
    assert filt_size % 2 == 1, "Can only use odd filter sizes."

    # Compute discontinuities
    offset = (filt_size - 1) / 2
    patches = 1.0 * im2col(depthMap, filt_size)
    mids = patches[:, :, offset, offset]
    mins = np.min(patches, axis=(2, 3))
    maxes = np.max(patches, axis=(2, 3))

    discont = np.maximum(np.abs(mins - mids),
                         np.abs(maxes - mids))
    mark = discont > thresh

    # Account for offsets
    final_mark = np.zeros((480, 640), dtype=np.uint16)
    final_mark[offset:offset + mark.shape[0],
               offset:offset + mark.shape[1]] = mark

    return depthMap * (1 - final_mark)

def registerDepthMap(unregisteredDepthMap,
                     rgbImage,
                     depthK,
                     rgbK,
                     H_RGBFromDepth):

    unregisteredHeight = unregisteredDepthMap.shape[0]
    unregisteredWidth = unregisteredDepthMap.shape[1]

    registeredHeight = rgbImage.shape[0]
    registeredWidth = rgbImage.shape[1]

    registeredDepthMap = np.zeros((registeredHeight, registeredWidth))

    xyzDepth = np.empty((4,1))
    xyzRGB = np.empty((4,1))

    # Ensure that the last value is 1 (homogeneous coordinates)
    xyzDepth[3] = 1

    invDepthFx = 1.0 / depthK[0,0]
    invDepthFy = 1.0 / depthK[1,1]
    depthCx = depthK[0,2]
    depthCy = depthK[1,2]

    rgbFx = rgbK[0,0]
    rgbFy = rgbK[1,1]
    rgbCx = rgbK[0,2]
    rgbCy = rgbK[1,2]

    undistorted = np.empty(2)
    for v in range(unregisteredHeight):
      for u in range(unregisteredWidth):

            depth = unregisteredDepthMap[v,u]
            if depth == 0:
                continue

            xyzDepth[0] = ((u - depthCx) * depth) * invDepthFx
            xyzDepth[1] = ((v - depthCy) * depth) * invDepthFy
            xyzDepth[2] = depth

            xyzRGB[0] = (H_RGBFromDepth[0,0] * xyzDepth[0] +
                         H_RGBFromDepth[0,1] * xyzDepth[1] +
                         H_RGBFromDepth[0,2] * xyzDepth[2] +
                         H_RGBFromDepth[0,3])
            xyzRGB[1] = (H_RGBFromDepth[1,0] * xyzDepth[0] +
                         H_RGBFromDepth[1,1] * xyzDepth[1] +
                         H_RGBFromDepth[1,2] * xyzDepth[2] +
                         H_RGBFromDepth[1,3])
            xyzRGB[2] = (H_RGBFromDepth[2,0] * xyzDepth[0] +
                         H_RGBFromDepth[2,1] * xyzDepth[1] +
                         H_RGBFromDepth[2,2] * xyzDepth[2] +
                         H_RGBFromDepth[2,3])

            invRGB_Z  = 1.0 / xyzRGB[2]
            undistorted[0] = (rgbFx * xyzRGB[0]) * invRGB_Z + rgbCx
            undistorted[1] = (rgbFy * xyzRGB[1]) * invRGB_Z + rgbCy

            uRGB = int(undistorted[0] + 0.5)
            vRGB = int(undistorted[1] + 0.5)

            if (uRGB < 0 or uRGB >= registeredWidth) or (vRGB < 0 or vRGB >= registeredHeight):
                continue

            registeredDepth = xyzRGB[2]
            if registeredDepth > registeredDepthMap[vRGB,uRGB]:
                registeredDepthMap[vRGB,uRGB] = registeredDepth

    return registeredDepthMap

def registeredDepthMapToPointCloud(depthMap, rgbImage, rgbK, organized=True):
    rgbCx = rgbK[0,2]
    rgbCy = rgbK[1,2]
    invRGBFx = 1.0/rgbK[0,0]
    invRGBFy = 1.0/rgbK[1,1]

    height = depthMap.shape[0]
    width = depthMap.shape[1]

    if organized:
      cloud = np.empty((height, width, 6), dtype=np.float)
    else:
      cloud = np.empty((1, height*width, 6), dtype=np.float)

    goodPointsCount = 0
    for v in range(height):
        for u in range(width):

            depth = depthMap[v,u]

            if organized:
              row = v
              col = u
            else:
              row = 0
              col = goodPointsCount

            if depth <= 0:
                if organized:
                    if depth <= 0:
                       cloud[row,col,0] = float('nan')
                       cloud[row,col,1] = float('nan')
                       cloud[row,col,2] = float('nan')
                       cloud[row,col,3] = 0
                       cloud[row,col,4] = 0
                       cloud[row,col,5] = 0
                continue

            cloud[row,col,0] = (u - rgbCx) * depth * invRGBFx
            cloud[row,col,1] = (v - rgbCy) * depth * invRGBFy
            cloud[row,col,2] = depth
            cloud[row,col,3] = rgbImage[v,u,0]
            cloud[row,col,4] = rgbImage[v,u,1]
            cloud[row,col,5] = rgbImage[v,u,2]
            if not organized:
              goodPointsCount += 1

    if not organized:
      cloud = cloud[:,:goodPointsCount,:]

    return cloud

def writePLY(filename, cloud, faces=[]):
    if len(cloud.shape) != 3:
        print("Expected pointCloud to have 3 dimensions. Got %d instead" % len(cloud.shape))
        return

    color = True if cloud.shape[2] == 6 else False
    num_points = cloud.shape[0]*cloud.shape[1]

    header_lines = [
        'ply',
        'format ascii 1.0',
        'element vertex %d' % num_points,
        'property float x',
        'property float y',
        'property float z',
        ]
    if color:
        header_lines.extend([
        'property uchar diffuse_red',
        'property uchar diffuse_green',
        'property uchar diffuse_blue',
        ])
    if faces != None:
        header_lines.extend([
        'element face %d' % len(faces),
        'property list uchar int vertex_indices'
        ])

    header_lines.extend([
      'end_header',
      ])

    f = open(filename, 'w+')
    f.write('\n'.join(header_lines))
    f.write('\n')

    lines = []
    for i in range(cloud.shape[0]):
        for j in range(cloud.shape[1]):
            if color:
                lines.append('%s %s %s %d %d %d' % tuple(cloud[i, j, :].tolist()))
            else:
                lines.append('%s %s %s' % tuple(cloud[i, j, :].tolist()))

    for face in faces:
        lines.append(('%d' + ' %d'*len(face)) % tuple([len(face)] + list(face)))

    f.write('\n'.join(lines) + '\n')
    f.close()

def writePCD(pointCloud, filename, ascii=False):
    if len(pointCloud.shape) != 3:
        print("Expected pointCloud to have 3 dimensions. Got %d instead" % len(pointCloud.shape))
        return
    with open(filename, 'w') as f:
        height = pointCloud.shape[0]
        width = pointCloud.shape[1]
        f.write("# .PCD v.7 - Point Cloud Data file format\n")
        f.write("VERSION .7\n")
        if pointCloud.shape[2] == 3:
            f.write("FIELDS x y z\n")
            f.write("SIZE 4 4 4\n")
            f.write("TYPE F F F\n")
            f.write("COUNT 1 1 1\n")
        else:
            f.write("FIELDS x y z rgb\n")
            f.write("SIZE 4 4 4 4\n")
            f.write("TYPE F F F F\n")
            f.write("COUNT 1 1 1 1\n")
        f.write("WIDTH %d\n" % width)
        f.write("HEIGHT %d\n" % height)
        f.write("VIEWPOINT 0 0 0 1 0 0 0\n")
        f.write("POINTS %d\n" % (height * width))
        if ascii:
          f.write("DATA ascii\n")
          for row in range(height):
            for col in range(width):
                if pointCloud.shape[2] == 3:
                    f.write("%f %f %f\n" % tuple(pointCloud[row, col, :]))
                else:
                    f.write("%f %f %f" % tuple(pointCloud[row, col, :3]))
                    r = int(pointCloud[row, col, 3])
                    g = int(pointCloud[row, col, 4])
                    b = int(pointCloud[row, col, 5])
                    rgb_int = (r << 16) | (g << 8) | b
                    packed = pack('i', rgb_int)
                    rgb = unpack('f', packed)[0]
                    f.write(" %.12e\n" % rgb)
        else:
          f.write("DATA binary\n")
          if pointCloud.shape[2] == 6:
              # These are written as bgr because rgb is interpreted as a single
              # little-endian float.
              dt = np.dtype([('x', np.float32),
                             ('y', np.float32),
                             ('z', np.float32),
                             ('b', np.uint8),
                             ('g', np.uint8),
                             ('r', np.uint8),
                             ('I', np.uint8)])
              pointCloud_tmp = np.zeros((height*width, 1), dtype=dt)
              for i, k in enumerate(['x', 'y', 'z', 'r', 'g', 'b']):
                  pointCloud_tmp[k] = pointCloud[:, :, i].reshape((height*width, 1))
              pointCloud_tmp.tofile(f)
          else:
              dt = np.dtype([('x', np.float32),
                             ('y', np.float32),
                             ('z', np.float32),
                             ('I', np.uint8)])
              pointCloud_tmp = np.zeros((height*width, 1), dtype=dt)
              for i, k in enumerate(['x', 'y', 'z']):
                  pointCloud_tmp[k] = pointCloud[:, :, i].reshape((height*width, 1))
              pointCloud_tmp.tofile(f)

def getRGBFromDepthTransform(calibration, camera, referenceCamera):
    irKey = "H_{0}_ir_from_{1}".format(camera, referenceCamera)
    rgbKey = "H_{0}_from_{1}".format(camera, referenceCamera)

    rgbFromRef = calibration[rgbKey][:]
    irFromRef = calibration[irKey][:]

    return np.dot(rgbFromRef, np.linalg.inv(irFromRef))

def get_data(ycb_data_folder, target_object, viewpoint_camera, viewpoint_angle):
    referenceCamera = "NP5"
    basename = "{0}_{1}".format(viewpoint_camera, viewpoint_angle)
    depthFilename = os.path.join(ycb_data_folder+ target_object, basename + ".h5")
    rgbFilename = os.path.join(ycb_data_folder + target_object, basename + ".jpg")

    calibrationFilename = os.path.join(ycb_data_folder+target_object, "calibration.h5")
    calibration = h5.File(calibrationFilename)

    rgbImage = imread(rgbFilename)
    depthK = calibration["{0}_ir_K".format(viewpoint_camera)][:]
    rgbK = calibration["{0}_rgb_K".format(viewpoint_camera)][:]
    depthScale = np.array(calibration["{0}_ir_depth_scale".format(viewpoint_camera)]) * .0001 # 100um to meters
    H_RGBFromDepth = getRGBFromDepthTransform(calibration, viewpoint_camera, referenceCamera)

    unregisteredDepthMap = h5.File(depthFilename)["depth"][:]
    unregisteredDepthMap = filterDiscontinuities(unregisteredDepthMap) * depthScale

    registeredDepthMap = registerDepthMap(unregisteredDepthMap,
                                          rgbImage,
                                          depthK,
                                          rgbK,
                                          H_RGBFromDepth)

    pointCloud = registeredDepthMapToPointCloud(registeredDepthMap, rgbImage, rgbK)
    return pointCloud
"""

PyCall.py"get_data"("/home/nishadg/Downloads/003_cracker_box_berkeley_rgbd/","003_cracker_box","NP3",15)

# +
import HDF5
filename = "/home/nishadg/Downloads/003_cracker_box_berkeley_rgbd/003_cracker_box/NP2_0.h5"
d = HDF5.h5open(filename, "r+")
depth = Float64.(HDF5.read(d["depth"]))
filename = "/home/nishadg/Downloads/003_cracker_box_berkeley_rgbd/003_cracker_box/calibration.h5"
c = HDF5.h5open(filename, "r+");
K = HDF5.read(c["NP2_ir_K"])
depthScale = HDF5.read(c["NP2_ir_depth_scale"]) * .0001 # 100um to meters

depth .*= depthScale;
depth = collect(transpose(depth))
camera = GL.CameraIntrinsics(size(depth)[2], size(depth)[1], K[1,1], K[2,2],K[3,1],K[3,2], 0.01, 100.0)
cloud = GL.flatten_point_cloud(GL.depth_image_to_point_cloud(depth, camera))
MeshCatViz.viz(cloud)

# -
import Plots
Plots.heatmap(depth)

import Pkg; Pkg.add("Plots")

YCB_DIR = "/home/nishadg/mcs/ThreeDVision.jl/data/ycbv2"
world_scaling_factor = 100.0
id_to_cloud, id_to_shift, id_to_box  = T.load_ycbv_models_adjusted(YCB_DIR, world_scaling_factor);
all_ids = sort(collect(keys(id_to_cloud)));

IDX =100
@show T.get_ycb_scene_frame_id_from_idx(YCB_DIR,IDX)
gt_poses, ids, rgb_image, gt_depth_image, cam_pose, original_camera = T.load_ycbv_scene_adjusted(
    YCB_DIR, IDX, world_scaling_factor, id_to_shift
);
camera = T.scale_down_camera(original_camera, 4)
@show camera
img = I.colorview(I.Gray, gt_depth_image ./ maximum(gt_depth_image))
I.colorview(I.RGB, permutedims(Float64.(rgb_image)./255.0, (3,1,2)))

# +
resolution = 0.5

renderer = GL.setup_renderer(camera, GL.DepthMode())
for id in all_ids
    cloud = id_to_cloud[id]
    v,n,f = GL.mesh_from_voxelized_cloud(GL.voxelize(cloud, resolution), resolution)
    GL.load_object!(renderer, v, f)
end
depth_image = GL.gl_render(
    renderer, ids, gt_poses, T.IDENTITY_POSE)
depth_image[depth_image .== 50000.0] .= 200.0

function get_cloud_from_ids_and_poses(ids, poses, camera_pose)
    depth_image = GL.gl_render(renderer, ids, poses, camera_pose)
    cloud = GL.flatten_point_cloud(GL.depth_image_to_point_cloud(depth_image, camera))
    cloud = cloud[:,cloud[3,:] .< (camera.far - 1.0)]
    if size(cloud)[2] == 0
        cloud = zeros(3,1)
    end
    cloud
end

img = I.colorview(I.Gray, depth_image ./ maximum(depth_image))

# +

hypers = T.Hyperparams(slack_dir_conc=300.0, slack_offset_var=0.5, p_outlier=0.01,
    noise=0.2, resolution=0.5, parent_face_mixture_prob=0.99, floating_position_bounds=(-1000.0, 1000.0, -1000.0,1000.0,-1000.0,1000.0))

params = T.SceneModelParameters(
    boxes=vcat([id_to_box[id] for id in ids],[S.Box(40.0,40.0,0.1)]),
    get_cloud_from_poses_and_idx=
    (poses, idx, p) -> get_cloud_from_ids_and_poses(ids, poses[1:end-1], cam_pose),
    camera_pose=cam_pose,
    hyperparams=hypers, N=1)

table_pose = Pose(0.0, 0.0, -0.1)
num_obj = length(params.boxes)
g = T.graph_with_edges(num_obj, [])
constraints = Gen.choicemap((T.structure_addr()) => g, T.floating_pose_addr(num_obj) => table_pose)
for i=1:num_obj-1
   constraints[T.floating_pose_addr(i)] = T.get_c_relative_to_a(cam_pose, gt_poses[i])
end

obs_cloud = GL.flatten_point_cloud(GL.depth_image_to_point_cloud(gt_depth_image, original_camera))
obs_cloud = GL.move_points_to_frame_b(obs_cloud, cam_pose)
obs_cloud =obs_cloud[:, 1.2 .< obs_cloud[3, :] .< 40.0];
obs_cloud =obs_cloud[:, -25.0 .< obs_cloud[1, :] .< 25.0];
obs_cloud =obs_cloud[:, -25.0 .< obs_cloud[2, :] .< 25.0];
constraints[T.obs_addr()] = T.voxelize(obs_cloud, params.hyperparams.resolution)
MeshCatViz.reset_visualizer()
MeshCatViz.viz(constraints[T.obs_addr()])

trace, _ = generate(T.scene, (params,), constraints);
# -

MeshCatViz.reset_visualizer()
MeshCatViz.viz(T.get_obs_cloud(trace); color=I.colorant"black", channel_name=:obs_cloud)
MeshCatViz.viz(T.get_gen_cloud(trace); color=I.colorant"red", channel_name=:gen_cloud)

traces = [trace];

for i in 1:length(ids)
    if !T.isFloating(T.get_structure(trace),i)
        continue
    end
    pose_addr = T.floating_pose_addr(i)
    for _ in 1:50
        trace, acc = T.drift_move(trace, pose_addr, 0.5, 10.0)
        if acc traces = push!(traces, trace) end

        trace, acc = T.drift_move(trace, pose_addr, 1.0, 100.0)
        if acc traces = push!(traces, trace) end
        trace, acc = T.drift_move(trace, pose_addr, 1.5, 1000.0)
        if acc traces = push!(traces, trace) end
        trace, acc = T.drift_move(trace, pose_addr, 0.1, 1000.0)
        if acc traces = push!(traces, trace) end
        trace, acc = T.drift_move(trace, pose_addr, 0.1, 100.0)
        if acc traces = push!(traces, trace) end
        trace, acc = T.drift_move(trace, pose_addr, 0.5, 5000.0)
        if acc traces = push!(traces, trace) end
        trace, acc = T.pose_flip_move(trace, pose_addr, 1, 1000.0)
        if acc traces = push!(traces, trace) end
        trace, acc = T.pose_flip_move(trace, pose_addr, 2, 1000.0)
        if acc traces = push!(traces, trace) end
        trace, acc = T.pose_flip_move(trace, pose_addr, 3, 1000.0)
        if acc traces = push!(traces, trace) end
    end
end

CRAZY_POSE = Pose(-100.0 * ones(3), T.IDENTITY_ORIENTATION)
for i in 1:length(ids)
    if !T.isFloating(T.get_structure(trace),i)
        continue
    end
    addr = T.floating_pose_addr(i)

    t, = Gen.update(trace, Gen.choicemap(T.floating_pose_addr(i) => CRAZY_POSE))
    valid = T.get_unexplained_obs_cloud(t)

    p = T.icp_object_pose(trace[addr], trace[:obs][:,valid],
                          p-> get_cloud_from_ids_and_poses(
                            [ids[i]], [p], cam_pose), false)
    for _ in 1:3
        trace, acc = T.pose_mixture_move(trace, addr, [trace[addr], p], [0.5, 0.5], 0.001, 5000.0)
        if acc traces = push!(traces, trace) end
    end
end

MeshCatViz.reset_visualizer()
MeshCatViz.viz(T.get_obs_cloud(trace); color=I.colorant"black", channel_name=:obs_cloud)
MeshCatViz.viz(T.get_gen_cloud(trace); color=I.colorant"red", channel_name=:gen_cloud)

depth_image = GL.gl_render(
    renderer, ids, T.get_poses(trace)[1:end-1], cam_pose)
depth_image[depth_image .== 50000.0] .= 200.0
img = I.colorview(I.Gray, depth_image ./ maximum(depth_image))

# +
renderer_texture = GL.setup_renderer(original_camera, GL.TextureMode())
obj_paths = T.load_ycb_model_obj_file_paths(YCB_DIR)
texture_paths = T.load_ycb_model_texture_file_paths(YCB_DIR)
for id in all_ids
    v,n,f,t = renderer_texture.gl_instance.load_obj_parameters(
        obj_paths[id]
    )
    v = v * world_scaling_factor
    v .-= id_to_shift[id]'

    GL.load_object!(renderer_texture, v, n, f, t,
        texture_paths[id]
    )
end

# -

rgb_image, depth_image = GL.gl_render(
    renderer_texture, ids, T.get_poses(trace)[1:end-1], cam_pose)
I.colorview(I.RGBA, permutedims(rgb_image,(3,1,2)))


text:
# print all significant SNPs from logistic GWAS into latex table format
using CSV
df = CSV.read("ukb.final.logistic.pval.txt")
function run()
    tot = 0
    for i in 1:size(df, 1)
        if df[i, 6] < 5e-8
            tot += 1
            # snpid, chrom, pos, p-value
            # println(df[i, 3], " & ", df[i, 1], " & ", df[i, 2], " & ",
            #     round(df[i, 6], sigdigits=3), "\\\\")
        end
    end
    tot
end


text:
using Unitful, Statistics, StaticArrays, Interpolations, StructArrays

const Point = SVector{2, Float64}
point(x::Missing) = missing
point(x::AbstractVector{Float64}) = Point(x[1], x[2])
points(x::Matrix{Float64}) = point.(eachrow(x))
point(x::Instantaneous{Matrix{Float64}})= point(vec(x.data))
# point(x) = point(x[])
struct Track
    coords::Vector{Point}
    distance::Vector{Float64}
    direction::Vector{Float64}
    tp::Int
    Δt::Float64
end
function Track(x::Prolonged{Matrix{Float64}})
    xyt = !issorted(x.data[:,3]) ? sortslices(x.data, dims=1, lt=(x,y)->isless(x[3],y[3])) : x.data
    t = xyt[:,3]
    Δ = diff(t)
    δ = minimum(Δ)
    xy = points(xyt)
    n = length(xy)
    if !all(x -> isapprox(x, δ, atol = 1/30), Δ)
        itp = interpolate((t, ), xy, Gridded(Linear()))
        t_lin = range(t[1], t[end], length = n)
        δ = step(t_lin)
        xy .= itp.(t_lin)
    end
    _xy, v, a, m = smoothed(xy)
    xy .= Point.(eachcol(_xy))
    tp = findfirst(x -> x > 0.5, m)
    if tp ≡ nothing || tp ≥ n - 1
        tp = n
    end
    Track(xy, v, a, tp, δ)
end
struct TimedPoint
    xy::Point
    t::Float64
end
const PointCollection = StructVector{TimedPoint}
pointcollection(x::Missing, t₀) = StructVector{TimedPoint}(undef, 0)
pointcollection(x, t₀) = StructVector(TimedPoint(Point(i[1], i[2]), i[3] - t₀) for i in eachrow(x.data))
function getdata(x)
    feeder = point(x[:feeder])
    track = Track(x[:track])
    pellet = pointcollection(get(x, :pellet, missing), x[:track].data[1,3])
    feeder, track, pellet
end
abstract type DungMethod end
struct ClosedNest <: DungMethod
    feeder::Point
    nest::Point
    track::Track
    pellet::PointCollection
end
function ClosedNest(x)
    feeder, track, pellet = getdata(x.data)
    nest = point(x.data[:nest])
    ClosedNest(feeder, nest, track, pellet)
end
struct Transfer <: DungMethod
    feeder::Point
    track::Track
    pellet::PointCollection
    south::Point
    north::Point
    nest2feeder::Float64
    azimuth::Float64
end
function Transfer(x)
    feeder, track, pellet = getdata(x.data)
    south = point(x.data[:south])
    north = point(x.data[:north])
    d, u = split(string(x.metadata.setup[:nest2feeder]))
    nest2feeder = Float64(ustrip(uconvert(Unitful.cm, parse(Int, d)*getfield(Unitful, Symbol(u)))))
    d, u = split(string(x.metadata.setup[:azimuth]))
    azimuth = Float64(ustrip(uconvert(Unitful.rad, parse(Int, d)*getfield(Unitful, Symbol(u)))))
    Transfer(feeder, track, pellet, south, north, nest2feeder, azimuth)
end
struct TransferNest <: DungMethod
    feeder::Point
    track::Track
    pellet::PointCollection
    south::Point
    north::Point
    nest2feeder::Float64
    azimuth::Float64
    originalnest::Point
end
function TransferNest(x)
    y = Transfer(x)
    originalnest = point(x.data[:originalnest])
    TransferNest(y.feeder, y.track, y.pellet, y.south, y.north, y.nest2feeder, y.azimuth, originalnest)
end
Transfer(x::TransferNest) = Transfer(x.feeder, x.track, x.pellet, x.south, x.north, x.nest2feeder, x.azimuth)
struct TransferNestBelen <: DungMethod
    feeder::Point
    track::Track
    pellet::PointCollection
    southbefore::Point
    northbefore::Point
    feederbefore::Point
    nestbefore::Point
    south::Point
    north::Point
end
function TransferNestBelen(x)
    feeder, track, pellet = getdata(x.data)
    southbefore = point(x.data[:southbefore])
    northbefore = point(x.data[:northbefore])
    feederbefore = point(x.data[:feederbefore])
    nestbefore = point(x.data[:nestbefore])
    south = point(x.data[:south])
    north = point(x.data[:north])
    TransferNestBelen(feeder, track, pellet, rightdowninitial, leftdowninitial, rightupinitial, leftupinitial, rightdownfinal, leftdownfinal, rightupfinal, leftupfinal)
end
struct DawaySandpaper <: DungMethod
    feeder::Point
    nest::Point
    track::Track
    pellet::PointCollection
    rightdowninitial::Point
    leftdowninitial::Point
    rightupinitial::Point
    leftupinitial::Point
    rightdownfinal::Point
    leftdownfinal::Point
    rightupfinal::Point
    leftupfinal::Point
end
function DawaySandpaper(x)
    feeder, track, pellet = getdata(x.data)
    nest = point(x.data[:nest])
    rightdowninitial = point(x.data[:rightdowninitial])
    leftdowninitial = point(x.data[:leftdowninitial])
    rightupinitial = point(x.data[:rightupinitial])
    leftupinitial = point(x.data[:leftupinitial])
    rightdownfinal = point(x.data[:rightdownfinal])
    leftdownfinal = point(x.data[:leftdownfinal])
    rightupfinal = point(x.data[:rightupfinal])
    leftupfinal = point(x.data[:leftupfinal])
    DawaySandpaper(feeder, nest, track, pellet, rightdowninitial, leftdowninitial, rightupinitial, leftupinitial, rightdownfinal, leftdownfinal, rightupfinal, leftupfinal)
end
struct Daway <: DungMethod
    feeder::Point
    nest::Point
    track::Track
    pellet::PointCollection
    initialfeeder::Point
end
function Daway(x)
    feeder, track, pellet = getdata(x.data)
    nest = point(x.data[:nest])
    initialfeeder = point(x.data[:initialfeeder])
    Daway(feeder, nest, track, pellet, initialfeeder)
end

######################### DungMethod methods ###########

DungMethod(x, displace_direction::Missing, transfer::Missing, nest_coverage::Missing) = error("unidentified experimental setup")
DungMethod(x, displace_direction::Missing, transfer::Missing, nest_coverage) = ClosedNest(x)
function DungMethod(x, displace_direction::Missing, transfer, nest_coverage)
    if x.metadata.setup[:person] == "belen"
        TransferNestBelen(x)
    else
        if transfer == "back"
            TransferNest(x)
        else
            Transfer(x)
        end
    end
end
DungMethod(x, displace_direction, transfer::Missing, nest_coverage) = x.metadata.setup[:person] == "belen" ? DawaySandpaper(x) : Daway(x)

#=function DungMethod(displace_direction, person, transfer, nest_coverage)
if !ismissing(x.displace_direction)
if x.person == "belen"
DawaySandpaper(x)
else
Daway(x)
end
elseif !ismissing(x.transfer)
if x.person == "belen"
TransferNestBelen(x)
else
if x.transfer == "back"
TransferNest(x)
else
Transfer(x)
end
end
elseif x.nest_coverage == "closed"
ClosedNest(x)
else
error("unidentified experimental setup")
end
end=#

######################### Common methods ###########

mutable struct Common{N}
    feeder::Point
    nest::Point
    track::Track
    pellet::PointCollection
    originalnest::N
end
Common(x::ClosedNest) = Common(x.feeder, x.nest, x.track, x.pellet, x.nest)
nest(x::Common) = x.nest
turning(x::Common) = x.track.homing[end]
originalnest(x::Common) = x.originalnest

function Common(x::TransferNestBelen)
    v = x.northbefore - x.southbefore
    u = x.nestbefore - x.feederbefore
    azimuth = atan(v[2], v[1]) - atan(u[2], u[1])
    nest2feeder = norm(x.nestbefore - x.feederbefore)
    v = x.north - x.south
    α = atan(v[2], v[1]) + azimuth
    u = Point(cos(α), sin(α))
    nest = x.feeder + u*nest2feeder
    Common(x.feeder, nest, x.track, x.pellet, missing)
end

function Common(x::Transfer)
    v = x.north - x.south
    α = atan(v[2], v[1]) + x.azimuth - π
    u = Point(cos(α), sin(α))
    nest = x.feeder + u*x.nest2feeder
    Common(x.feeder, nest, x.track, x.pellet, missing)
end

function Common(x::TransferNest)
    y = Common(Transfer(x))
    Common(y.feeder, y.nest, y.track, y.pellet, x.originalnest)
end

function Common(x::DawaySandpaper)
    originalnest = x.nest
    initial = mean(getproperty(x, k) for k in [:rightdowninitial, :leftdowninitial, :rightupinitial, :leftupinitial])
    final = mean(getproperty(x, k) for k in [:rightdownfinal, :leftdownfinal, :rightupfinal, :leftupfinal])
    v = final - initial
    nest = originalnest + v
    _feeder = x.feeder
    feeder = _feeder + v
    Common(feeder, nest, x.track, x.pellet, originalnest)
end

function Common(x::Daway)
    originalnest = x.nest
    v = x.feeder - x.initialfeeder
    nest = originalnest + v
    Common(x.feeder, nest, x.track, x.pellet, originalnest)
end

######################### END ######################

Common(x) = Common(DungMethod(x, get(x.metadata.setup, :displace_direction, missing), get(x.metadata.setup, :transfer, missing), get(x.metadata.setup, :nest_coverage, missing)))

#=x = rand(100)
i = 1:5
z = StructArray((x_cm = view(x, i), y_cm = view(x, 2 .+ i), t = view(x, 3 .+ i)))
z.t .-= z.t[1]=#


text:
# POPular Utils KIt

include("ArrSlices.jl")
const ENDLISTTOKEN = "popuki.endlisttoken"

function AE_filter(line::Ref{UInt}, ns::Namespace, predicate::AetherlangObject, collection::AetherlangObject)::AetherlangObject
    arr = []
    for i in 1:length(collection.dataref)
        aetherlang_call(predicate.dataref, ns, AetherlangObject[AetherlangObject(collection.dataref[i])], line).type[] !== AETH_BUILTIN_TYPES.void && push!(arr, collection.dataref[i])
    end
    AetherlangObject(arr)
end
function AE_map(line::Ref{UInt}, ns::Namespace, f::AetherlangObject, collection::AetherlangObject)::AetherlangObject
    arr = []
    for i in 1:length(collection.dataref)
        push!(arr, aetherlang_call(f.dataref, ns, AetherlangObject[AetherlangObject(collection.dataref[i])], line).dataref)
    end
    AetherlangObject(arr)
end
function AE_reduce(line::Ref{UInt}, ns::Namespace, f::AetherlangObject, collection::AetherlangObject)::AetherlangObject
    x = collection.dataref[1]
    for i in 1:length(collection.dataref)
        x = aetherlang_call(f.dataref, ns, AetherlangObject[AetherlangObject(x), AetherlangObject(collection.dataref[i])], line).dataref
    end
    AetherlangObject(x)
end
function AE_startlist(line::Ref{UInt}, ns::Namespace, args...)::AetherlangObject
    if args[end].type[] !== ENDLISTTOKEN
        throw(AetherlangError("No matching \"]\" found in popuki.AE_startlist expression"))
    end
    AetherlangObject([args[i].dataref for i in 1:length(args)-1])
end
function AE_arrslice(line::Ref{UInt}, ns::Namespace, arr::AetherlangObject, ind1::AetherlangObject, ind2::AetherlangObject)::AetherlangObject
    # :end token works in ind2
    if ind2.dataref == "end"
        ind2 = AetherlangObject(length(arr.dataref))
    end
    if ind1.dataref == 1 && ind2.dataref == length(arr.dataref)
        return arr
    end
    AetherlangObject(ArrSlice(arr.dataref, ind1.dataref, ind2.dataref), Ref(AETH_BUILTIN_TYPES.collection))
end
function AE_arrslicerestore(line::Ref{UInt}, ns::Namespace, arrs::AetherlangObject)::AetherlangObject
    if typeof(arrs.dataref) == ArrSlice
        return AetherlangObject(arrs.dataref.arr)
    end
    arrs
end
function AE_arrslicebounds(line::Ref{UInt}, ns::Namespace, arrs::AetherlangObject)::AetherlangObject
    if typeof(arrs.dataref) != ArrSlice
        return AetherlangObject()
    end
    AetherlangObject([arrs.dataref.beg, arrs.dataref.fin])
end
function AE_eleq(line::Ref{UInt}, ns::Namespace, coll1::AetherlangObject, coll2::AetherlangObject)::AetherlangObject
    """can compare phantoms, arrays, and slices element-wise regardless of their actual type"""
    if length(coll1.dataref) != length(coll2.dataref)
        return AetherlangObject()
    end
    l::Int = length(coll1.dataref)
    for i in 1:l
        if coll1.dataref[i] != coll2.dataref[i]
            return AetherlangObject()
        end
    end
    AetherlangObject(1)
end

popuki_namespace_modify = Namespace(
    "]" => AetherlangObject{String}("end", Ref(ENDLISTTOKEN)),
    "[" => AetherlangObject(AE_startlist),
    "[]" => AetherlangObject(Bool[]),
    "filter" => AetherlangObject(AE_filter),
    "map" => AetherlangObject(AE_map),
    "reduce" => AetherlangObject(AE_reduce),
    "arrslice" => AetherlangObject(AE_arrslice),
    "arrslicerestore" => AetherlangObject(AE_arrslicerestore),
    "arrslicebounds" => AetherlangObject(AE_arrslicebounds),
    "elementsequal?" => AetherlangObject(AE_eleq)
)
popuki_namespace_modify["eleq?"] = popuki_namespace_modify["elementsequal?"]

max_stars_repo_name: Aetherlang/aetherlang

max_issues_repo_name: gojakuch/aetherlang

max_forks_repo_name: gojakuch/aetherlang


text:
pair_array=[[1,2],[3,4],[5,6]]
for pair in pair_array
    odd,even = pair
    println("odd = $odd, even = $even")
end

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
"""
    RangeQuery

Range query is a type of histogram query corresponding to an
intervals over the domain. i -> j (1 <= i <= j <= domain).
"""
struct RangeQuery <: Query
    interval::Interval
    domain::Int # the domain size
end

"""
    RangeQueriesRangeQueries

Range queries are type of histogram query corresponding to
intervals over the domain.
"""
struct RangeQueries <: Queries
    intervals::Array{Interval, 1}
    domain::Int # the domain size
end

"""
    SeriesRangeQueriesRangeQueries

Range queries are type of histogram query corresponding to all possible
intervals over the domain that start with the first element.
"""
struct SeriesRangeQueries <: Queries
    domain::Int # the domain size
end

"""
    get(interval::Interval, domain::Int64)

Return 0/1 indicator vector for the interval in the size of the domain.
"""
function get_query_vector_from_interval(interval::Interval, domain::Int64)
    query_vector = zeros(domain)
    start_index, end_index = interval
    query_vector[start_index : end_index] .= 1
    query_vector
end

"""
    get(query::RangeQuery)

Return 0/1 indicator vector of interval [j:k] of the query. (query.interval = (j,k))
"""
function get(query::RangeQuery)
    HistogramQuery(get_query_vector_from_interval(query.interval, query.domain))
end

"""
    get(queries::RangeQueries, i::Int)

Return 0/1 indicator vector of interval [j:k] of the i'th range query (queries[i] = (j,k)).
"""
function get(queries::RangeQueries, i::Int)
    HistogramQuery(get_query_vector_from_interval(queries.intervals[i], queries.domain))
end

"""
    get(queries::RangeQueries, i::Int)

Return 0/1 indicator vector of interval [1:i].
"""
function get(queries::SeriesRangeQueries, i::Int)
    HistogramQuery(vcat(ones(i), zeros(queries.domain-i)))#-1.0*ones(queries.domain-i)))
end

"""
   evaluate(query::RangeQuery, h::Histogram)

Evaluate the query on the given histogram.
"""
function evaluate(query::RangeQuery, h::Histogram)
    @assert query.domain == length(h.weights)
    answer = evaluate(get(query), h)
    answer
end

"""
   evaluate(queries::RangeQueries, h::Histogram)

Evaluate all the given range queries on the given histogram.
"""
function evaluate(queries::RangeQueries, h::Histogram)
    @assert queries.domain == length(h.weights)
    num_queries = length(queries.intervals)
    answers = Array{Float64}(undef, num_queries)
    @simd for j=1:num_queries
        @inbounds answers[j] = evaluate(get(queries, j), h)
    end
    answers
end

"""
   evaluate(queries::RangeQueries, h::Histogram)

Evaluate all range queries on the given histogram.
"""
function evaluate(queries::SeriesRangeQueries, h::Histogram)
    @assert queries.domain == length(h.weights)
    answers = zeros(Float64, queries.domain)
    running_sum = 0
    @simd for j=1:queries.domain
        @inbounds running_sum += h.weights[j]
        @inbounds answers[j] = running_sum
    end
    answers
end

function queriesMatrix(queries::SeriesRangeQueries)

    query_matrix = zeros(queries.domain, queries.domain)
    @simd for i=1:queries.domain
        @simd for j=1:i
            @inbounds query_matrix[i,j] = 1
        end
    end
    query_matrix
end


text:
using Test
using GoogleCloud
import GoogleCloud.api

@testset "test api functions" begin
    @test "/this/is/it" == GoogleCloud.api.path_replace(
                                "/{foo}/{bar}/{baz}",
                                ["this", "is", "it"])

    @test api.path_tokens("/{foo}/{bar}/x/{baz}") == ["{foo}", "{bar}", "{baz}"]
end


text:
"""
Created on Thu Apr 30 2020

@author: Yoann Pradat

    CentraleSupelec
    MICS laboratory
    9 rue Juliot Curie, Gif-Sur-Yvette, 91190 France

Rank selection method in NMF based on average stability of latent factor across perturbed versions of the input matrix.
matrix V.
Alexandrov et al, Deciphering Signatures of Mutational Processes Operative in Human Cancer. 2013, Cell Reports.
"""

"""
    _get_perturbed_matrix_multi(V, rng)

Draws a perturbed version of the input matrix V using multinomial distributions.
"""
function _get_perturbed_matrix_multi(V::Matrix{T}, rng::AbstractRNG) where T <: Real
    #### normalize each column to sum 1
    V_norm = V .* repeat(sum(V, dims=1) .^ -1, size(V, 1), 1)

    #### draw from multinomial
    #### WARNING: incorrect for input V that does not contain integer entries
    V_pert = zeros(Int64, size(V))

    for j = 1:size(V_pert, 2)
        dist = Multinomial(convert(Integer, sum(V[:,j])), V_norm[:, j])
        V_pert[:, j] = rand(rng, dist)
    end

    convert.(T, V_pert)
end

"""
    _get_perturbed_matrix_noise(V, dist, rng[, force_positivity])

Draws a perturbed version of the input matrix V adding random noise drawn from the specified distribution.
"""
function _get_perturbed_matrix_noise(V::Matrix{T}, dist::Distribution, rng::AbstractRNG; force_positivity::Bool=true) where T <: Real
    V_nois = rand(rng, dist, size(V))
    V_sign = sign.(rand(rng, size(V)...) .- 0.5)
    V_pert = V + V_sign .* V_nois

    force_positivity ? max.(V_pert, 0) : V_pert
end

"""
    _get_perturbed_matrix(V, n_iter, pert_meth, pert_dist, rng)

Draws `n_iter` perturbed versions of the input matrix V adding random noise or sampling from the specified distribution.
"""
function _get_perturbed_matrix(V::Matrix{T}, pert_meth::Symbol, pert_dist::Union{Nothing, Distribution}, rng::AbstractRNG) where T <: Real
    if pert_meth == :multinomial
        V_pert = _get_perturbed_matrix_multi(V, rng)
    elseif pert_meth == :noise
        V_pert = _get_perturbed_matrix_noise(V, pert_dist, rng)
    end
    V_pert
end


text:
# Load necessary packages
using OrdinaryDiffEq
using Trixi

# Create equations
advectionvelocity = (1.0, 1.0)
equations = LinearScalarAdvectionEquation2D(advectionvelocity)

# Create DGSEM solver for polynomial degree = 3
solver = DGSEM(3, flux_lax_friedrichs)

# Create a uniformely refined mesh with periodic boundaries
coordinates_min = (-1, -1) # minimum coordinates
coordinates_max = ( 1,  1) # maximum coordinates
mesh_static = TreeMesh(coordinates_min, coordinates_max,
                       initial_refinement_level=4,
                       n_cells_max=30_000)

# Create semidiscretization with all spatial discretization-related components
semi = SemidiscretizationHyperbolic(mesh_static, equations,
                                    initial_condition_convergence_test,
                                    solver)

# Create ODE problem from semidiscretization with time span from 0.0 to 1.0
ode = semidiscretize(semi, (0.0, 1.0))

# Evolve ODE problem in time using OrdinaryDiffEq's `solve` method
sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false), dt=2.5e-2);

nothing # To suppress output when include'ing this file in the REPL (otherwise `sol` is printed)