text:
struct KLocalLiouvillian{H,T,Vs,A<:AbsLinearSuperOperator{T},B,C} <: AbsLinearSuperOperator{T}
    hilb::H
    sites::Vs
    HnH_l::A
    HnH_r::B
    LLdag::C
end

function KLocalLiouvillian(HnH, Lops)
    hilb = basis(HnH)
    T = eltype(HnH)
    HnH_l = T(-1.0im) * KLocalOperatorTensor(HnH, KLocalIdentity(T, hilb))
    HnH_r = T(1.0im) * KLocalOperatorTensor(KLocalIdentity(T, hilb), HnH')

    LLdag_list = [KLocalOperatorTensor(L, conj(L)) for L=Lops]
    LLdag = isempty(LLdag_list) ? LocalOperator(SuperOpSpace(basis(HnH))) : sum(LLdag_list)

    return KLocalLiouvillian(SuperOpSpace(basis(HnH)),[], HnH_l, HnH_r, LLdag)
end

"""
    liouvillian(Ĥ, ops::Vector)

Constructs the LocalOperator representation of a liouvillian super-operator with
coherent part given by the hamilotnian `Ĥ`, and `ops` as the list of jump operators.
"""
function QuantumOpticsBase.liouvillian(H::AbsLinearOperator, Lops::AbstractVector)
    HnH = duplicate(H)
    for L=Lops
        HnH += -0.5im * L'*L
    end
    return KLocalLiouvillian(HnH, Lops)
end

function QuantumOpticsBase.liouvillian(Lops::AbstractVector)
    hilb = basis(first(Lops))
    liouvillian(LocalOperator(hilb), Lops)
end
sites(op::KLocalLiouvillian) = op.sites
QuantumOpticsBase.basis(op::KLocalLiouvillian) = op.hilb

conn_type(op::KLocalLiouvillian) = conn_type(typeof(op))
conn_type(::Type{KLocalLiouvillian{H,T,Vs,A,B,C}}) where {H,T,Vs,A,B,C}  =
    SuperOpConnection{conn_type(A), conn_type(B), conn_type(C)}

function accumulate_connections!(acc::AbstractAccumulator, op::KLocalLiouvillian, v)
    accumulate_connections!(acc, op.HnH_l, v)
    accumulate_connections!(acc, op.HnH_r, v)
    accumulate_connections!(acc, op.LLdag, v)

    return acc
end

function row_valdiff!(conn::SuperOpConnection, op::KLocalLiouvillian, v::ADoubleState; init=true)
    row_valdiff!(conn.op_conn_l_id, op.HnH_l, v; init=init)
    row_valdiff!(conn.op_conn_r_id, op.HnH_r, v; init=init)
    row_valdiff!(conn.op_conn_l_r,  op.LLdag, v; init=init)

    return conn
end

function map_connections(fun::Function, op::KLocalLiouvillian, v::ADoubleState)
    map_connections(fun, op.HnH_l, v)
    map_connections(fun, op.HnH_r, v)
    map_connections(fun, op.LLdag, v)
    return nothing
end

Base.eltype(::T) where {T<:KLocalLiouvillian} = eltype(T)
Base.eltype(t::Type{<:KLocalLiouvillian{H,T,A,B,C}}) where {H,T,A,B,C} =
    eltype(A)

function Base.show(io::IO, m::MIME"text/plain", op::KLocalLiouvillian)
    T    = eltype(op.HnH_l)
    h    = basis(op)

    print(io, "KLocalLiouvillian($T)\n\t- Hilb: $h")
end

function Base.show(io::IO, op::KLocalLiouvillian)
    T    = eltype(op.HnH_l)
    h    = basis(op)

    print(io, "KLocalLiouvillian($T) on Hilb: $h")
end


text:
module TDD

# Write your package code here.

#1

export DCneighbors
export ICneighbors
export AllComponents

#note to self I am using graph_input as argument for different functions which is fine since none of them are scope global but keep in mind...

function DCneighbors(graph_input, node::Int64) #Does type declaration for object input improve performance? for example `graph_input::Array{Array{Int64,1},1}`
                                               #The graph object can be of type Array{Array{Int64,1},1} or Array{Array{Any,1},1} if there are blank nodes [] in a directed graph
                                               #and Array{Array{Any,1},1} is apparently not a supertype of Array{Array{Int64,1},1} which will throw method error
                                               #any way around this? Would it give a performance boost?

    try
     Set(vcat(node, graph_input[node])) #Sets are unordered
    catch e
        if isa(e, BoundsError)
            error("node does not exist")
        else

        end
    end

end

#4 Adjacency Matrix Format
function DCneighbors(graph_input::Array{Bool,2}, node::Int64)

    print("test")
    #for this implementation, all the information about connectedness should be in the row of interest in Bool

end
    #turn node to row

function ICneighbors(graph_input, node::Int64) #this is way convoluted implementation. An exemplary example solution of iteration and recursion examples would be helpful

    #if node does not exist return error
    try
        graph_input[node]
    catch e
        error("node does not exist")
    end

    #list connected vertices of each vertex in graph
    AllNodes = []
    for i in 1:length(graph_input)
        push!(AllNodes, DCneighbors(graph_input, i::Int64)) #this returns a Set
    end

    #store connected vertices of all vertices that directly connect to `node`
    ConnectedToNode = []
    for i in 1:length(graph_input)
        if in(node, AllNodes[i]) == true
            union!(ConnectedToNode, AllNodes[i]) #ConnectedToNode is now an array cuz ConnectedToNode is an array
                                                 #note: a more logical implementation would be to just find all i that have connection to 'node'. But this works...
                                                 #wondering which way would perform faster...
        else
        end
    end

    ConnectedToNode_Array = []
    push!(ConnectedToNode_Array, Set([node])) #this starts ConnectedToNode_Array[1] so that ConnectedToNode pushes to ConnectedToNode_Array[2], so that
                                              #nodeset_input[z-1] does not eval nodeset_input[0]
                                              #note that this becomes a Set because pushing a Set into an empty array Any[]
    push!(ConnectedToNode_Array, Set( convert(Array{Int64,1}, ConnectedToNode) ) ) #https://stackoverflow.com/questions/35482527/how-do-i-change-the-data-type-of-a-julia-array-from-any-to-float64

    #now can start recursive function

    z_input = 2

    function recursor(graph_input_input, AllNodes_input, nodeset_input, z) #the reason I'm doing this again is I don't know how to get around reinitializing
                                                                         #the ConnectedToNode object with every run
                                                                         #so this way ConnectedToNode_Array is initialized outside the function...
                                                                         #also when I highlight this function to evaluate in REPL it seems to evaluate like 5 times...why? Seen this with other scripts too

        if ( nodeset_input[z] == nodeset_input[z - 1] )  &  ( length(nodeset_input) == (2 + length(graph_input_input)) )
                                                    #earlier had trouble with the in() test always returning false so recursor was missing union! steps. Now fixed, see below in.()
                                                    #At this point it worked for most cases except node 4,5,6,7,8 of train_graph
                                                    #Fixed by adding an & statement to if exit condition to fix problem with recursor starting again at node and triggering early exit. Force more repetitions bypasses early exit
            return ConnectedToNode_Array[z] #return ConnectedToNode_Array to see stack of all connected nodes for debugging

        else

            for i in 1:length(graph_input_input)
                if any( in.(nodeset_input[z], Ref(AllNodes_input[i])) ) == true #See docs for in() and any()
                    union!(nodeset_input[z], AllNodes_input[i]) #Set union with Set gives Set
                else
                end
            end

            push!(nodeset_input, nodeset_input[z]) #this already gives a Set since both those are Set
            recursor(graph_input_input, AllNodes_input, nodeset_input, z+1)

        end
    end

    #call
    recursor(graph_input, AllNodes, ConnectedToNode_Array, z_input)

end #function ICneighbors end

function AllComponents(graph_input)

    output = []

    for i in 1:length(graph_input)

    push!(output, ICneighbors(graph_input, i::Int64) )

    end #for loop

    return Set(unique(output))

end #function AllComponents

#Problem 4 generalize to adjacency matrix format

end #module end

ext: jl


text:
module ClimateSatellite

# Main file for the ClimateSatellite module that downloads climate-related
# variables from various satellite instruments

## Modules Used
using Dates
using DelimitedFiles
using FTPClient
using GeoRegions
using HDF5
using Logging
using NCDatasets
using Printf
using Statistics

## Exporting the following functions:
export
        clisatroot, clisatlonlat, clisatinfo!,
        clisatdownload, clisatanalysis, clisatcompile,
        clisatrawfol, clisatrawname, clisatrawread,
        clisatanafol, clisatananame, clisatanaread,
        clisatcmpfol, clisatcmpname, clisatcmpread,
        clisatrawregion, clisatrawpoint, clisatrawgrid,
        clisatanaregion, clisatanapoint, clisatanagrid,
        clisatsubregion,
        gpmdwn, gpmlonlat,
        mimicdwn, mimiclonlat

## Including other files in the module
include("frontend.jl")
include("backend.jl")
include("raw.jl")
include("analysis.jl")
include("compile.jl")
include("extract.jl")
include("gpm.jl")
include("mimic.jl")

end # module


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 1addcc10-5304-11eb-06d7-e1c69bd44c54
using Plots 	# plotting backend, takes a sec for first compile

# ╔═╡ 54f41f50-530c-11eb-1a8e-0b08dda08164
html"
Neutron Diffusion
"

# ╔═╡ 7d8eb0c0-5329-11eb-0741-dfb9627d9ea9
html"Note: the plotting backend, Plots.jl, takes a sec to compile; it's normal to take >30 seconds on first compile! 😳 Be patient! It's good!"

# ╔═╡ 69894990-530c-11eb-225c-312d8636850a
html"
First N, the number of grid points:
"

# ╔═╡ dc83a9b0-52b0-11eb-26f5-0d90caf487e5
@bind N html"<input type='range' min='20' max='200' step='10' value='100'>"

# ╔═╡ 76919200-52b2-11eb-3da1-6b26f1a52ec5
N

# ╔═╡ 7e14fdee-530c-11eb-2b26-c390d4d1997a
html"
τ, the time step for our discretized equation:
"

# ╔═╡ 86734650-52b2-11eb-34ef-e926ddb83096
@bind τ html"<input type='range' min='0.000000000001' max='0.0000000006' step='0.000000000005' value='0.0000000001'>"

# ╔═╡ 89a94400-52b2-11eb-1354-2da5b990bc13
τ

# ╔═╡ 9a0fece0-530c-11eb-15bd-29ebb4adb032
html"
Constants:
"

# ╔═╡ b08dc730-52b2-11eb-1fc9-2d30b44c56a0
begin
# Diffusion constant, D, and creation rate, C, for U-235
D = 10^5 	# m²/s
C = 10^8 	# 1/s
L = 1.0  	# m (length of our 1D space)
h = L/(N-1) # grid spacing
end

# ╔═╡ b539f410-5308-11eb-16eb-5354e7f4bae2
html"
Check stability of FTCS algorithm
"

# ╔═╡ 1003abb0-5305-11eb-2ac3-775355caece9
τ <= h^2/(2D)

# ╔═╡ 115a22b0-52ff-11eb-288b-f1f22602c21b
html"
The reaction will grow exponentially if the critical length, Lc, is less than the length of our 1D geometry.
"

# ╔═╡ fcf86f10-52ff-11eb-30b7-bb29ca63f98d
(L_c = π*√(D/C)) < L

# ╔═╡ 65f6b782-530d-11eb-1105-c1ede0ba670a
html"
Finally we set up & execute our main FTCS algorithm:
"

# ╔═╡ d7270d13-212f-42ca-8a9e-a914b5b6dbd1
html"This function will advance state of the neutron density array using a Forward in Time, Centered in Space algorithm (basically averaging neighboring points to solve the diffusion equation):"

# ╔═╡ 895801fb-e8ca-4290-9cb2-e0b0ee59a265
function step!(n)
	# Boundary Conditions are n[1] = n[end] = 0
	n[2:end-1] = n[2:end-1] .+
		( D*τ/h^2 ) * ( n[3:end] .+ n[1:end-2] .- 2n[2:end-1] ) .+
		C*τ*n[2:end-1]
end

# ╔═╡ 446c3d60-52b3-11eb-3f7a-777fb650bbde
begin
maxsteps::Int64 = 300
maxplots::Int64 = 50
plotstep::Int64 = maxsteps ÷ maxplots
end

# ╔═╡ 35da6fc0-5302-11eb-1a76-050089f82a58
begin
ρₙ = zeros(N) 	# initialize neutron density=0 everywhere
ρₙ[N÷2] = 1/h 	# IC is a delta spike in the middle
# Arrays which will contain data used for plotting later:
xplot = -L/2:h:L/2
nplot = zeros(N, maxplots)
tplot = zeros(maxplots)
iplot = 1 		# index used in loop below
for i ∈ 1:maxsteps
	step!(ρₙ)
	if i%plotstep == 1 	# plotting occasionally
		nplot[:, iplot] = copy(ρₙ)
		tplot[iplot] = i*τ
		iplot += 1
	end
end
end

# ╔═╡ 373273c0-5304-11eb-2ecf-5d731708e057
contourf(tplot, xplot, nplot, title="Neutron Diffusion",
		 xlabel="Time", ylabel="Position", c=:thermal)

# ╔═╡ fc5f0250-530c-11eb-15d9-bb6c5204bde9
surface(tplot, xplot, nplot, c=:thermal, title="Neutron Diffusion",
	    xlabel="Time", ylabel="Position", zlabel="Neutron Density", yticks=[])

# ╔═╡ aa5c5bd0-530a-11eb-3b31-298cca0ffbf1
@gif for i ∈ 1:maxplots
	plot(xplot, nplot[:, i], lw=3, legend=false, ylims=(0, nplot[N÷2, end]),
	     title="Neutron Diffusion Over $(τ*maxsteps) Seconds",
		 xlabel="Position", ylabel="Neutron Density")
end

# ╔═╡ Cell order:
# ╟─54f41f50-530c-11eb-1a8e-0b08dda08164
# ╟─7d8eb0c0-5329-11eb-0741-dfb9627d9ea9
# ╟─69894990-530c-11eb-225c-312d8636850a
# ╟─dc83a9b0-52b0-11eb-26f5-0d90caf487e5
# ╠═76919200-52b2-11eb-3da1-6b26f1a52ec5
# ╟─7e14fdee-530c-11eb-2b26-c390d4d1997a
# ╟─86734650-52b2-11eb-34ef-e926ddb83096
# ╠═89a94400-52b2-11eb-1354-2da5b990bc13
# ╟─9a0fece0-530c-11eb-15bd-29ebb4adb032
# ╠═b08dc730-52b2-11eb-1fc9-2d30b44c56a0
# ╟─b539f410-5308-11eb-16eb-5354e7f4bae2
# ╠═1003abb0-5305-11eb-2ac3-775355caece9
# ╟─115a22b0-52ff-11eb-288b-f1f22602c21b
# ╠═fcf86f10-52ff-11eb-30b7-bb29ca63f98d
# ╟─65f6b782-530d-11eb-1105-c1ede0ba670a
# ╟─d7270d13-212f-42ca-8a9e-a914b5b6dbd1
# ╠═895801fb-e8ca-4290-9cb2-e0b0ee59a265
# ╠═446c3d60-52b3-11eb-3f7a-777fb650bbde
# ╠═35da6fc0-5302-11eb-1a76-050089f82a58
# ╠═1addcc10-5304-11eb-06d7-e1c69bd44c54
# ╠═373273c0-5304-11eb-2ecf-5d731708e057
# ╠═fc5f0250-530c-11eb-15d9-bb6c5204bde9
# ╠═aa5c5bd0-530a-11eb-3b31-298cca0ffbf1

max_stars_repo_name: benharmonics/GeneralPhysics

max_issues_repo_name: benharmonics/GeneralPhysics

max_forks_repo_name: benharmonics/GeneralPhysics


text:
# # [Ecosystem Models](@id ecosystem_example)
#
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/examples/Ecosystem.ipynb)

using AlgebraicDynamics
using AlgebraicDynamics.DWDDynam
using AlgebraicDynamics.UWDDynam

using Catlab.CategoricalAlgebra
using Catlab.WiringDiagrams
using Catlab.Graphics
using Catlab.Programs

using LabelledArrays
using OrdinaryDiffEq
using Plots, Plots.PlotMeasures

# There will be 16 parameters in to the total ecosystem. With four distinct types:
# - parameters $\alpha$ represent the rate at which a species population grows
# - parameters $\beta$ represent the rate at which a population of prey declines in a predation interaction
# - parameters $\gamma$ represent the rate at which a population of predators grows in a predation interaction
# - parameters $\delta$ represent the rate at with a species population declines

params = LVector(αr=0.3, βrf=0.015, γrf=0.015, δf=0.7,
                 βrh=0.01, γrh=0.01, δh=0.5,
                 γfishh=0.001, βfishh=0.003,
                 αfish=0.35, βfishF=0.015, γfishF=0.015,
                 δF=0.7, βFs=0.017, γFs=0.017, δs=0.35)

# ## Land Ecosystem

# ### Rabbits and foxes

# A standard Lotka Volterra predator-prey model is the composition of three primitive resource sharers:

# 1. a model of rabbit growth: this resource sharer has dynamics $\dot r(t) = \alpha r(t)$ and one port which exposes the rabbit population.
# 2. a model of rabbit/fox predation: this resource sharer has dynamics $$\dot r(t) = -\beta r(t) f(t), \dot f(t) = \gamma r(t)f(t)$$ and two ports which expose the rabbit and fox populations respectively.
# 3. a model of fox population decline: this resource sharer has dynamics $\dot f(t) = -\delta f(t)$ and one port which exposes the fox population.

# However, there are not two independent rabbit populations -- one that grows and one that gets eaten by foxes. Likewise, there are not two independent fox populations -- one that declines and one that feasts on rabbits. To capture these interactions between the trio of resource sharers, we compose them by identifying the exposed rabbit populations and identifying the exposed fox populations.
# The syntax for this undirected composition is defined by an undirected wiring diagram.

## Define the primitive systems
dotr(u,p,t) = p.αr*u
dotrf(u,p,t) = [-p.βrf*u[1]*u[2], p.γrf*u[1]*u[2]]
dotf(u,p,t) = -p.δf*u

rabbit_growth = ContinuousResourceSharer{Float64}(1, dotr)
rabbitfox_predation = ContinuousResourceSharer{Float64}(2, dotrf)
fox_decline = ContinuousResourceSharer{Float64}(1, dotf)

## Define the composition pattern
rabbitfox_pattern = @relation (rabbits, foxes) begin
    rabbit_growth(rabbits)
    rabbitfox_predation(rabbits,foxes)
    fox_decline(foxes)
end

## Compose
submodels = Dict(
    :rabbit_growth => rabbit_growth,
    :rabbitfox_predation => rabbitfox_predation,
    :fox_decline => fox_decline
)
rabbitfox_system = oapply(rabbitfox_pattern, submodels)

# Previously, when we derived the Lotka-Volterra model via [undirected composition](https://algebraicjulia.github.io/AlgebraicDynamics.jl/dev/examples/Lotka-Volterra/#Undirected-composition), we by-hand defined the undirected wiring diagram that implements the composition pattern. In contrast, here we implement the same composition pattern as before but this time using the [`@relation` macro](https://algebraicjulia.github.io/Catlab.jl/stable/apis/programs/#Catlab.Programs.RelationalPrograms.@relation-Tuple). This strategy simplifies the definition and explicitly names the boxes and variables. We  visualize the composition pattern below.

to_graphviz(rabbitfox_pattern, box_labels = :name, junction_labels = :variable, edge_attrs=Dict(:len => ".75"))

# We can now construct an `ODEProblem` from the resource sharer `rabbitfox_system` and plot the solution.
u0 = [10.0, 100.0]
tspan = (0.0, 100.0)

prob = ODEProblem(rabbitfox_system, u0, tspan, params)
sol = solve(prob, Tsit5())

plot(sol, rabbitfox_system,
    lw=2,
    title = "Lotka-Volterra Predator-Prey Model",
    xlabel = "time", ylabel = "population size")

# ### Rabbits, foxes, and hawks
# Suppose we now have a three species ecosystem containing rabbits, foxes, and hawks. Foxes and hawks both prey upon rabbits but do not interact with each other. This ecosystem consists of five primitive systems which share variables.
# 1. rabbit growth:  $\dot r(t) = \alpha r(t)$
# 2. rabbit/fox predation:  $\dot r(t) = -\beta r(t) f(t), \dot f(t) = \delta r(t)f(t)$
# 3. fox decline:  $\dot f(t) = -\gamma f(t)$
# 4. rabbit/hawk predation: $\dot r(t) = -\beta' r(t)h(t), \dot h(t) = \delta' r(t)h(t)$
# 5. hawk decline:  $\dot h(t) = -\gamma' h(t)$

# This means the desired composition pattern has five boxes and many ports and wires to keep track of. Instead of implementing this composition pattern by hand, we construct it as a pushout.

## Define the composition pattern for rabbit growth
rabbit_pattern = @relation (rabbits,) -> rabbit_growth(rabbits)

## Define the composition pattern for the rabbit/hawk Lotka Volterra model
rabbithawk_pattern = @relation (rabbits, hawks) begin
    rabbit_growth(rabbits)
    rabbithawk_predation(rabbits,hawks)
    hawk_decline(hawks)
end

## Define transformations between the composition patterns
rabbitfox_transform  = ACSetTransformation((Box=[1], Junction=[1], Port=[1], OuterPort=[1]), rabbit_pattern, rabbitfox_pattern)
rabbithawk_transform = ACSetTransformation((Box=[1], Junction=[1], Port=[1], OuterPort=[1]), rabbit_pattern, rabbithawk_pattern)

## Take the pushout to define the composition pattern for the rabbit, fox, hawk system
rabbitfoxhawk_pattern = ob(pushout(rabbitfox_transform, rabbithawk_transform))

## Visualize the compsition pattern
to_graphviz(rabbitfoxhawk_pattern, box_labels = :name, junction_labels = :variable, edge_attrs=Dict(:len => ".9"))

#-
## Define the additional primitive systems
dotrh(u, p, t) = [-p.βrh*u[1]*u[2], p.γrh*u[1]*u[2]]
doth(u, p, t)  = -p.δh*u

rabbithawk_predation = ContinuousResourceSharer{Float64}(2, dotrh)
hawk_decline         = ContinuousResourceSharer{Float64}(1, doth)

## Compose
merge!(submodels, Dict(
    :rabbithawk_predation => rabbithawk_predation,
    :hawk_decline => hawk_decline
))

land_system = oapply(rabbitfoxhawk_pattern, submodels)

## Solve and plot
u0 = [10.0, 100.0, 50.0]
tspan = (0.0, 100.0)

prob = ODEProblem(land_system, u0, tspan, params)
sol = solve(prob, Tsit5())

plot(sol, land_system,
    lw=2,
    title = "Land Ecosystem",
    xlabel = "time", ylabel = "population size")

# Unfortunately, the hawks are going extinct in this model. We'll have to give hawks something else to eat!
#-
# ## Ocean Ecosystem

# Consider a ocean ecosystem containing three species —- little fish, big fish, and sharks -— with two predation interactions —- sharks eat big fish and big fish eat little fish.

# This ecosystem can be modeled as the composition of 3 machines:
# 1. Evolution of the little fish population:  this machine has one exogenous variable which represents a population of predators $h(t)$ that hunt little fish. This machine has one output which emits the little fish population. The dynamics of this machine is the driven ODE $$\dot f(t) = \alpha f(t) - \beta f(t)h(t)$$
# 2. Evolution of the big fish population:  this machine has two exogenous variables which represent a population of prey $e(t)$ that are eaten by big fish and a population of predators $h(t)$ which hunt big fish. This machine has one output which emits the big fish population. The dynamics of this machine is the drive ODE $$\dot F(t) = \gamma F(t)e(t) - \delta F(t) - \beta'F(t)h(t)$$
# 3. Evolution of the shark population:  this machine has one exogenous variable which represents a population of prey $e(t)$ that are eaten by sharks. This machine has one output which emits the shark population. The dynamics of this machine is the driven ODE $$\dot s(t) = \gamma's(t)e(t) - \delta's(t)$$

## Define the primitive systems
dotfish(f, x, p, t) = [p.αfish*f[1] - p.βfishF*x[1]*f[1]]
dotFISH(F, x, p, t) = [p.γfishF*x[1]*F[1] - p.δF*F[1] - p.βFs*x[2]*F[1]]
dotsharks(s, x, p, t) = [p.γFs*s[1]*x[1]-p.δs*s[1]]

fish   = ContinuousMachine{Float64}(1,1,1, dotfish,   (f,p,t)->f)
FISH   = ContinuousMachine{Float64}(2,1,1, dotFISH,   (F,p,t)->F)
sharks = ContinuousMachine{Float64}(1,1,1, dotsharks, (s,p,t)->s)

# We compose these machines by (1) sending the output of the big fish machine as the input to both the little fish and shark machines and (2) sending the output of the little fish and shark machines as the inputs to the big fish machine.
# The syntax for this directed composition is given by a directed wiring diagram.

## Define the composition pattern
ocean_pattern = WiringDiagram([], [:fish, :Fish, :shark])
fish_box = add_box!(ocean_pattern, Box(:fish, [:pop], [:pop]))
Fish_box = add_box!(ocean_pattern, Box(:Fish, [:pop, :pop], [:pop]))
shark_box = add_box!(ocean_pattern, Box(:shark, [:pop], [:pop]))

output_box = output_id(ocean_pattern)
add_wires!(ocean_pattern, Pair[
    (fish_box, 1)  => (Fish_box, 1),
    (shark_box, 1) => (Fish_box, 2),
    (Fish_box, 1)  => (fish_box, 1),
    (Fish_box, 1)  => (shark_box, 1),
    (fish_box, 1)  => (output_box, 1),
    (Fish_box, 1)  => (output_box, 2),
    (shark_box, 1) => (output_box, 3)
])

## Visualize the composition pattern
to_graphviz(ocean_pattern, orientation=TopToBottom)

#-

## Compose
ocean_system = oapply(ocean_pattern, [fish, FISH, sharks])

## Solve and plot
u0 = [100.0, 10, 2.0]
tspan = (0.0, 100.0)

prob = ODEProblem(ocean_system, u0, tspan, params)
sol = solve(prob, Tsit5())

plot(sol, ocean_system,
    lw=2, title = "Ocean Ecosystem",
    xlabel = "time", ylabel = "population size"
)

# ## Total ecosystem
# ### Another layer of composition

# We will introduce a final predation interaction -- hawks eat little fish --  which will combine the land and ocean ecosystems.
# The composition will be as resource shareres so the first thing we will do is use the dynamics of the machine `ocean_system` to define the dynamics of a resource sharer. We will also define a resource sharer that models hawk/little fish predation.

## Define the additional primitive systems
ocean_system_rs = ContinuousResourceSharer{Float64}(3, (u,p,t) -> eval_dynamics(ocean_system, u, [], p))

dothf(u,p,t) = [p.γfishh*u[1]*u[2], -p.βfishh*u[1]*u[2]]
fishhawk_predation = ContinuousResourceSharer{Float64}(2, dothf)

## Define the composition pattern
eco_pattern = @relation (rabbits, foxes, hawks, littlefish, BigFish, sharks)  begin
    turf(rabbits,foxes,hawks)
    air(hawks, littlefish)
    surf(littlefish, BigFish, sharks)
end

## Visualize the composition pattern
to_graphviz(eco_pattern, box_labels = :name, junction_labels = :variable, edge_attrs=Dict(:len => ".75"))

#-
## Compose

ecosystem = oapply(eco_pattern, Dict(
    :turf => land_system,
    :air => fishhawk_predation,
    :surf => ocean_system_rs)
)

# We can now plot the evolution of the total ecosystem.

## Solve and plot
u0 = [100.0, 50.0, 20.0, 100, 10, 2.0]
tspan = (0.0, 100.0)

prob = ODEProblem(ecosystem, u0, tspan, params)
sol = solve(prob, Tsit5())
plot(sol, ecosystem,
    lw=2, title = "Total Ecosystem",
    xlabel = "time", ylabel = "population size"
)

# Let's zoom in on a narrower time-window.
tspan = (0.0, 30.0)

prob = ODEProblem(ecosystem, u0, tspan, params)
sol = solve(prob, Tsit5())
plot(sol, ecosystem,
    lw=2, title = "Total Ecosystem",
    xlabel = "time", ylabel = "population size"
)
# As a sanity check we can define the rates for the hawk/little fish predation to be 0. This decouples the land and ocean ecosystems. As expected, the plot shows the original evolution of the land ecosystem overlayed with the original evolution of the ocean ecosystem. This shows that they two ecosystems now evolve independently.

tspan = (0.0, 100.0)
params.βfishh = 0; params.γfishh = 0

prob = ODEProblem(ecosystem, u0, tspan, params)
sol = solve(prob, Tsit5())
plot(sol, ecosystem,
    lw=2, title = "Decoupled Ecosystem",
    xlabel = "time", ylabel = "population size"
)


text:
function card_dict{C<:NastranCard}(::Type{C},deck::NastranDeck)
    Dict{Int64,C}(map(filter(c -> typeof(c)<:C,deck)) do card
        (card.id,card)
    end)
end

# type Grid
#     id::Int64
#     csys_id::Int64
#     xyz::XYZ
# end

# type Coord
#     id::Int64
#     csys::CoordinateSystem
# end

# Grid(card::GRID) = Grid(card.id,card.csys_id,XYZ(card.x,card.y,card.z))
# Grid(card::SPOINT) = Grid(card.id,0,zero(XYZ))

# immutable CoordSet
#     coords::Vector{Coord}
#     local_xyzs::Vector{Grid}
#     global_xyzs::Vector{Grid}
# end

# sort_by_id(x::Int64) = x
# sort_by_id(x::Grid) = x.id
# sort_by_id(x::Coord) = x.id

# function get_coord(coordset::CoordSet,id)
#     i = searchsorted(coordset.coords,id,by=sort_by_id)
#     @assert length(i) == 1
#     coordset.coords[first(i)]
# end

# function get_global_xyz(coordset::CoordSet,id)
#     i = searchsorted(coordset.global_xyzs,id,by=sort_by_id)
#     @assert length(i) == 1
#     coordset.global_xyzs[first(i)]
# end

# function CoordinateSystem(global_xyzs::Vector{Grid},card::CORD1R)
#     iA = searchsorted(global_xyzs,card.gridA,by=sort_by_id)
#     @assert length(iA) == 1
#     gridA = global_xyzs[first(iA)]

#     iB = searchsortedfirst(global_xyzs,card.gridB,by=sort_by_id)
#     @assert length(iB) == 1
#     gridB = global_xyzs[first(iB)]

#     iC = searchsortedfirst(global_xyzs,card.gridC,by=sort_by_id)
#     @assert length(iC) == 1
#     gridC = global_xyzs[first(iC)]

#     if gridA.csys_id == 0 && gridB.csys_id == 0 && gridC.csys_id == 0
#         Coord(card.id,CoordinateSystem(gridA.xyz,gridB.xyz,gridC.xyz))
#     end
# end

# function CoordinateSystem(global_xyzs::Vector{Grid},card::CORD2R)
#     if card.csys_id == 0
#         pA = XYZ(card.gridAx,card.gridAy,card.gridAz)
#         pB = XYZ(card.gridBx,card.gridBy,card.gridBz)
#         pC = XYZ(card.gridCx,card.gridCy,card.gridCz)
#         Coord(card.id,CoordinateSystem(pA,pB,pC))
#     else
#         throw(MethodError("Not Implemented"))
#     end
# end

# function CoordSet(grid_cards::Dict{Int64,GridCard},csys_cards::Dict{Int64,CoordinateSystemCard})
#     grid_ids = sort(collect(keys(grid_cards)))
#     local_xyzs::Vector{Grid} = map(grid_ids) do id
#         Grid(grid_cards[id])
#     end
#     # global_xyzs starts the same as local_xyzs. As the coordinates are converted the csys_id is
#     # updated to be 0 and the global coordinate is changed
#     global_xyzs = copy(local_xyzs)
#     csys_ids = sort(collect(keys(csys_cards)))
#     converted = false
#     # coords is initiated as an empty vector that is filled in as coordinate systems are defined
#     coords = Vector{Coord}(length(csys_ids))
#     # Loop through all the coordinate systems and grids repeatedly until all the coordinate
#     # systems are defined and therefore the global locations of all the grids are found
#     while !converted
#         # Set converted to true at the top of every loop. converted gets set to false when
#         # it either fails to define a coordinate system or convert a grid to global coordinates
#         converted = true
#         for (i,id) in enumerate(csys_ids)
#             if isdefined(coords,i)
#                 continue
#             else
#                 csys = CoordinateSystem(global_xyzs,csys_cards[id])
#                 if csys != nothing
#                     coords[i] = csys
#                 else
#                     converted = false
#                 end
#             end
#         end
#         for grid in global_xyzs
#             if grid.csys_id == 0
#                 continue
#             else
#                 i_csys = searchsortedfirst(csys_ids,grid.csys_id)
#                 if isdefined(coords,i_csys)
#                     coord = coords[i_csys]
#                     grid.csys_id = 0
#                     grid.xyz = to_global(coord.csys,grid.xyz)
#                 else
#                     converted = false
#                 end
#             end
#         end
#     end
#     CoordSet(coords,local_xyzs,global_xyzs)
# end

type NastranModel
    grids::Dict{Int64,GridCard}
    elements::Dict{Int64,ElementCard}
    properties::Dict{Int64,PropertyCard}
    materials::Dict{Int64,MaterialCard}
    constraints::Dict{Int64,ConstraintCard}
    coordinate_systems::Dict{Int64,CoordinateSystemCard}
    sets::Dict{Int64,SetCard}
    # coordset::CoordSet
end

function show(io::IO,model::NastranModel)
    print(io,"NastranModel(")
    print(io,length(model.grids),",")
    print(io,length(model.elements),",")
    print(io,length(model.properties),",")
    print(io,length(model.materials),",")
    print(io,length(model.constraints),",")
    print(io,length(model.coordinate_systems),")")
    print(io,length(model.sets),",")
end

function CoordinateSystem(model::NastranModel,card::CORD1R)
    if card.gridA in keys(model.grids) && card.gridB in keys(model.grids) && card.gridC in keys(model.grids)
        xyzA = get_global_xyz(model,card.gridA)
        xyzB = get_global_xyz(model,card.gridB)
        xyzC = get_global_xyz(model,card.gridC)
        CoordinateSystem(xyzA,xyzB,xyzC)
    else
        throw(ExceptionError("Not all grids found"))
    end
end

function CoordinateSystem(model::NastranModel,card::CORD2R)
    if card.csys_id == 0
        pA = XYZ(card.gridAx,card.gridAy,card.gridAz)
        pB = XYZ(card.gridBx,card.gridBy,card.gridBz)
        pC = XYZ(card.gridCx,card.gridCy,card.gridCz)
        CoordinateSystem(pA,pB,pC)
    else
        throw(MethodError("Not Implemented"))
    end
end

function CoordinateSystem(model::NastranModel,id::Int64)
    card = model.coordinate_systems[id]
    CoordinateSystem(model,card)
end

function get_global_xyz(model::NastranModel,id::Int64)
    grid = model.grids[id]
    if grid.csys_id == 0
        XYZ(grid.x,grid.y,grid.z)
    else
        coord = CoordinateSystem(model,grid.csys_id)
        local_xyz = XYZ(grid.x,grid.y,grid.z)
        to_global(coord,local_xyz)
    end
end

function NastranModel(filename::AbstractString)
    deck = NastranDeck(filename)
    grid_cards = card_dict(GridCard,deck)
    element_cards = card_dict(ElementCard,deck)
    property_cards = card_dict(PropertyCard,deck)
    material_cards = card_dict(MaterialCard,deck)
    constraint_cards = card_dict(ConstraintCard,deck)
    set_cards = card_dict(SetCard,deck)
    coordinate_system_cards = card_dict(CoordinateSystemCard,deck)
    # coordset = CoordSet(grid_cards,coordinate_system_cards)
    NastranModel(grid_cards,
                 element_cards,
                 property_cards,
                 material_cards,
                 constraint_cards,
                 coordinate_system_cards,
                 set_cards)
end


text:
"""Object life cycle framework

This is interface for handling the life cycle context of an object.
"""
module LifeCycle
    "Abstract life cycle object"
    abstract Object

    @enum(State,
        INITIALIZED,# Initialized but not yet started.
        STARTING,   # In the process of starting.
        STARTED,    # Has started.
        STOPPING,   # Stopping is in progress.
        STOPPED,    # Has stopped.
        INVALID)    # Something happened

    doc"""Life cycle states enumeration"""
    State
end

"Starts a life cycle"
start{T<:LifeCycle.Object}(lco::T) = throw(AssertionError("Function 'start' is not implemented for type $(typeof(lco))"))

"Stops a life cycle"
stop{T<:LifeCycle.Object}(lco::T) = throw(AssertionError("Function 'stop' is not implemented for type $(typeof(lco))"))

"Returns a life cycle state"
state{T<:LifeCycle.Object}(lco::T) = isdefined(lco, :state) ? lco.state : throw(AssertionError("Field 'state' is not defined in type $(typeof(lco))"))

"Sets a life cycle state"
state!{T<:LifeCycle.Object}(lco::T, st::LifeCycle.State) = isdefined(lco, :state) ? (lco.state = st) : throw(AssertionError("Field 'state' is not defined in type $(typeof(lco))"))

""" Abstract log message

Message implementations that can be logged. Messages can act as wrappers
around objects so that user can have control over converting objects to strings
when necessary without requiring complicated formatters and as a way to manipulate
the message based on information available at runtime such as the locale of the system.
"""
abstract Message
typealias MESSAGE Nullable{Message}

""" Returns message format as string. """
format(msg::Message)     = throw(AssertionError("Function 'format' is not implemented"))

""" Returns message formatted as string. """
formatted(msg::Message)  = throw(AssertionError("Function 'formatted' is not implemented"))

""" Returns message parameters, if any. """
parameters(msg::Message) = throw(AssertionError("Function 'parameters' is not implemented"))

""" Log event

Provides contextual information about a logged message.
"""
abstract Event

""" Returns the fully qualified module name (FQMN) of the caller of the logging API """
fqmn(evnt::Event) = isdefined(evnt, :fqmn) ? evnt.fqmn : throw(AssertionError("Define field or method 'fqmn' in type $(typeof(evnt))"))

""" Gets the level. """
level(evnt::Event) = isdefined(evnt, :level) ? evnt.level : throw(AssertionError("Define field or method 'level' in type $(typeof(evnt))"))

""" Gets the logger name. """
logger(evnt::Event) = isdefined(evnt, :logger) ? evnt.logger : throw(AssertionError("Define field or method 'logger' in type $(typeof(evnt))"))

""" Gets the Marker associated with the event. """
marker(evnt::Event) = isdefined(evnt, :marker) ? evnt.marker : throw(AssertionError("Define field or method 'marker' in type $(typeof(evnt))"))

""" Gets the message associated with the event. """
message(evnt::Event) = isdefined(evnt, :message) ? evnt.message : throw(AssertionError("Define field or method 'message' in type $(typeof(evnt))"))

""" Gets event timestamp. """
timestamp(evnt::Event) = isdefined(evnt, :timestamp) ? evnt.timestamp : throw(AssertionError("Define field or method 'timestamp' in type $(typeof(evnt))"))

""" Abstract layout

Lays out a `Event` in different formats.
"""
abstract Layout
typealias LAYOUT Nullable{Layout}

""" Returns the header as byte array for the layout format. """
header(lyt::Layout) = throw(AssertionError("Function 'header' is not implemented for type $(typeof(lyt))"))

""" Returns the footer as byte array for the layout format. """
footer(lyt::Layout) = throw(AssertionError("Function 'footer' is not implemented for type $(typeof(lyt))"))

"""
    serialize(lyt::Layout, evnt::Event) -> Vector{UInt8}

Formats the even suitable for display into byte array.
"""
serialize(lyt::Layout, evnt::Event) = throw(AssertionError("Function 'serialize' is not implemented for type $(typeof(lyt))"))
serialize(lyt::LAYOUT, evnt::Event) = !isnull(lyt) ? serialize(get(lyt), evnt) : convert(Vector{UInt8}, "No layout present!")

""" Returns the content type output by this layout """
contenttype(lyt::Layout) = throw(
    AssertionError("""\n
        Every layout should provide a content type.
        Implement function `contenttype` for type $(typeof(lyt))
        that returns a proper content type (e.g. "text/plain")
    """))

""" Abstract type for layouts that result in a String """
abstract StringLayout <: Layout
contenttype(lyt::StringLayout) = "text/plain"

"""
    string(lyt::Layout, evnt::Event) -> AbstractString

Formats the even suitable for display into string.
"""
string(lyt::StringLayout, evnt::Event) = throw(AssertionError("Function 'string' is not implemented for type $(typeof(lyt))"))
string(lyt::LAYOUT, evnt::Event) = !isnull(lyt) ? string(get(lyt), evnt) : "No layout present!"

""" Abstract event filter """
abstract Filter <: LifeCycle.Object
typealias FILTER Nullable{Filter}

""" Supports event filtering

All filterable types should have `state::LifeCycle.State` and `filter::FILTER` fields.
"""
abstract Filterable <: LifeCycle.Object

module FilterResult
    @enum ResultType ACCEPT NEUTRAL DENY
end

""" Returns `FILTER` object """
filter(fltrbl::Filterable) = isdefined(fltrbl, :filter) ? fltrbl.filter : throw(AssertionError("Define field 'filter' in type $(typeof(fltrbl))"))
"""Determines if the event should be filtered.
Returns `true` if the event should be filtered, `false` otherwise.
"""
isfiltered(flt::FILTER, evnt::Event) = !isnull(flt) && filter(get(flt), evnt) == FilterResult.DENY
isfiltered(fltrbl::Filterable, evnt::Event) = isfiltered(filter(fltrbl), evnt)

""" Abstract appender

Any appender implementation must have two fields:

- `name`::AbstractString - appender name for reference
- `layout`::Nullable{`Layout`} - layout object for output modification

In addition to basic fields, every appender should have a method `append!`:

- append!(apnd::`Appender`, evnt::`Event`) - appends event

*Note:* All derived types should have a constructor which accepts `Dict{Symbol,Any}` object as a parameter.
Passed dictionary object should contain various initialization parameters. One of the parameters should have a key `:layout` with `Layout` object as its value.
"""
abstract Appender <: Filterable
typealias APPENDER Nullable{Appender}
typealias APPENDERS Dict{AbstractString, Appender}

""" Returns name of the appender """
name(apnd::Appender) = isdefined(apnd, :name) ? apnd.name : throw(AssertionError("Define field 'name' in type $(typeof(apnd))"))

""" Returns layout of the appender """
layout(apnd::Appender) = isdefined(apnd, :layout) ? apnd.layout : throw(AssertionError("Define field 'layout' in type $(typeof(apnd))"))

""" Adds event to the appender """
append!(apnd::Appender, evnt::Event) = throw(AssertionError("Function 'append!' is not implemented for type $(typeof(apnd))"))

""" Ignore exception in appender.

Some appenders need to propagate exceptions back to the application.
When `ignoreexceptions` returns `false`, the exception will be allow to percolate.
"""
ignoreexceptions(apnd::Appender) = true

""" Abstract configuration

Any configuration type implementation must have two fields:

- `name`::AbstractString, the configuration name for a reference
- `source`::AbstractString, a source of the configuration
- `state`::LifeCycle.State, a life cycle state
- `root`::LoggerConfig, the root logger configuration

Any configuration must implement following set of methods:

- logger(Configuration, AbstractString) -> LoggerConfig
- appender(Configuration, AbstractString) -> Appender
- loggers(Configuration) -> Dict{AbstractString, LoggerConfig}
- appenders(Configuration) -> Dict{AbstractString, Appender}
"""
abstract Configuration <: Filterable
typealias CONFIGURATION Nullable{Configuration}

show(io::IO, cfg::Configuration) = print(io, """Configuration($(cfg.name), $(isempty(cfg.source) ? "" : cfg.source*", ")$(cfg.state))""")

"Configuration name for a reference"
name(cfg::Configuration) = isdefined(cfg, :name) ? cfg.name : throw(AssertionError("Define field 'name' in type $(typeof(cfg))"))

"Returns the source of this configuration"
source(cfg::Configuration) = isdefined(cfg, :source) ? cfg.source : throw(AssertionError("Define field 'source' in type $(typeof(cfg))"))

"Returns the root logger configuration"
root(cfg::Configuration) = isdefined(cfg, :root) ? cfg.root : throw(AssertionError("Define field 'root' in type $(typeof(cfg))"))

"Return a list of `Logger`s from the configuration"
loggers(cfg::Configuration) = throw(AssertionError("Function 'loggers' is not implemented for type $(typeof(cfg))"))

"Returns  `Appender`  with the specified `name`"
appender(cfg::Configuration, name::AbstractString) = throw(AssertionError("Function 'appender' is not implemented for type $(typeof(cfg))"))

"Returns a list of `Appender`s from the configuration"
appenders(cfg::Configuration) = throw(AssertionError("Function 'appenders' is not implemented for type $(typeof(cfg))"))

"""
Abstract Logger Type

All implementations of this type must have two methods:

- `name`, returns a logger name as string is needed to be implemented
- `message`, returns a message construction type (factory)
"""
abstract AbstractLogger

typealias LOGGERS Dict{AbstractString, AbstractLogger}

"Returns a factory type that generates messages"
message(lgr::AbstractLogger) = throw(AssertionError("Function 'message' is not implemented for type $(typeof(lgr))"))

"Returns a name of the logger"
name(lgr::AbstractLogger) = throw(AssertionError("Function 'name' is not implemented for type $(typeof(lgr))"))

"Returns a status level of the logger"
level(lgr::AbstractLogger) = throw(AssertionError("Function 'level' is not implemented for type $(typeof(lgr))"))

"Sets a status level of the logger"
level!(lgr::AbstractLogger, lvl::Level.EventLevel) = throw(AssertionError("Function 'level!' is not implemented for type $(typeof(lgr))"))

# Aliases
typealias MARKER Nullable{Symbol}
typealias FACTORY Nullable{DataType}
typealias NAME Nullable{AbstractString}
typealias PROPERTIES Dict{AbstractString, AbstractString}


text:
using .Distances

@adjoint function (::SqEuclidean)(x::AbstractVector, y::AbstractVector)
  δ = x .- y
  function sqeuclidean(Δ::Real)
    x̄ = (2 * Δ) .* δ
    return x̄, -x̄
  end
  return sum(abs2, δ), sqeuclidean
end

@adjoint function colwise(s::SqEuclidean, x::AbstractMatrix, y::AbstractMatrix)
  return colwise(s, x, y), function (Δ::AbstractVector)
    x̄ = 2 .* Δ' .* (x .- y)
    return nothing, x̄, -x̄
  end
end

@adjoint function pairwise(s::SqEuclidean, x::AbstractMatrix, y::AbstractMatrix; dims::Int=2)
  if dims==1
    return pairwise(s, x, y; dims=1), ∇pairwise(s, transpose(x), transpose(y), transpose)
  else
    return pairwise(s, x, y; dims=dims), ∇pairwise(s, x, y, identity)
  end
end

∇pairwise(s, x, y, f) =
  function(Δ)
    x̄ = 2 .* (x * Diagonal(vec(sum(Δ; dims=2))) .- y * transpose(Δ))
    ȳ = 2 .* (y * Diagonal(vec(sum(Δ; dims=1))) .- x * Δ)
    return (nothing, f(x̄), f(ȳ))
  end

@adjoint function pairwise(s::SqEuclidean, x::AbstractMatrix; dims::Int=2)
  if dims==1
    return pairwise(s, x; dims=1), ∇pairwise(s, transpose(x), transpose)
  else
    return pairwise(s, x; dims=dims), ∇pairwise(s, x, identity)
  end
end

∇pairwise(s, x, f) =
  function(Δ)
    d1 = Diagonal(vec(sum(Δ; dims=1)))
    d2 = Diagonal(vec(sum(Δ; dims=2)))
    return (nothing, x * (2 .* (d1 .+ d2 .- Δ .- transpose(Δ))) |> f)
  end

@adjoint function (::Euclidean)(x::AbstractVector, y::AbstractVector)
  D = x.-y
  δ = sqrt(sum(abs2, D))
  function euclidean(Δ::Real)
    x̄ = (Δ / δ) .* D
    return x̄, -x̄
  end
  return δ, euclidean
end

@adjoint function colwise(s::Euclidean, x::AbstractMatrix, y::AbstractMatrix)
  d = colwise(s, x, y)
  return d, function (Δ::AbstractVector)
    x̄ = (Δ ./ d)' .* (x .- y)
    return nothing, x̄, -x̄
  end
end

@adjoint function pairwise(::Euclidean, X::AbstractMatrix, Y::AbstractMatrix; dims=2)
  D, back = pullback((X, Y) -> pairwise(SqEuclidean(), X, Y; dims = dims), X, Y)
  D .= sqrt.(D)
  return D, Δ -> (nothing, back(Δ ./ (2 .* D))...)
end

@adjoint function pairwise(::Euclidean, X::AbstractMatrix; dims=2)
  D, back = pullback(X -> pairwise(SqEuclidean(), X; dims = dims), X)
  D .= sqrt.(D)
  return D, function(Δ)
    Δ = Δ ./ (2 .* D)
    Δ[diagind(Δ)] .= 0
    return (nothing, first(back(Δ)))
  end
end


text:
@inline initial_condition_generator(rng,n,model) = (φ=[3], X=[5.0], Y=[0.0])

const tau = first_exit_y(0.0,model.dq.model.b+1.0)

const x = range(0,1.6,length=10_001)


text:
mult(x, y) = x * y
typeof(mult) #> Function
m = mult
m(6, 6) #> 36

plustwo = function (x)
              x + 2
          end
plustwo(3) #> 5

3+4 #> 7
+(3,4) #> 7

function numerical_derivative(f, x, dx=0.01)
  derivative = (f(x+dx) - f(x-dx))/(2*dx)
  return derivative
end

f = x -> 2x^2 + 30x + 9
println(numerical_derivative(f, 1, 0.001)) #> 33.99999999999537

function derivative(f)
    return function(x)
        # pick a small value for h
        h = x == 0 ? sqrt(eps(Float64)) : sqrt(eps(Float64)) * x
        xph = x + h
        dx = xph - x
        f1 = f(xph) # evaluate f at x + h
        f0 = f(x) # evaluate f at x
        return (f1 - f0) / dx  # divide by h
    end
end

# closure:
function counter()
    n = 0
    () -> n += 1, () -> n = 0
end

(addOne, reset) = counter()
addOne()
addOne()
addOne()
reset()

# currying:
function add(x)
    return function f(y)
        return x + y
    end
end
add(1)(2) #> 3
add(x) = f(y) = x + y
add(x) = y -> x + y

# nested:
function afun(x)
    z = x * 2
    function b(z)
        z += 1
    end
    b(z)
end

d = 5
afun(d) #> 11

# recursive:
sum(n) =  n > 1 ? sum(n-1) + n : n
sum(100) #> 5050
fib(n) = n < 2 ? n : fib(n-1) + fib(n-2)
println(fib(25)) #> 75025

max_stars_repo_path: Module 1/Chapter03/first_class.jl

max_issues_repo_path: Module 1/Chapter03/first_class.jl

max_forks_repo_path: Module 1/Chapter03/first_class.jl


text:
#=
Contains functions for processing velocity histograms, smoothing them,
aligning them, and more.
Takes advantage of the Integer nature of velocity!
=#
using DSP, Statistics, MusicManipulations

##############################
# Histograms
##############################
"""
    velhist(notes::Notes)

Get a histogram of the velocity.
"""
function velhist(notes::Notes{N}) where N
    mv = maximum(n.velocity for n in notes) + 50 # some buffer
    h = zeros(Int, mv+1)
    for n in notes
        h[n.velocity+1] += 1
    end
    return h
end

function velhist(v::Vector{<:Real})
    h = zeros(Int, maximum(v)+1)
    for vi in v
        h[vi+1] += 1
    end
    return h
end

"""
    velhists(notes)
Get a dictionary of velocity histograms,
each for a specific pitch of the `notes`.
(can also accept a dictionary directly)
"""
function velhists(notes::Notes)
    sep = separatepitches(notes)
    velhists(sep)
end
function velhists(sep::Dict)
    histos = Dict{UInt8, Vector{Float64}}()
    for k in keys(sep)
        histos[k] = velhist(sep[k])
    end
    return histos
end

"""
    smoothhist(h, w = PEAKW)

Smoothen out a histogram via convolution with Gaussian of window size `w`.
"""
smoothhist(notes::Notes, w = PEAKW) = smoothhist(velhist(notes), w)
function smoothhist(hist::AbstractVector, w = PEAKW)
    #prepare
    normal1 = maximum(hist)
    gaußk = exp.(-range(-0.3, 0.3, length = w).^2)

    #smoothen
    histc = DSP.conv(hist,gaußk)
    histc = DSP.conv(histc,gaußk)
    histc = DSP.conv(histc,gaußk)
    histc = DSP.conv(histc,gaußk)

    #normalize, cut and clean
    normal2 = maximum(histc)
    histc = histc./(normal2/normal1)
    histc = histc[(2w-1):end-(2w-2)]
    for i = 1:length(histc)
        if histc[i] ≤ 1
            histc[i] = 0
        end
    end
    return histc
end

"good value for peak detection & alignment"
const PEAKW = 7

smoothhists(s::Dict{K, <: Notes}) where {K} = smoothhists!(velhists(s))
function smoothhists!(vhists::Dict, w = PEAKW)
    for k in keys(vhists)
        vhists[k] = smoothhist(vhists[k], w)
    end
    return vhists
end

##############################
# Velocity peaks
##############################
"""
    find_peaks(h, thres = 2)
Find the peak indices of the given (expected smooth) histogram `h`.
Only consider peaks that have at least `thres` height.
"""
function find_peaks(h, thres = 2)
    peaks = findall(i -> h[i-1] < h[i] > h[i+1], 2:length(h)-1)
    del = findall(i -> h[i] < thres, peaks)
    deleteat!(peaks, del)
    return peaks .+ 1 # because the search range starts from 2
end

"""
    aling_peaks(notes; thres = 2, w = PEAKW, method = :scale)
Separate the `notes` by pitch, and then for each pitch calculate the
velocity histogram. Align these histograms based on the peak of
highest velocity (the "forte" peak), and return the aligned histograms
as a dictionary `pitch => histogram`.

The `method` indicates whether the `:scale` the histograms or `:shift` them
to match.
"""
function aling_peaks(notes; w = PEAKW, thres = 2, method = :scale)
    sepa = align_notes(notes; w = w, thres = thres, method = method)
    vhists = velhists(sepa)
    return vhists
end

function align_notes(notes; w = PEAKW, thres = 2, method = :scale)
    sepa = separatepitches(notes)
    vhists = velhists(sepa)
    smoothhists!(vhists, w)
    allpeaks = Dict((k, find_peaks(vhists[k], thres)) for k in keys(vhists))
    center = maximum(maximum(p) for (k, p) in allpeaks if length(p) != 0)
    # modify the velocity for each notes
    for (k, n) in sepa
        length(allpeaks[k]) == 0 && continue
        forte = maximum(allpeaks[k])
        align!(n, center, forte, method)
    end
    return sepa
end

function align!(notes::Notes, center, forte, method)
    center == forte && return

    for i in 1:length(notes)
        n = notes[i]
        if method == :shift
            n.velocity += center - forte
        elseif method == :scale
            n.velocity = round(Int, n.velocity * (center/forte))
        end
    end
    return
end

#######################
# Aligned drum notes
#######################
"""
    aligned_drums_notes()
Return a vector of all the drum playalong notes processed in the following way:
1. Fake hihat notes are removed.
2. Remove drum notes that we decide to not analyze due to their nature.
3. Shift time of drums to 0 of grid.
4. Channel notes into the 5 main channels (Lennart's thesis).
5. Align the forte peaks of each channel.
"""
function aligned_drums_notes()
    playdir = datadir("exp_raw", "playalongs")
    fdir = readdir(playdir)
    total = []
    for (i, file) in enumerate(fdir)
        midi = readMIDIFile(joinpath(playdir, file))
        notes = Drums.getnotes_td50(midi, 2)
        notes = drums_cleanup(notes)
        channeled = channel_direct(notes)
        anotes = align_notes(channeled; thres = 3)
        # all good but for some files we need manual adjustments
        if i == 3
            for n in anotes[2] # snare
                v = n.velocity
                n.velocity = v > 45 ? v - 45 : 0
            end
        end
        push!(total, anotes)
    end
    return total
end

max_stars_repo_name: Datseris/VelocityMTD

max_issues_repo_name: Datseris/VelocityMTD

max_forks_repo_name: Datseris/VelocityMTD


text:
# functions for random walk SMC
using StatsBase
include("smc_utils.jl")

# include("Utils.jl")

# abstract type EdgeList <: Array{Int64,2} end # stores edges in a graph
# abstract type VertexMap <: Vector end
# abstract type ParticlePath <: Vector end # stores retained smc particle path

struct ParticleState
  # structure for keeping track of state information for particle states

  n_edges::Vector{Int64} # how many edges have been inserted
  edge_list::Array{Int64,2} # edges in order of insertion (pre-allocate)
  edge_idx_list::Array{Int64,1} # same as above, but indexes into data edge_list
  nbd_list::Array{Int64,2} # vector of sparse vectors, each of which stores vertex neighborhoods
  n_vertices::Vector{Int64} # number of edges in current graph
  new_vertex::Vector{Bool} # boolean indicating whether the most recent edge inserted a new vertex
  log_w::Vector{Float64} # unnormalized log-importance weight corresponding to edge_list proposals
  ancestor::Vector{Int64} # index of ancestor particle
  degrees::Vector{Float64} # degree vector
  vertex_map::Vector{Int64}  # vector where `v[j]` maps vertex `j` in the current graph to `v[j]` of the observed (final) graph
  vertex_unmap::Vector{Int64} # vector where v[j] maps vertex j in the final graph to v[j] in the current graph
  edge_queue::Vector{Bool} # vector of categories indicating which edges in the final graph are currently eligible to be inserted (=true)
  has_eigensystem::Vector{Bool} # indicates whether eigenvalue decomposition has been computed yet
  eig_vals::Vector{Float64} # eigenvalues of normalized Laplacian
  eig_vecs::Array{Float64,2} # eigenvectors of normalized Laplacian (transposed to optimize computations so that rows are eigenvectors)
  # eligible_vertices::Vector{S} # vector of vertex indices indicating which vertices in the final graph haven't been inserted yet
  # eigen_system::Base.LinAlg.Eigen  # eigen system of normalized Laplacian of graph

end

struct SamplerState
  # structure for keeping track of overall sampler quantities

  data_elist::Array{Int64,2} # observed graph, fixed, edge list indices
  # data_elist_vals::Vector{Int64} # observed graph, fixed, edge list values (e.g., 1 for a simple graph)
  degree_data::Vector{Int64}
  ne_data::Vector{Int64} # number of edges in `G_data`, fixed
  nv_data::Vector{Int64} # number of edges in `G_data`, fixed
  size_bias::Vector{Bool} # whether root vertex is sampled from size-biased distribution, fixed
  a_α::Vector{Float64} # hyper-parameter for alpha, fixed
  b_α::Vector{Float64} # hyper-parameter for alpha, fixed
  a_λ::Vector{Float64} # hyper-parameter for lambda, fixed
  b_λ::Vector{Float64} # hyper-parameter for lambda, fixed
  α::Vector{Float64} # current value of alpha, updates with sampler
  λ::Vector{Float64} # current value of lambda, updates with sampler
  B::Vector{Int64} # current vector of coin flips, updates with sampler
  K::Vector{Int64} # current vector of r.w. lengths, updates with sampler
  k_range::Vector{Int64} # where to truncate support of conditional distribution of K
  b_range::Vector{Int64} # {0,1}
  particle_path::Vector{Int64} # index of current graph sequence in particle_container

end

function new_sampler_state(elist_data::Array{Int64,2},size_bias::Bool,α::Float64,λ::Float64)
"""
Create new sampler state for SMC
"""
# not needed for SMC
pp = zeros(Int64,n_edges_data)
B_coins = rand(Bernoulli(α),n_edges_data)
K_walks = rand(Poisson(λ),n_edges_data)
B_coins[1] = zero(Int64)
K_walks[1] = zero(Int64)
k_trunc = 100*one(Int64)
a_α = 1.0
b_α = 1.0
a_λ = 1.0
b_λ = 1.0

return new_sampler_state(elist_data,sb,a_α,b_α,a_λ,b_λ,α,λ,B_coins,K_walks,k_trunc,pp)

end

function new_sampler_state(elist_data::Array{Int64,2},
                            size_bias::Bool,
                            a_α::Float64,b_α::Float64,
                            a_λ::Float64,b_λ::Float64,
                            α::Float64,λ::Float64,
                            B::Vector{Int64},K::Vector{Int64},k_trunc::Int64,
                            particle_path::Array{Int64,1})
"""
  Construct new sampler state object from given data and parameters
"""
  BB = copy(B)
  KK = copy(K)
  BB[1] == zero(Int64) ? nothing : BB[1] = zero(Int64)
  KK[1] == zero(Int64) ? nothing : KK[1] = zero(Int64)
  pp = copy(particle_path)
  return SamplerState(elist_data, # edge list
                      getDegrees(elist_data), # degrees
                      [size(elist_data,1)], # number of edges
                      [maximum(elist_data)], # number of vertices
                      [size_bias],
                      [a_α],[b_α],[a_λ],[b_λ], # hyper-parameters
                      [α],[λ], # parameters
                      BB, # B vector
                      KK, # K vector
                      collect(0:(k_trunc)), # k_range
                      [zero(Int64),one(Int64)], # b_range
                      pp) # particle_path

end

function initialize_blank_particle_state(n_edges::Int64,n_edges_max::Int64,deg_max::Int64,n_vertices_max::Int64)

  ps = ParticleState([n_edges],
                      zeros(Int64,n_edges,2), # edge_list
                      zeros(Int64,n_edges), # edge_idx_list
                      zeros(Int64,deg_max,n_vertices_max), # nbd_list
                      zeros(Int64,1), # n_vertices
                      [false], # new_vertex
                      zeros(Float64,1), # log_w
                      zeros(Int64,1), # ancestor
                      zeros(Int64,n_vertices_max), # degrees
                      zeros(Int64,n_vertices_max), # vertex_map
                      zeros(Int64,n_vertices_max), # vertex_unmap
                      falses(n_edges_max), # edge_queue
                      zeros(Bool,1), # has_eigensystem
                      zeros(Float64,n_edges+1), # eig_vals
                      zeros(Float64,n_edges+1,n_edges+1) # eig_vecs
                      )
                      # [spzeros(deg_max) for i=1:n_vertices_max], # nbd_list
                      # spzeros(Float64,n_vertices_max), # eig_vals
                      # spzeros(Float64,n_vertices_max,n_vertices_max) # eig_vecs
  return ps
end

function reset_particle_state!(particle::ParticleState)
"""
  Function to reset the fields of a particle state
"""
  # leave .n_edges unchanged
  particle.edge_list[:] = zero(Int64)
  particle.edge_idx_list[:] = zero(Int64)
  # for i = 1:length(particle.nbd_list)
  #   particle.nbd_list[i][:] = zero(Int64)
  #   dropzeros!(particle.nbd_list[i])
  # end
  particle.nbd_list[:] = zero(Int64)
  particle.n_vertices[:] = zero(Int64)
  particle.new_vertex[:] = false
  particle.log_w[:] = zero(Int64)
  particle.ancestor[:] = zero(Int64)
  particle.degrees[:] = zero(Int64)
  particle.vertex_map[:] = zero(Int64)
  particle.vertex_unmap[:] = zero(Int64)
  particle.edge_queue[:] = false
  particle.has_eigensystem[:] = false
  particle.eig_vals[:] = zero(Float64)
  # dropzeros!(particle.eig_vals)
  particle.eig_vecs[:] = zero(Float64)
  # dropzeros!(particle.eig_vecs)
end

function reset_particles!(particle_container::Array{Array{ParticleState,1},1})

  np = length(particle_container)
  ne = length(particle_container[1])

  for p = 1:np
    for t = 1:ne
      reset_particle_state!(particle_container[p][t])
    end
  end

end

function reset_particles!(particle_container::Array{Array{ParticleState,1},1},particle_path::Vector{Int64})

  np = length(particle_container)::Int64
  ne = length(particle_container[1])::Int64

  for p = 1:np
    for t = 1:ne
      p == particle_path[t] ? nothing : reset_particle_state!(particle_container[p][t])
    end
  end

end

function copyParticlelState!(particle_target::ParticleState,particle_source::ParticleState)

  particle_target.n_edges[:] = particle_source.n_edges[:]
  particle_target.edge_list[:] = particle_source.edge_list[:]
  particle_target.edge_idx_list[:] = particle_source.edge_idx_list[:]
  particle_target.nbd_list[:] = particle_source.nbd_list[:]
  particle_target.n_vertices[:] = particle_source.n_vertices[:]
  particle_target.new_vertex[:] = particle_source.new_vertex[:]
  particle_target.log_w[:] = particle_source.log_w[:]
  particle_target.ancestor[:] = particle_source.ancestor[:]
  particle_target.degrees[:] = particle_source.degrees[:]
  particle_target.vertex_map[:] = particle_source.vertex_map[:]
  particle_target.vertex_unmap[:] = particle_source.vertex_unmap[:]
  particle_target.edge_queue[:] = particle_source.edge_queue[:]
  particle_target.has_eigensystem[:] = particle_source.has_eigensystem[:]
  particle_target.eig_vals[:] = particle_source.eig_vals[:]
  particle_target.eig_vecs[:] = particle_source.eig_vecs[:]

end

function rw_csmc!(particle_container::Array{Array{ParticleState,1},1},
                  s_state::SamplerState,
                  n_particles::Int64,
                  L::Array{Float64,2},
                  W::Array{Float64,2},
                  eig_pgf::Array{Float64,1},
                  ancestors::Array{Int64,2},
                  ancestors_ed::Array{Int64,2},
                  unq::Array{Int64,2},
                  n_edges_in_queue::Array{Int64,1},
                  edge_samples_idx::Array{Int64,1},
                  edge_samples::Array{Int64,1},
                  log_w::Array{Float64,1},
                  free_edge_samples_idx::Array{Int64,1},
                  free_edge_samples::Array{Int64,1}
                )

  T = s_state.ne_data[1]::Int64
  nv_max = s_state.nv_data[1]::Int64

  # reset particle_container except for conditioned path
  # ***** might be unnecessary *****
  reset_particles!(particle_container,s_state.particle_path)
  resetArray!(L)
  resetArray!(W)
  resetArray!(eig_pgf)

  resetArray!(ancestors)
  resetArray!(ancestors_ed)
  resetArray!(unq)
  resetArray!(n_edges_in_queue)
  resetArray!(edge_samples_idx)
  resetArray!(edge_samples)
  resetArray!(log_w)
  resetArray!(free_edge_samples_idx)
  resetArray!(free_edge_samples)

    # pre-allocated arrays
  # ancestors = zeros(Int64,T,n_particles)
  # ancestors_ed = zeros(Int64,T,n_particles)
  #
  # n_edges_in_queue = zeros(Int64,n_particles)
  # edge_samples_idx = zeros(Int64,n_particles)
  # edge_samples = zeros(Int64,n_particles)
  # log_w = zeros(Float64,n_particles)
  # free_edge_samples_idx = zeros(Int64,n_particles-1)
  # free_edge_samples = zeros(Int64,n_particles-1)

  # sample t=1
  t = one(Int64)
  # freeParticles!(free_particles,t,T,s_state.particle_path)
  edge_idx = sample(1:T,n_particles-1)::Array{Int64,1}
  fillParticles!(edge_samples,edge_idx,t,s_state.particle_path,particle_container)

  log_weights = -log(n_particles)::Float64
  updateParticles!(particle_container,t,s_state,edge_samples,zeros(Int64,1,1),-log(T)*ones(Float64,n_particles))
  ancestors[1,:] .= 1:n_particles
  ancestors_ed[1,1:(s_state.particle_path[t]-1)] .= edge_idx[1:(s_state.particle_path[t]-1)]
  ancestors_ed[1,s_state.particle_path[t]] = particle_container[s_state.particle_path[t]][t].edge_idx_list[t]
  ancestors_ed[1,(s_state.particle_path[t]+1):end] .= edge_idx[s_state.particle_path[t]:end]

  uniqueParticles!(unq,t,edge_samples,ancestors)

  for t = 2:(T-2)

    # generate exhaustive proposal set
    for p = 1:n_particles
      n_edges_in_queue[p] = sum(particle_container[p][t-1].edge_queue)
    end
    cumu_eq = cumsum(n_edges_in_queue)::Array{Int64,1}
    total_eq = cumu_eq[end]::Int64
    edge_logp = zeros(Float64,total_eq)
    lse_w = zeros(Float64,total_eq)
    edge_proposals = zeros(Int64,total_eq)

    # freeParticles!(free_particles,t,n_particles,s_state.particle_path)
    for p = 1:n_particles
      p==1 ? (idx = 1:cumu_eq[p]) : (idx = (cumu_eq[p-1]+1):cumu_eq[p])
      dup = findfirst( (unq[t-1,p] .== unq[t-1,1:p]) )::Int64
      if dup==p
        edge_proposals[idx] = find(particle_container[p][t-1].edge_queue)
        edge_logp[idx] = generateProposalProbsRW!(L,W,eig_pgf,particle_container[p][t-1].n_vertices[1],particle_container[p][t-1],s_state)
        resetArray!(L)
        resetArray!(W)
        resetArray!(eig_pgf)
      else # same ancestor edge
        dup==1 ? (dup_idx = 1:cumu_eq[dup]) : (dup_idx = (cumu_eq[dup-1]+1):cumu_eq[dup])
        # println("t " * string(t) * ", " * string(p) * " idx " * string(idx) * ", " * string(dup) * " idx " * string(dup_idx))
        edge_proposals[idx] .= edge_proposals[dup_idx]
        edge_logp[idx] .= edge_logp[dup_idx]
      end
    end

    # sample edges for propogation
    logSumExpWeights!(lse_w,edge_logp)
    stratifiedResample!(free_edge_samples_idx,lse_w,n_particles-1) # index of edge in s_state.data_elist
    fillParticleIdx!(edge_samples_idx,free_edge_samples_idx,t,edge_proposals,cumu_eq,s_state.particle_path,particle_container) # fill in edge sample indices

    # free_edge_samples[:] = edge_proposals[free_edge_samples_idx]
    # fillParticles!(edge_samples,free_edge_samples,t,s_state.particle_path,particle_container) # fill in edge samples
    edge_samples[:] = edge_proposals[edge_samples_idx]
    log_w[:] = edge_logp[edge_samples_idx]

    # keep track of ancestors
    getAncestors!(ancestors,ancestors_ed,t,edge_samples_idx,cumu_eq)
    # println( "step " * string(t) * ", assigned ancestor is " * string(ancestors[t,s_state.particle_path[t]]) *
            # ", path ancestor is " * string(s_state.particle_path[t-1]) )
    assert( ancestors[t,s_state.particle_path[t]]==s_state.particle_path[t-1] )

    uniqueParticles!(unq,t,edge_samples,ancestors)

    # update particles
    updateParticles!(particle_container,t,s_state,edge_samples,ancestors,log_w)

  end

  t = (T-1)::Int64
  # for step T-1, need to look ahead to last step for weights
  # generate exhaustive proposal set
  for p = 1:n_particles
    n_edges_in_queue[p] = sum(particle_container[p][t-1].edge_queue)
  end
  cumu_eq = cumsum(n_edges_in_queue)
  total_eq = cumu_eq[end]
  assert(total_eq <= 2*n_particles)
  edge_logp = zeros(Float64,total_eq)
  # log_w = zeros(Float64,total_eq)
  lse_w = zeros(Float64,total_eq)
  edge_proposals = zeros(Int64,total_eq)

  last_edge_proposals = zeros(Int64,total_eq)
  new_vertex = zeros(Bool,2)
  temp_edge_list = zeros(Int64,t,2)
  # deg_max = convert(Int64,maximum(s_state.degrees))::Int64
  tmp_nbd_list = zeros(Int64,size(particle_container[1][1].nbd_list,1),nv_max)
  temp_edge_queue = zeros(Bool,T)
  temp_vertex_unmap = zeros(Int64,nv_max)
  edges_to_add = zeros(Bool,T)
  temp_degrees = zeros(Float64,nv_max)
  temp_nv = zeros(Int64,1)

  for p = 1:n_particles
    p==1 ? (idx = 1:cumu_eq[p]) : (idx = (cumu_eq[p-1]+1):cumu_eq[p])
    dup = findfirst( (ancestors_ed[t-1,p] .== ancestors_ed[t-1,:]) )::Int64
    if dup==p
      edge_proposals[idx] = find(particle_container[p][t-1].edge_queue)
      edge_logp[idx] = generateProposalProbsRW!(L,W,eig_pgf,particle_container[p][t-1].n_vertices[1],particle_container[p][t-1],s_state)
      resetArray!(L)
      resetArray!(W)
      resetArray!(eig_pgf)
    else # same ancestor edge
      dup==1 ? (dup_idx = 1:cumu_eq[dup]) : (dup_idx = (cumu_eq[dup-1]+1):cumu_eq[dup])
      edge_proposals[idx] .= edge_proposals[dup_idx]
      edge_logp[idx] .= edge_logp[dup_idx]
    end
    # look ahead to last step for proper weighting
    for i = 1:length(idx)
      temp_edge_list[1:(t-1),:] .= particle_container[p][t-1].edge_list
      temp_edge_list[t,:] .= s_state.data_elist[edge_proposals[idx[i]],:]

      temp_edge_queue .= particle_container[p][t-1].edge_queue
      temp_vertex_unmap .= particle_container[p][t-1].vertex_unmap
      new_vertex[:] = (particle_container[p][t-1].vertex_unmap[s_state.data_elist[[edge_proposals[idx[i]]],:]] .== 0)
      temp_nv[:] = particle_container[p][t-1].n_vertices[:] + sum(new_vertex)
      if sum(new_vertex) > 0
        vtx = s_state.data_elist[[edge_proposals[idx[i]]],find(new_vertex)][1]::Int64
        temp_vertex_unmap[vtx] = sum(temp_vertex_unmap .> 0) + 1
        edges_to_add[:] = [ in(vtx,s_state.data_elist[m,:]) for m=1:T ]
        temp_edge_queue[edges_to_add] = true
      end
      temp_edge_queue[edge_proposals[idx[i]]] = false
      temp_degrees .= particle_container[p][t-1].degrees
      temp_degrees[temp_vertex_unmap[s_state.data_elist[[edge_proposals[idx[i]]],1]]] += one(Float64)
      temp_degrees[temp_vertex_unmap[s_state.data_elist[[edge_proposals[idx[i]]],2]]] += one(Float64)

      tmp_nbd_list[:] = particle_container[p][t-1].nbd_list
      tmp_nbd_list[convert(Int64,temp_degrees[temp_vertex_unmap[temp_edge_list[t,1]]]),temp_vertex_unmap[temp_edge_list[t,1]]] = temp_edge_list[t,2]
      tmp_nbd_list[convert(Int64,temp_degrees[temp_vertex_unmap[temp_edge_list[t,2]]]),temp_vertex_unmap[temp_edge_list[t,2]]] = temp_edge_list[t,1]

      # sum(temp_edge_queue)==1 ? nothing : error("uh oh, " * sum(temp_edge_queue))
      assert(sum(temp_edge_queue)==1)
      edge_logp[idx[i]] += generateProposalProbsRW!(L,W,eig_pgf,temp_nv[1],temp_edge_list,tmp_nbd_list,temp_vertex_unmap,t,temp_degrees,temp_edge_queue,s_state)[1]::Float64
      resetArray!(L)
      resetArray!(W)
      eig_pgf[:] = zero(Float64)
      last_edge_proposals[idx[i]] = find(temp_edge_queue)[1]
    end

  end

    # sample edges for propogation
    logSumExpWeights!(lse_w,edge_logp)
    stratifiedResample!(free_edge_samples_idx,lse_w,n_particles-1) # index of edge in s_state.data_elist
    fillParticleIdx!(edge_samples_idx,free_edge_samples_idx,t,edge_proposals,cumu_eq,s_state.particle_path,particle_container) # fill in edge sample indices
    # free_edge_samples[:] = edge_proposals[free_edge_samples_idx]
    # fillParticles!(edge_samples,free_edge_samples,t,s_state.particle_path,particle_container) # fill in edge samples
    edge_samples[:] = edge_proposals[edge_samples_idx]

    log_w[:] = edge_logp[edge_samples_idx]

    # keep track of ancestors
    getAncestors!(ancestors,T-1,edge_samples_idx,cumu_eq)
    ancestors[T,:] = 1:n_particles # last step is deterministic

    # update particles for steps T-1 and T
    updateParticles!(particle_container,T-1,s_state,edge_samples,ancestors,log_w)
    edge_samples[:] = last_edge_proposals[edge_samples_idx]
    updateParticles!(particle_container,T,s_state,edge_samples,ancestors,zeros(Float64,n_particles))

end

function SMC(g::Array{Int64,2},n_particles::Int64,α::Float64,λ::Float64,sb::Bool)

  # initialize empty particle container
  n_edges_data = size(g,1)
  nv_max = maximum(g)::Int64
  deg_max = convert(Int64,maximum(getDegrees(g)))
  particle_container = Array{Array{ParticleState,1},1}(n_particles)
  for p in 1:n_particles
    particle_container[p] = [initialize_blank_particle_state(t,n_edges_data,deg_max,nv_max) for t in 1:n_edges_data]
  end

  # initialize SamplerState
  s_state = new_sampler_state(g,sb,α,λ)

  # run smc
  rw_smc!(particle_container,n_particles,s_state)
  log_p = marginalLogLikelihodEstimate(particle_container,n_edges_data)

  return particle_container,log_p

end

function ParticleGibbs(g::Array{Int64,2},n_particles::Int64,
                        a_α::Float64,b_α::Float64,a_λ::Float64,b_λ::Float64,
                        α_start::Float64,λ_start::Float64,sb::Bool,k_trunc::Int64,
                        n_mcmc_iter::Int64,n_burn::Int64,n_collect::Int64,n_print::Int64)

  # initialize empty particle container
  n_edges_data = size(g,1)
  nv_max = maximum(g)::Int64
  deg_max = convert(Int64,maximum(getDegrees(g)))
  particle_container = Array{Array{ParticleState,1},1}(n_particles)
  for p in 1:n_particles
    particle_container[p] = [initialize_blank_particle_state(t,n_edges_data,deg_max,nv_max) for t in 1:n_edges_data]
  end

  pp = zeros(Int64,n_edges_data)
  B_coins = rand(Bernoulli(α_start),n_edges_data)
  K_walks = rand(Poisson(λ_start),n_edges_data)
  B_coins[1] = zero(Int64)
  K_walks[1] = zero(Int64)

  # initialize SamplerState
  s_state = new_sampler_state(g,sb,a_α,b_α,a_λ,b_λ,α_start,λ_start,B_coins,K_walks,k_trunc,pp)

  # initialize sampler
  rw_smc!(particle_container,n_particles,s_state)
  p_idx = sample(1:n_particles,1)[1]
  getParticlePath!(s_state.particle_path,particle_container,p_idx)

  # pre-allocate arrays for computation in rw_csmc!
  L = zeros(Float64,nv_max,nv_max)
  W = zeros(Float64,nv_max,nv_max)
  eig_pgf = zeros(Float64,nv_max)
  ancestors = zeros(Int64,n_edges_data,n_particles)
  ancestors_ed = zeros(Int64,n_edges_data,n_particles)
  unq = zeros(Int64,n_edges_data,n_particles)
  n_edges_in_queue = zeros(Int64,n_particles)
  edge_samples_idx = zeros(Int64,n_particles)
  edge_samples = zeros(Int64,n_particles)
  log_w = zeros(Float64,n_particles)
  free_edge_samples_idx = zeros(Int64,n_particles-1)
  free_edge_samples = zeros(Int64,n_particles-1)
  lp_b = zeros(Float64,2)
  p_b = zeros(Float64,2)
  lp_k = zeros(Float64,k_trunc+1)
  p_k = zeros(Float64,k_trunc+1)

  # pre-allocate sample collection arrays
  n_samples = convert(Int64,floor((n_mcmc_iter - n_burn)/n_collect))
  alpha_samples = zeros(Float64,n_samples)
  lambda_samples = zeros(Float64,n_samples)
  particle_trajectory_samples = zeros(Int64,n_samples,n_edges_data)
  edge_sequence_samples = zeros(Int64,n_samples,n_edges_data)
  log_marginal_samples = zeros(Float64,n_samples)

  n_s = zero(Int64)
  t_elapsed = zero(Float64)
  tic();
  for s = 1:n_mcmc_iter
      # update edge sequence
      rw_csmc!(particle_container,s_state,n_particles,
              L,W,eig_pgf,
              ancestors,
              ancestors_ed,
              unq,
              n_edges_in_queue,
              edge_samples_idx,
              edge_samples,
              log_w,
              free_edge_samples_idx,
              free_edge_samples)

      p_idx = sample(1:n_particles,1)[1]
      getParticlePath!(s_state.particle_path,particle_container,p_idx)

      # update B, K, α, λ
      updateBandK!(s_state,particle_container,
                  L,eig_pgf,lp_b,p_b,lp_k,p_k)
      updateAlphaAndLambda!(s_state)

      if s > n_burn && mod((s - n_burn),n_collect)==0
          # collect samples
          n_s += one(Int64)
          log_marginal_samples[n_s] = marginalLogLikelihodEstimate(particle_container,n_edges_data)
          saveSamples!(alpha_samples,lambda_samples,particle_trajectory_samples,edge_sequence_samples,
                      s_state,particle_container,p_idx,n_s)

      end

      if mod(s,n_print)==0
          t_elapsed += toq();
          println( "Finished with " * string(s) * " / " * string(n_mcmc_iter) * " iterations. Elapsed time: $(t_elapsed) sec." )
          tic();
      end

  end

  samples = Dict("alpha" => alpha_samples,
                  "lambda" => lambda_samples,
                  "edge_sequence" => edge_sequence_samples)

  return samples
end

function rw_smc!(particle_container::Array{Array{ParticleState,1},1},n_particles::Int64,s_state::SamplerState)

  T = s_state.ne_data[1]::Int64
  nv_max = s_state.nv_data[1]::Int64

  # reset particle_container
  # ***** might be unnecessary *****
  reset_particles!(particle_container)

  # pre-allocate arrays for edge probability computation
  L = zeros(Float64,nv_max,nv_max)
  W = zeros(Float64,nv_max,nv_max)
  eig_pgf = zeros(Float64,nv_max)

  # sample t=1
  t = one(Int64)
  edge_idx = sample(1:T,n_particles)::Array{Int64,1}
  log_weights = -log(n_particles)::Float64
  updateParticles!(particle_container,t,s_state,edge_idx,zeros(Int64,1,1),-log(T)*ones(Float64,n_particles))

  ancestors = zeros(Int64,T,n_particles)
  ancestors[1,:] .= 1:n_particles
  ancestors_ed = zeros(Int64,T,n_particles)
  ancestors_ed[1,:] .= edge_idx

  unq = zeros(Int64,T,n_particles)
  # unq_pos = zeros(Bool,T,n_particles)
  uniqueParticles!(unq,t,edge_idx,ancestors)

  # pre-allocated arrays
  n_edges_in_queue = zeros(Int64,n_particles)
  edge_samples_idx = zeros(Int64,n_particles)
  edge_samples = zeros(Int64,n_particles)
  log_w = zeros(Float64,n_particles)

  for t in 2:(T-2)

    # generate exhaustive proposal set
    for p = 1:n_particles
      n_edges_in_queue[p] = sum(particle_container[p][t-1].edge_queue)
    end
    cumu_eq = cumsum(n_edges_in_queue)::Array{Int64,1}
    total_eq = cumu_eq[end]::Int64
    edge_logp = zeros(Float64,total_eq)
    lse_w = zeros(Float64,total_eq)
    edge_proposals = zeros(Int64,total_eq)
    # println(total_eq)

    for p = 1:n_particles
      p==1 ? (idx = 1:cumu_eq[p]) : (idx = (cumu_eq[p-1]+1):cumu_eq[p])
      dup = findfirst( (unq[t-1,p] .== unq[t-1,1:p]) )::Int64
      if dup==p
        edge_proposals[idx] = find(particle_container[p][t-1].edge_queue)
        edge_logp[idx] = generateProposalProbsRW!(L,W,eig_pgf,particle_container[p][t-1].n_vertices[1],particle_container[p][t-1],s_state)
        resetArray!(L)
        resetArray!(W)
        resetArray!(eig_pgf)
      else # same ancestor edge
        dup==1 ? (dup_idx = 1:cumu_eq[dup]) : (dup_idx = (cumu_eq[dup-1]+1):cumu_eq[dup])
        edge_proposals[idx] .= edge_proposals[dup_idx]
        edge_logp[idx] .= edge_logp[dup_idx]
      end
    end

    # down-sample edges
    logSumExpWeights!(lse_w,edge_logp)
    stratifiedResample!(edge_samples_idx,lse_w,n_particles)
    edge_samples[:] = edge_proposals[edge_samples_idx] # index of edge in s_state.data_elist

    log_w[:] = edge_logp[edge_samples_idx]

    # keep track of ancestors
    # ancestors[t,:] = getAncestors(edge_samples_idx,cumu_eq)
    getAncestors!(ancestors,ancestors_ed,t,edge_samples_idx,cumu_eq)

    # keep track of unique particles
    uniqueParticles!(unq,t,edge_samples,ancestors)

    # update particles
    updateParticles!(particle_container,t,s_state,edge_samples,ancestors,log_w)

  end

  t = (T-1)::Int64
  # for step T-1, need to look ahead to last step for weights
  # generate exhaustive proposal set
  # n_edges_in_queue = zeros(Int64,n_particles)
  for p = 1:n_particles
    n_edges_in_queue[p] = sum(particle_container[p][t-1].edge_queue)
  end
  cumu_eq = cumsum(n_edges_in_queue)
  total_eq = cumu_eq[end]
  assert(total_eq <= 2*n_particles)
  edge_logp = zeros(Float64,total_eq)
  # log_w = zeros(Float64,total_eq)
  lse_w = zeros(Float64,total_eq)
  edge_proposals = zeros(Int64,total_eq)

  last_edge_proposals = zeros(Int64,total_eq)
  new_vertex = zeros(Bool,2)
  temp_edge_list = zeros(Int64,t,2)
  # deg_max = convert(Int64,maximum(s_state.degrees))::Int64
  tmp_nbd_list = zeros(Int64,size(particle_container[1][1].nbd_list,1),nv_max)
  temp_edge_queue = zeros(Bool,T)
  temp_vertex_unmap = zeros(Int64,nv_max)
  edges_to_add = zeros(Bool,T)
  temp_degrees = zeros(Float64,nv_max)
  temp_nv = zeros(Int64,1)

  for p = 1:n_particles
    p==1 ? (idx = 1:cumu_eq[p]) : (idx = (cumu_eq[p-1]+1):cumu_eq[p])
    edge_proposals[idx] = find(particle_container[p][t-1].edge_queue)
    edge_logp[idx] = generateProposalProbsRW!(L,W,eig_pgf,particle_container[p][t-1].n_vertices[1],particle_container[p][t-1],s_state)
    resetArray!(L)
    resetArray!(W)
    eig_pgf[:] = zero(Float64)
    # look ahead to last step for proper weighting
    for i = 1:length(idx)
      temp_edge_list[1:(t-1),:] .= particle_container[p][t-1].edge_list
      temp_edge_list[t,:] .= s_state.data_elist[edge_proposals[idx[i]],:]

      temp_edge_queue .= particle_container[p][t-1].edge_queue
      temp_vertex_unmap .= particle_container[p][t-1].vertex_unmap
      new_vertex[:] = (particle_container[p][t-1].vertex_unmap[s_state.data_elist[[edge_proposals[idx[i]]],:]] .== 0)
      temp_nv[:] = particle_container[p][t-1].n_vertices[:] + sum(new_vertex)
      if sum(new_vertex) > 0
        vtx = s_state.data_elist[[edge_proposals[idx[i]]],find(new_vertex)][1]::Int64
        temp_vertex_unmap[vtx] = sum(temp_vertex_unmap .> 0) + 1
        edges_to_add[:] = [ in(vtx,s_state.data_elist[m,:]) for m=1:T ]
        temp_edge_queue[edges_to_add] = true
      end
      temp_edge_queue[edge_proposals[idx[i]]] = false
      temp_degrees .= particle_container[p][t-1].degrees
      temp_degrees[temp_vertex_unmap[s_state.data_elist[[edge_proposals[idx[i]]],1]]] += one(Float64)
      temp_degrees[temp_vertex_unmap[s_state.data_elist[[edge_proposals[idx[i]]],2]]] += one(Float64)

      tmp_nbd_list[:] = particle_container[p][t-1].nbd_list[:]
      tmp_nbd_list[convert(Int64,temp_degrees[temp_vertex_unmap[temp_edge_list[t,1]]]),temp_vertex_unmap[temp_edge_list[t,1]]] = temp_edge_list[t,2]
      tmp_nbd_list[convert(Int64,temp_degrees[temp_vertex_unmap[temp_edge_list[t,2]]]),temp_vertex_unmap[temp_edge_list[t,2]]] = temp_edge_list[t,1]

      # sum(temp_edge_queue)==1 ? nothing : error("uh oh, " * sum(temp_edge_queue))
      assert(sum(temp_edge_queue)==1)
      edge_logp[idx[i]] += generateProposalProbsRW!(L,W,eig_pgf,temp_nv[1],temp_edge_list,tmp_nbd_list,temp_vertex_unmap,t,temp_degrees,temp_edge_queue,s_state)[1]
      resetArray!(L)
      resetArray!(W)
      eig_pgf[:] = zero(Float64)
      last_edge_proposals[idx[i]] = find(temp_edge_queue)[1]
    end

  end

  # sample edges for propogation
  logSumExpWeights!(lse_w,edge_logp)
  stratifiedResample!(edge_samples_idx,lse_w,n_particles) # index of edge in s_state.data_elist
  edge_samples[:] = edge_proposals[edge_samples_idx]

  log_w[:] = edge_logp[edge_samples_idx]

  # keep track of ancestors
  getAncestors!(ancestors,T-1,edge_samples_idx,cumu_eq)
  ancestors[T,:] = 1:n_particles # last step is deterministic

  # update particles for steps T-1 and T
  updateParticles!(particle_container,T-1,s_state,edge_samples,ancestors,log_w)
  edge_samples[:] = last_edge_proposals[edge_samples_idx]
  updateParticles!(particle_container,T,s_state,edge_samples,ancestors,zeros(Float64,n_particles))

end

function getAncestors(es::Array{Int64,1},ceq::Array{Int64,1})::Array{Int64,1}
  anc = zeros(Int64,length(es))
  for i = 1:length(es)
    anc[i] = findfirst( ceq .>= es[i] )
  end
  return anc
end

function getAncestors!(ancestors::Array{Int64,2},t::Int64,es::Array{Int64,1},ceq::Array{Int64,1})

  for i = 1:length(es)
    ancestors[t,i] = findfirst( ceq .>= es[i] )
  end

end

function getAncestors!(ancestors::Array{Int64,2},ancestors_ed::Array{Int64,2},t::Int64,es::Array{Int64,1},ceq::Array{Int64,1})

  for i = 1:length(es)
    ancestors[t,i] = findfirst( ceq .>= es[i] )
    # ancestors[t,i]==1 ? st_ed = zero(Int64) : st_ed = ceq[ancestors[t,i]-1]
    ancestors_ed[t,i] = es[i]
  end

end

function getAncestors!(ancestors::Array{Int64,2},ancestors_ed::Array{Int64,2},t::Int64,es_idx::Array{Int64,1},ceq::Array{Int64,1},
                        ed_p::Array{Int64,1},particle_path::Array{Int64,1},particle_container::Array{Array{ParticleState,1},1})

  for i = 1:(particle_path[t]-1)
    ancestors[t,i] = findfirst( ceq .>= es_idx[i])
    ancestors[t,i]==1 ? st_ed = zero(Int64) : st_ed = ceq[ancestors[t,i]-1]
    ancestors_ed[t,i] = es_idx[i] - st_ed
  end

  i = particle_path[t]::Int64
  ancestors[t,i] = particle_path[t-1]
  ancestors[t,i]==1 ? st_ed = zero(Int64) : st_ed = ceq[ancestors[t,i]-1]
  ancestors_ed[t,i] = es_idx[i] - st_ed

  for i = (particle_path[t]+1):length(es)
    ancestors[t,i] = findfirst( ceq .>= es_idx[i])
    ancestors[t,i]==1 ? st_ed = zero(Int64) : st_ed = ceq[ancestors[t,i]-1]
    ancestors_ed[t,i] = es_idx[i] - st_ed
  end

end

function updateParticles!(particle_container::Array{Array{ParticleState,1},1},t::Int64,s_state::SamplerState,
            sampled_edges::Vector{Int64},ancestors::Array{Int64,2},log_w::Vector{Float64})
"""
  Updates particle_container at step t to append newly sampled edges
"""
  T = length(particle_container[1])::Int64

  if t==1
    for p = 1:length(particle_container)
      particle_container[p][1].edge_list[1,:] = s_state.data_elist[sampled_edges[p],:] # add new edge
      particle_container[p][1].edge_idx_list[1] = sampled_edges[p]
      particle_container[p][1].nbd_list[1,1] = 2 # neighborhood of vertex 1
      particle_container[p][1].nbd_list[1,2] = 1 # neighborhood of vertex 1
      particle_container[p][1].n_vertices[:] = 2 # update number of vertices
      particle_container[p][1].new_vertex[:] = true # update new_vertex flag
      particle_container[p][1].log_w[:] = log_w[p]
      particle_container[p][1].degrees[1:2] = one(Float64) # update degrees
      particle_container[p][1].vertex_map[1:2] = s_state.data_elist[sampled_edges[p],:] # update vertex map
      particle_container[p][1].vertex_unmap[s_state.data_elist[sampled_edges[p],:]] = [1,2] # update vertex unmap
      # construct initial edge queue
      eq = [ ( in(s_state.data_elist[m,1],s_state.data_elist[sampled_edges[p],:]) || in(s_state.data_elist[m,2],s_state.data_elist[sampled_edges[p],:]) ) for m in 1:T ]::Vector{Bool}
      eq[sampled_edges[p]] = false
      particle_container[p][1].edge_queue .= eq
      particle_container[p][1].ancestor[:] = p # ancestor for book-keeping
    end

  else
    new_vertex = falses(2)
    edge = zeros(Int64,2)
    edge_unmap = zeros(Int64,2)
    edges_to_add = zeros(Bool,T)
    for p = 1:n_particles
      edge .= s_state.data_elist[sampled_edges[p],:]
      # update edge list
      particle_container[p][t].edge_list[1:(t-1),:] .= particle_container[ancestors[t,p]][t-1].edge_list # copy previous edges
      particle_container[p][t].edge_list[t,:] .= edge # add new edge

      # update edge index list
      particle_container[p][t].edge_idx_list[1:(t-1)] = particle_container[ancestors[t,p]][t-1].edge_idx_list
      particle_container[p][t].edge_idx_list[t] = sampled_edges[p]

      # update vertex_map
      particle_container[p][t].vertex_map .= particle_container[ancestors[t,p]][t-1].vertex_map # copy previous vertex_map
      new_vertex[:] = [ particle_container[ancestors[t,p]][t-1].vertex_unmap[edge[1]]==0;
                      particle_container[ancestors[t,p]][t-1].vertex_unmap[edge[2]]==0 ]
      # new_vertex = [!in(s_state.data_elist[sampled_edges[p],1],particle_container[p][t].vertex_map);
      #                 !in(s_state.data_elist[sampled_edges[p],2],particle_container[p][t].vertex_map) )]
      assert(0 <= sum(new_vertex) <= 1)
      particle_container[p][t].new_vertex[:] = (new_vertex[1] || new_vertex[2])
      particle_container[p][t].n_vertices[:] = sum(particle_container[p][t].vertex_map .> 0) + sum(new_vertex)
      any(new_vertex) ? particle_container[p][t].vertex_map[particle_container[p][t].n_vertices] = edge[find(new_vertex)] : nothing

      # update log_w
      particle_container[p][t].log_w[:] = log_w[p]

      # update ancestor
      particle_container[p][t].ancestor[:] = ancestors[t,p]

      # update vertex_unmap
      particle_container[p][t].vertex_unmap .= particle_container[ancestors[t,p]][t-1].vertex_unmap
      any(new_vertex) ? new_v = find(new_vertex)[1]::Int64 : nothing
      any(new_vertex) ? particle_container[p][t].vertex_unmap[particle_container[p][t].edge_list[t,new_v]] = particle_container[p][t].n_vertices[1] : nothing

      # update degrees
      edge_unmap .= particle_container[p][t].vertex_unmap[edge]
      particle_container[p][t].degrees .= particle_container[ancestors[t,p]][t-1].degrees
      particle_container[p][t].degrees[edge_unmap[1]] += one(Float64)
      particle_container[p][t].degrees[edge_unmap[2]] += one(Float64)

      # update neighborhoods
      particle_container[p][t].nbd_list[:] = particle_container[ancestors[t,p]][t-1].nbd_list[:]
      particle_container[p][t].nbd_list[convert(Int64,particle_container[p][t].degrees[edge_unmap[1]]),edge_unmap[1]] = edge_unmap[2]
      particle_container[p][t].nbd_list[convert(Int64,particle_container[p][t].degrees[edge_unmap[2]]),edge_unmap[2]] = edge_unmap[1]

      # update edge queue
      particle_container[p][t].edge_queue .= particle_container[ancestors[t,p]][t-1].edge_queue
      if new_vertex[1] || new_vertex[2] # add new edges to the queue
        # mark edges
        vtx = s_state.data_elist[sampled_edges[p],new_v][1]::Int64
        edges_to_add[:] = [ in(vtx,s_state.data_elist[m,:]) for m=1:T ]::Vector{Bool}
        particle_container[p][t].edge_queue[edges_to_add] = true

      end
      particle_container[p][t].edge_queue[sampled_edges[p]] = false # remove new edge from queue

    end

  end

end

function generateProposalProbsRW!(L::Array{Float64,2},W::Array{Float64,2},eig_pgf::Array{Float64,1},nv::Int64,p_state::ParticleState,s_state::SamplerState)::Array{Float64,1}

  !p_state.has_eigensystem[1] ? updateEigenSystem!(L,nv,p_state) : nothing

  return generateProposalProbsRW!(L,W,eig_pgf,p_state.eig_vals,p_state.eig_vecs,nv,p_state.edge_list,p_state.nbd_list,p_state.vertex_unmap,p_state.n_edges[1],p_state.degrees[1:p_state.n_vertices[1]],p_state.edge_queue,s_state)
end

function generateProposalProbsRW!(L::Array{Float64,2},W::Array{Float64,2},eig_pgf::Array{Float64,1},
        nv::Int64,edge_list::Array{Int64,2},nbd_list::Union{Array{Array{Int64,1},1},Array{Int64,2},Array{SparseVector{Int64,Int64},1}},vertex_unmap::Vector{Int64},
        n_edges::Int64,degrees::Vector{Float64},
        edge_queue::Vector{Bool},s_state::SamplerState)::Array{Float64,1}

    # need to generate eigensystem
    denseNormalizedLaplacian!(L,vertex_unmap[edge_list],degrees,nv)
    eig_vals,eig_vecs = generateEigenSystem!(L,nv)

    return generateProposalProbsRW!(L,W,eig_pgf,eig_vals,eig_vecs,nv,edge_list,nbd_list,vertex_unmap,n_edges,degrees,edge_queue,s_state)

end

function generateProposalProbsRW!(L::Array{Float64},W::Array{Float64,2},eig_pgf::Array{Float64,1},
        eig_vals::Union{Array{Float64,1},SparseVector{Float64,Int64}},
        eig_vecs::Union{Array{Float64,2},SparseMatrixCSC{Float64,Int64}},
        nv::Int64,edge_list::Array{Int64,2},nbd_list::Union{Array{Array{Int64,1},1},Array{Int64,2},Array{SparseVector{Int64,Int64},1}},vertex_unmap::Vector{Int64},
        n_edges::Int64,degrees::Vector{Float64},
        edge_queue::Vector{Bool},s_state::SamplerState)::Array{Float64,1}
"""
  Returns the log probability under the prior of the entire set of eligible edges,
  as indexed in the enumeration of data edges. Operates on L and W in-place.
"""

  t_step = (n_edges + 1)::Int64
  t_remaining = (s_state.ne_data[1] - n_edges)::Int64

  # calculate predictive probabiliy of a new vertex
  P_I_1 = ((s_state.a_α[1] + sum(s_state.B) - s_state.B[t_step])/(s_state.a_α[1] + s_state.b_α[1] + s_state.ne_data[1] - 1))::Float64

  # calculate parameters of predictive distribution for random walk length
  a_lambdaPrime = (s_state.a_λ[1] + sum(s_state.K) - s_state.K[t_step])::Float64
  b_lambdaPrime = (1/(s_state.b_λ[1] + s_state.ne_data[1]))::Float64

  # calculate r.w. probs (negative binomial predictive distribution)
  # eig_pgf = zeros(Float64,nv)
  nbPred!(eig_pgf,nv,a_lambdaPrime,b_lambdaPrime,eig_vals)
  # eig_pgf[:] = (((1.0 - b_lambdaPrime)^(a_lambdaPrime)).*( (1.0 .- esys[:values]).*((1.0 .- b_lambdaPrime.*(1.0 .- esys[:values])).^(-a_lambdaPrime)) ))::Array{Float64,1}
  # any(isinf.(eig_pgf)) ? eig_pgf[ isinf.(eig_pgf) ] = zero(Float64) : nothing # moore-penrose pseudoinverse
  # randomWalkProbs!(W,nv,degrees,eig_pgf,eig_vecs)
  # W[1:nv,1:nv] = Diagonal(degrees.^(-1/2)) * esys[:vectors] * Diagonal(eig_pgf) * esys[:vectors]' * Diagonal(degrees.^(1/2))
  # rwProbs = diag(degrees.^(-1/2)) * U * diag(eigenValues) * U' * diag(degrees.^(1/2));

  # account for initial vertex choice
  # if s_state.size_bias[1] # size-biased selection
  #   for m = 1:nv
  #     s_state.size_bias[1] ? W[1:nv,m] .= W[1:nv,m].*(degrees[1:nv]./(2*n_edges)) : nothing
  #   end
  # else # uniform selection
  #   W[1:nv,1:nv] .= W[1:nv,1:nv]./nv
  # end
  # elMax!(W[1:nv,1:nv], 0.0) # for numerical stability
  # elMin!(W[1:nv,1:nv], 1.0) # for numerical stability
  # W[1:nv,1:nv] .= W[1:nv,1:nv]./sum(W[1:nv,1:nv]) # for numerical stability

  # calculate log-probabilities
  edge_proposal_idx = find(edge_queue)::Array{Int64,1}
  n_eligible = sum(edge_queue)::Int64
  log_p = zeros(Float64,n_eligible)
  ed = zeros(Int64,1,2)
  p_tmp = zero(Float64)
  log_labelprob = log(s_state.nv_data[1] - nv)::Float64
  w = zeros(Float64,1,2)
  frwp = zeros(Float64,1)
  nbd = zeros(Int64,size(nbd_list,1))
  rt = zeros(Int64,1,2)

  if isempty(nbd_list)
    edgelist2adj!(L,vertex_unmap[edge_list]) # L is used for adjacency matrix
  end
  for n = 1:n_eligible
    resetArray!(nbd)
    resetArray!(ed)
    resetArray!(w)
    resetArray!(rt)
    ed[:] = vertex_unmap[s_state.data_elist[edge_proposal_idx[n],:]]
    if ed[1] > 0 && ed[2] > 0 # vertices not yet inserted have index 0 in current particle state
      randomWalkProbs!(w,ed,nv,degrees,eig_pgf,eig_vecs)
      s_state.size_bias[1] ? (w .*= (degrees[ed]./(2*n_edges))) : (w[:] *= 1/nv)
      log_p[n] = (log(1 - P_I_1) + log(sum(w)))::Float64
    else
      root_vtx = findfirst(ed .> 0)::Int64
      rt[:] = root_vtx
      # nbd[1] = root_vtx
      if(!isempty(nbd_list))
        nbd[:] = nbd_list[:,root_vtx]
        # nbd[:] = [root_vtx; nbd_list[:,root_vtx]]
      else
        nb = find( L[:,root_vtx] .> 0 ) # here, L is the adjacency matrix
        nbd[1:length(nb)] = nb
      end
      randomWalkProbs!(w,rt,nv,degrees,eig_pgf,eig_vecs) # prob of ending where it starts
      frwp[1] = (randomWalkProbs(root_vtx,nbd,nv,degrees,eig_pgf,eig_vecs) + w[1])::Float64 # prob of fruitless r.w.
      s_state.size_bias[1] ? (frwp[:] *= degrees[root_vtx]/(2*n_edges)) : (frwp[:] *= 1/nv)
      # isempty(nbd_list) ? p_tmp = fruitlessRWProb(W,L,root_vtx,nv)::Float64 : p_tmp = fruitlessRWProb(W,nbd_list,root_vtx)::Float64
      log_p[n] = (log(P_I_1 + (1 - P_I_1)*frwp[1] ) - log_labelprob)::Float64 # last term accounts for the "random permutation" applied to obtain the observed labels
    end
  end

  return log_p
end

function marginalLogLikelihodEstimate(particle_container::Array{Array{ParticleState,1},1},T::Int64)
"""
  Returns log of (unbiased) marginal likelihood estimate up to and including
  step T of input particle set
"""
  n_particles = length(particle_container)
  log_p = -T*log(n_particles)
  for t = 1:T
    tmp_w = zero(Float64)
    for p = 1:n_particles
      tmp_w += exp(particle_container[p][t].log_w[1])
      # println(exp(particle_container[p][t].log_w[1]))
    end
    log_p += log(tmp_w)
  end
  return log_p
end

function getParticlePath!(particle_path::Array{Int64,1},particle_container::Array{Array{ParticleState,1},1},p_idx::Int64)

  T = length(particle_path)
  particle_path[T] = p_idx
  for t in (T-1):-1:1
    particle_path[t] = particle_container[particle_path[t+1]][t+1].ancestor[1]
  end

end

function fillParticles!(edge_samples::Array{Int64,1},free_edge_samples::Array{Int64,1},t::Int64,particle_path::Array{Int64,1},particle_container::Array{Array{ParticleState,1},1})

  for i = 1:(particle_path[t]-1)
    edge_samples[i] = free_edge_samples[i]
  end
  edge_samples[particle_path[t]] = particle_container[particle_path[t]][t].edge_idx_list[t]
  for i = (particle_path[t]+1):length(edge_samples)
    edge_samples[i] = free_edge_samples[i-1]
  end

end

function fillParticleIdx!(es_idx::Array{Int64,1},free_es_idx::Array{Int64,1}, t::Int64,edge_proposals::Array{Int64,1},cumu_eq::Array{Int64,1}, particle_path::Array{Int64,1},particle_container::Array{Array{ParticleState,1},1})

  for i = 1:(particle_path[t]-1)
    es_idx[i] = free_es_idx[i]
  end

  p_anc = particle_path[t-1]::Int64
  p_anc==1 ? st_ed = one(Int64) : st_ed = (cumu_eq[p_anc-1])::Int64
  es_idx[particle_path[t]] = st_ed + findfirst( particle_container[particle_path[t]][t].edge_idx_list[t] .== edge_proposals[(st_ed+1):cumu_eq[p_anc]])

  for i = (particle_path[t]+1):length(es_idx)
    es_idx[i] = free_es_idx[i-1]
  end

end

function freeParticles!(free_particles::Array{Int64,1},t::Int64,T::Int64,particle_path::Array{Int64,1})

  tt = zero(Int64)
  for i = 1:(particle_path[t]-1)
    free_particles[i] = i
  end
  for i = particle_path[t]:(T-1)
    free_particles[i] = i+1
  end

end