text:
import Graphs: in_edges, in_degree, in_neighbors, out_edges, out_degree, out_neighbors
import Base: isempty

type BayesianNetwork <: AbstractGraph{BayesianNode, BayesianEdge}
    nodes::Array{BayesianNode,1}
    edges::Array{BayesianEdge,1}

    finclist::Array{Array{BayesianEdge,1},1}   #forward incidence list
    binclist::Array{Array{BayesianEdge,1},1}   #backward incidence list

    cpds::Dict{CPD,Any} ## Add types if needed
    checked_cpds::Dict{CPD,Uint8} # Keep a list of having or not: 0 = not processed, 1 = being processed, 2 = not derivable, 3 = derivable

    function BayesianNetwork{V <: BayesianNode}(_nodes::Array{V,1}, _edges::Array{BayesianEdge,1})
        _nodes = convert(Array{BayesianNode,1}, _nodes)
        l_nodes = length(_nodes)
        _cpds=Dict{CPD, Any}()
        _checked_cpds=Dict{CPD,Uint8}()
        if l_nodes > 0
            map(x -> assign_index(x[2],x[1]), enumerate(_nodes))
            for _n in _nodes
                if has_pd(_n)
                    if length(probabilities(_n.pd)) > 0
                        _cpds[P(_n.label)] = _n.pd
                        _checked_cpds[P(_n.label)] = 0
                    end
                end
            end
        end
        if length(_edges) > 0
            map(x -> assign_index(x[2],x[1]), enumerate(_edges))
        end
        b = new(_nodes, Array(BayesianEdge,0), multivecs(BayesianEdge, l_nodes), multivecs(BayesianEdge, l_nodes), _cpds, _checked_cpds)
        for edge in _edges
            add_edge!(b, edge)
        end
        b
    end
end

BayesianNetwork(n::Array{None,1}, e::Array{None,1}) = BayesianNetwork(Array(BayesianNode,0),Array(BayesianEdge,0))
BayesianNetwork{V <: BayesianNode}(n::Array{V,1}) = BayesianNetwork(n,Array(BayesianEdge,0))
BayesianNetwork{V <: BayesianNode}(n::Array{V,1}, e::Array{None,1}) = BayesianNetwork(n,Array(BayesianEdge,0))
BayesianNetwork() = BayesianNetwork([], [])

nodes(g::BayesianNetwork) = g.nodes
num_nodes(g::BayesianNetwork) = length(nodes(g))

edges(g::BayesianNetwork) = g.edges
num_edges(g::BayesianNetwork) = length(edges(g))

cpds(g::BayesianNetwork) = g.cpds
num_cpds(g::BayesianNetwork) = length(cpds(g))

checked_cpds(g::BayesianNetwork) = g.checked_cpds
num_checked_cpds(g::BayesianNetwork) = length(calculated_cpds(g))

function reset_checked_cpds!(g::BayesianNetwork)
    for i in collect(keys(checked_cpds(g)))
        checked_cpds(g)[i] = 0
    end
end

function set_cpd_derivable!(g::BayesianNetwork, cpd::CPD, value::Bool)
    v = value ? 3 : 2
    checked_cpds(g)[cpd] = v
end

function cpd_derivable(g::BayesianNetwork, cpd::CPD)
    if (check_checked_cpd(g, cpd))
        checked_cpds(g)[cpd] == 3
    else
        false
    end
end

function start_processing_cpd!(g::BayesianNetwork, cpd::CPD)
    checked_cpds(g)[cpd] = 1
end

function cpd_is_processing(g::BayesianNetwork, cpd::CPD)
    if check_checked_cpd(g, cpd)
        checked_cpds(g)[cpd] == 1
    else
        false
    end
end

function cpd_has_been_processed(g::BayesianNetwork, cpd::CPD)
    if check_checked_cpd(g, cpd)
        checked_cpds(g)[cpd] > 1
    else
        false
    end
end

function check_checked_cpd(g::BayesianNetwork, cpd::CPD)
    try
        getindex(checked_cpds(g), cpd)
        true
    catch
        false
    end
end

function assign_index(g_elem, i::Int)
    if g_elem.index != 0
        throw("Attempting to reassign node")
    end
    g_elem.index=i
end

function add_node!{T <: BayesianNode}(g::BayesianNetwork, n::T)
    assign_index(n, num_nodes(g) + 1)

    push!(g.nodes, n)

    # Create entries for vertex in backward and forward incidence lists
    push!(g.finclist, Int[])
    push!(g.binclist, Int[])

    ##Update what is put into the cpds dictionary when decided
    if typeof(n) == DBayesianNode && has_pd(n)
        add_probability!(g, P(n.label), n.pd)
        ##Check just made in case someone construct a graph in an non-obvious manner, will probably not be run
        if length(in_edges(n, g)) > 0
            add_probability!(g, P(n.label|[edge.source.label for edge in in_edges(n, g)]))
        end
    else
        add_probability!(g, P(n.label))
    end

    n
end

add_edge!{S <: BayesianNode, T <: BayesianNode}(g::BayesianNetwork, s::S, t::T) = add_edge!(g, BayesianEdge(s,t))

function add_edge!(g::BayesianNetwork, e::BayesianEdge)
    assign_index(e, num_edges(g) + 1)

    u = source(e)
    v = target(e)

    if !nodes_in_network(g, [u, v])
        throw("Attempting to add an edge to a network where the nodes are not present")
    end

    ui = node_index(u)::Int
    vi = node_index(v)::Int

    # Add edge to network and to forward list for source and to backward list for target
    push!(g.edges, e)
    push!(g.finclist[ui], e)
    push!(g.binclist[vi], e)
    add_cpd_for_edge!(g,u,v)

    e
end

function add_cpd_for_edge!(g::BayesianNetwork, s::BayesianNode,t::BayesianNode)
    if typeof(s) == DBayesianNode && has_pd(s) && has_pd(t)
        # g.cpds[CPD(t.label, [edge.source.label for edge in in_edges(t,g)])] = null
    else
        # g.cpds[CPD(t.label, s.label)] = null
    end
end

function nodes_in_network{V <: BayesianNode}(g::BayesianNetwork, ns::Array{V,1})
    for node in ns
        if !node_in_network(g,node)
            return false
        end
    end
    true
end

function node_in_network{V <: BayesianNode}(g::BayesianNetwork, n::V)
    if node_index(n) > length(nodes(g)) || node_index(n) == 0
        false
    else
        nodes(g)[node_index(n)] == n
    end
end

function symbols_in_network(g::BayesianNetwork, ss::Array{Symbol,1})
    for s in ss
        if find_node_by_symbol(g,s) == false
            return false
        end
    end
    true
end

function find_node_by_symbol(g::BayesianNetwork, s::Symbol)
    for node in g.nodes
        if node.label == s
            return node
        end
    end
    false
end

function in_edges{V <: BayesianNode}(n::V, g::BayesianNetwork)
    if !node_in_network(g, n)
        Array{BayesianEdge,1}[]
    else
        g.binclist[node_index(n)]
    end
end

in_degree{V <: BayesianNode}(n::V, g::BayesianNetwork) = length(in_edges(n, g))
in_neighbors{V <: BayesianNode}(n::V, g::BayesianNetwork) = map(e -> e.source, in_edges(n, g))

function out_edges{V <: BayesianNode}(n::V, g::BayesianNetwork)
    if !node_in_network(g, n)
        Array{BayesianEdge,1}[]
    else
        g.finclist[node_index(n)]
    end
end

out_degree{V <: BayesianNode}(n::V, g::BayesianNetwork) = length(out_edges(n, g))
out_neighbors{V <: BayesianNode}(n::V, g::BayesianNetwork) = map(e -> e.target, out_edges(n, g))

add_probability!{V <: PDistribution}(g::BayesianNetwork, cpd::CPD, pd::V) = _add_probability!(g, cpd, pd)
add_probability!(g::BayesianNetwork, cpd::CPD) = _add_probability!(g, cpd, nothing)

function _add_probability!(g::BayesianNetwork, cpd::CPD, pd)
    if symbols_in_network(g, distribution(cpd)) && symbols_in_network(g, conditionals(cpd))
        g.cpds[cpd] = pd
        checked_cpds(g)[cpd] = 0
        true
    else
        false
    end
end

multivecs{T}(::Type{T}, n::Int) = [T[] for _ =1:n]

#################################################
#
#  query
#
################################################

function legal_configuration(bn::BayesianNetwork, cpd::CPD)
    nodes = Set(cpd.conditionals)
    size = 0
    checkedNodes = Set()
    uncheckedNodes = Set()
    for label in cpd.distribution
        uncheckedNodes = union(Set(gather_nodes_by_edge_out(out_edges(find_node_by_symbol(bn,label), bn))), gather_nodes_by_edge_in(in_edges(find_node_by_symbol(bn,label), bn)))
    end
    while 0 < length(uncheckedNodes)
        for node in uncheckedNodes
            inNodes = gather_nodes_by_edge_in(in_edges(find_node_by_symbol(bn,node.label), bn))
            outNodes = gather_nodes_by_edge_out(out_edges(find_node_by_symbol(bn,node.label), bn))
            union!(uncheckedNodes, setdiff(inNodes,checkedNodes))
            union!(uncheckedNodes, setdiff(outNodes,checkedNodes))
            push!(checkedNodes,node)
            delete!(uncheckedNodes,node)
        end
    end

    for node in checkedNodes
        if node.label in nodes
            delete!(nodes,node.label)
        end
    end
    if isempty(nodes)
        true
    else
        false
    end
end

gather_nodes_by_edge_in(edges::Array{BayesianEdge,1}) = Set(map(edge -> edge.source, edges))

gather_nodes_by_edge_out(edges::Array{BayesianEdge,1}) = Set(map(edge -> edge.target, edges))

function check_requirements(g::BayesianNetwork, cpd::CPD)
    reset_checked_cpds!(g)
    _check_requirements(g, cpd)
end

function _check_requirements(g::BayesianNetwork, cpd::CPD)
    if cpd_has_been_processed(g, cpd)
        cpd_derivable(g, cpd)
    else
        if check_cpd(g, cpd)
            set_cpd_derivable!(g, cpd, true)
            true
        elseif cpd_is_processing(g, cpd)
            false
        else
            start_processing_cpd!(g, cpd)

            val = false
            conds = copy(conditionals(cpd))
            l_conds = length(conds)
            dist = copy(distribution(cpd))
            l_dist = length(dist)
            # If there are conditionals, then try bayes rule
            if l_conds > 0 && check_bayes_rule(g, cpd)
                val = true
            # If there are more than one distribution and there are conditionals, then try to split it
            elseif l_dist > 1 && l_conds > 0 && check_split_distribution(g, cpd)
                val = true
            # If there are more than one distribution but no conditionals
            elseif l_dist > 1 && l_conds == 0 && check_joint_probability_distribution(g, cpd)
                val = true
            end
            set_cpd_derivable!(g, cpd, val)
            val
        end
    end
end

function check_split_distribution(g::BayesianNetwork, cpd::CPD)
    conds = copy(conditionals(cpd))
    dist = copy(distribution(cpd))

    new_dists = 0
    if length(dist) > 1
        new_dists = splice!(dist, 2:length(dist))
    else
        throw("Length of distribution must be larger than 1.")
    end
    _check_requirements(g, P(new_dists|conds)) && _check_requirements(g, P(dist|conds))
end

function check_bayes_rule(g::BayesianNetwork, cpd::CPD)
    conds = copy(conditionals(cpd))
    dist = copy(distribution(cpd))

    opposite = P(conds|dist)
    if _check_requirements(g, opposite)
        to_check = [P(conds), P(dist)]
        for i in to_check
            if !_check_requirements(g, i)
                return false
            end
        end
        true
    else
        false
    end
end

function check_joint_probability_distribution(bn::BayesianNetwork, cpd::CPD)
    for d in distribution(cpd)
        node = find_node_by_symbol(bn, d)
        parent_labels = convert(Array{Symbol,1}, map(p -> p.label, in_neighbors(node, bn)))
        if !_check_requirements(bn, P(d|parent_labels))
            return false
        end
    end
    true
end

function check_cpd(g::BayesianNetwork, cpd::CPD)
    try
        getindex(cpds(g), cpd)
        true
    catch
        false
    end
end

has_cached_result(bn::BayesianNetwork, cpd::CPD) = cpd in collect(keys(cpds(bn))) ? true : false;

function query(bn::BayesianNetwork, cpd::CPD)
    if is_naive_bayes_model(bn, cpd)
        query_naive_bayes(bn,cpd)
    else
        throw("Form of cpd not supported by query at the moment")
    end
end

function is_naive_bayes_model(bn::BayesianNetwork, cpd::CPD)
    if typeof(find_node_by_symbol(bn,distribution(cpd)[1])) == DBayesianNode
        for label in conditionals(cpd)
            if !(typeof(find_node_by_symbol(bn, label)) == CBayesianNode)
                return false
            end
        end
    else
        false
    end
    true
end

function query_naive_bayes(bn::BayesianNetwork, cpd::CPD)
    if !legal_configuration(bn,cpd)
        throw("Invaled configuration")
    elseif has_cached_result(bn,cpd)
        cpds(bn)[cpd]
    elseif check_requirements(bn, cpd)
        calculate_probabilities(bn,cpd)
    end
end

function calculate_probabilities(bn::BayesianNetwork, cpd::CPD)
    res = Dict()
    temp = create_probabilities(cpd)
    for p in temp
        has_cached_result(bn, p) ? calculate_probability!(res, cpds(bn)[p], cpd.parameters[p.distribution[1]]) : nothing;
    end
    for key in keys(res)
        res[key] = get_probability(bn, find_node_by_symbol(bn,cpd.distribution[1]), key)*prod(res[key])
    end
    res
end

get_probability(bn::BayesianNetwork, node::BayesianNode, state) = has_cached_result(bn, P(node.label)) ? cpds(bn)[P(node.label)][state][1] : throw("Node mismatch");

create_probabilities(cpd::CPD) = map(cond -> CPD(cond,cpd.distribution), cpd.conditionals)

function calculate_probability!(dict, dist::ProbabilityDensityDistribution, v)
    for state in states(dist)
        if !(state in collect(keys(dict)))
            dict[state] = Array(Float64,0)
        end
    end

    for i = 1:length(states(dist))
        push!(dict[states(dist)[i]],pdfs(dist)[i](v))
    end
end


text:
# ==================================================================================================================================
# Module https://github.com/mkretlow/MOID.jl : Mike Kretlow [astrodynamics.de], Start 2018, MIT "Expat" License
# ==================================================================================================================================

module MOID

using Libdl
using Printf

export moidF

# Load shared library

const lib = find_library(["libmoid"],["deps", joinpath(dirname(@__FILE__), "..", "deps")])

if isempty(lib)
	error("Could not find shared library")
end

@info("Using shared library $lib")

"""
   `MOID = moidF(sma1, ecc1, peri1, node1, incl1, sma2, ecc2, peri2, node2, incl2)`

Calculating Accurate Minimum Orbit Intersection Distances (MOID) by calling original Fortran routine

## Args

   sma1, ecc1, peri1, node1, incl1, sma2, ecc2, peri2, node2, incl2 [angles in deg !]

## Return
   MOID in units like sma (usually AU)

## Reference

T.Wisniowski and H.Rickman "A Fast, Geometric Method for Calculating Accurate Minimum Orbit Intersection Distances (MOIDs)"
published in 2013 in Acta Astronomica.

See also: http://moid.cbk.waw.pl/orbity/default/index

"""
function moidF(sma1::Float64, ecc1::Float64, peri1::Float64, node1::Float64, incl1::Float64,
                     sma2::Float64, ecc2::Float64, peri2::Float64, node2::Float64, incl2::Float64)

   ccall((:moid_,lib),Float64,(Ref{Float64},Ref{Float64},Ref{Float64},Ref{Float64},Ref{Float64},
                               Ref{Float64},Ref{Float64},Ref{Float64},Ref{Float64},Ref{Float64}),
                               sma1, ecc1, peri1, node1, incl1, sma2, ecc2, peri2, node2, incl2)
end

end  # module

#==================================================================================================================================#

ext: jl


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
        el
    end
end

# ╔═╡ 6c1969e0-02f5-11eb-3fa2-09931a63b1ac
begin
    using PoreMatMod, PlutoUI, Bio3DView
    HOME = joinpath(homedir(), ".PoreMatModGO")
	rc[:paths][:crystals] = HOME
	rc[:paths][:moieties] = HOME
	if !isdir(HOME)
		mkdir(HOME)
	end
	cd(HOME)
md"""
# 💠 PoreMatModGO 🚀

This notebook interactively substitutes moieties within a crystal using a modified implementation of Ullmann's algorithm to perform substructure searches and applying singular value decomposition to align fragments of the generated materials. Read the docs [here](SimonEnsemble.github.io/PoreMatMod.jl).

See the original publication on PoreMatMod.jl here: [(article)](https://doi.org/10.33774/chemrxiv-2021-vx5r3) [(GitHub)](https://github.com/SimonEnsemble/PoreMatMod.jl)
"""
end

# ╔═╡ 5dc43a20-10b8-11eb-26dc-7fb98e9aeb1a
md"""
Adrian Henle, [Simon Ensemble](http://simonensemble.github.io), 2021

$(Resource("https://simonensemble.github.io/osu_logo.jpg", :width => 250))
"""

# ╔═╡ 90696d20-10b7-11eb-20b5-6174faeaf613
@bind load_inputs Button("Reset")

# ╔═╡ 50269ffe-02ef-11eb-0614-f11975d991fe
begin load_inputs
    # input fields: query, replacement, xtal
    md"""
    ##### Input Files

    Parent Crystal $(@bind parent_crystal FilePicker())

    Search Moiety $(@bind search_moiety FilePicker())

    Replace Moiety $(@bind replace_moiety FilePicker())
    """
end

# ╔═╡ 33b1fb50-0f73-11eb-2ab2-9d2cb6c5a533
# write file input strings to files in temp directory
begin
    # dict for tracking load status of inputs
    isloaded = Dict([:replacement => false, :query => false, :parent => false])
    # replacement loader
    if !isnothing(replace_moiety)
        write("replacement.xyz", replace_moiety["data"])
        replacement = moiety("replacement.xyz")
        isloaded[:replacement] = true
    end
    # query loader
    if !isnothing(search_moiety)
        write("query.xyz", search_moiety["data"])
        query = moiety("query.xyz")
        isloaded[:query] = true
    end
    # xtal loader
    if !isnothing(parent_crystal)
        write("parent.cif", parent_crystal["data"])
        xtal = Crystal("parent.cif", check_overlap=false)
        Xtals.strip_numbers_from_atom_labels!(xtal)
        infer_bonds!(xtal, true)
        isloaded[:parent] = true
    end
    # run search and display terminal message
    if isloaded[:query] && isloaded[:parent]
        search = query ∈ xtal
        with_terminal() do
            @info "Search Results" isomorphisms=nb_isomorphisms(search) locations=nb_locations(search)
        end
    end
end

# ╔═╡ 415e9210-0f71-11eb-15c8-e7484b5be309
# choose replacement type
if all(values(isloaded))
    md"""
    ### Find/Replace Options

    Mode $(@bind replace_mode Select(["", "random replacement at each location", "random replacement at n random locations", "random replacement at specific locations", "specific replacements"]))
    """
end

# ╔═╡ 3997c4d0-0f75-11eb-2976-c161879b8d0c
# options populated w/ conditional logic based on mode selection
if all(values(isloaded))
	local output = nothing
	x = ["$(x)" for x in 1:nb_locations(search)]
	if replace_mode == "random replacement at each location"
		output = nothing
	elseif replace_mode == "random replacement at n random locations"
		output = md"Number of locations $(@bind nb_loc Slider(1:nb_locations(search)))"
	elseif replace_mode == "random replacement at specific locations"
		output = md"Locations $(@bind loc MultiSelect(x))"
	elseif replace_mode == "specific replacements"
		output = md"""
	Locations $(@bind loc MultiSelect(x))

	Orientations $(@bind ori TextField())
	"""
	end
	output
end

# ╔═╡ 69edca20-0f94-11eb-13ba-334438ca2406
if all(values(isloaded))
	new_xtal_flag = true
	if replace_mode == "random replacement at each location"
		new_xtal = substructure_replace(search, replacement)
	elseif replace_mode == "random replacement at n random locations" && nb_loc > 0
		new_xtal = substructure_replace(search, replacement, nb_loc=nb_loc)
	elseif replace_mode == "random replacement at specific locations" && loc ≠ []
		new_xtal = substructure_replace(search, replacement, loc=[parse(Int, x) for x in loc])
	elseif replace_mode == "specific replacements"
		if loc ≠ [] && ori ≠ "" && length(loc) == length(split(ori, ","))
			new_xtal = substructure_replace(search, replacement,
				loc=[parse(Int, x) for x in loc],
				ori=[parse(Int, x) for x in split(ori, ",")])
		else
			new_xtal_flag = false
		end
	else
		new_xtal_flag = false
	end
	if new_xtal_flag
		with_terminal() do
			if replace_mode == "random replacement at each location"
				@info replace_mode new_xtal
			elseif replace_mode == "random replacement at n random locations"
				@info replace_mode nb_loc new_xtal
			elseif replace_mode == "random replacement at specific locations"
				@info replace_mode loc new_xtal
			elseif replace_mode == "specific replacements"
				@info replace_mode loc ori new_xtal
			end
		end
	end
end;

# ╔═╡ 5918f770-103d-11eb-0537-81036bd3e675
if all(values(isloaded)) && new_xtal_flag
	write_cif(new_xtal, "crystal.cif")
	write_xyz(new_xtal, "atoms.xyz")
	write_vtk(new_xtal.box, "unit_cell.vtk")
	write_bond_information(new_xtal, "bonds.vtk")
	no_pb = deepcopy(new_xtal)
	drop_cross_pb_bonds!(no_pb)
	write_mol2(new_xtal, filename="crystal.mol2")
	write_mol2(no_pb, filename="view.mol2")
	viewfile("view.mol2", "mol2", vtkcell="unit_cell.vtk", axes=Axes(4, 0.25))
end

# ╔═╡ 31832e30-1054-11eb-24ed-219fd3e236a1
if all(values(isloaded)) && new_xtal_flag
    download_cif = DownloadButton(read("crystal.cif"), "crystal.cif")
    download_box = DownloadButton(read("unit_cell.vtk"), "unit_cell.vtk")
    download_xyz = DownloadButton(read("atoms.xyz"), "atoms.xyz")
    download_bonds = DownloadButton(read("bonds.vtk"), "bonds.vtk")
	download_mol2 = DownloadButton(read("crystal.mol2"), "crystal.mol2")
md"""
### Output Files
Complete Crystal $download_mol2 $download_cif

Components $download_xyz $download_bonds $download_box
"""
end

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
Bio3DView = "99c8bb3a-9d13-5280-9740-b4880ed9c598"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
PoreMatMod = "2de0d7f0-0963-4438-8bc8-7e7ffe3dc69a"

[compat]
Bio3DView = "~0.1.3"
PlutoUI = "~0.7.14"
PoreMatMod = "~0.2.0"
"""

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

[[ANSIColoredPrinters]]
git-tree-sha1 = "574baf8110975760d391c710b6341da1afa48d8c"
uuid = "a4c015fc-c6ff-483c-b24f-f7ea428134e9"
version = "0.0.1"

[[AbstractPlutoDingetjes]]
deps = ["Pkg"]
git-tree-sha1 = "abb72771fd8895a7ebd83d5632dc4b989b022b5b"
uuid = "6e696c72-6542-2067-7265-42206c756150"
version = "1.1.2"

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

[[ArnoldiMethod]]
deps = ["LinearAlgebra", "Random", "StaticArrays"]
git-tree-sha1 = "62e51b39331de8911e4a7ff6f5aaf38a5f4cc0ae"
uuid = "ec485272-7323-5ecc-a04f-4719b315124d"
version = "0.2.0"

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

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

[[BinaryProvider]]
deps = ["Libdl", "Logging", "SHA"]
git-tree-sha1 = "ecdec412a9abc8db54c0efc5548c64dfce072058"
uuid = "b99e7846-7c00-51b0-8f62-c81ae34c0232"
version = "0.5.10"

[[Bio3DView]]
deps = ["Requires"]
git-tree-sha1 = "7f472efd9b6af772307dd017f9deeff2a243754f"
uuid = "99c8bb3a-9d13-5280-9740-b4880ed9c598"
version = "0.1.3"

[[CSV]]
deps = ["Dates", "Mmap", "Parsers", "PooledArrays", "SentinelArrays", "Tables", "Unicode"]
git-tree-sha1 = "b83aa3f513be680454437a0eee21001607e5d983"
uuid = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
version = "0.8.5"

[[CategoricalArrays]]
deps = ["DataAPI", "Future", "JSON", "Missings", "Printf", "Statistics", "StructTypes", "Unicode"]
git-tree-sha1 = "18d7f3e82c1a80dd38c16453b8fd3f0a7db92f23"
uuid = "324d7699-5711-5eae-9e2f-1d82baa6b597"
version = "0.9.7"

[[ChainRulesCore]]
deps = ["Compat", "LinearAlgebra", "SparseArrays"]
git-tree-sha1 = "f885e7e7c124f8c92650d61b9477b9ac2ee607dd"
uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
version = "1.11.1"

[[ChangesOfVariables]]
deps = ["LinearAlgebra", "Test"]
git-tree-sha1 = "9a1d594397670492219635b35a3d830b04730d62"
uuid = "9e997f8a-9a97-42d5-a9f1-ce6bfc15e2c0"
version = "0.1.1"

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

[[Conda]]
deps = ["JSON", "VersionParsing"]
git-tree-sha1 = "299304989a5e6473d985212c28928899c74e9421"
uuid = "8f4d0f93-b110-5947-807f-2305c1781a2d"
version = "1.5.2"

[[Configurations]]
deps = ["ExproniconLite", "OrderedCollections", "TOML"]
git-tree-sha1 = "79e812c535bb9780ba00f3acba526bde5652eb13"
uuid = "5218b696-f38b-4ac9-8b61-a12ec717816d"
version = "0.16.6"

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

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

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

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

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

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

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

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

[[DocStringExtensions]]
deps = ["LibGit2"]
git-tree-sha1 = "b19534d1895d702889b219c382a6e18010797f0b"
uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
version = "0.8.6"

[[Documenter]]
deps = ["ANSIColoredPrinters", "Base64", "Dates", "DocStringExtensions", "IOCapture", "InteractiveUtils", "JSON", "LibGit2", "Logging", "Markdown", "REPL", "Test", "Unicode"]
git-tree-sha1 = "f425293f7e0acaf9144de6d731772de156676233"
uuid = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
version = "0.27.10"

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

[[ExproniconLite]]
git-tree-sha1 = "8b08cc88844e4d01db5a2405a08e9178e19e479e"
uuid = "55351af7-c7e9-48d6-89ff-24e801d99491"
version = "0.6.13"

[[FIGlet]]
deps = ["BinaryProvider", "Pkg"]
git-tree-sha1 = "bfc6b52f75b4720581e3e49ae786da6764e65b6a"
uuid = "3064a664-84fe-4d92-92c7-ed492f3d8fae"
version = "0.2.1"

[[FileIO]]
deps = ["Pkg", "Requires", "UUIDs"]
git-tree-sha1 = "2db648b6712831ecb333eae76dbfd1c156ca13bb"
uuid = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
version = "1.11.2"

[[FileWatching]]
uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee"

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

[[FromFile]]
git-tree-sha1 = "81e918d0ed5978fcdacd06b7c64c0c5074c4d55a"
uuid = "ff7dd447-1dcb-4ce3-b8ac-22a812192de7"
version = "0.1.2"

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

[[FuzzyCompletions]]
deps = ["REPL"]
git-tree-sha1 = "2cc2791b324e8ed387a91d7226d17be754e9de61"
uuid = "fb4132e2-a121-4a70-b8a1-d5b831dcdcc2"
version = "0.4.3"

[[GitHubActions]]
deps = ["JSON", "Logging"]
git-tree-sha1 = "56e01ec63d13e1cf015d9ff586156eae3cc7cd6f"
uuid = "6b79fd1a-b13a-48ab-b6b0-aaee1fee41df"
version = "0.1.4"

[[Graphs]]
deps = ["ArnoldiMethod", "DataStructures", "Distributed", "Inflate", "LinearAlgebra", "Random", "SharedArrays", "SimpleTraits", "SparseArrays", "Statistics"]
git-tree-sha1 = "92243c07e786ea3458532e199eb3feee0e7e08eb"
uuid = "86223c79-3864-5bf0-83f7-82e725a168b6"
version = "1.4.1"

[[HTTP]]
deps = ["Base64", "Dates", "IniFile", "Logging", "MbedTLS", "NetworkOptions", "Sockets", "URIs"]
git-tree-sha1 = "0fa77022fe4b511826b39c894c90daf5fce3334a"
uuid = "cd3eb016-35fb-5094-929b-558a96fad6f3"
version = "0.9.17"

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

[[HypertextLiteral]]
git-tree-sha1 = "2b078b5a615c6c0396c77810d92ee8c6f470d238"
uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
version = "0.9.3"

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

[[Inflate]]
git-tree-sha1 = "f5fc07d4e706b84f72d54eedcc1c13d92fb0871c"
uuid = "d25df0c9-e2be-5dd7-82c8-3ad0b3e990b9"
version = "0.1.2"

[[IniFile]]
deps = ["Test"]
git-tree-sha1 = "098e4d2c533924c921f9f9847274f2ad89e018b8"
uuid = "83e8ac13-25f8-5344-8a64-a9f2b223428f"
version = "0.5.0"

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

[[InverseFunctions]]
deps = ["Test"]
git-tree-sha1 = "a7254c0acd8e62f1ac75ad24d5db43f5f19f3c65"
uuid = "3587e190-3f89-42d0-90ee-14403ec27112"
version = "0.1.2"

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

[[IrrationalConstants]]
git-tree-sha1 = "7fd44fd4ff43fc60815f8e764c0f352b83c49151"
uuid = "92d709cd-6900-40b7-9082-c6be49f344b6"
version = "0.1.1"

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

[[JLD2]]
deps = ["DataStructures", "FileIO", "MacroTools", "Mmap", "Pkg", "Printf", "Reexport", "TranscodingStreams", "UUIDs"]
git-tree-sha1 = "46b7834ec8165c541b0b5d1c8ba63ec940723ffb"
uuid = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
version = "0.4.15"

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

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

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

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

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

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

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

[[LogExpFunctions]]
deps = ["ChainRulesCore", "ChangesOfVariables", "DocStringExtensions", "InverseFunctions", "IrrationalConstants", "LinearAlgebra"]
git-tree-sha1 = "be9eef9f9d78cecb6f262f3c10da151a6c5ab827"
uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
version = "0.3.5"

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

[[MacroTools]]
deps = ["Markdown", "Random"]
git-tree-sha1 = "3d3e902b31198a27340d0bf00d6ac452866021cf"
uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
version = "0.5.9"

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

[[MbedTLS]]
deps = ["Dates", "MbedTLS_jll", "Random", "Sockets"]
git-tree-sha1 = "1c38e51c3d08ef2278062ebceade0e46cefc96fe"
uuid = "739be429-bea8-5141-9913-cc70e7f3736d"
version = "1.0.3"

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

[[MetaGraphs]]
deps = ["Graphs", "JLD2", "Random"]
git-tree-sha1 = "2af69ff3c024d13bde52b34a2a7d6887d4e7b438"
uuid = "626554b9-1ddb-594c-aa3c-2596fe9399a5"
version = "0.7.1"

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

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

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

[[MsgPack]]
deps = ["Serialization"]
git-tree-sha1 = "a8cbf066b54d793b9a48c5daa5d586cf2b5bd43d"
uuid = "99f44e22-a591-53d1-9472-aa23ef4bd671"
version = "1.1.0"

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

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

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

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

[[Pluto]]
deps = ["Base64", "Configurations", "Dates", "Distributed", "FileWatching", "FuzzyCompletions", "HTTP", "InteractiveUtils", "Logging", "Markdown", "MsgPack", "Pkg", "REPL", "Sockets", "TableIOInterface", "Tables", "UUIDs"]
git-tree-sha1 = "a5b3fee95de0c0a324bab53a03911395936d15d9"
uuid = "c3e4b0f8-55cb-11ea-2926-15256bba5781"
version = "0.17.2"

[[PlutoSliderServer]]
deps = ["Base64", "Configurations", "Distributed", "FromFile", "GitHubActions", "HTTP", "Logging", "Pkg", "Pluto", "SHA", "Sockets", "TOML", "UUIDs"]
git-tree-sha1 = "ed9660bb2c9eee9d389601bd80a10cee3dd64f0b"
uuid = "2fc8631c-6f24-4c5b-bca7-cbb509c42db4"
version = "0.2.7"

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

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

[[PoreMatMod]]
deps = ["Bio3DView", "CSV", "DataFrames", "FIGlet", "Graphs", "LinearAlgebra", "MetaGraphs", "PlutoSliderServer", "PlutoUI", "Reexport", "StatsBase", "Xtals"]
git-tree-sha1 = "e4b40bedba7a4aadefe93b52b809de10539bc915"
uuid = "2de0d7f0-0963-4438-8bc8-7e7ffe3dc69a"
version = "0.2.6"

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

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

[[PyCall]]
deps = ["Conda", "Dates", "Libdl", "LinearAlgebra", "MacroTools", "Serialization", "VersionParsing"]
git-tree-sha1 = "4ba3651d33ef76e24fef6a598b63ffd1c5e1cd17"
uuid = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0"
version = "1.92.5"

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

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

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

[[Requires]]
deps = ["UUIDs"]
git-tree-sha1 = "4036a3bd08ac7e968e27c203d45f5fff15020621"
uuid = "ae029012-a4dd-5104-9daa-d747884805df"
version = "1.1.3"

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

[[SentinelArrays]]
deps = ["Dates", "Random"]
git-tree-sha1 = "f45b34656397a1f6e729901dc9ef679610bd12b5"
uuid = "91c51154-3ec4-41a3-a24f-3f23e20d615c"
version = "1.3.8"

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

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

[[SimpleTraits]]
deps = ["InteractiveUtils", "MacroTools"]
git-tree-sha1 = "5d7e3f4e11935503d3ecaf7186eac40602e7d231"
uuid = "699a6c99-e7fa-54fc-8d76-47d257e15c1d"
version = "0.9.4"

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

[[SortingAlgorithms]]
deps = ["DataStructures", "Random", "Test"]
git-tree-sha1 = "03f5898c9959f8115e30bc7226ada7d0df554ddd"
uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c"
version = "0.3.1"

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

[[StaticArrays]]
deps = ["LinearAlgebra", "Random", "Statistics"]
git-tree-sha1 = "3c76dde64d03699e074ac02eb2e8ba8254d428da"
uuid = "90137ffa-7385-5640-81b9-e52037218182"
version = "1.2.13"

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

[[StatsAPI]]
git-tree-sha1 = "0f2aa8e32d511f758a2ce49208181f7733a0936a"
uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0"
version = "1.1.0"

[[StatsBase]]
deps = ["DataAPI", "DataStructures", "LinearAlgebra", "LogExpFunctions", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics", "StatsAPI"]
git-tree-sha1 = "2bb0cb32026a66037360606510fca5984ccc6b75"
uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
version = "0.33.13"

[[StructTypes]]
deps = ["Dates", "UUIDs"]
git-tree-sha1 = "d24a825a95a6d98c385001212dc9020d609f2d4f"
uuid = "856f2bd8-1eba-4b0a-8007-ebc267875bd4"
version = "1.8.1"

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

[[TableIOInterface]]
git-tree-sha1 = "9a0d3ab8afd14f33a35af7391491ff3104401a35"
uuid = "d1efa939-5518-4425-949f-ab857e148477"
version = "0.1.6"

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

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

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

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

[[TranscodingStreams]]
deps = ["Random", "Test"]
git-tree-sha1 = "216b95ea110b5972db65aa90f88d8d89dcb8851c"
uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa"
version = "0.9.6"

[[URIs]]
git-tree-sha1 = "97bbe755a53fe859669cd907f2d96aee8d2c1355"
uuid = "5c2747f8-b7ea-4ff2-ba2e-563bfd36b1d4"
version = "1.3.0"

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

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

[[VersionParsing]]
git-tree-sha1 = "e575cf85535c7c3292b4d89d89cc29e8c3098e47"
uuid = "81def892-9a0e-5fdd-b105-ffc91e053289"
version = "1.2.1"

[[Xtals]]
deps = ["Bio3DView", "CSV", "DataFrames", "Documenter", "FIGlet", "Graphs", "JLD2", "LinearAlgebra", "Logging", "MetaGraphs", "Printf", "PyCall", "UUIDs"]
git-tree-sha1 = "3147503cd35c4f2b3744fe36301c7de3efee98c5"
uuid = "ede5f01d-793e-4c47-9885-c447d1f18d6d"
version = "0.3.9"

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

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

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

# ╔═╡ Cell order:
# ╟─6c1969e0-02f5-11eb-3fa2-09931a63b1ac
# ╟─50269ffe-02ef-11eb-0614-f11975d991fe
# ╟─33b1fb50-0f73-11eb-2ab2-9d2cb6c5a533
# ╟─415e9210-0f71-11eb-15c8-e7484b5be309
# ╟─3997c4d0-0f75-11eb-2976-c161879b8d0c
# ╟─69edca20-0f94-11eb-13ba-334438ca2406
# ╟─5918f770-103d-11eb-0537-81036bd3e675
# ╟─31832e30-1054-11eb-24ed-219fd3e236a1
# ╟─5dc43a20-10b8-11eb-26dc-7fb98e9aeb1a
# ╟─90696d20-10b7-11eb-20b5-6174faeaf613
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
using DataFrames, DataArrays

token  = getRepositoryToken(mPulseAPITenant, mPulseAPIToken)
domain = getRepositoryDomain(token, appName="mPulse Demo")
appKey = domain["attributes"]["appKey"]

# Summary
# We need to add a random string to the appKey to cache bust API results from the previous test
@test_throws mPulseAPIAuthException mPulseAPI.getSummaryTimers("invalid-token", appKey * "-" * base(16, round(Int, rand()*100000), 5))
@test_throws mPulseAPIRequestException mPulseAPI.getSummaryTimers(token, appKey, filters = Dict("invalid-token" => "foo"))

summary = mPulseAPI.getSummaryTimers(token, appKey)

@test length(summary) == 5

hasDemoData = true

if summary["n"] != nothing
    @test isa(summary["n"], Int)
    @test isa(summary["median"], Int)
    @test isa(summary["p95"], Int)
    @test isa(summary["p98"], Int)
    @test isa(summary["moe"], Float64)
else
    warn("mPulse Demo app has no demo data")
    hasDemoData = false
end

function testDimensionTable(method, first_symbol, first_friendly)
    local x = "--"
    try
        x = getfield(mPulseAPI, method)(token, appKey)

        @test size(x, 2) == 5
        @test names(x) == [first_symbol, :t_done_median, :t_done_moe, :t_done_count, :t_done_total_pc]

        x = getfield(mPulseAPI, method)(token, appKey, friendly_names=true)
        @test size(x, 2) == 5
        @test names(x) == [Symbol(first_friendly), Symbol("Median Time (ms)"), Symbol("MoE (ms)"), Symbol("Measurements"), Symbol("% of total")]
    catch ex
        warn("mPulseAPI.$method")
        show(x)
        println()
        if :response ∈ fieldnames(ex) && isa(ex.response, HttpCommon.Response)
            show(ex.response)
            println()
        end
        rethrow(ex)
    end
end

# Page Groups
testDimensionTable(:getPageGroupTimers, :page_group, "Page Group")

# Browsers
testDimensionTable(:getBrowserTimers, :user_agent, "User Agent")

# ABTests
testDimensionTable(:getABTestTimers, :test_name, "Test Name")

# Geo
testDimensionTable(:getGeoTimers, :country, "Country")

# MetricsByDimension

for dimension in ["browser", "page_group", "country", "bw_block", "ab_test"]
    local metrics = []
    try
        metrics = mPulseAPI.getMetricsByDimension(token, appKey, dimension)

        @test size(metrics, 2) == 1 + length(domain["custom_metrics"])
        @test sort(names(metrics)) == sort(map(Symbol, [ dimension; collect(keys(domain["custom_metrics"])) ]))

        hasDemoData && @test size(metrics, 1) > 0
    catch ex
        warn("mPulseAPI.getMetricsByDimension:$dimension")
        show(metrics)
        println()
        if :response ∈ fieldnames(ex) && isa(ex.response, HttpCommon.Response)
            show(ex.response)
            println()
        end
        rethrow(ex)
    end
end

@test_throws mPulseAPIRequestException mPulseAPI.getMetricsByDimension(token, appKey, "some-fake-dimension")

# TimersMetrics
tm = mPulseAPI.getTimersMetrics(token, appKey)

fixed_cols = [:Beacons, :PageLoad] # mPulse now only returns Beacons & PageLoad by default
# fixed_cols = [:Beacons, :PageLoad, :Sessions, :BounceRate, :DNS, :TCP, :SSL, :FirstByte, :DomLoad, :DomReady, :FirstLastByte]
varia_cols = Symbol[] #map(Symbol, [collect(keys(domain["custom_timers"])); collect(keys(domain["custom_metrics"]))])

@test size(tm, 2) == length(fixed_cols ∪ varia_cols)

@test sort(names(tm)) == sort(fixed_cols ∪ varia_cols)

hasDemoData && @test size(tm, 1) == 1441

if hasDemoData
    tm = mPulseAPI.getTimersMetrics(token, appKey, filters=Dict("page-group" => ["Account", "Search"]))

    fixed_cols = [:Beacons, :PageLoad] # mPulse now only returns Beacons & PageLoad by default
    # fixed_cols = [:Beacons, :PageLoad, :Sessions, :BounceRate, :DNS, :TCP, :SSL, :FirstByte, :DomLoad, :DomReady, :FirstLastByte]
    varia_cols = Symbol[] #map(Symbol, [collect(keys(domain["custom_timers"])); collect(keys(domain["custom_metrics"]))])

    @test size(tm, 2) >= length(fixed_cols)

    @test size(tm, 1) == 1441
end

# Histogram

hgm = mPulseAPI.getHistogram(token, appKey)

@test length(hgm) == 4

if hasDemoData
    @test isa(hgm["median"], Int)
    @test isa(hgm["p95"], Int)
    @test isa(hgm["p98"], Int)
end

@test isa(hgm["buckets"], DataFrame)

@test size(hgm["buckets"], 2) == 3
@test names(hgm["buckets"]) == [:bucket_start, :bucket_end, :element_count]
for col in names(hgm["buckets"])
    @test isa(hgm["buckets"][col], DataArrays.DataArray{Real,1})
end

# Sessions/Metrics OverPageLoadTime
metric_frames = []
metrics = [
    (:getSessionsOverPageLoadTime, :Sessions),
    (:getMetricOverPageLoadTime, :BounceRate)
] ∪ map(m -> (:getMetricOverPageLoadTime, Symbol(m), m), collect(keys(domain["custom_metrics"])))
for tuple in metrics
    local x = nothing
    try
        if length(tuple) == 3
            x = getfield(mPulseAPI, tuple[1])(token, appKey, metric=tuple[3])
        else
            x = getfield(mPulseAPI, tuple[1])(token, appKey)
        end

        push!(metric_frames, x)

        @test size(x, 2) == 2
        @test names(x) == [:t_done, tuple[2]]

        hasDemoData && @test size(x, 1) > 0
    catch ex
        warn("mPulseAPI.$(tuple[1])", length(tuple) == 3 ? ":$(tuple[3])" : "")
        show(x)
        println()
        if :response ∈ fieldnames(ex) && isa(ex.response, HttpCommon.Response)
            show(ex.response)
            println()
        end

        rethrow(ex)
    end
end

# MergeMetrics
merged_frame = mPulseAPI.mergeMetrics(metric_frames...)

@test size(merged_frame, 2) == 3 + length(domain["custom_metrics"])
@test names(merged_frame) == [:t_done, :Sessions, :BounceRate] ∪ map(Symbol, collect(keys(domain["custom_metrics"])))
hasDemoData && @test size(merged_frame, 1) > 0

# TimerByMinute
for timer in mPulseAPI.supported_timers ∪ collect(keys(domain["custom_timers"]))
    local tbm = []
    try
        tbm = mPulseAPI.getTimerByMinute(token, appKey, timer=timer)

        if timer == "ResourceTimer" || !hasDemoData
            @test size(tbm) == (0, 3)                                   # ResourceTimer has no data
        else
            @test size(tbm) == (1440, 3) || size(tbm) == (1439, 3)      # The mPulse API will sometimes not return the last minute of the day
        end
        @test names(tbm) == [:timestamp, Symbol(timer), :moe]
    catch ex
        warn("mPulseAPI.getTimerByMinute($timer)")
        show(tbm)
        println()
        if :response ∈ fieldnames(ex) && isa(ex.response, HttpCommon.Response)
            show(ex.response)
            println()
        end

        rethrow(ex)
    end
end

@test_throws mPulseAPIRequestException mPulseAPI.getTimerByMinute(token, appKey, timer="Some-Bad-Timer")


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

"""
    BoundedSearch(method, nmax)

A method for searching at most `nmax` neighbors using `method`.
"""
struct BoundedSearch{M<:NeighborSearchMethod} <: BoundedNeighborSearchMethod
  method::M
  nmax::Int
end

object(method::BoundedSearch) = object(method.method)

maxneighbors(method::BoundedSearch) = method.nmax

function search!(neighbors, xₒ::AbstractVector,
                 method::BoundedSearch; mask=nothing)
  inds = search(xₒ, method.method)
  nmax = method.nmax

  if isnothing(mask)
    nneigh = min(length(inds), nmax)
    @inbounds for i in 1:nneigh
      neighbors[i] = inds[i]
    end
  else
    nneigh = 0
    @inbounds for ind in inds
      if mask[ind]
        nneigh += 1
        neighbors[nneigh] = ind
      end
      nneigh == nmax && break
    end
  end

  nneigh
end


text:
module DDL
using HTTP, ZipFile, Tar
end


text:
"""

Exported names are
$(EXPORTS)
"""
module Geometry

using Test
using DocStringExtensions
using FillArrays

using LinearAlgebra: ⋅

using Gridap.Helpers
using Gridap.Arrays
using Gridap.Fields
using Gridap.ReferenceFEs
using Gridap.TensorValues
using Gridap.Io

using Gridap.ReferenceFEs: _num_faces
using Gridap.ReferenceFEs: _num_facets
using Gridap.ReferenceFEs: _num_vertices
using Gridap.ReferenceFEs: _num_edges
using Gridap.ReferenceFEs: _get_facedims
using Gridap.ReferenceFEs: _get_offsets
using Gridap.ReferenceFEs: _get_offset
using Gridap.ReferenceFEs: _find_unique_with_indices

import Gridap.Arrays: array_cache
import Gridap.Arrays: getindex!
import Gridap.Arrays: get_array
import Gridap.Arrays: lazy_append
import Gridap.Arrays: lazy_map

import Gridap.Arrays: return_cache
import Gridap.Arrays: evaluate!
import Gridap.Arrays: get_children

import Gridap.Io: to_dict
import Gridap.Io: from_dict

import Gridap.ReferenceFEs: ReferenceFE
import Gridap.ReferenceFEs: get_node_coordinates
import Gridap.ReferenceFEs: num_nodes
import Gridap.ReferenceFEs: is_first_order
import Gridap.ReferenceFEs: get_faces
import Gridap.ReferenceFEs: get_face_vertices
import Gridap.ReferenceFEs: get_dimranges
import Gridap.ReferenceFEs: get_dimrange
import Gridap.ReferenceFEs: num_faces
import Gridap.ReferenceFEs: num_facets
import Gridap.ReferenceFEs: num_edges
import Gridap.ReferenceFEs: num_vertices
import Gridap.ReferenceFEs: get_facedims
import Gridap.ReferenceFEs: get_offsets
import Gridap.ReferenceFEs: get_offset
import Gridap.ReferenceFEs: get_vertex_coordinates
import Gridap.ReferenceFEs: get_face_nodes
import Gridap.ReferenceFEs: get_face_own_nodes
import Gridap.ReferenceFEs: get_polytope
import Gridap.ReferenceFEs: get_vertex_node
import Gridap.ReferenceFEs: is_simplex
import Gridap.ReferenceFEs: is_n_cube
import Gridap.ReferenceFEs: get_reffaces
import Gridap.ReferenceFEs: get_face_type
import Gridap.ReferenceFEs: num_dims
import Gridap.ReferenceFEs: num_cell_dims
import Gridap.ReferenceFEs: num_point_dims
import Gridap.ReferenceFEs: simplexify
import Gridap.ReferenceFEs: get_facet_normal
import Gridap.ReferenceFEs: Quadrature

export GridTopology
export num_cells
export get_polytopes
export get_cell_type
export get_cell_vertices
export get_reffaces_offsets
export compute_reffaces
export compute_cell_faces
export compute_face_vertices
export compute_isboundary_face
export get_cell_permutations
export compute_cell_permutations
export test_grid_topology
export get_cell_faces
export get_isboundary_face
export OrientationStyle
export Oriented
export NonOriented
export RegularityStyle
export Regular
export Irregular
export is_oriented
export is_regular
export expand_cell_data
export compress_cell_data
export compress_contributions
export compress_ids

export UnstructuredGridTopology

export Triangulation
export get_reffes
export get_cell_coordinates
export get_cell_ref_coordinates
export get_cell_reffe
export get_cell_shapefuns
export get_facet_normal
export test_triangulation
export get_cell_map

export get_background_model
export get_active_model
export get_grid
export get_glue

export Grid
export get_cell_node_ids
export test_grid
export compute_linear_grid
export compute_reference_grid

export GridPortion
export UnstructuredGrid

export CartesianGrid
export CartesianDescriptor
export get_cartesian_descriptor

export FaceLabeling
export num_tags
export num_entities
export get_face_entity
export get_tag_entities
export get_tag_name
export get_tag_from_name
export get_tags_from_names
export add_tag!
export add_tag_from_tags!
export get_face_mask
export get_face_tag
export get_face_tag_index
export get_cell_entity

export DiscreteModel
export DiscreteModelFromFile
export get_grid
export get_grid_topology
export get_face_labeling
export test_discrete_model
export compute_face_nodes
export compute_face_own_nodes
export compute_vertex_node
export get_node_face_owner
export compute_node_face_owner
export get_triangulation

export UnstructuredDiscreteModel
export CartesianDiscreteModel

export GenericTriangulation
export BoundaryTriangulation
export DiscreteModelPortion

export Interior
export Boundary
export Skeleton
export Interface

export SkeletonPair
export SkeletonTriangulation
export InterfaceTriangulation
#export get_left_boundary
#export get_right_boundary

export get_parent_model

export AppendedTriangulation

export GridMock

export extend
export pos_neg_data
export is_change_possible
export best_target
export FaceToFaceGlue
export BodyFittedTriangulation
export move_contributions

include("Grids.jl")

include("GridMocks.jl")

include("UnstructuredGrids.jl")

include("CartesianGrids.jl")

include("GridTopologies.jl")

include("GridTopologyMocks.jl")

include("UnstructuredGridTopologies.jl")

include("FaceLabelings.jl")

include("DiscreteModels.jl")

include("DiscreteModelMocks.jl")

include("UnstructuredDiscreteModels.jl")

include("CartesianDiscreteModels.jl")

include("GridPortions.jl")

include("DiscreteModelPortions.jl")

include("Triangulations.jl")

include("BoundaryTriangulations.jl")

include("SkeletonTriangulations.jl")

include("AppendedTriangulations.jl")

include("CompressedCellArrays.jl")

end # module


text:
using LinearAlgebra: cross, norm
import AstroTime: seconds

export kepler, period

function period(a, μ)
    2π * sqrt(abs(a)^3/μ) * seconds
end

function kepler(μ, r₀, v₀, Δt, numiter=50, rtol=sqrt(eps()))
    if abs(Δt) < rtol
        return r₀, v₀
    end

    rm = norm(r₀)
    rdv = r₀ ⋅ v₀
    vdv = v₀ ⋅ v₀
    α = -vdv / μ + 2 /rm

    # Elliptic orbit
    if α > rtol
        # For α == 1 the first guess will be too close to converge
        if α ≈ 1
            χ₁ = √μ * Δt * α * 0.97
        else
            χ₁ = √μ * Δt * α
        end
    # Parabolic orbit
    elseif abs(α) < rtol
        h = cross(r₀, v₀)
        hm = norm(h)
        p = hm^2 / μ
        s = 0.5 * (π/2 - atan(3 * sqrt(μ / p^3) * Δt))
        w = atan(tan(s)^(1/3))
        χ₁ = √p * (2 * cot(2 * w))
        # TODO: Check if this is necessary.
        #= α = 0.0 =#
    # Hyperbolic orbit
    else
        a = 1/α
        χ₁ = (sign(Δt) * sqrt(-a)
            * log(-2 * μ * Δt
            / (a * (rdv + sign(Δt) * sqrt(-μ * α) * (1 - rm * alpha)))))
    end

    counter = 0
    converged = false
    χ = 0.0
    c₂ = 0.0
    c₃ = 0.0
    r = 0.0
    ψ = 0.0
    while counter < numiter
        counter += 1
        χ = χ₁
        χ² = χ^2
        ψ = χ² * α
        c₂ = c2(ψ)
        c₃ = c3(ψ)
        r = χ²*c₂ + rdv / √μ * χ * (1 - ψ * c₃) + rm * (1 - ψ * c₂)
        δt = χ^3 * c₃ + rdv / √μ * χ² * c₂ + rm * χ * (1 - ψ * c₃)
        χ₁ = χ + (Δt * √μ - δt) / r
        if abs(χ - χ₁) < rtol
            converged = true
            break
        end
    end

    if !converged
        error("Not converged.")
    end

    f = 1 - χ^2 / rm * c₂
    g = Δt - χ^3 / √μ * c₃
    fdot = √μ / (r * rm) * χ * (ψ * c₃ - 1)
    gdot = 1 - χ^2 / r * c₂
    @assert f * gdot - fdot * g ≈ 1

    return f * r₀ + g * v₀, fdot * r₀ + gdot * v₀
end


text:
@testset "1.1.4 Improper" begin
    include("1.1.4.2 (c x)^m (a x^j+b x^n)^p.jl")
    include("1.1.4.3 (e x)^m (a x^j+b x^k)^p (c+d x^n)^q.jl")
end

max_stars_repo_path: test/1 Algebraic functions/1.1 Binomial products/1.1.4 Improper/.jl


text:
export Mesh

struct Mesh

    xmin :: Float64
    xmax :: Float64
    nx   :: Int
    dx   :: Float64
    ymin :: Float64
    ymax :: Float64
    ny   :: Int
    dy   :: Float64

    function Mesh( xmin, xmax, nx, ymin, ymax, ny )

        dx = (xmax - xmin) / nx
        dy = (ymax - ymin) / ny

        new( xmin, xmax, nx, dx, ymin, ymax, ny, dy )

    end

end

export MeshFields

struct MeshFields

    mesh :: Mesh

    e :: Array{Float64,3}
    ρ :: Array{Float64,2}

    function MeshFields( mesh :: Mesh )

	     nx, ny = mesh.nx, mesh.ny

         e = zeros(Float64, (2, nx+1, ny+1))
         ρ = zeros(Float64, (nx+1, ny+1))

         new( mesh, e, ρ)

    end

end


text:
#=##############################################################################
# DESCRIPTION
Test rotor model simulating an isolated 9.4in rotor in hover. The rotor and
configuration matches the 3D-printed propeller described in Ning, Z.,
*Experimental investigations on the aerodynamic and aeroacoustic characteristics
of small UAS propellers*, Sec. 5.2. This rotor roughly resembles a DJI Phantom
II rotor.

# AUTHORSHIP
  * Author    : Eduardo J. Alvarez
  * Email     : Edo.AlvarezR@gmail.com
  * Created   : Dec 2019
  * License   : MIT
=###############################################################################

# ------------ MODULES ---------------------------------------------------------
# Load simulation engine
import FLOWUnsteady
uns = FLOWUnsteady
vlm = uns.vlm

import GeometricTools
gt = GeometricTools

using PyPlot

# ------------ GLOBAL VARIABLES ------------------------------------------------
# Default path where to save data
extdrive_path = "/media/edoalvar/MyExtDrive/simulationdata7/"
# extdrive_path = "temps/"

# ------------ DRIVERS ---------------------------------------------------------
function run_singlerotor_hover(; xfoil=true, prompt=true)

    J = 0.00                # Advance ratio Vinf/(nD)
    angle = 0.0             # (deg) angle of freestream (0 == climb, 90==forward flight)

    singlerotor(;   xfoil=xfoil,
                    VehicleType=uns.VLMVehicle,
                    J=J,
                    DVinf=[cos(pi/180*angle), sin(pi/180*angle), 0],
                    save_path=extdrive_path*"singlerotor_hover_test00/",
                    prompt=prompt)
end

function run_singlerotor_forwardflight(; xfoil=true, prompt=true)

    J = 0.15                # Advance ratio Vinf/(nD)
    angle = 60.0            # (deg) angle of freestream (0 == climb, ~90==forward flight)

    singlerotor(;   xfoil=xfoil,
                    VehicleType=uns.VLMVehicle,
                    J=J,
                    DVinf=[cos(pi/180*angle), sin(pi/180*angle), 0],
                    nrevs=2,
                    nsteps_per_rev=120,
                    save_path=extdrive_path*"singlerotor_fflight_test01/",
                    prompt=prompt)
end

# ------------------------------------------------------------------------------

function singlerotor(;  xfoil       = true,             # Whether to run XFOIL
                        VehicleType = uns.VLMVehicle,   # Vehicle type
                        J           = 0.0,              # Advance ratio
                        DVinf       = [1.0, 0, 0],      # Freestream direction
                        nrevs       = 6,                # Number of revolutions
                        nsteps_per_rev = 72,            # Time steps per revolution
                        shed_unsteady = true,
                        lambda      = 2.125,
                        n           = 10,
                        overwrite_overwrite_sigma = nothing,
                        # OUTPUT OPTIONS
                        save_path   = nothing,
                        run_name    = "singlerotor",
                        prompt      = true,
                        verbose     = true,
                        v_lvl       = 0,
                        rotor_file = "DJI-II.csv",           # Rotor geometry
                        optargs...)

    # TODO: Wake removal ?

    # ------------ PARAMETERS --------------------------------------------------

    # Rotor geometry
    data_path = uns.def_data_path       # Path to rotor database
    pitch = 0.0                         # (deg) collective pitch of blades
    # n = 50                              # Number of blade elements
    # n = 10
    CW = false                          # Clock-wise rotation
    # xfoil = false                     # Whether to run XFOIL

    # Read radius of this rotor and number of blades
    R, B = uns.read_rotor(rotor_file; data_path=data_path)[[1,3]]

    # Simulation parameters
    RPM = 81*60                         # RPM
    # J = 0.00001                       # Advance ratio Vinf/(nD)
    rho = 1.225                         # (kg/m^3) air density
    mu = 1.81e-5                        # (kg/ms) air dynamic viscosity
    ReD = 2*pi*RPM/60*R * rho/mu * 2*R  # Diameter-based Reynolds number

    magVinf = J*RPM/60*(2*R)
    Vinf(X,t) = magVinf*DVinf           # (m/s) freestream velocity

    # Solver parameters
    # nrevs = 6                         # Number of revolutions in simulation
    # nsteps_per_rev = 72                 # Time steps per revolution
    p_per_step = 2                      # Sheds per time step
    ttot = nrevs/(RPM/60)               # (s) total simulation time
    nsteps = nrevs*nsteps_per_rev       # Number of time steps
    # lambda = 2.125                      # Core overlap
    overwrite_sigma = overwrite_overwrite_sigma != nothing ? overwrite_overwrite_sigma :
                                        lambda * 2*pi*R/(nsteps_per_rev*p_per_step) # Smoothing core size

    surf_sigma = R/10                   # Smoothing radius of lifting surface
    vlm_sigma = surf_sigma              # Smoothing radius of VLM
    # shed_unsteady = true                # Shed particles from unsteady loading

    max_particles = ((2*n+1)*B)*nrevs*nsteps_per_rev*p_per_step # Max particles for memory pre-allocation
    plot_disc = true                    # Plot blade discretization for debugging

    # ------------ SIMULATION SETUP --------------------------------------------
    # Generate rotor
    rotor = uns.generate_rotor(rotor_file; pitch=pitch,
                                            n=n, CW=CW, ReD=ReD,
                                            verbose=verbose, xfoil=xfoil,
                                            data_path=data_path,
                                            plot_disc=plot_disc)
    # ----- VEHICLE DEFINITION
    # System of all FLOWVLM objects
    system = vlm.WingSystem()
    vlm.addwing(system, run_name, rotor)

    # Systems of rotors
    rotors = vlm.Rotor[rotor]   # Defining this rotor as its own system
    rotor_systems = (rotors,)

    # Wake-shedding system (doesn't include the rotor if quasi-steady vehicle)
    wake_system = vlm.WingSystem()

    if VehicleType != uns.QVLMVehicle
        vlm.addwing(wake_system, run_name, rotor)
    else
        # Mute colinear warnings. This is needed since the quasi-steady solver
        #   will probe induced velocities at the lifting line of the blade
        uns.vlm.VLMSolver._mute_warning(true)
    end

    # FVS's Vehicle object
    vehicle = VehicleType(   system;
                                rotor_systems=rotor_systems,
                                wake_system=wake_system
                             )

    # ----- MANEUVER DEFINITION
    RPM_fun(t) = 1.0                # RPM (normalized by reference RPM) as a
                                    # function of normalized time

    angle = ()                      # Angle of each tilting system (none in this case)
    sysRPM = (RPM_fun, )              # RPM of each rotor system
    Vvehicle(t) = zeros(3)          # Translational velocity of vehicle over Vcruise
    anglevehicle(t) = zeros(3)      # (deg) angle of the vehicle

    # FVS's Maneuver object
    maneuver = uns.KinematicManeuver(angle, sysRPM, Vvehicle, anglevehicle)

    # Plot maneuver path and controls
    uns.plot_maneuver(maneuver; vis_nsteps=nsteps)

    # ----- SIMULATION DEFINITION
    RPMref = RPM
    Vref = 0.0
    simulation = uns.Simulation(vehicle, maneuver, Vref, RPMref, ttot)

    monitor = generate_monitor(J, rho, RPM, nsteps; save_path=save_path,
                                                            run_name=run_name)

    # ------------ RUN SIMULATION ----------------------------------------------
    pfield = uns.run_simulation(simulation, nsteps;
                                      # SIMULATION OPTIONS
                                      Vinf=Vinf,
                                      # SOLVERS OPTIONS
                                      p_per_step=p_per_step,
                                      overwrite_sigma=overwrite_sigma,
                                      vlm_sigma=vlm_sigma,
                                      surf_sigma=surf_sigma,
                                      max_particles=max_particles,
                                      shed_unsteady=shed_unsteady,
                                      extra_runtime_function=monitor,
                                      # OUTPUT OPTIONS
                                      save_path=save_path,
                                      run_name=run_name,
                                      prompt=prompt,
                                      verbose=verbose, v_lvl=v_lvl,
                                      optargs...
                                      )
    return pfield, rotor
end

"""
Generate monitor for rotor performance parameters
"""
function generate_monitor(J, rho, RPM, nsteps; save_path=nothing,
                            run_name="singlerotor",
                            figname="monitor", disp_conv=true,
                            nsteps_savefig=10)

    fcalls = 0                  # Number of function calls

    colors="rgbcmy"^100
    stls = "o^*.px"^100

    # Name of convergence file
    if save_path!=nothing
        fname = joinpath(save_path, run_name*"_convergence.csv")
    end

    # Function for run_vpm! to call on each iteration
    function extra_runtime_function(sim::uns.Simulation{V, M, R},
                                    PFIELD::uns.vpm.ParticleField,
                                    T::Real, DT::Real
                                   ) where{V<:uns.AbstractVLMVehicle, M, R}

        rotors = vcat(sim.vehicle.rotor_systems...)
        angle = T*360*RPM/60

        # Call figure
        if disp_conv; fig = figure(figname, figsize=(7*3,5*2)); end;

        if fcalls==0
            # Format subplots
            if disp_conv
                subplot(231)
                title("Circulation Distribution")
                xlabel("Element index")
                ylabel(L"Circulation $\Gamma$ (m$^2$/s)")
                grid(true, color="0.8", linestyle="--")
                subplot(232)
                title("Plane-of-rotation Normal Force")
                xlabel("Element index")
                ylabel(L"Normal Force $N_p$ (N)")
                grid(true, color="0.8", linestyle="--")
                subplot(233)
                title("Plane-of-rotation Tangential Force")
                xlabel("Element index")
                ylabel(L"Tangential Force $T_p$ (N)")
                grid(true, color="0.8", linestyle="--")
                subplot(234)
                xlabel(L"Age $\psi$ ($^\circ$)")
                ylabel(L"Thrust $C_T$")
                grid(true, color="0.8", linestyle="--")
                subplot(235)
                xlabel(L"Age $\psi$ ($^\circ$)")
                ylabel(L"Torque $C_Q$")
                grid(true, color="0.8", linestyle="--")
                subplot(236)
                xlabel(L"Age $\psi$ ($^\circ$)")
                ylabel(L"Propulsive efficiency $\eta$")
                grid(true, color="0.8", linestyle="--")
            end

            # Convergence file header
            if save_path!=nothing
                f = open(fname, "w")
                print(f, "age (deg),T,DT")
                for (i, rotor) in enumerate(rotors)
                    print(f, ",RPM_$i,CT_$i,CQ_$i,eta_$i")
                end
                print(f, "\n")
                close(f)
            end
        end

        # Write rotor position and time on convergence file
        if save_path!=nothing
            f = open(fname, "a")
            print(f, angle, ",", T, ",", DT)
        end

        # Plot circulation and loads distributions
        if disp_conv

            cratio = PFIELD.nt/nsteps
            cratio = cratio > 1 ? 1 : cratio
            clr = fcalls==0 && false ? (0,0,0) : (1-cratio, 0, cratio)
            stl = fcalls==0 && false ? "o" : "-"
            alpha = fcalls==0 && false ? 1 : 0.5

            # Circulation distribution
            subplot(231)
            this_sol = []
            for rotor in rotors
                this_sol = vcat(this_sol, [vlm.get_blade(rotor, j).sol["Gamma"] for j in 1:rotor.B]...)
            end
            plot(1:size(this_sol,1), this_sol, stl, alpha=alpha, color=clr)

            # Np distribution
            subplot(232)
            this_sol = []
            for rotor in rotors
                this_sol = vcat(this_sol, rotor.sol["Np"]["field_data"]...)
            end
            plot(1:size(this_sol,1), this_sol, stl, alpha=alpha, color=clr)

            # Tp distribution
            subplot(233)
            this_sol = []
            for rotor in rotors
                this_sol = vcat(this_sol, rotor.sol["Tp"]["field_data"]...)
            end
            plot(1:size(this_sol,1), this_sol, stl, alpha=alpha, color=clr)
        end

        # Plot performance parameters
        for (i,rotor) in enumerate(rotors)
            CT, CQ = vlm.calc_thrust_torque_coeffs(rotor, rho)
            eta = J*CT/(2*pi*CQ)

            if disp_conv
                subplot(234)
                plot([angle], [CT], "$(stls[i])", alpha=alpha, color=clr)
                subplot(235)
                plot([angle], [CQ], "$(stls[i])", alpha=alpha, color=clr)
                subplot(236)
                plot([angle], [eta], "$(stls[i])", alpha=alpha, color=clr)
            end

            if save_path!=nothing
                print(f, ",", rotor.RPM, ",", CT, ",", CQ, ",", eta)
            end
        end

        if disp_conv
            # Save figure
            if fcalls%nsteps_savefig==0 && fcalls!=0 && save_path!=nothing
                savefig(joinpath(save_path, run_name*"_convergence.png"),
                                                            transparent=true)
            end
        end

        # Close convergence file
        if save_path!=nothing
            print(f, "\n")
            close(f)
        end

        fcalls += 1

        return false
    end

    return extra_runtime_function
end


text:
using Voronoi
using Base.Test

points = Point[]
segments = Segment[Segment(Point(0, 0), Point(0, 1)),
                   Segment(Point(0, 1), Point(1, 1)),
                   Segment(Point(1, 1), Point(1, 0)),
                   Segment(Point(1, 0), Point(0, 0))]

vertices, edges =  get_edges(points, segments)

@test length(vertices) == 5
@test length(filter(e -> e.is_primary == true, edges)) == 4

segments = Segment[Segment(Point(0, 0), Point(0, 2)),
                   Segment(Point(0, 2), Point(1, 2)),
                   Segment(Point(1, 2), Point(1, 0)),
                   Segment(Point(1, 0), Point(0, 0))]

vertices, edges =  get_edges(points, segments)

@test length(vertices) == 6
@test length(filter(e -> e.is_primary == true, edges)) == 5


text:
# This file is a part of PropDicts.jl, licensed under the MIT License (MIT).

__precompile__(true)

module PropDicts

import JSON

include("dictmerge.jl")
include("varsubst.jl")
include("propdict.jl")

end # module


text:
"""Poisson Likelihood"""
struct PoissonLikelihood{T<:Real} <: EventLikelihood{T}
    λ::LatentArray{T}
    c::LatentArray{Vector{T}}
    γ::LatentArray{Vector{T}}
    θ::LatentArray{Vector{T}}
    function PoissonLikelihood{T}(λ::AbstractVector{T}) where {T<:Real}
        new{T}(λ)
    end
    function PoissonLikelihood{T}(λ,c,γ,θ) where {T<:Real}
        new{T}(λ,c,γ,θ)
    end
end

function PoissonLikelihood(λ::T=1.0) where {T<:Real}
    PoissonLikelihood{T}([λ])
end

function init_likelihood(likelihood::PoissonLikelihood{T},inference::Inference{T},nLatent::Integer,nSamplesUsed::Int,nFeatures::Int) where T
    PoissonLikelihood{T}(
    [likelihood.λ[1] for _ in 1:nLatent],
    [zeros(T,nSamplesUsed) for _ in 1:nLatent],
    [zeros(T,nSamplesUsed) for _ in 1:nLatent],
    [zeros(T,nSamplesUsed) for _ in 1:nLatent])
end

function pdf(l::PoissonLikelihood,y::Real,f::Real)
    pdf(Poisson(l.λ[1]*logistic(f)),y) #WARNING not valid for multioutput
end

function Base.show(io::IO,model::PoissonLikelihood{T}) where T
    print(io,"Poisson Likelihood")
end

function compute_proba(l::PoissonLikelihood{T},μ::Vector{T},σ²::Vector{T}) where {T<:Real}
    N = length(μ)
    pred = zeros(T,N)
    for i in 1:N
        if σ²[i] <= 0.0
            pred[i] = l.λ[1]*logistic(μ[i]) #WARNING Not valid for multioutput
        else
            nodes = pred_nodes.*sqrt2.*sqrt.(σ²[i]).+μ[i]
            pred[i] =  dot(pred_weights,l.λ[1].*logistic.(nodes)) #WARNING not valid for multioutput
        end
    end
    return pred
end

## Local Updates ##

function local_updates!(model::VGP{T,<:PoissonLikelihood,<:AnalyticVI}) where {T}
    model.likelihood.c .= broadcast((μ,Σ)->sqrt.(abs2.(μ) + Σ) ,model.μ,diag.(model.Σ))
    model.likelihood.γ .= broadcast((c,μ,λ)->0.5*λ*exp.(-0.5*μ)./cosh.(0.5*c),model.likelihood.c,model.μ,model.likelihood.λ)
    model.likelihood.θ .= broadcast((y,γ,c)->(y+γ)./c.*tanh.(0.5*c),model.y,model.likelihood.γ,model.likelihood.c)
    model.likelihood.λ .= broadcast((y,μ)->sum(y)./sum(logistic.(μ)),model.y,model.μ)
end

function local_updates!(model::SVGP{T,<:PoissonLikelihood,<:AnalyticVI}) where {T}
    model.likelihood.c .= broadcast((κ,μ,Σ,K̃)->sqrt.(abs2.(κ*μ) + opt_diag(κ*Σ,κ) + K̃),model.κ,model.μ,model.Σ,model.K̃)
    model.likelihood.γ .= broadcast((c,κμ,λ)->0.5*λ*exp.(-0.5*κμ)./cosh.(0.5*c),model.likelihood.c,model.κ.*model.μ,model.likelihood.λ)
    model.likelihood.θ .= broadcast((y,γ,c)->(y+γ)./c.*tanh.(0.5*c),model.inference.y,model.likelihood.γ,model.likelihood.c)
    model.likelihood.λ .= broadcast((y,κμ)->sum(y)./sum(logistic.(κμ)),model.inference.y,model.κ.*model.μ)
end

## Global Updates ##

@inline ∇E_μ(model::AbstractGP{T,<:PoissonLikelihood,<:GibbsorVI}) where {T} = 0.5.*(model.inference.y.-model.likelihood.γ)
@inline ∇E_μ(model::AbstractGP{T,<:PoissonLikelihood,<:GibbsorVI},i::Int) where {T} = 0.5*(model.inference.y[i]-model.likelihood.γ[i])
@inline ∇E_Σ(model::AbstractGP{T,<:PoissonLikelihood,<:GibbsorVI}) where {T} = 0.5.*model.likelihood.θ
@inline ∇E_Σ(model::AbstractGP{T,<:PoissonLikelihood,<:GibbsorVI},i::Int) where {T} = 0.5*model.likelihood.θ[i]

## ELBO Section ##

function ELBO(model::AbstractGP{T,<:PoissonLikelihood,<:AnalyticVI}) where {T}
    return expecLogLikelihood(model) - GaussianKL(model) - PoissonKL(model) - PolyaGammaKL(model)
end

function expecLogLikelihood(model::VGP{T,<:PoissonLikelihood,<:AnalyticVI}) where {T}
    model.likelihood.c .= broadcast((μ,Σ)->sqrt.(abs2.(μ) + Σ) ,model.μ,diag.(model.Σ))
    tot = sum(broadcast((y,λ,γ)->sum(y*log(λ))-sum(lfactorial.(y))-log(2.0)*sum((y+γ)),model.y,model.likelihood.λ,model.likelihood.γ))
    tot += sum(broadcast((μ,y,γ,c,θ)->0.5*dot(μ,(y-γ))-0.5*dot(c.^2,θ),model.μ,model.inference.y,model.likelihood.γ,model.likelihood.c,model.likelihood.θ))
    return tot
end

function expecLogLikelihood(model::SVGP{T,<:PoissonLikelihood,<:AnalyticVI}) where {T}
    model.likelihood.c .= broadcast((κ,μ,Σ,K̃)->sqrt.(abs2.(κ*μ) + opt_diag(κ*Σ,κ) + K̃),model.κ,model.μ,model.Σ,model.K̃)
    tot = sum(broadcast((y,λ,γ)->sum(y*log(λ))-sum(lfactorial.(y))-log(2.0)*sum(y+γ),model.inference.y,model.likelihood.λ,model.likelihood.γ))
    tot += sum(broadcast((κμ,y,γ,c,θ)->0.5*dot(κμ,(y-γ))-0.5*dot(c.^2,θ),model.κ.*model.μ,model.inference.y,model.likelihood.γ,model.likelihood.c,model.likelihood.θ))
    return model.inference.ρ*tot
end

function PoissonKL(model::AbstractGP{T,<:PoissonLikelihood}) where {T}
    return NaN
    #TODO replace with correct expectations
    model.inference.ρ*sum(broadcast(PoissonKL,model.likelihood.γ,model.likelihood.λ))
end

function PolyaGammaKL(model::VGP{T,<:PoissonLikelihood}) where {T}
    sum(broadcast(PolyaGammaKL,model.y.+model.likelihood.γ,model.likelihood.c,model.likelihood.θ))
end

function PolyaGammaKL(model::SVGP{T,<:PoissonLikelihood}) where {T}
    model.inference.ρ*sum(broadcast(PolyaGammaKL,model.inference.y.+model.likelihood.γ,model.likelihood.c,model.likelihood.θ))
end


text:
const CFNAME_EXT = ".cache.jls"

## ----------------------------------------------------------------------------
function is_cfname(fname)
    fname = basename(fname)
    !isvalid_dfname(fname) && return false
    head, params, ext = parse_dfname(fname)
    return isempty(head) && length(params) == 1 &&
        haskey(params, "hash") && (ext == CFNAME_EXT)
end

## ----------------------------------------------------------------------------
function _cfname(arg, args...)
    _hash = hash(hash.((arg, args...)))
    dfname((;hash = _hash), CFNAME_EXT)
end

_cfname() = _cfname(threadid(), time(), rand())

cfname() = _cfname()
cfname(str::AbstractString) = is_cfname(str) ? str : _cfname(str)
function cfname(arg, args...)
    dir, args = _extract_dir(arg, args...)
    return joinpath(dir, _cfname(args...))
end

## ----------------------------------------------------------------------------
function _scache(dat, cfile;
        onerr::Function = (err) -> rethrow(err),
        print_fun::Function = global_conf(:PRINT_FUN),
        mkdir::Bool = global_conf(:MK_DIR),
        verbose::Bool = global_conf(:VERBOSE),
        msg::String = "",
    )

    try
        mkdir && mkpath(dirname(cfile))
        serialize(cfile, Dict(DATA_KEY => dat))
        verbose && _io_print(print_fun, "CACHE SAVED", msg, dat, cfile)
    catch err
        verbose && _io_error_print(print_fun, err, cfile)
        onerr(err)
    end
    return cfile
end

# defaults
function scache(dat, dirv::Vector{<:AbstractString}, args...; kwargs...)
    dir, args = _extract_dir(dirv, args...)
    cfile = joinpath(dir, cfname(args...))
    _scache(dat, cfile; kwargs...)
end
scache(dat, args...; kwargs...) =
    scache(dat, [global_conf(:CACHE_DIR)], args...; kwargs...)
scache(f::Function, args...; kwargs...) = scache(f(), args...; kwargs...)

## ----------------------------------------------------------------------------
function _lcache(f::Function, savecache::Bool,
        dirv::Vector{<:AbstractString}, arg, args...;
        onerr::Function = (err) -> rethrow(err),
        print_fun::Function = global_conf(:PRINT_FUN),
        mkdir::Bool = global_conf(:MK_DIR),
        verbose::Bool = global_conf(:VERBOSE),
        msg::AbstractString = "",
    )

    dir, args = _extract_dir(dirv, arg, args...)
    cfile = joinpath(dir, cfname(args...))
    try
        mkdir && mkpath(dirname(cfile))
        dat = isfile(cfile) ? deserialize(cfile)[DATA_KEY] : f()
        savecache && _scache(dat, cfile; verbose, onerr, print_fun)
        verbose && _io_print(print_fun, "CACHE LOADED", msg, dat, cfile)
        return dat
    catch err
        verbose && _io_error_print(print_fun, err, cfile)
        return onerr(err)
    end
end

lcache(f::Function, dirv::Vector{<:AbstractString}, arg, args...; kwargs...) =
    _lcache(f, true, dirv, arg, args...; kwargs...)

lcache(f::Function, arg, args...; kwargs...) =
    lcache(f::Function, [global_conf(:CACHE_DIR)], arg, args...; kwargs...)

lcache(dirv::Vector{<:AbstractString}, arg, args...; kwargs...) =
    _lcache(() -> nothing, false, dirv, arg, args...; kwargs...)

lcache(arg, args...; kwargs...) =
    lcache([global_conf(:CACHE_DIR)], arg, args...; kwargs...)

## ----------------------------------------------------------------------------
function delcache(dirv::Vector{<:AbstractString}, arg, args...;
        verbose::Bool = global_conf(:VERBOSE),
        print_fun::Function = global_conf(:PRINT_FUN)
    )
    dir, args = _extract_dir(dirv, arg, args...)
    path = isempty(args) ? dir : joinpath(dir, cfname(args...))
    if isfile(path)
        rm(path; force = true, recursive = true)
        verbose && print_fun(relpath(path), " deleted!!!")
        return path
    elseif isdir(path)
        tcaches = filter(is_cfname, readdir(path))
        for tc in tcaches
            tc = joinpath(path, tc)
            rm(tc, force = true)
            verbose && print_fun(relpath(tc), " deleted!!!")
        end
        return path
    end
end

delcache(args...; kwargs...) =
    delcache([global_conf(:CACHE_DIR)], args...; kwargs...)

## ----------------------------------------------------------------------------
function exist_cache(dirv::Vector{<:AbstractString}, arg, args...)
    dir, args = _extract_dir(dirv, arg, args...)
    cfile = joinpath(dir, cfname(args...))
    return isfile(cfile)
end
exist_cache(arg, args...) = exist_cache([global_conf(:CACHE_DIR)], arg, args...)

## ----------------------------------------------------------------------------
function backup_cachedir(;
        cache_dir::AbstractString = global_conf(:CACHE_DIR),
        backup_dir::AbstractString = string(cache_dir, "_backup")
    )
    tcaches = filter(is_cfname, readdir(cache_dir))
    !isdir(backup_dir) && mkpath(backup_dir)
    for file in tcaches
        src_file, dest_file = joinpath(cache_dir, file), joinpath(backup_dir, file)
        isfile(dest_file) && mtime(src_file) < mtime(dest_file) && continue
        cp(src_file, dest_file; force = true, follow_symlinks = true)
    end
    return backup_dir
end


text:
module test_test_warn_nowarn

using Test

@test_warn   "warning" println(stderr, :warning)

if VERSION >= v"1.8.0-DEV.363"
@test_warn   "warning" (@warn :warning)
end

@test_nowarn           (1 + 2)

end # module test_test_warn_nowarn


text:
using DataStructures

function readPlayer(f, N)
	p = Queue{Int}()
	for i=1:N
		k = parse(Int, readline(f))
		enqueue!(p, k)
	end
	p
end

function parsegame(N)
	f = open("22.txt")
	readline(f)
	p1 = readPlayer(f, N)
	readline(f)
	readline(f)
	p2 = readPlayer(f, N)
	(p1, p2)
end

# returns solution of part 1
function x()
	N = 25
	(p1, p2) = parsegame(N)

	while length(p1) > 0 && length(p2) > 0
		a = dequeue!(p1)
		b = dequeue!(p2)
		if a < b
			enqueue!(p2, b)
			enqueue!(p2, a)
		else
			enqueue!(p1, a)
			enqueue!(p1, b)
		end
	end
	x = 0
	for i=(2*N):-1:1
		x += dequeue!(p1) * i
	end
	x
end

function game(p1, p2)
	hist = Set{Tuple{Vector{Int}, Vector{Int}}}()
	while length(p1) > 0 && length(p2) > 0
		#println(p1)
		#println(p2)
		if (collect(p1), collect(p2)) in hist
			#println("repeat of ", p1, p2)
			return (true, p1, p2)
		end
		push!(hist, (collect(p1), collect(p2)))
		a = dequeue!(p1)
		b = dequeue!(p2)
		if a > length(p1) || b > length(p2)
			if a < b
				enqueue!(p2, b)
				enqueue!(p2, a)
			else
				enqueue!(p1, a)
				enqueue!(p1, b)
			end
		else
			q1 = copy(collect(p1))
			q2 = copy(collect(p2))
			r1 = Queue{Int}()
			r2 = Queue{Int}()
			for i=1:a
				enqueue!(r1, q1[i])
			end
			for i=1:b
				enqueue!(r2, q2[i])
			end
			p1won = game(r1, r2)[1]
			if p1won
				enqueue!(p1, a)
				enqueue!(p1, b)
			else
				enqueue!(p2, b)
				enqueue!(p2, a)
			end
		end
	end
	(length(p2) == 0, p1, p2)
end

# returns solution of part 2 (slow solution: 5.997827 seconds)
function y()
	N = 25
	(p1, p2) = parsegame(N)
	who, q1, q2 = game(p1, p2)
	x = 0
	if who
		for i=(2*N):-1:1
			x += dequeue!(q1) * i
		end
	else
		for i=(2*N):-1:1
			x += dequeue!(q2) * i
		end
	end
	x
end


text:
using BinDeps

@BinDeps.setup

deps = [sqlite3_lib = library_dependency("sqlite3_lib", aliases=["sqlite3","sqlite3-64","libsqlite3"])]

@BinDeps.if_install begin

provides(HomeBrew, {"sqlite" => sqlite3_lib})
provides(Yum, {"sqlite3-devel" => sqlite3_lib})
provides(AptGet, {"sqlite3" => sqlite3_lib})

provides(Binaries, {URI("http://drive.google.com/uc?export=download&id=0B6x_qYLl89S0TnAxS1Jtd3g3QkU") => sqlite3_lib}, os = :Windows)

@BinDeps.install

end


text:
doc"""
one layer of a multilayer perceptron (neural net)
ninᵢ := noutᵢ₋₁ + 1  (we add the bias term automatically, so there's one extra input)

This is a normalized layer, with inspiration from:
  `Batch Normalization: Accelerating Deep Network Training by Reducing Internal Covariate Shift`
    2015 Ioffe and Szegedy
  `http://arxiv.org/pdf/1502.03167.pdf`

forward value := a = f(Σ)
  Σ = yw + b
  y = xhat * β + α
  xhat = (x - μ) / σ

Note: that we solve for w, b, β, and α, and require intermediate sensitivities δΣ and δy
to compute all gradients.

note: w is a parameter for the case of tied weights (it can be a TransposeView!)
"""
type NormalizedLayer{A <: Mapping,
                     MATF <: AbstractMatrix{Float64},
                     GSTATE <: ParameterUpdaterState,
                     WGT <: Weight} <: NeuralNetLayer
  nin::Int
  nout::Int
  activation::A
  ploss::ParameterLoss
  p::Float64  # dropout retention probability

  # the state of the layer

  dwState::Matrix{GSTATE}
  dbState::Matrix{GSTATE}
  dβState::Matrix{GSTATE}
  dαState::Matrix{GSTATE}

  xvar::Vector{Variance{WGT}}  # nin x 1 -- online variances to calculate μ and σ for the normalization step
  w::MATF               # nout x nin -- weights connecting y --> Σ

  x::VecF               # nin x 1 -- input
  xhat::VecF            # nin x 1 -- xhat = standardize(x)
  α::VecF               # nin x 1 --
  y::VecF               # nin x 1 -- y = xhat .* β + α
  δy::VecF              # nin x 1 -- sensitivities for y:  δyᵢ := dC / dyᵢ   (calculated during backward pass)

  b::VecF               # nout x 1 -- bias terms
  Σ::VecF               # nout x 1 -- inner products
  a::VecF               # nout x 1 -- activation := f(Σ)
  δΣ::VecF               # nout x 1 -- sensitivities for Σ:  δΣᵢ := dC / dΣᵢ (calculated during backward pass)

  β::VecF               # nin x 1 -- updated by SGD step, used to compute y
  r::VecF               # nin x 1 -- vector of dropout retention... 0 if we drop this incoming weight, 1 if we keep it
  nextr::VecF           # nout x 1 -- retention of the nodes of this layer (as opposed to r
                        #             which applies to the incoming weights)
end

function NormalizedLayer(nin::Integer, nout::Integer, activation::Mapping,
                         updater::ParameterUpdater, p::Float64 = 1.0;
                         ploss::ParameterLoss = NoParameterLoss(),
                         weightInit::Function = _initialWeights,
                         wgt = ExponentialWeight(500))

  w = weightInit(nin, nout, activation)
  NormalizedLayer(nin, nout, activation, ploss, p,
                  ParameterUpdaterState(updater, nout, nin),
                  ParameterUpdaterState(updater, nout, 1),
                  ParameterUpdaterState(updater, nin, 1),
                  ParameterUpdaterState(updater, nin, 1),
                  [Variance(wgt) for i in 1:nin],
                  w,
                  [zeros(nin) for i in 1:5]...,
                  [zeros(nout) for i in 1:4]...,
                  ones(nin),
                  ones(nin),
                  ones(nout))
end

function Base.print{A,M,G}(io::IO, l::NormalizedLayer{A,M,G})
  print(io, "NormalizedLayer{$(l.nin)=>$(l.nout) $(l.activation) p=$(l.p) ")
  print(io, "‖δΣ‖₁=$(sumabs(l.δΣ)) ‖δy‖₁=$(sumabs(l.δy)) ")
  print(io, "$(M<:TransposeView ? "T" : "")}")
end
function Base.show(io::IO, l::NormalizedLayer)
  print(io, l)
end

function Base.display(::Base.REPL.REPLDisplay, ::MIME"text/plain", l::NormalizedLayer)
  println(l)
  @show l.r
  @show l.nextr
  @show l.x
  println(l.xvar)
  @show l.xhat
  @show l.y
  @show l.Σ
  @show l.a
  @show l.δy
  @show l.δΣ
  @show l.w
  @show l.b
  @show l.β
  @show l.α
  println()
end

# gemv! :: Σ += p * w * x  (note: 'T' would imply p * w' * x)
# note: we generate 2 different functions here, one for w::Matrix, and one for w::TransposeView
@generated function dosigmamult!(layer::NormalizedLayer, α::Float64)
  if layer.parameters[2] <: TransposeView
    return quote
      copy!(layer.Σ, layer.b)
      BLAS.gemv!('T', α, layer.w.mat, layer.y, 1.0, layer.Σ)
    end
  else
    return quote
      copy!(layer.Σ, layer.b)
      BLAS.gemv!('N', α, layer.w, layer.y, 1.0, layer.Σ)
    end
  end
end

# takes input vector, and computes Σⱼ = wⱼ'x + bⱼ  and  Oⱼ = A(Σⱼ)
function forward!(layer::NormalizedLayer, x::AVecF, istraining::Bool)

  # update r and get scale factor p
  if istraining
    # train... randomly drop out incoming nodes
    # note: I said incoming, since this layers weights are the weights connecting the previous layer to this one
    #       So on dropout, we are actually dropping out thr previous layer's nodes...
    for i in 1:layer.nin
      layer.r[i] = rand() <= layer.p ? 1.0 : 0.0
      fit!(layer.xvar[i], x[i])
    end
    p = 1.0
  else
    fill!(layer.r, 1.0)
    p = layer.p
  end

  # update x, xvar, xhat, y
  for i in 1:layer.nin
    layer.x[i] = x[i]

    # update xvar and standardize xi in one call
    layer.xhat[i] = standardize(layer.xvar[i], x[i])

    # we compute y = xhat * β + α, then multiply by r in case we're dropping out
    layer.y[i] = ((layer.xhat[i] * layer.β[i]) + layer.α[i]) * layer.r[i]
  end

  # update Σ
  dosigmamult!(layer, p)

  # forward!(layer.activation, layer.a, layer.Σ)     # activate
  value!(layer.a, layer.activation, layer.Σ)  # activate
end

# function updateSensitivities!(layer::NormalizedLayer, costmult::AVecF, multiplyDerivative::Bool)
#   updateδΣ!(layer, costmult, multiplyDerivative)
#   updateδy!(layer)
# end

function updateSensitivities!(layer::NormalizedLayer, mloss::ModelLoss, output::AVec, target::AVec)
    sensitivity!(layer.δΣ, layer.activation, mloss, layer.Σ, output, target)
    updateδy!(layer)
end

function updateSensitivities!(layer::NormalizedLayer, nextlayer::NormalizedLayer)
  updateδΣ!(layer, nextlayer)
  updateδy!(layer)
end

# # backward step for the final (output) layer
# # note: costmult is the amount to multiply against f'(Σ)... L2 case should be: (yhat-y)
# function updateδΣ!(layer::NormalizedLayer, costmult::AVecF, multiplyDerivative::Bool)
#   copy!(layer.δΣ, costmult)
#   if multiplyDerivative
#     for i in 1:layer.nout
#       layer.δΣ[i] *= backward(layer.activation, layer.Σ[i])
#     end
#   end
# end

# this is the backward step for a hidden layer
# notes: we are figuring out the effect of each node's activation value on the next sensitivities
# TODO: this is unstable... how can we stabilize correctly????
function updateδΣ!(layer::NormalizedLayer, nextlayer::NormalizedLayer)
  for i in 1:layer.nout
    σi = std(nextlayer.xvar[i]) |> if0then1
    βi = nextlayer.β[i]
    d = deriv(layer.activation, layer.Σ[i])
    layer.δΣ[i] = nextlayer.δy[i] * (βi / σi) * d
  end
end

function updateδy!(layer::NormalizedLayer)
  for i in 1:layer.nin
    δyi = 0.0
    for j in 1:layer.nout
      δyi += layer.δΣ[j] * layer.w[j,i]
    end
    layer.δy[i] = δyi
  end
end

doc"""
Update weights/bias one column at a time... skipping over the dropped out nodes.

gradients:
```
  bGrad = δΣ
  wGrad = δΣ * y
  αGrad = δy
  βGrad = δy * xhat
```
"""
function updateWeights!(layer::NormalizedLayer, updater::ParameterUpdater)

    # note: i refers to the output, j refers to the input
    # TODO: change gradient state to operate on one vector only... then have dbState, dwState, etc

    # update w and b
    # note: make sure we only update when the node is retained (from dropout)
    for i in 1:layer.nout
        if layer.nextr[i] > 0.0
            # layer.b[i] += Δij(updater, layer.dbState, layer.δΣ[i], 0.0, i, 1)
            layer.b[i] += param_change!(layer.dbState[i], updater, layer.ploss, layer.δΣ[i], 0.0)

            for j in 1:layer.nin
                if layer.r[j] > 0.0
                    # layer.w[i,j] += Δij(updater, layer.dwState, layer.δΣ[i] * layer.y[j], layer.w[i,j], i, j)
                    layer.w[i,j] += param_change!(layer.dwState[i,j], updater, layer.ploss, layer.δΣ[i] * layer.y[j], layer.w[i,j])
                end
            end
        end
    end

    # update β and α... only when the input is retained (from dropout)
    for i in 1:layer.nin
        if layer.r[i] > 0.0
            # layer.α[i] += Δij(updater, layer.dαState, layer.δy[i], 0.0, i, 1)
            # layer.β[i] += Δij(updater, layer.dβState, layer.δy[i] * layer.xhat[i], 0.0, i, 1)
            layer.α[i] += param_change!(layer.dαState[i], updater, layer.ploss, layer.δy[i], 0.0)
            layer.β[i] += param_change!(layer.dβState[i], updater, layer.ploss, layer.δy[i] * layer.xhat[i], 0.0)
        end
    end

end


text:
# automatically generated by the FlatBuffers compiler, do not modify

NamespaceA.NamespaceB.eval(quote

@enum EnumInNestedNS::Int8 begin
    EnumInNestedNSA = 0
    EnumInNestedNSB = 1
    EnumInNestedNSC = 2
end

end)

max_stars_repo_name: rjkat/flatbuffers-julia

max_issues_repo_name: rjkat/flatbuffers-julia

max_forks_repo_name: rjkat/flatbuffers-julia


text:
function TimeSpanMap(arg0::Object)
    return TimeSpanMap((Object,), arg0)
end

function add_valid_after(obj::TimeSpanMap, arg0::Object, arg1::AbsoluteDate)
    return jcall(obj, "addValidAfter", void, (Object, AbsoluteDate), arg0, arg1)
end

function add_valid_before(obj::TimeSpanMap, arg0::Object, arg1::AbsoluteDate)
    return jcall(obj, "addValidBefore", void, (Object, AbsoluteDate), arg0, arg1)
end

function extract_range(obj::TimeSpanMap, arg0::AbsoluteDate, arg1::AbsoluteDate)
    return jcall(obj, "extractRange", TimeSpanMap, (AbsoluteDate, AbsoluteDate), arg0, arg1)
end

function get(obj::TimeSpanMap, arg0::AbsoluteDate)
    return jcall(obj, "get", Object, (AbsoluteDate,), arg0)
end

function get_span(obj::TimeSpanMap, arg0::AbsoluteDate)
    return jcall(obj, "getSpan", TimeSpanMap_Span, (AbsoluteDate,), arg0)
end

function get_transitions(obj::TimeSpanMap)
    return jcall(obj, "getTransitions", NavigableSet, ())
end


text:
##
import LinearAlgebra as linalg
import Plots as plt

include("KalmanFilter.jl")
import .KalmanFilter as kf

##
function get_data()
    # Observations
    x = [4000.0 4260.0 4550.0 4860.0 5110.0] # position [m]
    v = [ 280.0  282.0  285.0  286.0  290.0] # velocity [m/s]
    # x = [4000.0 4260.0 4550.0 4860.0 5110.0 5220.0 5115.0 4995.0 4871.0] # position [m]
    # v = [ 280.0  282.0  285.0  286.0  290.0  295.0  292.0  287.0  284.0] # velocity [m/s]
    X = [x; v]

    # Initial values
    a = 2.0 # acceleration
    Δt = 1.0 # Time step
    ΔP = [20.0; 5.0] # Process errors in process covariance matrix
    Δ = [25.0; 6.0] # Observation errors
    return X, a, Δ, ΔP, Δt
end

## Define the state space model

function define_model(X, a, Δ, ΔP, Δt)
    m, n = size(X)
    Q = zeros(m, m) # Process noise covariance matrix
    F = [1.0 Δt; 0.0 1.0] # Transition matrix
    B = [Δt^2/2.0 Δt]' # Control-input matrix
    u = 2.0 # # Control vector, acceleration
    R = linalg.Diagonal([Δ[1]^2 0.0; 0.0 Δ[2]^2]) # Sensor noise covariance matrix
    H = kf.eye(m) # Measurement matrix
    w = zeros(m)
    v = zeros(m)
    # Initial predicted process covariance matrix
    P = linalg.Diagonal([ΔP[1]^2 ΔP[1]*ΔP[2]; ΔP[2]*ΔP[1] ΔP[2]^2])
    return P, Q, R, F, H, B, u, m, n
end

## Data input
function data_input(X, n)
    X = [ [X[1,i], X[2,i]] for i in 1:n ]
    Xp = [ X[1] ] # Initialize predicted state
    Xk = copy(Xp) # Initialized new state matrix
    return Xp, Xk, X
end

## Run Kalman Filter
function run_kf(X, P, F, H, Q, R, B, u; filter=kf.CholeskySqrt())
    n = size(X, 2)
    _, _, Xx = data_input(X, n)
    k = 1
    S = kf._innovation_covariance(P, H, R)
    K = kf._kalman_gain(P, H, S)
    skf = [ kf.LinearKF(Xx[k], Xx[k], P, K, S, Xx[1], filter) ]
    while k < n
        s = kf.linear_kalman(
            Xx[k+1], skf[k].xk, skf[k].P, F, H, Q, R;
            B=B, u=u, sqrt_filter=filter,
        )
        push!(skf, s)
        k += 1
    end
    return skf
end

## Smoothers
function run_smoothers_bts(Xk, F, Σ)
    k = length(Xk)
    sm = Vector{kf.RTSSmoother}(undef, k)
    sm[k] = kf.RTSSmoother(Xk[k], Σ[k])
    ksmoother = kf.RauchTungStriebel()
    while k > 1
        s = kf.kalman_smooth(ksmoother, sm[k].xk, Xk[k-1], Xk[k], Σ[k-1], Σ[k], F)
        k -= 1
        sm[k] = s
    end
    return sm
end

function run_smoothers_mbf(skf, F, H)
    m = size(skf[1].xk, 1)
    k = length(skf)
    sm = Vector{kf.MBFSmoother}(undef, k)
    ksmoother = kf.ModifiedBrysonFrazier()
    # Λ, λ = zeros(size(lsm.F)), zeros(size(lsm.F,1))
    SH = skf[k].S \ H'
    Λ, λ = SH*H, -SH*skf[k].y
    sm[k] = kf.MBFSmoother(skf[k].xk, skf[k].P, Λ, λ)
    while k > 1
        s = kf.kalman_smooth(
            ksmoother,
            skf[k-1].xk,
            skf[k-1].P,
            skf[k-1].K,
            F,
            H,
            sm[k].L,
            sm[k].l,
            skf[k-1].S,
            skf[k-1].y,
        )
        k -= 1
        sm[k] = s
    end
    return sm
end

## Plots
function plot_kf(Xp, Xk, Xx, Δ)
    # Xp, Xk, Xx = reduce(hcat, Xp), reduce(hcat, Xk), reduce(hcat, Xx)
    p1 = plt.plot(Xp[1,:], m=:o, label="Predicted", legend=:topleft, ylabel="Position [m]", xlabel="Time [s]")
    plt.plot!(p1, Xx[1,:], m=:s, label="Measured", ribbon=Δ[1])
    plt.plot!(p1, Xk[1,:], m=:p, label="KF")
    p2 = plt.plot(Xp[2,:], m=:o, label="Predicted", legend=:topleft, ylabel="Velocity [m/s]", xlabel="Time [s]")
    plt.plot!(p2, Xx[2,:], m=:s, label="Measured", ribbon=Δ[2])
    plt.plot!(p2, Xk[2,:], m=:p, label="KF")
    p = plt.plot(p1, p2, size=(700,350))
    return p
end

function plot_kf_smooth(Xp, Xk, Xx, Xks, Δ)
    # Xp, Xk = reduce(hcat, Xp), reduce(hcat, Xk)
    # Xx, Xks = reduce(hcat, Xx), reduce(hcat, Xks)
    p1 = plt.plot(Xp[1,:], m=:o, label="Predicted", legend=:topleft, ylabel="Position [m]", xlabel="Time [s]")
    plt.plot!(p1, Xx[1,:], m=:s, label="Measured", ribbon=Δ[1])
    plt.plot!(p1, Xk[1,:], m=:p, label="KF")
    plt.plot!(p1, Xks[1,:], m=:h, label="Smooth")
    p2 = plt.plot(Xp[2,:], m=:o, label="Predicted", legend=:topleft, ylabel="Velocity [m/s]", xlabel="Time [s]")
    plt.plot!(p2, Xx[2,:], m=:s, label="Measured", ribbon=Δ[2])
    plt.plot!(p2, Xk[2,:], m=:p, label="KF")
    plt.plot!(p2, Xks[2,:], m=:h, label="Smooth")
    p = plt.plot(p1, p2, size=(700,350))
    return p
end

## Run KFs
X, a, Δ, ΔP, Δt = get_data()
P, Q, R, F, H, B, u, m, n = define_model(X, a, Δ, ΔP, Δt)

## Simple Kalman filter
skf = run_kf(X, P, F, H, Q, R, B, u; filter=kf.CholeskyModSqrt())
Xp, Xk = kf.unpack(skf, [:xp, :xk])
p = plot_kf(Xp, Xk, X, Δ)

## Rauch Tung Striebel Smoother
Σ, Xk2m = kf.unpack_matrix(skf, [:P, :xk])
skfm = run_smoothers_bts(Xk2m, F, Σ)
Xk2m = kf.unpack(skfm, :xk)
p2m = plot_kf_smooth(Xp, Xk, X, Xk2m, Δ)

## Modified Bryson–Frazier Smoother
skfmb = run_smoothers_mbf(skf, F, H)
Xk2mb = kf.unpack(skfmb, :xk)
p2m = plot_kf_smooth(Xp, Xk, X, Xk2mb, Δ)