text:
Tags_Default = Dict(

# ==============================================================================
# SHADOW DOM stuff
# ==============================================================================
# h-scroll v-scroll: scrollbar-track scrollbar-thumb
"h-scroll" => Dict( "offset" => 0,
                    "color" => "darkgrey",
										"width" =>  12,
										"height" =>  12,
										"opacity" => 0.3,
                    "hover" => Dict("opacity" => 0.9 ),
										"nodes" => [ Dict(">" => "scrollbar-thumb", "color" => "grey", "height" =>  12,"width" =>  10 ) ]
									),
"v-scroll" => Dict( "offset" => 0,
                    "color" => "darkgrey",
										"width" =>  12,
										"height" =>  12,
										"opacity" => 0.3,
                    "hover" => Dict("opacity" => 0.9 ),
										"nodes" => [ Dict(">" => "scrollbar-thumb", "color" => "grey", "height" =>  12,"width" =>  10 ) ]
									),

"page" => Dict("color" => "white"),
"circle" => Dict("color" => "white"),
"window" => Dict("color" => "white"),
# ==============================================================================
# NORMAL DOM stuff
# ==============================================================================
# 	        Defines a hyperlink
"a" => Dict( "display" => "inline", "margin" => [0,0,4,4],
             "hover"   => Dict( "text" => Dict("color" => "blue" )),
			 "font"    => Dict( "size" => 13, "family" => "Sans", "weight" => "bold",
							 "color" => [0.23,0.2695,0.6757] )
	),
# 	    Defines an abbreviation or an acronym
"abbr" => Dict("display" => "inline", "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
# <address>	    Defines contact information for the author/owner of a document
"address" => Dict( "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
# <area>	    Defines an area inside an image-map
"area" => Dict( "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
#
	    Defines an article
"article" => Dict( "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
# 	    Defines content aside from the page content
"aside" => Dict( "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
# <audio>	    Defines sound content
"audio" => Dict( "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
# 	        Defines bold text
"b" => Dict("display" => "inline",  "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
# <base>	    Specifies the base URL/target for all relative URLs in a document
#"base" => Dict( "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
# <bdi>	        Isolates a part of text that might be formatted in a different direction from other text outside it
"bdi" => Dict( "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
# <bdo>	        Overrides the current text direction
"bdo" => Dict( "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
#
	Defines a section that is quoted from another source
"blockquote" => Dict("display" => "inline",  "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
# 	    Defines the document's body
"body" => Dict("display" => "block",  "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
#
	        Defines a single line break
"br" => Dict( "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
# <button>	    Defines a clickable button
"button" => Dict( "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
# <canvas>	    Used to draw graphics, on the fly, via scripting (usually JavaScript)
"canvas" => Dict( "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
#
	    Defines a table caption
"caption" => Dict( "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
# 	    Defines the title of a work
"cite" => Dict( "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
# 	    Defines a piece of computer code
"code" => Dict( "font" =>  Dict( "size" => 13, "family" => "Courier New", "color" => "black" )),
# <col>	        Specifies column properties for each column within a # <colgroup> element
"col" => Dict( "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
# <colgroup>	Specifies a group of one or more columns in a table f"" => Dict( "font" =>  ""),or formatting
"colgroup" => Dict( "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
# <datalist>	Specifies a list of pre-defined options for input controls
"datalist" => Dict( "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
#
	        Defines a description/value of a term in a description list
"dd" => Dict( "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
# 	        Defines text that has been deleted from a document
"del" => Dict( "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
# 	    Defines additional details that the user can view or hide
"details" => Dict( "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
# <dfn>	        Represents the defining instance of a term
"dfn" => Dict( "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
# <dialog>	    Defines a dialog box or window
"dialog" => Dict( "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
#
	        Defines a section in a document
"div" => Dict(  "display" => "block", "margin" => [0,0,0,0] ),
#
	        Defines a description list
"dl" => Dict( "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
#
	        Defines a term/name in a description list
"dt" => Dict( "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
# 	        Defines emphasized text
"em" => Dict( "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
#
	Defines a caption for a #
 element
"figcaption" => Dict( "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
#
	    Specifies self-contained... content
"figure" => Dict( "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
# 	    Defines a footer for a document or section
"footer" => Dict( "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
# <form>	    Defines an HTML form for user input
"form" => Dict( "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
#
 to #
	Defines HTML headings
"h1" => Dict( "display" => "block", "margin" => [0,0,9,9],
			  "font" => Dict( "size" => 32, "family" => "Sans",
			  "weight" => "bold", "color" => "black" )
							),
"h2" => Dict( "display" => "block", "margin" => [0,0,7,7],
			  "font" => Dict( "size" => 24, "family" => "Sans",
			  "weight" => "bold", "color" => "black" ) ),
"h3" => Dict( "display" => "block", "margin" => [0,0,6,6],
			  "font" => Dict( "size" => 18, "family" => "Sans",
			  "weight" => "bold", "color" => "black" ) ),
"h4" => Dict( "display" => "block", "margin" => [0,0,5,5],
			  "font" => Dict( "size" => 16, "family" => "Sans",
			  "weight" => "bold", "color" => "black" ) ),
"h5" => Dict( "display" => "block", "margin" => [0,0,4,4],
			  "font" => Dict( "size" => 13, "family" => "Sans",
			  "weight" => "bold", "color" => "black" ) ),
"h6" => Dict( "display" => "block", "margin" => [0,0,3,3],
			  "font" => Dict( "size" => 10, "family" => "Sans",
			  "weight" => "bold", "color" => "black" ) ),
# 	    Defines a header for a document or section
"header" => Dict( "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
#
	        Defines a thematic change in the content
"hr" => Dict( "display" => "block", "margin" => [20,20,20,20], "height" => 1, "color" => "darkgrey" ),
# 	        Defines a part of text in an alternate voice or mood
"i" => Dict( "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
# 	    Defines an inline frame
"iframe" => Dict( "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
# 	        Defines an image
"img" => Dict( "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
# 	    Defines a label for an # <input> element
"label" => Dict( "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
# <legend>	    Defines a caption for a <fieldset> element
"legend" => Dict( "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
#
	        Defines a list item"" => Dict( "font" =>  ""),
"li" => Dict( "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
# 	    Specifies the main content of a document
"main" => Dict( "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
# <map>	        Defines a client-side image-map
"map" => Dict( "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
# 	    Defines marked/highlighted text
"mark" => Dict( "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
# <menu>	    Defines a list/menu of commands
"menu" => Dict( "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
# <menuitem>	Defines a command/menu item that the user can invoke from a popup menu
"menuitem" => Dict( "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
# 	        Defines navigation links
"nav" => Dict( "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
#
	        Defines an ordered list
"ol" => Dict( "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
# <optgroup>	Defines a group of related options in a drop-down list
"optgroup" => Dict( "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
# <option>	    Defines an option i""n a drop-down list
"option" => Dict( "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
#
	        Defines a paragraph
"p" => Dict( "display" => "block", "margin" => [0, 16, 0, 16], "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
#
	        Defines preformatted text
"pre" => Dict( "display"=>"block", "font" =>  Dict( "size" => 13, "family" => "Courier New", "color" => "black" )),
# <progress>	Represents the progress of a task
"progress" => Dict( "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
# <q>	        Defines a short quotation
"q" => Dict( "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
# <s>	        Defines text that is no longer correct
"s" => Dict( "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
#
	    Defines a section in a document
"section" => Dict( "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
# <select>	    Defines a drop-down list
"select" => Dict( "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
# 	    Defines smaller text
"small" => Dict( "font" =>  Dict( "size" => 10, "family" => "Sans", "color" => "black" )),
# 	    Defines a section in a document
"span" => Dict( "display"=>"inline", "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
# 	    Defines important text
"strong" => Dict( "font" =>  Dict( "weight"=>"bold", "size" => 13, "family" => "Sans", "color" => "black" )),
# 	        Defines subscripted text
"sub" => Dict( "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
# 	    Defines a visible heading for a #  element
"summary" => Dict( "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
# 	        Defines superscripted text""
"sup" => Dict( "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
#
	    Defines a table
"table" => Dict( "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
#
	    Groups the body content in a table
"tbody" => Dict( "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
# 	        Defines a cell in a table
"td" => Dict( "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
# <textarea>	Defines a multiline input control (text area)
"textarea" => Dict( "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
#
	    Groups the footer content in a table
"tfoot" => Dict( "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
# 	        Defines a header cell in a table
"th" => Dict( "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
#
	    Groups the header content in a table
"thead" => Dict( "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
#
	        Defines a row in a table
"tr" => Dict( "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
# <track>	    Defines text tracks for media elements (# <video> and # <audio>)
"track" => Dict( "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
# 	        Defines text that should be stylistically om normal text
"u" => Dict( "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
#
	        Defines an unordered list
"ul" => Dict( "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
# <video>	    Defines a video or movie
"video" => Dict( "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
# <wbr>	        Defines a possible line-break
"wbr" => Dict( "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" ))

)
#= Other element types which may be errelevant in this project

# <fieldset>	Groups related elements in a form
"fieldset" => Dict( "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
# <embed>	    Defines a container for an external (non-HTML) application
"embed" => Dict( "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
# , Defines information about the document
#"head" => Dict( "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
# 	    Defines the root of an HTML document
#"html" => Dict( "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
# <input>	    Defines an input control
"input" => Dict( "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
# 	        Defines a text that has been inserted into a document
"ins" => Dict( "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
# 	        Defines keyboard input
"kbd" => Dict( "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
# <keygen>	    Defines a key-pair generator field (for forms)
"keygen" => Dict( "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
# <link>	    Defines the relationship between a document and an external resource (link to style sheets)
"link" => Dict( "font" =>  Dict( "weight "=> "normal", "size" => 13, "family" => "Sans", "color" => "black" )),
# <meter>	    Defines a scalar measurement within a known range (a gauge)
"meter" => Dict( "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
# <meta>	    Defines metadata about an HTML document
"meta" => Dict( "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
# 	Defines an alternate content for users that do not support client-side scripts
"noscript" => Dict( "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
# <object>	    Defines an embedded object
"object" => Dict( "font" =>  Dict( "size" => 13, "family" => "Sans", "color" => "black" )),
#


text:
function PrintEstimationResults(results, names)
    clabels = ["estimate", "st. err", "t-stat", "p-value"]
    prettyprint(results, clabels, names)
end

ext: jl

max_stars_repo_name: nilshg/Econometrics

max_issues_repo_name: nilshg/Econometrics

max_forks_repo_name: nilshg/Econometrics


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

module Variography

using Meshes
using GeoStatsBase

using Optim
using Distances
using InteractiveUtils: subtypes
using SpecialFunctions: gamma, besselk
using NearestNeighbors: MinkowskiMetric
using Transducers: Map, foldxt
using LinearAlgebra
using RecipesBase
using Parameters
using Setfield
using Printf

import Base: merge, +, *
import GeoStatsBase: fit

include("utils.jl")
include("empirical.jl")
include("partition.jl")
include("varioplane.jl")
include("theoretical.jl")
include("nested.jl")
include("pairwise.jl")
include("fitting.jl")

# plot recipes
include("plotrecipes/empirical.jl")
include("plotrecipes/varioplane.jl")
include("plotrecipes/theoretical.jl")

export
  # empirical variograms
  EmpiricalVariogram,
  EmpiricalVarioplane,
  DirectionalVariogram,
  PlanarVariogram,
  distance,

  # theoretical variograms
  Variogram,
  NuggetEffect,
  GaussianVariogram,
  ExponentialVariogram,
  MaternVariogram,
  SphericalVariogram,
  CubicVariogram,
  PentasphericalVariogram,
  PowerVariogram,
  SineHoleVariogram,
  NestedVariogram,
  isstationary,
  sill, nugget,
  structures,
  distance,
  pairwise,
  pairwise!,

  # fitting methods
  VariogramFitAlgo,
  WeightedLeastSquares

end


text:
module SMC

using Random
using Zygote
using StatsFuns
using Statistics
using StatsBase
using Distributions
using LinearAlgebra
using Printf
using Distributed
using Crayons
using Base.Iterators
using Dates

include("resampler.jl")
include("smc.jl")
include("kalman.jl")
include("smc2_inference.jl")

export @NT

module Models
import ..@NT
include("models/lineargaussian.jl")
include("models/lineargaussian_cached.jl")
include("models/lineargaussian_inference.jl")
end

end


text:
using DrWatson
@quickactivate "ExercisesRepo"
include(srcdir("style.jl"))
using DynamicalSystems, PyPlot, Random

# %% mutual information
using DynamicalSystems, PyPlot
using InformationMeasures, Random
MI = get_mutual_information #shortcut
lo, N = Systems.logistic(), 100
x = trajectory(lo, N-1)
y = trajectory(lo, N-1; Ttr = 1)
x .+= randn(N)/25; y .+= randn(N)/25
m = MI(x, y)
null = [MI(shuffle!(x), shuffle!(y))
        for _ in 1:10000]

# plot stuff
fig = figure(figsize = (0.9figx/2, 0.8figx/2))
using Statistics
μ, σ = mean(null), std(null)
plt.hist(null, 50, label = "null pdf")
axvline(μ, color = "C1", label = "\$\\mu\$", ls = "dashed")
axvline(μ-3σ, color = "C2", label = "\$\\mu \\pm 3\\sigma\$", ls = "dashed")
axvline(μ+3σ, color = "C2", ls = "dashed")
axvline(m, color = "C3", label = "\$m\$")

legend(framealpha = 1.0, loc = "upper right")
yticks([])
xlabel("mutual information (bits)")
tight_layout()
# fsave(joinpath(figdir, "mutualinfo"), fig)

# %% transfer entropy
using DynamicalSystems, CausalityTools, PyPlot, Random

function ulam(dx, x, p, t)
    f(x) = 2 - x^2; ε = p[1]; N = length(x)
    for i in 1:N
        dx[i] = f(ε*x[mod1(i-1, N)] + (1-ε)*x[i])
    end
end
ds = DiscreteDynamicalSystem(ulam, rand(100), [0.04])

methods = [RectangularBinning(r) for r in (0.01, 0.1, 0.4)]
εs = 0.0:0.01:1.0
tes = [zeros(length(εs), 2) for j in 1:length(methods)]

for (i, ε) in enumerate(εs), (j, meth) in enumerate(methods)
    set_parameter!(ds, 1, ε)
    tr = trajectory(ds, 10000; Ttr = 10000)
    X1 = tr[:, 1]; X2 = tr[:, 2]
    tes[j][i, 1] = transferentropy(X1, X2, meth, 1, 1, 1)
    tes[j][i, 2] = transferentropy(X2, X1, meth, 1, 1, 1)
end

fig = figure(figsize = (6figx/10, figx/3))
rs = (0.01, 0.1, 0.4)

for j in 1:length(methods)
    plot(εs, tes[j][:, 1], color = "C$(j-1)", label = "\$r=$(rs[j])\$")
    plot(εs, tes[j][:, 2], color = "C$(j-1)", ls = "dashed")
end
xlabel("coupling strength \$\\epsilon\$")
ylabel("transfer entropy")
legend(loc = "center")
tight_layout()
# fsave(joinpath(figdir, "transfer"), fig)

# %% Surrogate application (From Lancaster)
using DynamicalSystems, TimeseriesSurrogates,
      PyPlot, Statistics, Random

fig = figure()
ax1 = fig.add_subplot(1,2,1)
ax2 = fig.add_subplot(1,4,3)
ax3 = fig.add_subplot(1,4,4)
axs = [ax1, ax2, ax3]

ro = Systems.roessler(ones(3); a = 0.165, b = 0.2, c = 10.0)

tr = trajectory(ro, 500; dt = 0.1, Ttr = 500)
x = tr[:, 1]
x ./= std(x)
x .+= randn(length(x))*std(x)*0.1
axs[1].plot(x .+ 2,  label = "Rössler", lw = 1.2)

# generate autoregressive process
Random.seed!(77163)
η = randn(5000)
s = ones(5000)
for n in 4:5000
    s[n] = 1.625s[n-1] - 0.284s[n-2] - 0.355s[n-3] + η[n] - 0.96η[n-1]
end
s ./= std(s)

# # equivalent ARFIMA (doesn't work)
# using ARFIMA
# φ = SVector(1.625, -0.284, -0.355)
# θ = SVector(0.96)
# Random.seed!(77163)
# s2 = arma(5000, 1.0, φ, θ)
# s2 ./= std(s2)
# plot(s)
# plot(s2)
# xlim(0, 1000)

extrema(s)
axs[1].plot(s .- 2,label = "ARMA", lw = 1.0)
axs[1].set_xlim(0, 1000)
axs[1].set_ylim(-5, 5)
axs[1].set_yticks(-4:2:4)
# axs[1].set_xticklabels([])
leg = axs[1].legend(bbox_to_anchor=(0., 1.02, 1., .102), loc="lower left",
           ncol=2, mode="expand", borderaxespad=0.)
# leg = axs[1].legend(ncol=1, loc = "center right", fontsize = 24)
for line in leg.get_lines()
    line.set_linewidth(4.0)
end
# axs[1].set_ylabel("timeseries")

# Do the surrogate calculation
εro, εma = std.((x, s))./4
algs = [RandomFourier(), AAFT(), IAAFT()]
names = ["FT", "AAFT", "IAAFT"]
sgx = [surrogenerator(x, m) for m in algs]
sgs = [surrogenerator(s, m) for m in algs]
τx = estimate_delay(x, "ac_zero")
τs = estimate_delay(s, "ac_zero")
Cx = takens_best_estimate(embed(x, 3, τx), εro)
Cs = takens_best_estimate(embed(s, 4, τs), εma)
A = length(algs)

xboxes = []
sboxes = []
for i in 1:A
    sx, ss = sgx[i], sgs[i]
    bx, bs = [], []
    for j in 1:100
        X = embed(sx(), 3, τx)
        S = embed(ss(), 4, τs)
        push!(bx, takens_best_estimate(X, εro))
        push!(bs, takens_best_estimate(S, εma))
    end
    push!(xboxes, bx)
    push!(sboxes, bs)
end

for (j, boxes) in enumerate((xboxes, sboxes))
    c = "C$(j-1)"
    ax = axs[j+1]
    for (i, b) in enumerate(boxes)
        ax.boxplot(b; positions = [i],
        patch_artist=true, boxprops=Dict(:facecolor=>c, :color=>c),
        medianprops=Dict(:color=>"w"), flierprops=Dict(:markeredgecolor=>c))
    end
end

l = [1-0.2, A+0.2]
axs[2].plot(l, fill(Cx, 2), color = "C0", ls = "dashed")
# axs[2].plot(l, fill(2.75, 2), color = "C0", ls = "dashed")
axs[3].plot(l, fill(Cs, 2), color = "C1", ls = "dashed")
# axs[2].set_ylim(2.7, 3)
axs[2].set_yticks(2.6:0.2:3.0)
axs[2].set_title("\$D_C\$, Rössler")
axs[3].set_title("\$D_C\$, ARMA")
axs[3].set_yticks(3.7:0.1:3.8)

for ax in axs[2:3]
    ax.grid(false; axis = "x")
    ax.set_xticks(1:A)
    ax.set_xticklabels(names, rotation = 0, size = 20)
end
fig.subplots_adjust(top = 0.85, left = 0.08, bottom = 0.1, right = 0.98, wspace = 0.3)
# fsave(joinpath(figdir, "surrogates"), fig)

# %% Convergent Cross Mapping
using DynamicalSystems, Neighborhood,
Statistics, LinearAlgebra
function ccm(x, y, d, τ)
    Mx = embed(x, d, τ); theiler = Theiler(0); tree = KDTree(Mx)
    idxs = bulkisearch(tree, Mx.data, NeighborNumber(d+1), theiler)
    ỹ = copy(y)
    for i in 1:length(Mx)
        J = idxs[i]
        xᵢ = Mx[i]; n1 = norm(xᵢ - Mx[J[1]])
        w = [exp(-norm(xᵢ - Mx[j])/n1) for j in J]
        w ./= sum(w)
        ỹ[i] = sum(w[k]*y[j] for (k, j) in enumerate(J))
    end
    return cor(y, ỹ)
end

ds = Systems.lorenz()
tr = trajectory(ds, 5000; dt = 0.1)
x, y = columns(tr)
Ns = 5000:5000:length(tr)
ρs = [ccm(x[1:N], y[1:N], 3, 5) for N in Ns]

ccm(x, y, 3, 5)

plot(Ns, ρs)

max_stars_repo_name: moritz-g/ExercisesRepo

max_issues_repo_name: moritz-g/ExercisesRepo

max_forks_repo_name: moritz-g/ExercisesRepo


text:
function plot_chains(chain; options...)
    p1 = plot(chain, xaxis=font(5), yaxis=font(5), seriestype=(:traceplot),
      grid=false, size=(275,125), titlefont=font(5); options...)
    p2 = plot(chain, xaxis=font(5), yaxis=font(5), seriestype=(:autocorplot),
      grid=false, size=(275,125), titlefont=font(5); options...)
    p3 = plot(chain, xaxis=font(5), yaxis=font(5), seriestype=(:mixeddensity),
      grid=false, size=(275,125), titlefont=font(5); options...)
    p = plot(p1, p2, p3, layout=(3,1), size=(300,300); options...)
    return p
end


text:
@enum NeutrinoFlavour begin
  Electron = 1
  Muon = 2
  Tau = 3
end

struct OscillationParameters{T}
    mixing_angles::SparseMatrixCSC{T,<:Integer}
    mass_squared_diff::SparseMatrixCSC{T,<:Integer}
    cp_phases::SparseMatrixCSC{T,<:Integer}
    dim::Int64
    OscillationParameters(dim::Int64) = begin
        new{ComplexF64}(
                spzeros(dim, dim),
                spzeros(dim, dim),
                spzeros(dim, dim),
                dim)
    end
end

function _generate_ordered_index_pairs(n::Integer)
    indices = Vector{Pair{Int64,Int64}}(undef, mixingangles(n))
    a = 1
    for i in 1:n
        for j in 1:i-1
            indices[a] = Pair(j,i)
            a += 1
        end
    end
    indices
end

"""
$(SIGNATURES)

Set a mixing angle of an oscillation parameters struct

# Arguments
- `osc::OscillationParameters`: Oscillation parameters
- `indices::Pair{<:Integer, <:Integer}`: The indices of the mixing angle
- `value<:Real` The value which should be applied to the oscillation parameters

"""
function mixingangle!(osc::OscillationParameters, indices::Pair{T, T}, value::S) where {T <: Integer, S <: Real}
    fromidx = first(indices)
    toidx = last(indices)
    if fromidx < toidx
        osc.mixing_angles[fromidx, toidx] = value
    else
        osc.mixing_angles[toidx, fromidx] = value
    end
end

"""
$(SIGNATURES)

Set a mixing angle of an oscillation parameters struct

# Arguments
- `osc::OscillationParameters`: Oscillation parameters
- `args::Tuple{Pair{<:Integer, <:Integer}, <:Real}`: The indices of the mixing angle

"""
function mixingangle!(osc::OscillationParameters, (args::Tuple{Pair{T, T}, S})...) where {T <: Integer, S<: Real}
    for a in args
        mixingangle!(osc, first(a), last(a))
    end
end

const setθ! = mixingangle!

function _mass_matrix_fully_determined(osc::OscillationParameters)
    I, J, _ = findnz(osc.mass_squared_diff)
    set_elements = collect(zip(I, J))
    indices = Set([first.(set_elements)..., last.(set_elements)...])
    (length(indices) >= osc.dim) & (length(set_elements) >= (osc.dim - 1))
end

function _mass_matrix_overdetermined(osc::OscillationParameters)
    I, J, _ = findnz(osc.mass_squared_diff)
    set_elements = collect(zip(I, J))
    indices = Set([first.(set_elements)..., last.(set_elements)...])
    length(set_elements) >= length(indices)
end

function Base.isvalid(osc::OscillationParameters)
    return _mass_matrix_fully_determined(osc)
end

function _completed_mass_matrix(osc::OscillationParameters)
    tmp = Matrix(osc.mass_squared_diff)
    tmp = tmp - transpose(tmp)
    if _mass_matrix_fully_determined(osc)
        I, J, _ = findnz(osc.mass_squared_diff)
        given_idx = collect(zip(I, J))
        wanted_idx = filter(x->(x[1] < x[2]) & (x ∉ given_idx), collect(Iterators.product(1:osc.dim, 1:osc.dim)))
        graph = SimpleDiGraph(map(x->x!=0.0, tmp))
        for (from, to) in wanted_idx
            path = a_star(graph, from, to)
            for edge in path
                tmp[from, to] += tmp[edge.src, edge.dst]
            end
        end
    else
        error("Mass squared differences not fully determined!")
    end
    UpperTriangular(tmp)
end

"""
$(SIGNATURES)

Set a mass squared difference of an oscillation parameters struct

# Arguments
- `osc::OscillationParameters`: Oscillation parameters
- `indices::Pair{<:Integer, <:Integer}`: The indices of the mass squared difference
- `value` The value which should be applied to the oscillation parameters

"""
function masssquareddiff!(osc::OscillationParameters, indices::Pair{T, T}, value::S) where {T <: Integer, S <: Number}
    fromidx = first(indices)
    toidx = last(indices)
    if fromidx < toidx
        osc.mass_squared_diff[fromidx, toidx] = value
    elseif fromidx == toidx
        error("Mass squared difference with equal index cannot be modified.")
    else
        osc.mass_squared_diff[toidx, fromidx] = -value
    end
    if _mass_matrix_overdetermined(osc)
        @warn "Mass squared difference fields (partially) overdetermined!"
    end
end

"""
$(SIGNATURES)

Set a mass squared difference of an oscillation parameters struct

# Arguments
- `osc::OscillationParameters`: Oscillation parameters
- `args::Tuple{Pair{<:Integer, <:Integer}, <:Number}`: Indices and values of the mass squared difference

"""
function masssquareddiff!(osc::OscillationParameters, (args::Tuple{Pair{<:Integer, <:Integer}, <:Number})...)
    for a in args
        masssquareddiff!(osc, first(a), last(a))
    end
end

const setΔm²! = masssquareddiff!

"""
$(SIGNATURES)

Set a CP phase of an oscillation parameters struct

# Arguments
- `osc::OscillationParameters`: Oscillation parameters
- `indices::Pair{<:Integer, <:Integer}`: The indices of the mass difference
- `value` The value which should be applied to the oscillation parameters

"""
function cpphase!(osc::OscillationParameters, indices::Pair{T, T}, value::S) where {T <: Integer, S <: Real}
    fromidx = first(indices)
    toidx = last(indices)
    if fromidx < toidx
        osc.cp_phases[fromidx, toidx] = value
    else
        osc.cp_phases[toidx, fromidx] = value
    end
end

"""
$(SIGNATURES)

Set a CP phase of an oscillation parameters struct

# Arguments
- `osc::OscillationParameters`: Oscillation parameters
- `args::Tuple{Pair{<:Integer, <:Integer}, <:Number}`: Indices and values of the CP phase

"""
function cpphase!(osc::OscillationParameters, (args::Tuple{Pair{T, T}, S})...) where {T <: Integer, S <: Real}
    for a in args
        cpphase!(osc, first(a), last(a))
    end
end

const setδ! = cpphase!

function PMNSMatrix(osc_params::OscillationParameters; anti=false)
"""
$(SIGNATURES)

Create rotation matrix (PMNS) based on the given oscillation parameters

# Arguments
- `osc_params::OscillationParameters`: Oscillation parameters
- `anti`: Is anti neutrino

"""
    dim = size(osc_params.mixing_angles)[1]
    pmns = Matrix{ComplexF64}(1.0I, dim, dim)
    indices = _generate_ordered_index_pairs(dim)
    for (i, j) in indices
        rot = sparse((1.0+0im)I, dim, dim)
        mixing_angle = osc_params.mixing_angles[i, j]
        c, s = cos(mixing_angle), sin(mixing_angle)
        rot[i, i] = c
        rot[j, j] = c
        rot[i, j] = s
        rot[j, i] = -s
        if CartesianIndex(i, j) in findall(!iszero, osc_params.cp_phases)
            cp_phase = osc_params.cp_phases[i, j]
            cp_term = exp(-1im * cp_phase)
            if anti
                cp_term = conj(cp_term)
            end
            rot[i, j] *= cp_term
            rot[j, i] *= conj(cp_term)
        end
        pmns = rot * pmns
    end
    pmns
end

function Hamiltonian(osc_params::OscillationParameters)
"""
$(SIGNATURES)

Create modified hamiltonian matrix consisting of the squared mass differences
based on the given oscillation parameters

# Arguments
- `osc_params::OscillationParameters`: Oscillation parameters

"""
    Hamiltonian(osc_params, zeros(Float64, osc_params.dim))
end

function Hamiltonian(osc_params::OscillationParameters, lambda)
"""
$(SIGNATURES)

Create modified hamiltonian matrix consisting of the squared mass differences
based on the given oscillation parameters

# Arguments
- `osc_params::OscillationParameters`:  Oscillation parameters
- `lambda`:                             Decay parameters for each mass eigenstate

"""
    full_mass_squared_matrix = _completed_mass_matrix(osc_params)
    H = zeros(ComplexF64, osc_params.dim)
    for i in 1:osc_params.dim
        for j in 1:osc_params.dim
            if i < j
                H[i] += full_mass_squared_matrix[i,j]
            elseif j < i
                H[i] -= full_mass_squared_matrix[j,i]
            end
        end
        H[i] += 1im * lambda[i]
    end
    H /= osc_params.dim
    H
end

function _oscprobampl(U, H, energy, baseline)
    H_diag = 2.5338653580781976 * Diagonal{ComplexF64}(H) * baseline / energy
    U * exp(1im * H_diag) * adjoint(U)
end

function _make_flavour_range(size::Integer)
    if size <= 3
        return NeutrinoFlavour.(1:size)
    else
        return [NeutrinoFlavour.(1:3)..., 4:size...]
    end
end

"""
$(SIGNATURES)

Calculate the transistion probabilities between the neutrino flavours

# Arguments
- `U`:          PMNS Matrix
- `H`:          Hamiltonian
- `energy`:     Energies [GeV]
- `baseline`:   Baselines [km]

"""
function oscprob(U, H, energy::Vector{T}, baseline::Vector{S}) where {T,S <: Real}
    s = (size(U)..., length(energy), length(baseline))
    combinations = collect(Iterators.product(energy, baseline))
    tmp = map(x->abs.(_oscprobampl(U, H, first(x), last(x))).^2, combinations)
    P = reshape(hcat(collect(Iterators.flatten(tmp))), s...)
    P = permutedims(P, (3,4,1,2))
    flavrange = _make_flavour_range(first(size(U)))
    AxisArray(P; Energy=energy, Baseline=baseline, InitFlav=flavrange, FinalFlav=flavrange)
end

const oscprob(U, H, energy::T, baseline::Vector{S}) where {S,T <: Real} = oscprob(U, H, [energy], baseline)
const oscprob(U, H, energy, baseline::T) where {T <: Real} = oscprob(U, H, energy, [baseline])

"""
$(SIGNATURES)

Calculate the transistion probabilities between the neutrino flavours

# Arguments
- `osc_params::OscillationParameters`:  Oscillation parameters
- `energy`:                             Energy [GeV]
- `baseline`:                           Baseline [km]
- `anti`:                               Is anti neutrino

"""
function oscprob(osc_params::OscillationParameters, energy, baseline; anti=false)
    H = Hamiltonian(osc_params)
    U = PMNSMatrix(osc_params; anti=anti)
    Pνν(U, H, energy, baseline)
end

const Pνν = oscprob

"""
$(SIGNATURES)

Returns the number of CP violating phases at given number of neutrino types

# Arguments
- `n`: number of neutrino types in the supposed model

# Examples
```julia-repl
julia> cpphases(3)
1
```
"""
function cpphases(n)
    n < 1 && return 0
    div((n - 1) * (n - 2), 2 )
end

"""
$(SIGNATURES)

Returns the number of mixing angles at given number of neutrino types

# Arguments
- `n`: number of neutrino types in the supposed model

# Examples
```julia-repl
julia> mixingangles(3)
3
```
"""
function mixingangles(n)
    n < 1 && return 0
    div(n * (n - 1), 2)
end


text:
show_version() = "_Version:_ $(today()).\n"

"""
    generate_samples(mean, sd, n)

Return `n` samples from a Gaussian distribution with `mean` and `sd`.

Using `StableRNG` and not `Random.seed` to have a stable seed even for different
Julia versions.
"""
function generate_samples(mean, sd, n)
    rng = StableRNG(1)
    dist = Normal(mean, sd)
    data = rand(rng, dist, n)
    round.(data; digits=2)
end

function commandos_data()
    n = 110
    group = fill("commandos", n)
    neuroticism = generate_samples(111.9, 16.7, n)
    extraversion = generate_samples(161.6, 12.8, n)
    DataFrame(; group, neuroticism, extraversion)
end

function civilians_data()
    n = 275
    group = fill("civilians", n)
    neuroticism = generate_samples(130.9, 37.2, n)
    extraversion = generate_samples(157.4, 33.1, n)
    DataFrame(; group, neuroticism, extraversion)
end

"""
    add_missing!(df)

Real data contains missing values, so we add a few.
At least the second value is also `missing`, so that it is visible for the reader.
"""
function add_missing!(df)
    rng = StableRNG(1)
    n = nrow(df)
    n_indices = round(Int, n/10)
    indices = [2; rand(rng, 1:n, n_indices)]
    allowmissing!(df, [:neuroticism, :extraversion])
    for i in indices
        # Missing can be in :neuroticism or :extraversion.
        # To make the filter a bit more difficult.
        col = rand(rng, [:neuroticism, :extraversion])
        df[i, col] = missing
    end
end

function dataset()
    df = vcat(commandos_data(), civilians_data())
    df[!, :participant] = 1:nrow(df)
    select!(df, :participant, :)
    add_missing!(df)
    # Using missingstring because it's non-trivial to filter empty cells.
    data_path = joinpath(Books.BUILD_DIR, "data.csv")
    CSV.write(data_path, df; missingstring="-1")
    df
end

function stacked_dataset()
    df = dataset()

    variable_name = :domain
    value_name = :score
    sdf = stack(df, [:neuroticism, :extraversion]; variable_name, value_name)
end

function first_and_last_few_rows(df)
    firsts = first(df, 4)
    lasts = last(df, 4)
    cols_and_dots = zip(names(df), fill(["..."], ncol(df)))
    pairs = [a => b for (a, b) in cols_and_dots]
    middle = DataFrame(pairs...)
    vcat(firsts, middle, lasts)
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 62c13856-92f9-4c51-927c-fdd52314aea4
begin
	using Pkg
	Pkg.activate()
	using Revise, Protos, BenchmarkTools
	using Protos.Specs
	import Protos.Serialization: _write_uleb, _write_key, writeproto, writefield, _read_key, readfield, jtype, wiretype
end

# ╔═╡ 81d7c2e7-6eeb-4d01-a183-1b6e669d0d7d
begin
	io3 = IOBuffer()
	_write_uleb(io3, 300)
	take!(io3)
end

# ╔═╡ e8744639-d84c-4a82-8211-5c64fcd6c416
begin
	first_bit_set = 0x80
	seven_bits_set = 0x7f
end

# ╔═╡ e5b8f286-ab63-4e33-ba40-74793d8fda94
# get the rightmost 7 bits as a UInt8
get_seven_bits(i::Integer) = UInt8(seven_bits_set & i)

# ╔═╡ 398dcf24-a03c-49be-a98b-54d6ad91b802
@assert ~zero(UInt16) |> get_seven_bits == 0x7f

# ╔═╡ 38c1160a-ba01-4a05-b226-6ca82ff7147b
# this method constructs a UIntXX from a string of 0's and 1's and
# spaces (which are skipped)
function frombits_str(s, ::Type{T}) where T <: Unsigned
	i = zero(T)
	for c in s
		if c == ' '
			continue
		end
		i = i << 1
		if c == '1'
			i = i | one(T)
		end
	end
	i
end

# ╔═╡ bbd33f3d-27d6-483e-b3d5-b298968a952b
# some string macros for diffent sizes of UInts
begin
	macro frombits8_str(s)
		frombits_str(s, UInt8)
	end
	macro frombits64_str(s)
		frombits_str(s, UInt64)
	end
end

# ╔═╡ 5455261f-0838-49cf-9327-1b76793f1930
target300 = UInt8[frombits8"1010 1100", frombits8"0000 0010"]

# ╔═╡ 5a0b5e21-8867-4aab-8a41-1a1bc0f2b103
frombits8"1 1111 1010"

# ╔═╡ 37a46a78-0e82-4543-8584-8cab5a104292
frombits64"1 1111 1010"

# ╔═╡ 2238d1bb-2a54-4b27-8956-4863323ebc6a
md"""
# Serialization

This is a demonstration of proto serialization.  We assume a couple things:

1. We geneate immutable structs where each field is `Union{Missing, SomeType}`. This lets us treat all fields as optional and make the default value `missing`.
2. We generate a `writeproto` method for each struct which also takes in `IO`... this is just a bunch of code blocks (one for each field) that writes a key and then does field-specific serialization.
"""

# ╔═╡ fd3de1a2-28e3-4ddb-a9b1-ba358d7ba5c5
# # TODO use generated functions here
# begin
# 	writefield(io::IO, idx::Int64, val, ::Val{:uint32}) = writefield(io, idx, val, WIRETYPES[:uint32]...)
# 	writefield(io::IO, idx::Int64, val, ::Val{:string}) = writefield(io, idx, val, WIRETYPES[:string]...)

# 	# fallback for all others
# 	# function writefield(io::IO, idx::Int64, val, ::Val{T}) where T
# 	# 	writefield(io, idx, val, WIRETYP_LENDELIM, writeproto, readproto, typeof(val))
# 	# end
# end

# ╔═╡ 60748857-3b38-48ef-8dc3-5da3bb06515a
# function writefield(io::IO, idx::Int64, val, wiretype, writefn, readfn, T)
# 	n = 0
# 	n += _write_key(io, idx, wiretype)
# 	n += writefn(io, val)
# 	n
# end

# ╔═╡ 777d3098-1fed-4d7e-b8c3-a5896ee30c3e
# suppose this is the generated struct for a message with 2 fields
struct XXX
	f1::Union{Missing, UInt32} # uint32 f1 = 1;
	f2::Union{Missing, String} # string f2 = 2;
end

# ╔═╡ 4cacfe92-665b-4a6d-b33c-d6d780e20a1a
# then this is the serialization
function Protos.Serialization.writeproto(io::IO, x::XXX)
	n = 0

	# this is a generated block for uint32 values.
	# we only need the field name, field number, and field type during generation
	# though maybe we inline the next function call and remove the need
	# for generated code??
	if !ismissing(x.f1)
		n += writefield(io, 1, x.f1, Val(:uint32))
	end

	# same as above but for string
	if !ismissing(x.f2)
		n += writefield(io, 2, x.f2, Val(:string))
	end

	n
end

# ╔═╡ 7a6e441b-9d3f-4bcd-8ead-5b72cff8cb68
x = XXX(150, missing)

# ╔═╡ 7b71d148-213e-4066-9d33-92cac8e17e01
y = XXX(missing, "testing")

# ╔═╡ 35245834-c4f2-446a-ba23-7866db85833f
writeproto(y)

# ╔═╡ cf186cb3-6395-4e8c-aa4e-3cc5ecd4d6d1
struct ZZZ
	x::Union{Missing, XXX} # XXX x = 3;
end

# ╔═╡ 83d54da3-ff87-41be-ae67-118684aad6b2
# TODO this should be generated code
function Protos.Serialization.writeproto(io::IO, z::ZZZ)
	n = 0

	# this is a generated block for a nested proto (i.e. non-primitive).
	# as you can see we will let dispatch handle the complexity of it...
	# this code is pretty generic.
	if !ismissing(z.x)
		n += _write_key(io, 3, Protos.Serialization.WIRETYP_LENDELIM)
		# write the nested proto to a temp buffer,
		# then write then length as a varint and then the raw bytes
		bytes = writeproto(z.x)
		n += _write_uleb(io, length(bytes))
		n += write(io, bytes)
	end

	n
end

# ╔═╡ 8015df0d-13d2-4912-814b-8bdc70904653
# should be: [1a 03 08 96 01]
writeproto(ZZZ(XXX(150, missing)))

# ╔═╡ 000447d3-45e6-4e16-a3e8-0f5b12fcf9ef
md"""
# Deserialization

Similar to serialization, we need to generate a constructor `MyMessage(io::IO)`... i.e. a method which takes an IO stream and constructs the julia object from it.
"""

# ╔═╡ 3127f08e-d7bf-4b50-a516-719634b44046
# # TODO use generated functions here
# begin
# 	readfield(io::IO, ::Val{:uint32}) = readfield(io, WIRETYPES[:uint32]...)
# 	readfield(io::IO, ::Val{:string}) = readfield(io, WIRETYPES[:string]...)

# 	# fallback for all others
# 	# function writefield(io::IO, idx::Int64, val, ::Val{T}) where T
# 	# 	writefield(io, idx, val, WIRETYP_LENDELIM, writeproto, readproto, typeof(val))
# 	# end
# end

# ╔═╡ e2c93b2e-bfa7-44cd-bd10-1a5dcaa876a0
# function readfield(io::IO, wiretype, writefn, readfn, ::Type{T}) where T
# 	readfn(io, T)
# end

# ╔═╡ 7962fd02-2d67-4d12-a220-c76f289f62eb
struct XXX2
	f1::Union{Missing, UInt32} # uint32 f1 = 1;
	f2::Union{Missing, String} # string f2 = 2;

	# this is the additional constructor
	function XXX2(io::IO)
		# the size of this should be the number of fields
		fields = Vector{Any}(missing, 2)

		while (!eof(io))
			idx, wiretype = _read_key(io)

			# first field
			if idx == 1 # <-- this is the field number
				#      |
				# but  V  is the ordering
				fields[1] = readfield(io, Val(:uint32))
			elseif idx == 2 # <-- this is the field number
				#      |
				# but  V  is the ordering
				fields[2] = readfield(io, Val(:string))

			# TODO need better handling around repeated fields, maps, etc

			else
				skip_field(io, wiretype)
			end
		end

		# construct the immutable object using the fields array
		new(fields...)
	end
end

# ╔═╡ 73faa3b0-eda9-4645-a4f3-917fcd93c756
# verify that we can deserialize a previously-serialized XXX into an XXX2
begin
	iob = IOBuffer()
	writeproto(iob, x)
	seekstart(iob)
	XXX2(iob).f1 |> Int
end

# ╔═╡ c03e80d4-22a1-4c9f-bfc1-b2fabf5d1e38
begin
	generate_struct(m::Message) = """
	Base.@kwdef struct $(m.name)
		$(join(generate_struct.(m.fields), "\n    "))
	end
	"""

	generate_struct(f::NormalField) = "$(f.name)::Union{Missing, $(jtype(f.t))} = missing"
end

# ╔═╡ 467db54e-2696-43c5-aad2-854bc83758ae
begin
	is_primitive(f::NormalField) = haskey(Protos.Serialization.WIRETYPES, Symbol(f.t))
end

# ╔═╡ 2645e0fe-0329-46c2-b6da-c341a0326bff
function generate_writefield(f::NormalField)
	if is_primitive(f)
		"""
			if !ismissing(m.$(f.name))
				n += writefield(io, $(f.num), m.$(f.name), Val(:$(f.t)))
			end
		"""
	else
		"""
			if !ismissing(m.$(f.name))
				n += _write_key(io, $(f.num), $(wiretype(f.t)))
				# write the nested proto to a temp buffer,
				# then write then length as a varint and then the raw bytes
				bytes = writeproto(m.$(f.name))
				n += _write_uleb(io, length(bytes))
				n += write(io, bytes)
			end
		"""
	end
end

# ╔═╡ b6f500d8-a98e-43d5-a0e8-9977353e5958
generate_writeproto(m::Message) = """
function Protos.Serialization.writeproto(io::IO, m::$(m.name))
	n = 0

$(join(generate_writefield.(m.fields), "\n"))

	n
end
"""

# ╔═╡ c6ac93fc-3247-4d61-a946-84bfad277e73
generate_readfield(f::NormalField, i) = """
		$(i==1 ? "" : "else")if idx == $(f.num)  # $(f.name)
			fields[$i] = readfield(io, Val(:$(Symbol(f.t))))
"""

# ╔═╡ bae37f75-a99a-4a1d-bb13-76a36d5b2c31
generate_readproto(m::Message) = """
function $(m.name)(io::IO)
	# the size of this should be the number of fields
	fields = Vector{Any}(missing, $(length(m.fields)))

	while (!eof(io))
		idx, wiretype = _read_key(io)

$(join([generate_readfield(f, i) for (i, f) in enumerate(m.fields)], "\n"))
		else
			skip_field(io, wiretype)
		end
	end

	# construct the immutable object using the fields array
	$(m.name)(fields...)
end
"""

# ╔═╡ 6a0a729f-1eaa-42a8-80b4-d5b550537809
generate_message_file(m::Message) = """
$(generate_struct(m))

$(generate_readproto(m))

$(generate_writeproto(m))
"""

# ╔═╡ 5b4fd6f6-887f-4aa3-871b-466a13713d5c
pf = ProtoFile("""
syntax = "proto3";

message Test1 {
	uint32 f1 = 1;
	string f2 = 2;
}

message Test2 {
	Test1 x = 1;
}
""")

# ╔═╡ 3fbebff1-05a2-43c1-a3e9-c26dd4d01b5f
display_as_code(text) = Markdown.parse("""
```julia
$text
```
""")

# ╔═╡ 63feb0fa-9c6f-47bb-a955-c1920595a2d7
generate_message_file(pf.messages[1]) |> display_as_code

# ╔═╡ a541625e-fe4d-4f29-8e4f-152ddbc41e1a
generate_message_file(pf.messages[2]) |> display_as_code

# ╔═╡ 53837856-9021-408a-87eb-f93612f533bf
begin

"""
docs
"""
Base.@kwdef struct Test1
    f1::Union{Missing, UInt32} = missing
    f2::Union{Missing, AbstractString} = missing
end

function Test1(io::IO)
    # the size of this should be the number of fields
    fields = Vector{Any}(missing, 2)

    while (!eof(io))
        idx, wiretype = _read_key(io)

        if idx == 1  # f1
            fields[1] = readfield(io, Val(:uint32))

        elseif idx == 2  # f2
            fields[2] = readfield(io, Val(:string))

        else
            skip_field(io, wiretype)
        end
    end

    # construct the immutable object using the fields array
    Test1(fields...)
end

function Protos.Serialization.writeproto(io::IO, m::Test1)
    n = 0

    if !ismissing(m.f1)
        n += writefield(io, 1, m.f1, Val(:uint32))
    end

    if !ismissing(m.f2)
        n += writefield(io, 2, m.f2, Val(:string))
    end

    n
end

end

# ╔═╡ e1a879dd-ecbc-4855-9566-295a4e08880b
begin

Base.@kwdef struct Test2
    x::Union{Missing, Test1} = missing
end

function Test2(io::IO)
    # the size of this should be the number of fields
    fields = Vector{Any}(missing, 1)

    while (!eof(io))
        idx, wiretype = _read_key(io)

        if idx == 1  # x
            fields[1] = readfield(io, Val(:Test1))

        else
            skip_field(io, wiretype)
        end
    end

    # construct the immutable object using the fields array
    Test2(fields...)
end

function Protos.Serialization.writeproto(io::IO, m::Test2)
    n = 0

    if !ismissing(m.x)
        n += _write_key(io, 1, 2)
        # write the nested proto to a temp buffer,
        # then write then length as a varint and then the raw bytes
        bytes = writeproto(m.x)
        n += _write_uleb(io, length(bytes))
        n += write(io, bytes)
    end

    n
end

end

# ╔═╡ 310a528b-5056-46e3-8525-4649d32e6ef0
# should be 0x08 0x96 0x01
Protos.Serialization.writeproto(x)

# ╔═╡ b7bb0428-a5c2-4991-919f-6f3ce8f64644
Test1(32, "hello")

# ╔═╡ 68803c8a-c847-4f31-89a4-5eaf0b736800
t1 = Test1(f2="world")

# ╔═╡ e30c0ed0-778b-4c0e-b569-60a842ef81bc
bytes = writeproto(t1)

# ╔═╡ b67b05fa-2d5c-4448-b972-5c363a6d5a6b
Test1(IOBuffer(bytes))

# ╔═╡ 7a920003-0888-493d-b7fb-73dcfa25b5b7
t2 = Test2(Test1(f2="hi"))

# ╔═╡ e1ea9818-e913-47ed-b949-afe0be53b46e
bytes2 = writeproto(t2)

# ╔═╡ d30170e9-309a-4cb9-87cd-5049b45e3c8a
Test2(IOBuffer(bytes2))

# ╔═╡ 813080f0-54f5-467d-9915-a64b3e639953
@benchmark Test2(Test1(f1=50))

# ╔═╡ 6548f526-4eb7-4b36-94ed-d75197c40cae
@benchmark writeproto(t1)

# ╔═╡ Cell order:
# ╠═62c13856-92f9-4c51-927c-fdd52314aea4
# ╠═81d7c2e7-6eeb-4d01-a183-1b6e669d0d7d
# ╠═e8744639-d84c-4a82-8211-5c64fcd6c416
# ╠═e5b8f286-ab63-4e33-ba40-74793d8fda94
# ╠═398dcf24-a03c-49be-a98b-54d6ad91b802
# ╠═5455261f-0838-49cf-9327-1b76793f1930
# ╟─38c1160a-ba01-4a05-b226-6ca82ff7147b
# ╟─bbd33f3d-27d6-483e-b3d5-b298968a952b
# ╠═5a0b5e21-8867-4aab-8a41-1a1bc0f2b103
# ╠═37a46a78-0e82-4543-8584-8cab5a104292
# ╟─2238d1bb-2a54-4b27-8956-4863323ebc6a
# ╠═fd3de1a2-28e3-4ddb-a9b1-ba358d7ba5c5
# ╠═60748857-3b38-48ef-8dc3-5da3bb06515a
# ╠═777d3098-1fed-4d7e-b8c3-a5896ee30c3e
# ╠═4cacfe92-665b-4a6d-b33c-d6d780e20a1a
# ╠═7a6e441b-9d3f-4bcd-8ead-5b72cff8cb68
# ╠═310a528b-5056-46e3-8525-4649d32e6ef0
# ╠═7b71d148-213e-4066-9d33-92cac8e17e01
# ╠═35245834-c4f2-446a-ba23-7866db85833f
# ╠═cf186cb3-6395-4e8c-aa4e-3cc5ecd4d6d1
# ╠═83d54da3-ff87-41be-ae67-118684aad6b2
# ╠═8015df0d-13d2-4912-814b-8bdc70904653
# ╟─000447d3-45e6-4e16-a3e8-0f5b12fcf9ef
# ╠═3127f08e-d7bf-4b50-a516-719634b44046
# ╠═e2c93b2e-bfa7-44cd-bd10-1a5dcaa876a0
# ╠═7962fd02-2d67-4d12-a220-c76f289f62eb
# ╠═73faa3b0-eda9-4645-a4f3-917fcd93c756
# ╠═c03e80d4-22a1-4c9f-bfc1-b2fabf5d1e38
# ╠═b6f500d8-a98e-43d5-a0e8-9977353e5958
# ╠═467db54e-2696-43c5-aad2-854bc83758ae
# ╠═2645e0fe-0329-46c2-b6da-c341a0326bff
# ╠═bae37f75-a99a-4a1d-bb13-76a36d5b2c31
# ╠═c6ac93fc-3247-4d61-a946-84bfad277e73
# ╠═6a0a729f-1eaa-42a8-80b4-d5b550537809
# ╠═5b4fd6f6-887f-4aa3-871b-466a13713d5c
# ╟─3fbebff1-05a2-43c1-a3e9-c26dd4d01b5f
# ╠═63feb0fa-9c6f-47bb-a955-c1920595a2d7
# ╠═a541625e-fe4d-4f29-8e4f-152ddbc41e1a
# ╠═53837856-9021-408a-87eb-f93612f533bf
# ╠═e1a879dd-ecbc-4855-9566-295a4e08880b
# ╠═b7bb0428-a5c2-4991-919f-6f3ce8f64644
# ╠═68803c8a-c847-4f31-89a4-5eaf0b736800
# ╠═e30c0ed0-778b-4c0e-b569-60a842ef81bc
# ╠═b67b05fa-2d5c-4448-b972-5c363a6d5a6b
# ╠═7a920003-0888-493d-b7fb-73dcfa25b5b7
# ╠═e1ea9818-e913-47ed-b949-afe0be53b46e
# ╠═d30170e9-309a-4cb9-87cd-5049b45e3c8a
# ╠═813080f0-54f5-467d-9915-a64b3e639953
# ╠═6548f526-4eb7-4b36-94ed-d75197c40cae


text:
function _update_manifest(ctx::Context, pkg::PackageSpec, hash::Union{SHA1, Nothing})
    env = ctx.env
    uuid, name, version, path, special_action, repo = pkg.uuid, pkg.name, pkg.version, pkg.path, pkg.special_action, pkg.repo
    hash === nothing && @assert (path != nothing || pkg.uuid in keys(ctx.stdlibs) || pkg.repo.url != nothing)
    # TODO I think ^ assertion is wrong, add-repo should have a hash
    entry = get!(env.manifest, uuid, Types.PackageEntry())
    entry.name = name
    is_stdlib = uuid in keys(ctx.stdlibs)
    if !is_stdlib
        entry.version = version
        entry.tree_hash = hash
        entry.path = path
        if special_action == PKGSPEC_DEVELOPED
            entry.pinned = false
            entry.repo.url = nothing
            entry.repo.rev = nothing
        elseif special_action == PKGSPEC_FREED
            if entry.pinned
                entry.pinned = false
            else
                entry.repo.url = nothing
                entry.repo.rev = nothing
            end
        elseif special_action == PKGSPEC_PINNED
            entry.pinned = true
        elseif special_action == PKGSPEC_REPO_ADDED
            entry.repo.url = repo.url
            entry.repo.rev = repo.rev
            path = find_installed(name, uuid, hash)
        end
        if entry.repo.url !== nothing
            path = find_installed(name, uuid, hash)
        end
    end

    empty!(entry.deps)
    if path !== nothing || is_stdlib
        if is_stdlib
            path = Types.stdlib_path(name)
        else
            path = joinpath(dirname(ctx.env.project_file), path)
        end

        deps = Dict{String,UUID}()

        # Check for deps in project file
        project_file = projectfile_path(path; strict=true)
        if nothing !== project_file
            project = read_project(project_file)
            deps = project.deps
        else
            # Check in REQUIRE file
            # Remove when packages uses Project files properly
            dep_pkgs = PackageSpec[]
            stdlib_deps = find_stdlib_deps(ctx, path)
            for (uuid, name) in stdlib_deps
                push!(dep_pkgs, PackageSpec(name, uuid))
            end
            reqfile = joinpath(path, "REQUIRE")
            if isfile(reqfile)
                for r in Pkg2.Reqs.read(reqfile)
                    r isa Pkg2.Reqs.Requirement || continue
                    push!(dep_pkgs, PackageSpec(name=r.package))
                end
                registry_resolve!(ctx, dep_pkgs)
                project_deps_resolve!(ctx, dep_pkgs)
                ensure_resolved(ctx, dep_pkgs; registry=true)
            end
            for dep_pkg in dep_pkgs
                dep_pkg.name == "julia" && continue
                deps[dep_pkg.name] = dep_pkg.uuid
            end
        end
        entry.deps = deps
    else
        for path in registered_paths(ctx, uuid)
            data = load_package_data(UUID, joinpath(path, "Deps.toml"), version)
            if data !== nothing
                entry.deps = data
                break
            end
        end
    end
    return
end

# Resolve a set of versions given package version specs
function _resolve_versions!(
    ctx::Context,
    pkgs::Vector{PackageSpec},
    target::Union{Nothing, String} = nothing,
)::Dict{UUID,VersionNumber}
    printpkgstyle(ctx, :Resolving, "package versions...")
    # anything not mentioned is fixed
    uuids = UUID[pkg.uuid for pkg in pkgs]
    uuid_to_name = Dict{UUID, String}(uuid => stdlib for (uuid, stdlib) in ctx.stdlibs)
    uuid_to_name[uuid_julia] = "julia"

    for (name::String, uuid::UUID) in get_deps(ctx, target)
        uuid_to_name[uuid] = name

        uuid_idx = findfirst(isequal(uuid), uuids)
        entry = manifest_info(ctx, uuid)
        if entry !== nothing && entry.version !== nothing # stdlibs might not have a version
            ver = VersionSpec(entry.version)
        else
            ver = VersionSpec()
        end
        if uuid_idx != nothing
            pkg = pkgs[uuid_idx]
            if entry !== nothing && pkg.special_action != PKGSPEC_FREED && entry.pinned
                # This is a pinned package, fix its version
                pkg.version = ver
            end
        else
            push!(pkgs, PackageSpec(name, uuid, ver))
        end
    end

    # construct data structures for resolver and call it
    # this also sets pkg.version for fixed packages
    fixed = _collect_fixed!(ctx, pkgs, uuid_to_name)

    # compatibility
    proj_compat = Types.project_compatibility(ctx, "julia")
    v = intersect(VERSION, proj_compat)
    if isempty(v)
        @warn "julia version requirement for project not satisfied" _module=nothing _file=nothing
    end

    for pkg in pkgs
        proj_compat = Types.project_compatibility(ctx, pkg.name)
        v = intersect(pkg.version, proj_compat)
        if isempty(v)
            pkgerror(string("empty intersection between $(pkg.name)@$(pkg.version) and project ",
                            "compatibility $(proj_compat)"))
        end
        # Work around not clobbering 0.x.y+ for checked out old type of packages
        if !(pkg.version isa VersionNumber)
            pkg.version = v
        end
    end

    reqs = Requires(pkg.uuid => VersionSpec(pkg.version) for pkg in pkgs if pkg.uuid ≠ uuid_julia)
    fixed[uuid_julia] = Fixed(VERSION)
    graph = deps_graph(ctx, uuid_to_name, reqs, fixed)

    simplify_graph!(graph)
    vers = resolve(graph)
    find_registered!(ctx, collect(keys(vers)))
    # update vector of package versions
    for pkg in pkgs
        # Fixed packages are not returned by resolve (they already have their version set)
        haskey(vers, pkg.uuid) && (pkg.version = vers[pkg.uuid])
    end
    uuids = UUID[pkg.uuid for pkg in pkgs]
    for (uuid, ver) in vers
        uuid in uuids && continue
        name = (uuid in keys(ctx.stdlibs)) ? ctx.stdlibs[uuid] : registered_name(ctx, uuid)
        push!(pkgs, PackageSpec(;name=name, uuid=uuid, version=ver))
    end
    return vers
end
# This also sets the .path field for fixed packages in `pkgs`
function _collect_fixed!(ctx::Context, pkgs::Vector{PackageSpec}, uuid_to_name::Dict{UUID, String})
    fixed_pkgs = PackageSpec[]
    fix_deps_map = Dict{UUID,Vector{PackageSpec}}()
    uuid_to_pkg = Dict{UUID,PackageSpec}()
    for pkg in pkgs
        local path
        entry = manifest_info(ctx, pkg.uuid)
        if pkg.special_action == PKGSPEC_FREED && !entry.pinned
            continue
        elseif pkg.special_action == PKGSPEC_DEVELOPED
            @assert pkg.path !== nothing
            path = pkg.path
        elseif pkg.special_action == PKGSPEC_REPO_ADDED
            @assert pkg.tree_hash !== nothing
            path = find_installed(pkg.name, pkg.uuid, pkg.tree_hash)
        elseif entry !== nothing && entry.path !== nothing
            path = pkg.path = entry.path
        elseif entry !== nothing && entry.repo.url !== nothing
            path = find_installed(pkg.name, pkg.uuid, entry.tree_hash)
            pkg.repo = entry.repo
            pkg.tree_hash = entry.tree_hash
        else
            continue
        end

        path = project_rel_path(ctx, path)
        if !isdir(path)
            pkgerror("path $(path) for package $(pkg.name) no longer exists. Remove the package or `develop` it at a new path")
        end

        uuid_to_pkg[pkg.uuid] = pkg
        uuid_to_name[pkg.uuid] = pkg.name
        found_project = collect_project!(ctx, pkg, path, fix_deps_map)
        if !found_project
            collect_require!(ctx, pkg, path, fix_deps_map)
        end
    end

    fixed = Dict{UUID,Fixed}()
    # Collect the dependencies for the fixed packages
    for (uuid, fixed_pkgs) in fix_deps_map
        fix_pkg = uuid_to_pkg[uuid]
        v = Dict{VersionNumber,Dict{UUID,VersionSpec}}()
        q = Dict{UUID, VersionSpec}()
        for deppkg in fixed_pkgs
            uuid_to_name[deppkg.uuid] = deppkg.name
            q[deppkg.uuid] = deppkg.version
        end
        fixed[uuid] = Fixed(fix_pkg.version, q)
    end
    return fixed
end

# install & update manifest
function apply_versions(ctx::Context, pkgs::Vector{PackageSpec}; mode=:add, update=true)::Vector{UUID}
    hashes, urls = _version_data!(ctx, pkgs)
    apply_versions(ctx, pkgs, hashes, urls; mode=mode, update=update)
end

function apply_versions(ctx::Context, pkgs::Vector{PackageSpec}, hashes::Dict{UUID,SHA1}, urls::Dict{UUID,Vector{String}}; mode=:add, update=true)
    probe_platform_engines!()
    new_versions = UUID[]

    pkgs_to_install = Tuple{PackageSpec, String}[]
    for pkg in pkgs
        !is_stdlib(ctx, pkg.uuid) || continue
        pkg.path === nothing || continue
        pkg.repo.url === nothing || continue
        path = find_installed(pkg.name, pkg.uuid, hashes[pkg.uuid])
        if !ispath(path)
            push!(pkgs_to_install, (pkg, path))
            push!(new_versions, pkg.uuid)
        end
    end

    widths = [textwidth(pkg.name) for (pkg, _) in pkgs_to_install]
    max_name = length(widths) == 0 ? 0 : maximum(widths)

    ########################################
    # Install from archives asynchronously #
    ########################################
    jobs = Channel(ctx.num_concurrent_downloads);
    results = Channel(ctx.num_concurrent_downloads);
    @async begin
        for pkg in pkgs_to_install
            put!(jobs, pkg)
        end
    end

    for i in 1:ctx.num_concurrent_downloads
        @async begin
            for (pkg, path) in jobs
                if ctx.preview
                    put!(results, (pkg, true, path))
                    continue
                end
                if ctx.use_libgit2_for_all_downloads
                    put!(results, (pkg, false, path))
                    continue
                end
                try
                    success = install_archive(urls[pkg.uuid], hashes[pkg.uuid], path)
                    if success && mode == :add
                        set_readonly(path) # In add mode, files should be read-only
                    end
                    if ctx.use_only_tarballs_for_downloads && !success
                        pkgerror("failed to get tarball from $(urls[pkg.uuid])")
                    end
                    put!(results, (pkg, success, path))
                catch err
                    put!(results, (pkg, err, catch_backtrace()))
                end
            end
        end
    end

    missed_packages = Tuple{PackageSpec, String}[]
    for i in 1:length(pkgs_to_install)
        pkg, exc_or_success, bt_or_path = take!(results)
        exc_or_success isa Exception && pkgerror("Error when installing package $(pkg.name):\n",
                                                 sprint(Base.showerror, exc_or_success, bt_or_path))
        success, path = exc_or_success, bt_or_path
        if success
            vstr = pkg.version != nothing ? "v$(pkg.version)" : "[$h]"
            printpkgstyle(ctx, :Installed, string(rpad(pkg.name * " ", max_name + 2, "─"), " ", vstr))
        else
            push!(missed_packages, (pkg, path))
        end
    end

    ##################################################
    # Use LibGit2 to download any remaining packages #
    ##################################################
    for (pkg, path) in missed_packages
        uuid = pkg.uuid
        if !ctx.preview
            install_git(ctx, pkg.uuid, pkg.name, hashes[uuid], urls[uuid], pkg.version::VersionNumber, path)
            if mode == :add
                set_readonly(path)
            end
        end
        vstr = pkg.version != nothing ? "v$(pkg.version)" : "[$h]"
        printpkgstyle(ctx, :Installed, string(rpad(pkg.name * " ", max_name + 2, "─"), " ", vstr))
    end

    ##########################################
    # Installation done, update the manifest #
    ##########################################
    if update
        for pkg in pkgs
            uuid = pkg.uuid
            if pkg.path !== nothing || uuid in keys(ctx.stdlibs)
                hash = nothing
            elseif pkg.repo.url !== nothing
                hash = pkg.tree_hash
            else
                hash = hashes[uuid]
            end
            _update_manifest(ctx, pkg, hash)
        end
    end

    prune_manifest(ctx)
    return new_versions
end

function _add_or_develop(ctx::Context, pkgs::Vector{PackageSpec}; new_git = UUID[], mode=:add)
    # copy added name/UUIDs into project
    for pkg in pkgs
        ctx.env.project.deps[pkg.name] = pkg.uuid
    end
    # if a package is in the project file and
    # the manifest version in the specified version set
    # then leave the package as is at the installed version
    for (name::String, uuid::UUID) in ctx.env.project.deps
        entry = manifest_info(ctx, uuid)
        entry !== nothing && entry.version !== nothing || continue
        for pkg in pkgs
            pkg.uuid == uuid && entry.version ∈ pkg.version || continue
            pkg.version = entry.version
        end
    end
    # resolve & apply package versions
    _resolve_versions!(ctx, pkgs)
    new_apply = apply_versions(ctx, pkgs; mode=mode)
    write_env(ctx) # write env before building
    build_versions(ctx, union(new_apply, new_git))
end
# Find repos and hashes for each package UUID & version
function _version_data!(ctx::Context, pkgs::Vector{PackageSpec})
    hashes = Dict{UUID,SHA1}()
    clones = Dict{UUID,Vector{String}}()
    for pkg in pkgs
        !is_stdlib(pkg.uuid) || continue
        pkg.repo.url === nothing || continue
        pkg.path === nothing || continue
        uuid = pkg.uuid
        ver = pkg.version::VersionNumber
        clones[uuid] = String[]
        for path in registered_paths(ctx, uuid)
            info = parse_toml(path, "Package.toml")
            repo = info["repo"]
            repo in clones[uuid] || push!(clones[uuid], repo)
            vers = load_versions(path; include_yanked = true)
            hash = get(vers, ver, nothing)
            hash !== nothing || continue
            if haskey(hashes, uuid)
                hash == hashes[uuid] || @warn "$uuid: hash mismatch for version $ver!"
            else
                hashes[uuid] = hash
            end
        end
        @assert haskey(hashes, uuid)
    end
    foreach(sort!, values(clones))
    return hashes, clones
end

function collect_target_deps!(
    ctx::Context,
    pkgs::Vector{PackageSpec},
    pkg::PackageSpec,
    target::String,
)
    # Find the path to the package
    if pkg.uuid in keys(ctx.stdlibs)
        path = Types.stdlib_path(pkg.name)
    elseif Types.is_project_uuid(ctx, pkg.uuid)
        path = dirname(ctx.env.project_file)
    else
        entry = manifest_info(ctx, pkg.uuid)
        path = (entry.path !== nothing) ?
            project_rel_path(ctx, entry.path) :
            find_installed(pkg.name, pkg.uuid, entry.tree_hash)
    end

    project_path = nothing
    for project_name in Base.project_names
        project_path_cand = joinpath(path, project_name)
        if isfile(project_path_cand)
            project_path = project_path_cand
            break
        end
    end
    project = nothing
    if project_path !== nothing
        project = read_package(project_path)
    end

    # Pkg2 compatibiity with test/REQUIRE
    has_project_test_target = false
    if project !== nothing && !isempty(project.targets)
        has_project_test_target = true
    end
    if target == "test" && !has_project_test_target
        pkg2_test_target_compatibility!(ctx, path, pkgs)
        return
    end

    # Collect target deps from Project
    if project !== nothing
        targets = project.targets
        haskey(targets, target) || return
        for pkg in targets[target]
            uuid = project.extras[pkg]
            push!(pkgs, PackageSpec(pkg, uuid))
        end
    end
    return
end

# When testing or building a dependency, we want that dependency to be able to load its own dependencies
# at top level. Therefore we would like to execute the build or testing of a dependency using its own Project file as
# the current environment. Being backwards compatible with REQUIRE file complicates the story a bit since these packages
# do not have any Project files.
function with_dependencies_loadable_at_toplevel(f, mainctx::Context, pkg::PackageSpec; might_need_to_resolve=false)
    # localctx is the context for the temporary environment we run the testing / building in
    localctx = deepcopy(mainctx)
    localctx.currently_running_target = true
    # If pkg or its dependencies are checked out, we will need to resolve
    # unless we already have resolved for the current environment, which the calleer indicates
    # with `might_need_to_resolve`
    need_to_resolve = false
    is_project = Types.is_project(localctx, pkg)

    target = nothing
    if pkg.special_action == PKGSPEC_TESTED
        target = "test"
    end

    # In order to fix dependencies at their current versions we need to
    # add them as explicit dependencies to the project, but we need to
    # remove them after the resolve to make them not-loadable.
    # See issue https://github.com/JuliaLang/Pkg.jl/issues/1144
    should_be_in_project = Set{UUID}()
    should_be_in_manifest = Set{UUID}()

    # Only put `pkg` and its deps + target deps (recursively) in the temp project
    collect_deps!(seen, pkg; depth) = begin
        # See issue https://github.com/JuliaLang/Pkg.jl/issues/1144
        # depth = 0 - the package itself (should be in project)
        # depth = 1 - direct dependencies/and or target dependencies (should be in project)
        # depth > 1 - indirect dependencies that should only be in the manifest in the end
        if depth <= 1
            push!(should_be_in_project, pkg.uuid)
        else
            push!(should_be_in_manifest, pkg.uuid)
        end
        pkg.uuid in seen && return
        push!(seen, pkg.uuid)
        entry = manifest_info(localctx, pkg.uuid)
        entry === nothing && return
        need_to_resolve |= (entry.path !== nothing)
        localctx.env.project.deps[pkg.name] = pkg.uuid
        for (name, uuid) in entry.deps
            collect_deps!(seen, PackageSpec(name, uuid); depth=depth+1)
        end
    end

    if is_project # testing the project itself
        # the project might have changes made to it so need to resolve
        need_to_resolve = true
        # Since we will create a temp environment in another place we need to extract the project
        # and put it in the Project as a normal `deps` entry and in the Manifest with a path.
        foreach(k->setfield!(localctx.env.project, k, nothing), (:name, :uuid, :version))
        localctx.env.pkg = nothing
        localctx.env.project.deps[pkg.name] = pkg.uuid
        localctx.env.manifest[pkg.uuid] = Types.PackageEntry(
            name=pkg.name,
            deps=get_deps(mainctx, target),
            path=dirname(localctx.env.project_file),
            version=pkg.version,
        )
    else
        # Only put `pkg` and its deps (recursively) in the temp project
        empty!(localctx.env.project.deps)
        localctx.env.project.deps[pkg.name] = pkg.uuid
    end
    # Only put `pkg` and its deps (recursively) in the temp project
    seen_uuids = Set{UUID}()
    collect_deps!(seen_uuids, pkg; depth=0)

    pkgs = PackageSpec[]
    if target !== nothing
        collect_target_deps!(localctx, pkgs, pkg, target)
        for dpkg in pkgs
            # Also put eventual deps of target deps in new manifest
            collect_deps!(seen_uuids, dpkg; depth=1)
        end
    end

    mktempdir() do tmpdir
        localctx.env.project_file = joinpath(tmpdir, "Project.toml")
        localctx.env.manifest_file = joinpath(tmpdir, "Manifest.toml")

        function rewrite_manifests(manifest)
            # Rewrite paths in Manifest since relative paths won't work here due to the temporary environment
            for (uuid, entry) in manifest
                if uuid in keys(localctx.stdlibs)
                    entry.path = Types.stdlib_path(entry.name)
                end
                if entry.path !== nothing
                    entry.path = project_rel_path(mainctx, entry.path)
                end
            end
        end

        rewrite_manifests(localctx.env.manifest)

        # Add target deps to deps (https://github.com/JuliaLang/Pkg.jl/issues/427)
        if !isempty(pkgs)
            target_deps = deepcopy(pkgs)
            _add_or_develop(localctx, pkgs)
            need_to_resolve = false # add resolves
            entry = manifest_info(localctx, pkg.uuid)
            for deppkg in target_deps
                entry.deps[deppkg.name] = deppkg.uuid
            end
        end

        # Might have added stdlibs in `add` above
        rewrite_manifests(localctx.env.manifest)

        local new
        will_resolve = might_need_to_resolve && need_to_resolve
        if will_resolve
            _resolve_versions!(localctx, pkgs)
            new = apply_versions(localctx, pkgs)
        else
            prune_manifest(localctx)
        end
        # Remove deps that we added to the project just to keep their versions fixed,
        # since we know all of the packages should be in the manifest this should be
        # a trivial modification of the project file only.
        # See issue https://github.com/JuliaLang/Pkg.jl/issues/1144
        not_loadable = setdiff(should_be_in_manifest, should_be_in_project)
        Pkg.API.rm(localctx, [PackageSpec(uuid = uuid) for uuid in not_loadable])

        write_env(localctx, display_diff = false)
        will_resolve && build_versions(localctx, new)

        sep = Sys.iswindows() ? ';' : ':'
        withenv("JULIA_LOAD_PATH" => "@$sep$tmpdir", "JULIA_PROJECT"=>nothing) do
            f(localctx)
        end
    end
end

function backwards_compatibility_for_test(
    ctx::Context, pkg::PackageSpec, testfile::String, pkgs_errored::Vector{String},
    coverage; julia_args=``, test_args=``
)
    printpkgstyle(ctx, :Testing, pkg.name)
    if ctx.preview
        @info("In preview mode, skipping tests for $(pkg.name)")
        return
    end
    code = """
        $(Base.load_path_setup_code(false))
        cd($(repr(dirname(testfile))))
        append!(empty!(ARGS), $(repr(test_args.exec)))
        include($(repr(testfile)))
        """
    cmd = ```
        $(Base.julia_cmd())
        --code-coverage=$(coverage ? "user" : "none")
        --color=$(Base.have_color ? "yes" : "no")
        --compiled-modules=$(Bool(Base.JLOptions().use_compiled_modules) ? "yes" : "no")
        --check-bounds=yes
        --inline=$(Bool(Base.JLOptions().can_inline) ? "yes" : "no")
        --startup-file=$(Base.JLOptions().startupfile == 1 ? "yes" : "no")
        --track-allocation=$(("none", "user", "all")[Base.JLOptions().malloc_log + 1])
        $(julia_args)
        --eval $code
    ```
    run_test = () -> begin
        try
            run(cmd)
            printpkgstyle(ctx, :Testing, pkg.name * " tests passed ")
        catch err
            push!(pkgs_errored, pkg.name)
        end
    end
    with_dependencies_loadable_at_toplevel(ctx, pkg; might_need_to_resolve=true) do localctx
        if !Types.is_project_uuid(ctx, pkg.uuid)
            Display.status(localctx, mode=PKGMODE_MANIFEST)
        end

        run_test()
    end
end

function backwards_compat_for_build(ctx::Context, pkg::PackageSpec, build_file::String, verbose::Bool,
                                    might_need_to_resolve::Bool, max_name::Int)
    log_file = splitext(build_file)[1] * ".log"
    printpkgstyle(ctx, :Building,
        rpad(pkg.name * " ", max_name + 1, "─") * "→ " * Types.pathrepr(log_file))
    code = """
        $(Base.load_path_setup_code(false))
        cd($(repr(dirname(build_file))))
        include($(repr(build_file)))
        """
    cmd = ```
        $(Base.julia_cmd()) -O0 --color=no --history-file=no
        --startup-file=$(Base.JLOptions().startupfile == 1 ? "yes" : "no")
        --compiled-modules=$(Bool(Base.JLOptions().use_compiled_modules) ? "yes" : "no")
        --eval $code
        ```
    run_build = () -> begin
        ok = open(log_file, "w") do log
            success(pipeline(cmd, stdout = verbose ? stdout : log, stderr = verbose ? stderr : log))
        end
        if !ok
            n_lines = isinteractive() ? 100 : 5000
            # TODO: Extract last n  lines more efficiently
            log_lines = readlines(log_file)
            log_show = join(log_lines[max(1, length(log_lines) - n_lines):end], '\n')
            full_log_at, last_lines =
            if length(log_lines) > n_lines
                "\n\nFull log at $log_file",
                ", showing the last $n_lines of log"
            else
                "", ""
            end
            @error "Error building `$(pkg.name)`$last_lines: \n$log_show$full_log_at"
        end
    end
    with_dependencies_loadable_at_toplevel(ctx, pkg;
                                           might_need_to_resolve=might_need_to_resolve) do localctx
        run_build()
    end
end


text:
function _default{M<:GLNormalAttributeMesh}(mesh::TOrSignal{M}, s::Style, data::Dict)
    @gen_defaults! data begin
        main = mesh
        boundingbox = const_lift(GLBoundingBox, mesh)
        shader      = GLVisualizeShader("fragment_output.frag", "util.vert", "attribute_mesh.vert", "standard.frag")
    end
end

function _default{M<:GLNormalMesh}(mesh::TOrSignal{M}, s::Style, data::Dict)
    @gen_defaults! data begin
        main = mesh
        color       = default(RGBA{Float32}, s)
        boundingbox = const_lift(GLBoundingBox, mesh)
        shader      = GLVisualizeShader("fragment_output.frag", "util.vert", "standard.vert", "standard.frag")
    end
end
function _default{M<:GLNormalVertexcolorMesh}(mesh::TOrSignal{M}, s::Style, data::Dict)
    @gen_defaults! data begin
        main = mesh
        boundingbox  = const_lift(GLBoundingBox, mesh)
        shader       = GLVisualizeShader("fragment_output.frag", "util.vert", "vertexcolor.vert", "standard.frag")
    end
end

function _default(mesh::GLNormalColorMesh, s::Style, data::Dict)
    data[:color] = decompose(RGBA{Float32}, mesh)
    _default(GLNormalMesh(mesh), s, data)
end

function _default{M<:GLPlainMesh}(main::TOrSignal{M}, ::style"grid", data::Dict)
    @gen_defaults! data begin
        primitive::GLPlainMesh = main
        color           = default(RGBA, s, 1)
        bg_colorc       = default(RGBA, s, 2)
        grid_thickness  = Vec3f0(2)
        gridsteps       = Vec3f0(5)
        boundingbox     = const_lift(GLBoundingBox, mesh)
        shader          = GLVisualizeShader("fragment_output.frag", "grid.vert", "grid.frag")
    end
end


text:
export AbstractPolicy

"""
A policy is a functional object to generate an action given a state.
"""
abstract type AbstractPolicy end

learner(π::AbstractPolicy) = π.learner

function update!(π::AbstractPolicy, buffer::AbstractTurnBuffer; kw...)
    transitions = extract_transitions(buffer, π)
    if !isnothing(transitions)
        update!(learner(π), transitions; kw...)
    end
end


text:
# Hamming Code
# ============

type HammingCode{k} <: AbstractCode{k}
    codewords::Vector{DNAKmer{k}}
    mindist::Int
end

function distance{k}(::Type{HammingCode{k}}, x::DNAKmer{k}, y::DNAKmer{k})
    return mismatches(x, y)
end

function genbarcodes{k}(::Type{HammingCode{k}}, mindist::Integer, n::Integer=typemax(Int))
    @assert k < 32
    codewords = DNAKmer{k}[]
    x = typemin(UInt64)
    while length(codewords) < n && x < 4^k
        kmer_x = convert(DNAKmer{k}, x)
        for kmer_y in codewords
            if distance(HammingCode{k}, kmer_x, kmer_y) < mindist
                @goto next
            end
        end
        push!(codewords, kmer_x)
        @label next
        x += 1
    end
    return HammingCode(codewords, mindist)
end

function demultiplex{k}(code::HammingCode{k}, seq::DNASequence)
    barcode = extractkmer(DNAKmer{k}, seq, 1)
    i_min = 0
    dist_min = typemax(Int)
    for (i, codeword) in enumerate(code.codewords)
        dist = distance(HammingCode{k}, codeword, barcode)
        if dist < dist_min
            i_min = i
            dist_min = dist
        end
    end
    return i_min
end

function extractkmer{k}(::Type{DNAKmer{k}}, seq::DNASequence, from::Integer)
    x = UInt64(0)
    for i in from:from+k-1
        x <<= 2
        nt = seq[i]
        if isambiguous(nt)
            x |= UInt64(DNA_A)
        else
            x |= UInt64(nt)
        end
    end
    return DNAKmer{k}(x)
end


text:
function g_ij(i::Int64,j::Int64,t)
    ##OverFlowError if i>20

    s=0 ##The sum of each term
    for k in j:i
        e=exp(-k*(k-1)*t/2)
        numerator=((2*k)-1)*((-1)^(k-j))
        denominator=factorial(j)*factorial(k-j)*(j+k-1)

        mult_term=1
        for y in 0:(k-1)
            mult_term*=(((j+y)*(i-y))/(i+y))
        end
        val=e*(numerator/denominator)*mult_term
        s+=val
    end
    return(s)
end

function nHistories(i::Int64,j::Int64)
    ##the number of histories where i lineages coalesce into j lineages

    mult_term=1
    (j==1) && (j=2) ## coalescing into 1 or 2 lineages have the same number of histories but binomial(k,2) returns 0 when k=1
    for k in (j+1):i
        mult_term*=binomial(k,2)
    end
    return(mult_term)
end

function nTops(n::Int64,labeled=true::Bool,rooted=true::Bool)
    if labeled ##Number of labeled tree topologies
        n==1 && return(1)
        return(Int(factorial(2*n-3)/((2^(n-2)*factorial(n-2)))))

    else ##number of unlabeled topologies i.e. tree shapes
        n==1 && return(1) ##base case for recursion
        odd=true
        (mod(n,2)==0) && (odd=false)

        if odd ## n is odd
            val=0
            for p in 1:((n-1)/2)
                val+=nTops(p,false)*nTops(n-p,false)
            end
        else ## n is even
            t2=nTops(n/2,false)
            term=t2*(t2+1)/2

            s=0
            for p in 1:((n/2)-1)
                s+=nTops(p,false)*nTops(n-p,false)
            end
            val=term+s
        end
        return(Int(val))
    end
end

function nLabeledShape(n::Int64,s::Int64)
    ##Compute the number of labeled topologies of a given tree shapes
    ##Where:
        #n= # of taxa
        #s= # of symetric nodes i.e. nodes where its two subtrees have the same shape
    return(factorial(n)/(2^s))
end

function nHistoriesTop(desc_nodes::Array{Int64,1})
    return(factorial(length(desc_nodes))*prod((x->(1/(x+1)).(desc_nodes))))
end

function unlabeledGenerate(n)
    ##we will create dictionaries to store the bits we want for unlabeled tree shapes
    ##Each key in the dictionary will point to a vector of the bits for each tree shape
        ##e.g. sym_tops[5] will give us a vector of how many points of symmetry each tree shape with 5 taxa has
    vect_tops=Dict{Int64, Array{Array{Int64,1},} }() ##this bit is the node vector representation of the tree shape. Each element in the vector represents an internal node and it denotes how many internal nodes it has as descendents
    sym_tops=Dict{Int64,Array{Int64,1}}() ##this bit stores the number of reflection points for a given tree shape
    prod_tops=Dict{Int64,Array{Float64,1}}() ##if x is the node vector representation of the tree shape, this bit stores the product of 1/(1+x)

    ##Initialize the first two sets
    vect_tops[1]=[[]]
    sym_tops[1]=[0]
    prod_tops[1]=[1]

    vect_tops[2]=[[0]]
    sym_tops[2]=[1]
    prod_tops[2]=[1]

    for i in 3:n
        vect_tops[i]=[]
        sym_tops[i]=[]
        prod_tops[i]=[]

        for j in 1:(floor(i/2))
            ##make new shapes by combining all the the j_th and (n-j)_th arrays

            bottom=j
            top=i-j

            for k in 1:(length(vect_tops[bottom]))
                for l in 1:(length(vect_tops[top]))
                   first=vect_tops[bottom][k]
                   second=vect_tops[top][l]

                   c=Int64[]
                   append!(c,first)
                   append!(c,second)
                   push!(c,length(c))
                   push!(vect_tops[i],c)

                   sym_new=(sym_tops[bottom][k]+sym_tops[top][l]) ##The number of symmetry points is the sum of the symmetry points of the two subtrees
                   (first==second) && (sym_new+=1) ## If the two shapes are the same then this makes another point of symmetry
                   push!(sym_tops[i],sym_new)

                   push!( prod_tops[i],prod( (x->(1/(x+1))).(c) ) ) ##TODO could make this more efficient by instead computing a new w from the c vector, just multiply the to previous products together
                end
            end
        end
    end
    return ((vect_tops,sym_tops,prod_tops));
end

function splitsHistories(splits::Array{Int64,1}, tops_mat)
    total=[0]

    spltsHisRec(splits,1,1,tops_mat,total)

    return total[1];

end

function spltsHisRec(splits::Array{Int64,1},ind_splits,running_val,tops_mat,total::Array{Int64,1} )
    split=splits[ind_splits]
    ##Get all the symmetries and products for the split we're interested in
    syms=tops_mat[2][splits[ind_splits]]
    prods=tops_mat[3][splits[ind_splits]]

    nshapes=length(syms) ##the number of treeshapes for a given number of taxa in the split

    for shape in 1:nshapes ##Do calculation for each shape
        val=running_val ##reset the value for each shape

        val*=factorial(split)/(2^(syms[shape]))
        val*=prods[shape]

        if ind_splits==length(splits) ##Base Case: we went thru each split and did the calculation
            val*=factorial(sum((x->x-1).(splits)))
            total[1]+=val
        else ##we still have more splits to do calculations
            spltsHisRec(splits,ind_splits+1,val,tops_mat,total)
        end
    end
end

function findLowNode(net::HybridNetwork) ##TODO. Make this generate a list of all the nodes to do in order instead of one at a time. Also, since net.node never changes until the end, we could have it store the node indices instead of pointers to the nodes themselves.
    numbers=(x->x.number).(net.node) ##store old numbers '
    nd_changed= net.nodes_changed ##store node ordering
    rankNodes!(net)

    ##get all hybrid nodes
    hybs=net.hybrid

    nds=Set() ##This will have all the nodes we need to decompose
    for hyb in hybs
        union!(nds,Set(preorder(net,hyb))) ##add the hyb node and all its descendents to the list
    end
    nds=collect(nds) ##turn the nodes into an array
    nds= nds[(x->(!x.leaf)).(nds)] ##get rid of tips

    nds=nds[sortperm((x->x.number).(nds))]
    ((x,y)->x.number=y).(net.node,numbers) ##reset node numbers to their old values
    net.nodes_changed=nd_changed ##reset node ordering

    return nds
end

function decompILS(network::HybridNetwork,decomp_node::Node,prb::Float64,tm,hyb_sorting::Dict{Int64, Dict{Int64,Set{String}}})

    ## Make a copy of the network that we can muck around with.
    nd_ind=findfirst(x->x==decomp_node,network.node)
    net=deepcopy(network)
    nd=net.node[nd_ind]

    par=getParents(nd)[1]
    nmes=(x->x.name).(getChildren(nd))

    c_edge=getChildrenEdges(nd)[1].length   ##get the edge length of the first child, ignore other child edge lengths, assuming ultrametricy they should all be the same anyways
    p_edge=getChildrenEdges(par)[1].length  ##get the edge length of the parent
    elen=c_edge+p_edge ##This will be the edge length we assign to our coalesced edges

    removeClade!(net,nd,false)
    par_ind=findfirst(x->x==par,net.node)

    subsets=collect(partitions(nmes)) ##we need the power set (excluding the empty set) of all names for decomposition

    ##We ultimately want ot report the decomposed trees and their respective probabilities after everything is all said and done
    decomposed_trees=HybridNetwork[]
    probs=Float64[]

    for subset in subsets
        ##start with a fresh tree for adding things back
        decomp_net=deepcopy(net)
        decomp_par=decomp_net.node[par_ind]

        new_node_num=maximum((x-> x.number).(net.node))
        new_edge_num=maximum((x-> x.number).(net.edge))

        nms_coal=String[] ##These are the names of the coalesced tips we're adding to the tree
        for nameset  in subset ##We need the individual sets within subset to concatenate the names
            if length(nameset)==1
                nm=nameset[1]
            else
                nm=join(nameset,"|")
            end
            push!(nms_coal,nm)
        end

        for nm in nms_coal #Sequentially add each coalesced tip back to the tree
            ##create a new node and edge for the tip
            ntip=Node(new_node_num,true)
            ntip.name=nm
            nedge=Edge(new_edge_num,elen)
            new_edge_num+=1
            new_node_num+=1

            ##make appropriate connections
            push!(decomp_par.edge,nedge)
            push!(nedge.node,ntip)
            push!(nedge.node,decomp_par)
            push!(ntip.edge,nedge)

            ##add the new nodes and edges to the network
            push!(decomp_net.node,ntip)
            push!(decomp_net.edge,nedge)
            push!(decomp_net.leaf,ntip)
        end

        ##update the attributes of the network
        num_added=length(nms_coal)
        decomp_net.numTaxa+=num_added
        decomp_net.numNodes+=num_added
        decomp_net.numEdges+=num_added
        preorder!(decomp_net) ##preorder things for good measure

        ##Compute the probability of the tree
        gij=g_ij(length(nmes),length(nms_coal),p_edge) ## the probability of i lineages coalescing into j lineages by time t
        hp=splitsHistories(length.(subset),tm) ##The number of coalescent histories that lead to these partitions
        ch=nHistories(length(nmes),length(nms_coal)) ##The total number of coalescent histories for i coalescing into j lineages

        treeprob=gij*(hp/ch)*prb ##TODO I'm concerned with underflow so make this deal with small numbers. perhaps switch to log probs
        treeprob==0 && error("Likely underflow of the probability")

        push!(decomposed_trees,decomp_net)
        push!(probs,treeprob)
    end
    ##Make copies of hyb_sorting for each decomposed tree
    hybsort= Dict{Int64, Dict{Int64,Set{String}}}[]
    for i in 1:length(probs)
        push!(hybsort,deepcopy(hyb_sorting))
    end
    return collect(zip(decomposed_trees,probs,hybsort))
end

function emptyHybSorting(net::HybridNetwork;type="name")
    if type == "name"
        hybsort= Dict{Int64, Dict{Int64,Set{String}}}()
        hybs=net.hybrid
        for hyb in hybs
            hybsort[hyb.number]=Dict{Int64,Set{String}}()
            for nd in getParents(hyb)
            hybsort[hyb.number][nd.number]=Set{String}()
            end
        end
    elseif type == "node"
        hybsort= Dict{Int64, Dict{Int64,Set{Node}}}()
        hybs=net.hybrid
        for hyb in hybs
            hybsort[hyb.number]=Dict{Int64,Set{Node}}()
            for nd in getParents(hyb)
            hybsort[hyb.number][nd.number]=Set{Node}()
            end
        end
    else
        error("we didn't get a valid type. Either 'node' or 'name' is accepted")
    end
    return hybsort
end

function decompHyb(network::HybridNetwork,nd::Node,prb::Float64,hyb_sorting::Dict{Int64, Dict{Int64,Set{String}}} )
    ## Make a copy of the network that we can muck around with.
    nd_ind=findfirst(x->x==nd,network.node)
    net=deepcopy(network)
    nd=net.node[nd_ind]

    hyb_num= nd.number

    lpar=getMajorParent(nd) ##left parent
    rpar=getMinorParent(nd) ##right parent
    lnum=lpar.number
    rnum=rpar.number
    lp_e=getMajorParentEdge(nd)
    rp_e=getMinorParentEdge(nd)

    c_elen=getChildrenEdges(nd)[1].length   ##get the edge length of the first child, ignore other child edge lengths, they should be the same anyway
    lp_elen=lp_e.length
    rp_elen=rp_e.length

    l_gamma=lp_e.gamma
    r_gamma=rp_e.gamma

    removeClade!(net,nd,false)

    lpar_ind=findfirst(x->x==lpar,net.node) ##The index of the left parent in net.node
    rpar_ind=findfirst(x->x==rpar,net.node) ##The index of the right parent in net.node

    ##create all combinations for how the lineages could sort between the two hybrid getParents
    nmes=(x->x.name).(getChildren(nd)) ##get each lineage
    lsort=collect(combinations(nmes)) ##Generate all possible ways that things could sort to the left
    push!(lsort,[]) ## the function above doesn't account for the fact that no species could sort left
    rsort=(x->setdiff(nmes,x)).(lsort) ## everything that doesn't sort left must go right. NOTE: we assume that no lineages have the same name here
    combs=collect(zip(lsort,rsort)) ##associate the proper sortings together to generate all ways that lineages can sort

    ##These will store the decomposed trees and their probabilities
    decomposed_trees=HybridNetwork[]
    probs=Float64[]
    hybsorts=Dict{Int64, Dict{Int64,Set{String}}}[]

    for i in combs
        ##start with a fresh slate for adding things
        hybsort=deepcopy(hyb_sorting)
        decomp_net=deepcopy(net)
        ldecomp_par=decomp_net.node[lpar_ind]
        rdecomp_par=decomp_net.node[rpar_ind]

        leftgoing=first(i)
        rightgoing=last(i)

        ####################
        ###Deal with left###
        ####################
        if length(leftgoing)>1 ##if more than 1 tip, we need to make a node for the mrca of all tips going left
            nedge=Edge(-1,lp_elen)
            ##Make appropriate connections
            push!(ldecomp_par.edge,nedge)
            push!(nedge.node,ldecomp_par)
            ldecomp_par=Node(-1,false,false)
            push!(ldecomp_par.edge,nedge)
            pushfirst!(nedge.node,ldecomp_par)
            ##add the new nodes and tips to the network
            push!(decomp_net.node,ldecomp_par)
            push!(decomp_net.edge,nedge)
            ##update attributes of the network
            decomp_net.numNodes+=1
            decomp_net.numEdges+=1
            r_elen=c_elen
        elseif  length(leftgoing)==0 ##if there is nothing going to the left we need to resolve the parent node
            deleteleaf!(decomp_net,findfirst(x->x===ldecomp_par,decomp_net.node);index=true,simplify=false)
        else ##only 1 tip going left
            r_elen=c_elen+rp_elen
        end
        for nm in leftgoing ##add tips going left to the network
            ntip=Node(-1,true)
            ntip.name=nm
            nedge=Edge(-1,lp_elen+c_elen)

             ##make appropriate connections
             push!(ldecomp_par.edge,nedge)
             push!(ntip.edge,nedge)
             push!(nedge.node,ntip)
             push!(nedge.node,ldecomp_par)

            ##add the new nodes and tips to the network
            push!(decomp_net.node,ntip)
            push!(decomp_net.edge,nedge)
            push!(decomp_net.leaf,ntip)
        end

        ####################
        ###Deal with right##
        ####################
        if length(rightgoing)>1 ##if more than 1 tip, we need to make a node for the mrca of all tips going right
            nedge=Edge(-1,rp_elen)
            ##Make appropriate connections
            push!(rdecomp_par.edge,nedge)
            push!(nedge.node,rdecomp_par)
            rdecomp_par=Node(-1,false,false)
            push!(rdecomp_par.edge,nedge)
            pushfirst!(nedge.node,rdecomp_par)
            ##add the new nodes and tips to the network
            push!(decomp_net.node,rdecomp_par)
            push!(decomp_net.edge,nedge)
            ##update attributes of the network
            decomp_net.numNodes+=1
            decomp_net.numEdges+=1
            r_elen=c_elen
        elseif  length(rightgoing)==0 ##if there is nothing going to the right we need to resolve the parent node
            deleteleaf!(decomp_net,findfirst(x->x===rdecomp_par,decomp_net.node);index=true,simplify=false)
        else #only 1 tip going right
            r_elen=c_elen+rp_elen
        end
        for nm in rightgoing ##add tips to the right parent
            ntip=Node(-1,true)
            ntip.name=nm
            nedge=Edge(-1,r_elen)

             ##make appropriate connections
             push!(rdecomp_par.edge,nedge)
             push!(ntip.edge,nedge)
             push!(nedge.node,ntip)
             push!(nedge.node,rdecomp_par)

            ##add the new nodes and tips to the network
            push!(decomp_net.node,ntip)
            push!(decomp_net.edge,nedge)
            push!(decomp_net.leaf,ntip)
        end
        ##the tree should now be deocmposed

        ##Update the attributes of the network
        num_added=length(rightgoing)+length(leftgoing)
        decomp_net.numNodes+=num_added
        decomp_net.numEdges+=num_added
        decomp_net.numTaxa+=num_added
        preorder!(decomp_net) ##preorder things for good measure

        ##record which lineages went left and right along the reticulation
        left_names=String[]
        right_names=String[]
        lnms= (x->split(x,"|")).(leftgoing)
        rnms= (x->split(x,"|")).(rightgoing)
        println(lnms)
        println(rnms)
        for nms in lnms
            append!(left_names,String.(nms))
        end
        for nms in rnms
            append!(right_names,String.(nms))
        end
        hybsort[hyb_num][lnum]=Set(left_names)
        hybsort[hyb_num][rnum]=Set(right_names)

        ##Compute the probability of the tree
        treeprob=(r_gamma^length(rightgoing))*(l_gamma^length(leftgoing))*prb ##TODO I'm concerned with underflow so make this deal with small numbers. perhaps switch to log probs
        treeprob==0 && error("Likely underflow of the probability")

        push!(decomposed_trees,decomp_net)
        push!(probs,treeprob)
        push!(hybsorts,hybsort)
    end
    return collect(zip(decomposed_trees,probs,hybsorts))
end

function networkDecomposition(net::HybridNetwork)

    hybsorting=emptyHybSorting(net)
    prob=1.0
    q= Array{Tuple{HybridNetwork,Array{Dict{Int64,Dict{Int64,Set{String}}},1},Float64}}() ##Queue of partially decomposed trees along with their hyb_sorting profiles and probabilities
    fully_decomposed =  Array{Tuple{HybridNetwork,Array{Dict{Int64,Dict{Int64,Set{String}}},1},Float64}}() ##Where we store fully decompposed trees
    push!(q,(net,hybsorting,1.0))

    ##determine how big of an unlabeledGenerate we need to run. it should be the size of max(numberLeaves(hyb_node) for all hyb_node in hyb_nodes)
    nleaves=Int64[]
    for hyb in net.hybrid
        e=getChildrenEdges(hyb)[1]
        length(descendants(e))
    end

    while (!isempty(q)) ##while there are still things in the Queue...

        current=popfirst!(q)

        curr_net=current[1]
        curr_hybsort=current[2]
        curr_prob=current[3]

        decomp_nds=findLowNode

    end

    nd=findLowNode(net)

end


text:
include("imports.jl")

path_to_output = ("/scratchdata/alui2/cytof/results/repfam/test-sim-6-7-2/" *
                  "maxtemp10-ntempts20-degree2-N500/output.bson")

output = BSON.load(path_to_output)

# 1. For a fixed state (theta), compute_marg_loglike for different temperatures.
tempers = output[:tempers]
# tempers = CytofRepFAM.MCMC.WSPT.gentempers(2, 20, degree=1)
lls = [CytofRepFAM.Model.compute_marg_loglike(output[:lastState].theta,
                                              output[:c].constants,
                                              output[:d].data,
                                              tau) for tau in tempers];
plt.plot(tempers, lls, marker="o")
plt.savefig("tmp.pdf", bbox_inches="tight")
plt.close()

# plt.plot(log.(tempers), lls - [ll[end] for ll in output[:lls]], marker="o")
# plt.savefig("tmp.pdf", bbox_inches="tight")
# plt.close()