text:
# *****************************************************************************
# Written by Ritchie Lee, ritchie.lee@sv.cmu.edu
# *****************************************************************************
# Copyright ã 2015, United States Government, as represented by the
# Administrator of the National Aeronautics and Space Administration. All
# rights reserved.  The Reinforcement Learning Encounter Simulator (RLES)
# platform is licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License. You
# may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0. Unless required by applicable
# law or agreed to in writing, software distributed under the License is
# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the License for the specific language
# governing permissions and limitations under the License.
# _____________________________________________________________________________
# Reinforcement Learning Encounter Simulator (RLES) includes the following
# third party software. The SISLES.jl package is licensed under the MIT Expat
# License: Copyright (c) 2014: Youngjun Kim.
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED
# "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
# NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
# PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
# *****************************************************************************

#TODO: Convert to new Problem format and move to /modules/AntProblem

module GrammarDef

export create_grammar, to_function
export World, Ant, copy

using GrammaticalEvolution
import GrammaticalEvolution.isless
import Base: +, -, copy

function create_grammar()
  @grammar grammar begin
    start = block
    command = turnleft | turnright | forward | if_statement | loop
    turnleft = Expr(:call, :turn_left, :ant)
    turnright = Expr(:call, :turn_right, :ant)
    forward = Expr(:call, :forward, :world, :ant)
    if_statement = Expr(:if, condition, block, block)
    condition = food_ahead
    food_ahead = Expr(:call, :food_ahead, :world, :ant)
    loop = Expr(:for, Expr(:(=), :i, Expr(:(:), 1, digit)), block)
    block[make_block] = (command)^(1:5)
    digit = 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9
  end
  return grammar
end

function to_function(code)
  @eval f(world::World, ant::Ant) = $code
  return f
end

const FOOD = convert(Int64, '1')
const GAP = convert(Int64, '0')
const EMPTY = 46 #convert(Int64, '.')

type World
  width::Int64
  height::Int64
  value::Array{Int64, 2}
  start::Array{Int64, 1}

  function World(width::Int64, height::Int64)
    value = zeros(width, height)
    return new(width, height, value)
  end

  function World(map::Array{Int64, 2}, start::Array)
    return new(size(map, 1), size(map, 2), map, start)
  end
end

copy(world::World) = World(copy(world.value), world.start)

immutable Direction
  value::Int64
end

+{T <: Number}(d::Direction, v::T) = Direction(d.value + v)
-{T <: Number}(d::Direction, v::T) = Direction(d.value - v)
isless(d1::Direction, d2::Direction) = isless(d1.value, d2.value)
isless(d1::Direction, v::Int64) = isless(d1.value, value)

const north = Direction(1)
const east = Direction(2)
const south = Direction(3)
const west = Direction(4)

type Ant
  x::Int64
  y::Int64
  direction::Direction
  eaten::Int64

  function Ant(x::Int64, y::Int64)
    return new(x, y, north, 0)
  end
end

function calc_forward(world::World, ant::Ant)
  nx = ant.x
  ny = ant.y

  if ant.direction === north
    ny -= 1
    if ny < 1 ny = world.height end
  elseif ant.direction === east
    nx += 1
    if nx > world.width nx = 1 end
  elseif ant.direction === south
    ny += 1
    if ny > world.height ny = 1 end
  else
    nx -= 1
    if nx < 1 nx = world.width end
  end

  return (nx, ny)
end

function food_ahead(world::World, ant::Ant)
  (x, y) = calc_forward(world, ant)
  return world.value[x, y] == FOOD
end

function forward(world::World, ant::Ant)
  (ant.x, ant.y) = calc_forward(world, ant)
  if world.value[ant.x, ant.y] == FOOD
    world.value[ant.x, ant.y] = EMPTY
    ant.eaten += 1
  elseif world.value[ant.x, ant.y] == EMPTY
    ant.eaten -= 1
  end
end

function turn_right(ant::Ant)
  #ant.direction = ant.direction + 1 > south ? north : ant.direction+1
  ndir = ant.direction + 1
  ant.direction = ndir.value > south.value ? north : ndir;
end

function turn_left(ant::Ant)
  #ant.direction = ant.direction - 1 < 1 ? south : ant.direction-1
  ndir = ant.direction - 1
  ant.direction = ndir.value < 1 ? south : ndir
end

make_block(lst::Array) = Expr(:block, lst...)
make_loop(values::Array) = Expr(:for, :(=), :i, Expr(:(:), 1, values[1]), values[2])
make_if(values::Array) = Expr(:if, values[1], values[2], values[3])

macro make_if(condition, true_block, false_block)
  Expr(:if, condition, true_block, false_block)
end

macro make_call(fn, args...)
  Expr(:call, fn, args...)
end

macro make_for(start, stop, block)
  Expr(:for, Expr(:(=), :i, Expr(:(:), start, stop)), block)
end

end #module


text:
module TodoValidator

using App, SearchLight, Validation

function not_empty{T<:AbstractModel}(field::Symbol, m::T, args::Any...) :: ValidationStatus
  if isempty(getfield(m, field))
    false, :not_empty, "should be not empty"
  else
    true, :not_empty, ""
  end
end

end

ext: jl

max_stars_repo_name: essenciary/genie-todo-mvc

max_issues_repo_name: essenciary/genie-todo-mvc

max_forks_repo_name: essenciary/genie-todo-mvc


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

using Markdown
using InteractiveUtils

# ╔═╡ 3cd880f0-6856-11eb-187f-e7fe0238d041
include("config.jl")

# ╔═╡ 78625b00-6856-11eb-0e08-bf8915eb2f5e
md"# A most **basic** test and headline"

# ╔═╡ 9f489670-6857-11eb-16f7-53ecb28dc3a5
md"A paragraph containing a link to [Julia](http://www.julialang.org)."

# ╔═╡ 71383360-6890-11eb-14eb-e39ac006624e
stri = "Da jeg var på vei til kiirken ∈ dag morges så kom jeg forbi en liten sjømann. En frisk og hyggelig liten sjømann som hilste meg."

# ╔═╡ fcdff590-686a-11eb-10ad-0f67e129deda
l1 = text(stri, -0.45 * WI, 0)

# ╔═╡ 240c9390-6892-11eb-3ba7-e95f27ba7e5a
l2 = sethue("yellow")

# ╔═╡ 32b3cc10-6892-11eb-38d2-99a597b4b53c
   l3 = text("1234567890123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890",
        -WI / 2, EM)

# ╔═╡ e3ca9330-6892-11eb-2b65-2f98c90f0c9f

# ╔═╡ 5b66cdae-6892-11eb-0a16-91fdc7e0c19e
l4 = circle(O, 20, :stroke)

# ╔═╡ 9796abbe-6892-11eb-2faa-a326479b3b7c
t = Turtle()

# ╔═╡ 9c961ed2-6892-11eb-1fdc-9d379272db73
l5 = Pencolor(t, "cyan")

# ╔═╡ a1d7f210-6892-11eb-1b33-5958ade28f77
l6 = Penwidth(t, 1.5)

# ╔═╡ b5048ba0-6892-11eb-1225-f32767f28867
n = 4

# ╔═╡ bbe038c0-6892-11eb-1333-a70977b4d71c
l7 = for i in 1:200
        Forward(t, n)
        Turn(t, 89.5)
        HueShift(t)
        n += 0.75
    end

# ╔═╡ 53105580-6894-11eb-3af3-3feadd3c43bf
begin
	l1, l2, l3, l4, l5, l6, l7
	this_fig = empty_figure(filename = joinpath(@__DIR__, "pluto_1.png"))
	background("green")
end

# ╔═╡ cad199a0-6892-11eb-2e31-316dc6c8e541
begin
	l1, l2, l3, l4, l5, l6, l7
	finish()
end

# ╔═╡ c4ab80e0-6892-11eb-19a5-cf50d9dc0bdb

# ╔═╡ Cell order:
# ╠═3cd880f0-6856-11eb-187f-e7fe0238d041
# ╟─78625b00-6856-11eb-0e08-bf8915eb2f5e
# ╠═9f489670-6857-11eb-16f7-53ecb28dc3a5
# ╠═53105580-6894-11eb-3af3-3feadd3c43bf
# ╠═71383360-6890-11eb-14eb-e39ac006624e
# ╠═fcdff590-686a-11eb-10ad-0f67e129deda
# ╠═240c9390-6892-11eb-3ba7-e95f27ba7e5a
# ╠═32b3cc10-6892-11eb-38d2-99a597b4b53c
# ╠═e3ca9330-6892-11eb-2b65-2f98c90f0c9f
# ╠═5b66cdae-6892-11eb-0a16-91fdc7e0c19e
# ╠═9796abbe-6892-11eb-2faa-a326479b3b7c
# ╠═9c961ed2-6892-11eb-1fdc-9d379272db73
# ╠═a1d7f210-6892-11eb-1b33-5958ade28f77
# ╠═b5048ba0-6892-11eb-1225-f32767f28867
# ╠═bbe038c0-6892-11eb-1333-a70977b4d71c
# ╠═cad199a0-6892-11eb-2e31-316dc6c8e541
# ╠═c4ab80e0-6892-11eb-19a5-cf50d9dc0bdb


text:
# This file was generated, do not modify it. # hide
res = y_hat .- y_wind[test];


text:
"""
Training and Validation subset selection through either:
- Fast and Unique Representative Subset selection
- Random Representative Subset selection
"""
function selectSubset(method::String, fileIndex::Int, trainSize::Int, validSize::Int, degreeThreshold::Float64)
  if method == "furs"
    sortedNodeIndices = sortIndicesDegree(fileData[fileIndex][2], degreeThreshold)
    trainSubset, sortedNodeIndices, trainComplete = furs(fileIndex, trainSize, sortedNodeIndices)
    # TODO: remove training subset from neighbours/degree of remaining nodes
    validSubset, sortedNodeIndices, validComplete = furs(fileIndex, validSize, sortedNodeIndices)
    return (trainSubset, trainComplete, validSubset, validComplete)
  elseif method == "randrs"
    # get all node indices above the degree threshold
    nodeIndices = find(x -> x > degreeThreshold, fileData[fileIndex][2])
    trainSubset, trainComplete = randrs(trainSize, nodeIndices)
    validSubset, validComplete = randrs(validSize, nodeIndices)
    return (trainSubset, trainComplete, validSubset, validComplete)
  else
    error("selectSubset: wrong method argument.")
  end
end

"""
Fast and Unique Representative Subset selection (FURS)
Return array of selected node indices and boolean value (true if equal to requested sampleSize)
"""
function furs(fileIndex::Int, sampleSize::Int, sortedNodeIndices::Vector{Int})
  selectedSize = 0
  selectedSubset = Array(Int, 0)
  tempIndices = Array(Int, 0)
  while selectedSize < sampleSize
    # check if still viable nodes left
    if isempty(sortedNodeIndices)
      if isempty(tempIndices)
        # FURS unable to finnish
        return (selectedSubset, sortedNodeIndices, false)
      else
        # add and sort the removed indices back into sortedNodeIndices
        sortedNodeIndices = tempIndices[sortIndicesDegree(fileData[fileIndex][2][tempIndices])]
        empty!(tempIndices)
      end
    else
      # get index of node with highest degree
      nodeIndex = splice!(sortedNodeIndices, 1)
      neighbourArray = fileData[fileIndex][3][:, nodeIndex]
      push!(selectedSubset, nodeIndex)
      selectedSize += 1
      # remove selected node's neighbours from sortedNodeIndices and store them in tempIndices
      for neighbour in neighbourArray
        # stop loop if no more candidate nodes
        if isempty(sortedNodeIndices)
          break
        end
        # stop loop when -1 padding is reached
        if neighbour == -1
          break
        end
        index = findfirst(fileData[fileIndex][1][sortedNodeIndices], neighbour)
        if index != 0
          push!(tempIndices, splice!(sortedNodeIndices, index))
        end
      end
    end
  end
  # add tempIndices back into sortedNodeIndices
  if !isempty(tempIndices)
    push!(sortedNodeIndices, tempIndices...)
    sortedNodeIndices = sortedNodeIndices[sortIndicesDegree(fileData[fileIndex][2][sortedNodeIndices])]
  end
  return (selectedSubset, sortedNodeIndices, true)
end

"""
Random representative subset selection
Return array of randomly selected node indices and boolean value (true if equal to requested sampleSize)
"""
function randrs(sampleSize::Int, nodeIndices::Vector{Int})
  selectedSize = 0
  selectedSubset = Array(Int, 0)
  while selectedSize < sampleSize
    if isempty(nodeIndices)
      return (selectedSubset, false)
    end
    # randomly sample nodeIndices
    sample = rand(1:length(nodeIndices))
    push!(selectedSubset, splice!(nodeIndices, sample))
    selectedSize += 1
  end
  return (selectedSubset, true)
end

"""
Node index sorting functions
"""
function sortIndicesDegree(degreeArray::Array{Float64, 2}, threshold::Float64)
  # remove all node indices with degree lower than threshold
  nodeIndices = find(x -> x >= threshold, degreeArray)
  # get sorted node indices based on descending degree
  return nodeIndices[sortperm(degreeArray[nodeIndices], rev=true)]
end

function sortIndicesDegree(degreeArray::Array{Float64, 1})
  # get sorted node indices based on descending degree
  return sortperm(degreeArray, rev=true)
end


text:
"""
    fit(events::Array{<:Number, 1}, maxT::Number, its::Int64)

# Arguments

* `events` The array of event times to fit the Hawkes process too.
* `maxT` The boundary time over which the events were observed.
* `its` Number of iterations to sample for.

# Notes

* All `events` must be unique

"""
function fit(events::Array{<:Number, 1}, maxT::Number, its::Int64)

    nEvents = length(events)
    eventDifferences = event_difference_list(events)
    maxTDifferences = maxT .- events

    bgSamples = zeros(Float64, its)
    kappaSamples = zeros(Float64, its)
    kernSamples = zeros(Float64, its)

    bgSamples[1] = rand()
    kappaSamples[1] = rand()
    kernSamples[1] = rand()

    kernSample = kernSamples[1]
    kernDistribution::Exponential{Float64} = Distributions.Exponential(1/kernSample)
    kernFunc(x::Number) = Distributions.pdf(kernDistribution, x)

    parentVectorSample, numBG, chEvents, shiftTimes = sample_parents(events, bgSamples[1], kappaSamples[1], kernFunc, eventDifferences)

    for i = 2:its

        bgSamples[i] = Distributions.rand(Distributions.Gamma(0.01 + numBG, 1/(0.01+1))) / maxT

        H_tilde::Float64 = sum(cdf.(kernDistribution, maxTDifferences))
        kappaSamples[i] = Distributions.rand(Distributions.Gamma(0.01 + sum(chEvents), 1/(0.01 + H_tilde)))

        kernSample = Distributions.rand(Distributions.Gamma(0.01 + length(shiftTimes), 1/(0.01 + sum(shiftTimes))))
        kernSamples[i] = kernSample
        kernDistribution = Distributions.Exponential(1/kernSamples[i])

        parentVectorSample, numBG, chEvents, shiftTimes = sample_parents(events, bgSamples[i], kappaSamples[i], kernFunc, eventDifferences)
    end
    bgSamples, kappaSamples, kernSamples
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "LAMMPS"
version = v"1.0.1" # Equivalent to 2020-10-29

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/lammps/lammps.git", "88fd96ec52f86dba4b222623f3a06632a32e42f1")
]

# Bash recipe for building across all platforms
# LAMMPS DPD packages do not work on all platforms
script = raw"""
cd $WORKSPACE/srcdir/lammps/
mkdir build && cd build/
cmake -C ../cmake/presets/most.cmake -C ../cmake/presets/nolib.cmake ../cmake -DCMAKE_INSTALL_PREFIX=${prefix} \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
    -DCMAKE_BUILD_TYPE=Release \
    -DBUILD_SHARED_LIBS=ON \
    -DLAMMPS_EXCEPTIONS=ON \
    -DPKG_MPI=ON \
	-DPKG_SNAP=ON \
    -DPKG_DPD-BASIC=OFF \
	-DPKG_DPD-MESO=OFF \
	-DPKG_DPD-REACT=OFF \
        -DPKG_USER-MESODPD=OFF \
        -DPKG_USER-DPD=OFF \
        -DPKG_USER-SDPD=OFF \
	-DPKG_DPD-SMOOTH=OFF

make -j${nproc}
make install

if [[ "${target}" == *mingw* ]]; then
    cp *.dll ${prefix}/bin/
fi
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
# platforms = supported_platforms(; experimental=true)
platforms = supported_platforms()
platforms = filter(p -> !(Sys.isfreebsd(p) || libc(p) == "musl"), platforms)

# We need this since currently MPItrampoline_jll has a dependency on gfortran
platforms = expand_gfortran_versions(platforms)
# libgfortran3 does not support `!GCC$ ATTRIBUTES NO_ARG_CHECK`. (We
# could in principle build without Fortran support there.)
platforms = filter(p -> libgfortran_version(p) ≠ v"3", platforms)

platforms = expand_cxxstring_abis(platforms)

# The products that we will ensure are always built
products = [
    LibraryProduct("liblammps", :liblammps),
    ExecutableProduct("lmp", :lmp),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="CompilerSupportLibraries_jll")),
    Dependency(PackageSpec(name="MPItrampoline_jll"), compat="2"),
    Dependency(PackageSpec(name="MicrosoftMPI_jll"))
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6", preferred_gcc_version=v"8")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
module TestProject

using ComoniconTypes

const CASTED_COMMANDS = Dict(
    "main" => Entry(;
        version=v"1.2.0",
        root=LeafCommand(;
            fn=sin,
            name="sin",
            args=Argument[
                Argument(;name="x", type=Float64),
            ]
        )
    )
)

end


text:
######################################################################
# trajectory_generation.jl
######################################################################

# helper files
include("gradients.jl")
include("scoring.jl")
include("projection.jl")
include("printing.jl")

# actual methods
include("pto.jl")
include("rrt.jl")
include("pto_multi.jl")
include("pto_linear.jl")
include("smc.jl")
include("smc_multi.jl")
include("max.jl")
include("kmeans.jl")
#include("rig.jl")

# some additional helpers
include("info.jl")

function check_convergence(es::Float64, es_crit::Float64, i::Int, max_iters::Int, dd::Float64, dd_crit::Float64, verbose::Bool, es_count::Int)
    not_finished = true
    if es < es_crit
        not_finished = false
        if verbose
            println("reached ergodic criterion...")
        end
    end
    if i > max_iters
        not_finished = false
        if verbose
            println("max iterations reached...")
        end
    end
    if abs(dd) < dd_crit
        not_finished = false
        if verbose
            println("reached directional derivative criterion...")
        end
    end
    if es_count > 50
        not_finished = false
        if verbose
            println("We've been stuck for 50 iterations...")
        end
    end
    return not_finished
end

# called if logging, not meant for general use
function save(outfile::IOStream, xd::VVF)
    n = length(xd[1])
    for xi in xd
        for i = 1:(n-1)
            wi = xi[i]
            write(outfile,"$(xi[i]),")
        end
        write(outfile,"$(xi[n])\n")
    end
end


text:
export CompositionSystem

"""
    CompositionSystem(composition::Composition, systems_constructor) <: AbstractSystem

A system representing the composition of polynomial maps.
"""
struct CompositionSystem{S1<:AbstractSystem,S2<:AbstractSystem} <: AbstractSystem
    # The system is g ∘ f
    g::S2 # Never a composition system
    f::S1 # Can be a CompositionSystem again

    g_has_parameters::Bool
    f_has_parameters::Bool
end

function CompositionSystem(
    C::Composition,
    system_constructor;
    variables = nothing,
    parameters = nothing,
    homvars = nothing,
)
    n = length(C.polys)
    f = CompositionSystem(
        C.polys[n-1],
        C.polys[n],
        system_constructor;
        input_variables = variables,
        parameters = parameters,
        homvars = homvars,
    )
    for k = (n - 2):-1:1
        f = CompositionSystem(
            C.polys[k],
            f,
            system_constructor;
            parameters = parameters,
            homvars = homvars,
        )
    end
    f
end

function CompositionSystem(
    g::Vector{<:MP.AbstractPolynomialLike},
    f::Union{CompositionSystem,Vector{<:MP.AbstractPolynomialLike}},
    system_constructor;
    input_variables = nothing,
    parameters = nothing,
    homvars = nothing,
)

    vars = variables(g; parameters = parameters)
    homvars_to_end!(vars, homvars)
    G = system_constructor(g, variables = vars, parameters = parameters)

    if isa(f, CompositionSystem)
        CompositionSystem(G, f, hasparameters(g, parameters), hasparameters(f))
    else
        F = system_constructor(f, variables = input_variables, parameters = parameters)
        CompositionSystem(G, F, hasparameters(g, parameters), hasparameters(f, parameters))
    end
end

hasparameters(C::CompositionSystem) = C.g_has_parameters || C.f_has_parameters

homvars_to_end!(variables, ::Nothing) = variables
homvars_to_end!(vars, homvars::MP.AbstractVariable) = homvars_to_end!(vars, (homvars,))
function homvars_to_end!(vars, homvars)
    deleteat!(vars, findall(v -> v ∈ homvars, vars))
    append!(vars, homvars)
    vars
end

struct CompositionSystemCache{
    C1<:AbstractSystemCache,
    C2<:AbstractSystemCache,
    T,
} <: AbstractSystemCache
    cache_f::C1
    cache_g::C2

    eval_f::Vector{T}
    J_f::Matrix{T}
    J_g::Matrix{T}
    Jp_f::Union{Nothing,Matrix{T}}
    Jp_g::Union{Nothing,Matrix{T}}
end

function cache(C::CompositionSystem, x, p = nothing)
    f, g = C.f, C.g
    c_f = p === nothing ? cache(f, x) : cache(f, x, p)
    eval_f = p === nothing ? evaluate(f, x, c_f) : evaluate(f, x, p, c_f)
    c_g = p === nothing ? cache(g, eval_f) : cache(g, eval_f, p)
    J_f = p === nothing ? jacobian(f, x, c_f) : jacobian(f, x, p, c_f)
    J_g = p === nothing ? jacobian(g, eval_f, c_g) : jacobian(g, eval_f, p, c_g)
    # need to bring eval_f, J_f, J_g to the same element type
    T = promote_type(eltype(eval_f), eltype(J_f), eltype(J_g))
    eval_f_T = convert(Vector{T}, eval_f)
    J_f_T = convert(Matrix{T}, J_f)
    J_g_T = convert(Matrix{T}, J_g)
    if p === nothing
        Jp_f = Jp_g = nothing
    else
        Jp_f = similar(J_f_T, size(J_f_T, 1), length(p))
        Jp_g = similar(J_g_T, size(J_g_T, 1), length(p))
    end

    CompositionSystemCache(c_f, c_g, eval_f_T, J_f_T, J_g_T, Jp_f, Jp_g)
end

Base.size(C::CompositionSystem) = (length(C.g), size(C.f)[2])

function evaluate!(u, C::CompositionSystem, x, c::CompositionSystemCache)
    evaluate!(u, C.g, evaluate!(c.eval_f, C.f, x, c.cache_f), c.cache_g)
end
function evaluate!(u, C::CompositionSystem, x, p, c::CompositionSystemCache)
    evaluate!(u, C.g, evaluate!(c.eval_f, C.f, x, p, c.cache_f), p, c.cache_g)
end
function evaluate(C::CompositionSystem, x, c::CompositionSystemCache)
    evaluate(C.g, evaluate!(c.eval_f, C.f, x, c.cache_f), c.cache_g)
end
function evaluate(C::CompositionSystem, x, p, c::CompositionSystemCache)
    evaluate(C.g, evaluate!(c.eval_f, C.f, x, p, c.cache_f), p, c.cache_g)
end
function jacobian!(U, C::CompositionSystem, x, c::CompositionSystemCache)
    # chain rule
    evaluate_and_jacobian!(c.eval_f, c.J_f, C.f, x, c.cache_f)
    jacobian!(c.J_g, C.g, c.eval_f, c.cache_g)
    LinearAlgebra.mul!(U, c.J_g, c.J_f)
end
function jacobian!(U, C::CompositionSystem, x, p, c::CompositionSystemCache)
    # chain rule
    evaluate_and_jacobian!(c.eval_f, c.J_f, C.f, x, p, c.cache_f)
    jacobian!(c.J_g, C.g, c.eval_f, p, c.cache_g)
    LinearAlgebra.mul!(U, c.J_g, c.J_f)
end
function jacobian(C::CompositionSystem, x, c::CompositionSystemCache)
    jacobian!(similar(c.J_g, size(c.J_g, 1), size(c.J_f, 2)), C, x, c)
end
function jacobian(C::CompositionSystem, x, p, c::CompositionSystemCache)
    jacobian!(similar(c.J_g, size(c.J_g, 1), size(c.J_f, 2)), C, x, p, c)
end
function evaluate_and_jacobian!(u, U, C::CompositionSystem, x, c::CompositionSystemCache)
    evaluate_and_jacobian!(c.eval_f, c.J_f, C.f, x, c.cache_f)
    evaluate_and_jacobian!(u, c.J_g, C.g, c.eval_f, c.cache_g)
    LinearAlgebra.mul!(U, c.J_g, c.J_f)
    nothing
end
function evaluate_and_jacobian!(u, U, C::CompositionSystem, x, p, c::CompositionSystemCache)
    evaluate_and_jacobian!(c.eval_f, c.J_f, C.f, x, p, c.cache_f)
    evaluate_and_jacobian!(u, c.J_g, C.g, c.eval_f, p, c.cache_g)
    LinearAlgebra.mul!(U, c.J_g, c.J_f)
    nothing
end

function differentiate_parameters!(U, C::CompositionSystem, x, p, c::CompositionSystemCache)
    @assert c.Jp_f !== nothing && c.Jp_g !== nothing

    # d/(dp)(g(f(x, p), p)) = g_y(f(x, p), p)f_p(x, p) + g_p(f(x, p), p)
    evaluate!(c.eval_f, C.f, x, p, c.cache_f)
    # See whether we can simplify things
    if C.g_has_parameters && C.f_has_parameters
        differentiate_parameters!(c.Jp_g, C.g, c.eval_f, p, c.cache_g)
        differentiate_parameters!(c.Jp_f, C.f, x, p, c.cache_f)
        jacobian!(c.J_g, C.g, c.eval_f, p, c.cache_g)
        LinearAlgebra.mul!(U, c.J_g, c.Jp_f)
        U .= U .+ c.Jp_g
    elseif C.g_has_parameters # && !C.f_has_parameters
        differentiate_parameters!(c.Jp_g, C.g, c.eval_f, p, c.cache_g)
        U .= c.Jp_g
    elseif C.f_has_parameters  # !C.g_has_parameters &&
        differentiate_parameters!(c.Jp_f, C.f, x, p, c.cache_f)
        jacobian!(c.J_g, C.g, c.eval_f, p, c.cache_g)
        LinearAlgebra.mul!(U, c.J_g, c.Jp_f)
    else
        U .= zero(eltype(U))
    end

    U
end

function differentiate_parameters(F::CompositionSystem, x, p, c::CompositionSystemCache)
    U = similar(c.J_g, size(c.J_g, 1), size(c.Jp_f, 2))
    differentiate_parameters!(U, F, x, p, c)
end


text:
module Day10
export parse_input
export find_min_area

import PyPlot

function parse_input(input::String)
    lines = split(strip(input), '\n')
    position = reshape(Int[], 0, 2)
    velocity = reshape(Int[], 0, 2)
    for line in lines
       m = match(r"position=<( *-?\d+), *(-?\d+)> velocity=< *(-?\d+), *(-?\d+)>", line)
       data = map(x -> parse(Int, x), m.captures)
       position = [position; [data[1] data[2]]]
       velocity = [velocity; [data[3] data[4]]]
    end

    return position, velocity
end

function find_min_area(positions, velocities)
    current_area::Int64 = 2^62
    last_area::Int64 = 2^62
    t = -1
    while current_area <= last_area
        positions .+= velocities
        t += 1
        width = maximum(positions[:, 1]) - minimum(positions[:, 1])
        height = maximum(positions[:, 2]) - minimum(positions[:, 2])
        last_area = current_area
        current_area = width * height
    end

    println(t)
    positions .-= velocities

    PyPlot.axis("equal")
    PyPlot.scatter(positions[:, 1], -positions[:, 2])
    PyPlot.show()
end

end


text:
# LibGit2.jl "diff" example - shows how to use the diff API
#
# Written by the LibGit2.jl contributors, based on the libgit2 examples.
#
# To the extent possible under law, the author(s) have dedicated all copyright
# and related and neighboring rights to this software to the public domain
# worldwide. This software is distributed without any warranty.
#
# You should have received a copy of the CC0 Public Domain Dedication along
# with this software. If not, see
# <http://creativecommons.org/publicdomain/zero/1.0/>.

using LibGit2
using ArgParse

# This example demonstrates the use of the LibGit2.jl diff APIs
# emulating a number of core Git `diff` command line options.

const COLORS = {
    :reset => "\033[m",
    :bold =>  "\033[1m",
    :red =>   "\033[31m",
    :green => "\033[32m",
    :cyan =>  "\033[36m"
}

function parse_commandline()
    settings = ArgParseSettings(autofix_names=true)
    @add_arg_table settings begin
        "treeish1"
            help = "oid of the first treeish"

        "treeish2"
            help = "oid of the second treeish"

        "path"
            help = "path of the repository"

        "-p", "-u", "--patch"
            action = :store_true

        "--cached"
            action = :store_true

        "--name-only"
            action = :store_true

        "--name-status"
            action = :store_true

        "--raw"
            action = :store_true

        "--color"
            action = :store_true

        "-R"
            action = :store_true
            dest_name = "reverse"

        "-a", "--text"
            action = :store_true
            dest_name = "force_text"

        "--ignore-space-at-eol"
            action = :store_true
            dest_name = "ignore_whitespace_eol"

        "-b", "--ignore-space-change"
            action = :store_true
            dest_name = "ignore_whitespace_change"

        "-w", "--ignore-all-space"
            action = :store_true
            dest_name = "ignore_whitespace"

        "--ignored" # not an option of git-diff
            action = :store_true
            dest_name = "include_ignored"

        "--untracked" # not an option of git-diff
            action = :store_true
            dest_name = "include_untracked"

        "--patience"
            action = :store_true

        "--minimal"
            action = :store_true

        "--numstat"
            action = :store_true

        "--shortstat"
            action = :store_true
        #=
        "-M", "--find-renames"
          arg_type = Integer

        "-C", "--find-copies"
          arg_type = Integer

        "--find-copies-harder"
          action = :store_true
          dest_name = "copies_from_unmodified"

        "-B", "--break-rewrites"
          action = :store_true
          dest_name = "rewrites"
        =#

        "-U", "--unified"
            arg_type = Integer
            dest_name = "context_lines"

        "--inter-hunk-context"
            arg_type = Integer
            dest_name = "interhunk_lines"
        #=
        "--old", "--src-prefix"
          arg_type = String
          dest_name = "old_prefix"

        "--dst-prefix"
          arg_type = String
          dest_name = "new_prefix"
        =#
    end

    diff_flags = (:raw, :patch, :name_only, :name_status)
    diff_format = :patch

    args = Dict()
    for (k,v) in parse_args(settings)
        is(v, nothing) && continue
        args[symbol(k)] = v
        if k in diff_flags
            diff_flags = v
        end
    end
    return args, diff_format
end

function diff_print_shortstat(d::GitDiff)
    s = LibGit2.stat(d)
    if s.files > 0
        print(" $(s.files) ")
        @printf("%s", 1 == s.files ? "file changed" : "files changed")
        if s.adds > 0
            print(", $(s.adds) ")
            @printf("%s", 1 == s.adds ? "insertion(+)" : "insertions(+)")
        end
        if s.dels > 0
            print(", $(s.dels) ")
            @printf("%s", 1 == s.dels ? "deletion(-)" : "deletions(-)")
        end
        println()
    end
end

function diff_print_numstat(d::GitDiff)
    p = LibGit2.patches(d)
    is(p, nothing) && return
    for p in patches
        s = LibGit2.stat(p)
        d = LibGit2.delta(p)
        @printf("%d\t%d\t%s\n", s.adds, s.dels, d.new_file.path)
    end
end

function color_printer(d)
    color, last_color = :reset, nothing
    for p in LibGit2.patches(d)
        for h in LibGit2.hunks(p)
            # color only text between `@@`
            idx = search(h.header,"@@", 3)
            print(COLORS[:cyan],  h.header[1:idx[end]],
                  COLORS[:reset], h.header[idx[end]+1:end])

            #= FIXME I miss the header like which should be highlited in bold:
            ```
            diff --git a/src/diff.jl b/src/diff.jl
            index c27da2b..a01abcd 100644
            --- a/src/diff.jl
            +++ b/src/diff.jl
            ```
            to recover the file I have to pass through delta
            d = delta(hj.patch)
            d.new_file.path
            =#

            for l in lines(h)
                lo = l.line_origin
                if lo == :addition || lo == :eof_newline_added
                    color = :green
                elseif lo == :deletion || lo == :eof_newline_removed
                    color = :red
                #TODO: elseif lo == :file_header
                else
                    color = :reset
                end
                if color != last_color
                    if last_color == :bold || color == :bold
                        print(COLORS[:reset]) # FIXME
                    end
                    print(COLORS[color])
                    last_color = color
                end
                lo == :context && print(" ")
                lo == :addition && print("+")
                lo == :deletion && print("-")
                print(l.content)
            end
        end
    end
    return color
end

function main()
    o, diff_format = parse_commandline()
    repo = haskey(o, :path) ? LibGit2.GitRepo(o[:path]) :
                              LibGit2.repo_discover()
    if haskey(o, :treeish1) && haskey(o, :treeish2)
        println("diff trees")
        d = LibGit2.diff(repo, o[:treeish1], o[:treeish2], o)
    elseif haskey(o, :treeish1) && o[:cached]
        println("diff tree to index")
        d = LibGit2.diff(repo, o[:treeish1], GitIndex(repo), o)
    elseif haskey(o, :treeish1)
        println("diff tree to working directory")
        d = LibGit2.diff_workdir(repo, o[:treeish1], o)
    elseif haskey(o, :cached)
        println("diff head to index")
        d = LibGit2.diff(repo, "HEAD", GitIndex(repo), o)
    else
        println("diff index to working directory")
        d = LibGit2.diff_workdir(repo, o)
    end
    if o[:numstat]
        diff_print_numstat(d)
    elseif o[:shortstat]
        diff_print_shortstat(d)
    elseif o[:color]
        color_printer(d)
    else
        println(LibGit2.patch(d, format=diff_format))
    end
end

main()


text:
function mdl_Zload!(voltage_r,
                    voltage_i,
                    current_r,
                    current_i,
                    device::PSY.PowerLoad,
                    sys::DynamicSystem)

    #Get parameters
    P = device.activepower
    Q = device.reactivepower

    #Compute impedance and load
    #Z_load =  1.0/(P-Q*1im) #in pu
    #I = -(voltage_r[1] + voltage_i[1]*1im)/Z_load #in pu flowing out
    #current_r[1] += real(I)
    #current_i[1] += imag(I)

    #Compute current given constant power
    #I = (P-Q*1im)/(voltage_r[1] - voltage_i[1]*1im)
    #current_r[1] += real(I)
    #current_i[1] += imag(I)

    #Update current
    current_r[1] += -(voltage_r[1]*P + voltage_i[1]*Q) #in system pu flowing out
    current_i[1] += -(voltage_i[1]*P - voltage_r[1]*Q) #in system pu flowing out

    return
end


text:
const __bodyfunction__ = Dict{Method,Any}()

# Find keyword "body functions" (the function that contains the body
# as written by the developer, called after all missing keyword-arguments
# have been assigned values), in a manner that doesn't depend on
# gensymmed names.
# `mnokw` is the method that gets called when you invoke it without
# supplying any keywords.
function __lookup_kwbody__(mnokw::Method)
    function getsym(arg)
        isa(arg, Symbol) && return arg
        @assert isa(arg, GlobalRef)
        return arg.name
    end

    f = get(__bodyfunction__, mnokw, nothing)
    if f === nothing
        fmod = mnokw.module
        # The lowered code for `mnokw` should look like
        #   %1 = mkw(kwvalues..., #self#, args...)
        #        return %1
        # where `mkw` is the name of the "active" keyword body-function.
        ast = Base.uncompressed_ast(mnokw)
        if isa(ast, Core.CodeInfo) && length(ast.code) >= 2
            callexpr = ast.code[end-1]
            if isa(callexpr, Expr) && callexpr.head == :call
                fsym = callexpr.args[1]
                if isa(fsym, Symbol)
                    f = getfield(fmod, fsym)
                elseif isa(fsym, GlobalRef)
                    if fsym.mod === Core && fsym.name === :_apply
                        f = getfield(mnokw.module, getsym(callexpr.args[2]))
                    elseif fsym.mod === Core && fsym.name === :_apply_iterate
                        f = getfield(mnokw.module, getsym(callexpr.args[3]))
                    else
                        f = getfield(fsym.mod, fsym.name)
                    end
                else
                    f = missing
                end
            else
                f = missing
            end
        else
            f = missing
        end
        __bodyfunction__[mnokw] = f
    end
    return f
end

function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    isdefined(Base, Symbol("#630#631")) && precompile(Tuple{getfield(Base, Symbol("#630#631"))})
    isdefined(Base.Broadcast, Symbol("#31#32")) && precompile(Tuple{getfield(Base.Broadcast, Symbol("#31#32")),ArchiveSource})
    isdefined(Base.Broadcast, Symbol("#31#32")) && precompile(Tuple{getfield(Base.Broadcast, Symbol("#31#32")),DirectorySource})
    isdefined(BinaryBuilder.Auditor, Symbol("#110#113")) && precompile(Tuple{getfield(BinaryBuilder.Auditor, Symbol("#110#113")),String})
    isdefined(BinaryBuilder.Auditor, Symbol("#32#45")) && precompile(Tuple{getfield(BinaryBuilder.Auditor, Symbol("#32#45")),String})
    isdefined(BinaryBuilder.Auditor, Symbol("#33#46")) && precompile(Tuple{getfield(BinaryBuilder.Auditor, Symbol("#33#46")),String})
    isdefined(BinaryBuilder.Auditor, Symbol("#33#46")) && precompile(Tuple{getfield(BinaryBuilder.Auditor, Symbol("#33#46")),SubString{String}})
    isdefined(BinaryBuilder.Auditor, Symbol("#92#98")) && precompile(Tuple{getfield(BinaryBuilder.Auditor, Symbol("#92#98")),String})
    isdefined(BinaryBuilder.Auditor, Symbol("#check_set#87")) && precompile(Tuple{getfield(BinaryBuilder.Auditor, Symbol("#check_set#87")),String,Array{String,1}})
    isdefined(BinaryBuilder.BinaryBuilderBase, Symbol("#109#114")) && precompile(Tuple{getfield(BinaryBuilder.BinaryBuilderBase, Symbol("#109#114")),String})
    isdefined(BinaryBuilder.BinaryBuilderBase, Symbol("#110#115")) && precompile(Tuple{getfield(BinaryBuilder.BinaryBuilderBase, Symbol("#110#115")),String})
    isdefined(BinaryBuilder.BinaryBuilderBase, Symbol("#111#116")) && precompile(Tuple{getfield(BinaryBuilder.BinaryBuilderBase, Symbol("#111#116")),String})
    isdefined(BinaryBuilder.BinaryBuilderBase, Symbol("#171#172")) && precompile(Tuple{getfield(BinaryBuilder.BinaryBuilderBase, Symbol("#171#172"))})
    isdefined(BinaryBuilder.BinaryBuilderBase, Symbol("#GOOS#94")) && precompile(Tuple{getfield(BinaryBuilder.BinaryBuilderBase, Symbol("#GOOS#94")),Linux})
    isdefined(BinaryBuilder.BinaryBuilderBase, Symbol("#base_gcc_flags#76")) && precompile(Tuple{getfield(BinaryBuilder.BinaryBuilderBase, Symbol("#base_gcc_flags#76")),Linux})
    isdefined(BinaryBuilder.BinaryBuilderBase, Symbol("#cargo#101")) && precompile(Tuple{getfield(BinaryBuilder.BinaryBuilderBase, Symbol("#cargo#101")),IOStream,Linux})
    isdefined(BinaryBuilder.BinaryBuilderBase, Symbol("#cc#91")) && precompile(Tuple{getfield(BinaryBuilder.BinaryBuilderBase, Symbol("#cc#91")),IOStream,Linux})
    isdefined(BinaryBuilder.BinaryBuilderBase, Symbol("#clang#88")) && precompile(Tuple{getfield(BinaryBuilder.BinaryBuilderBase, Symbol("#clang#88")),IOStream,Linux})
    isdefined(BinaryBuilder.BinaryBuilderBase, Symbol("#clang_compile_flags#81")) && precompile(Tuple{getfield(BinaryBuilder.BinaryBuilderBase, Symbol("#clang_compile_flags#81")),Linux})
    isdefined(BinaryBuilder.BinaryBuilderBase, Symbol("#clang_flags#80")) && precompile(Tuple{getfield(BinaryBuilder.BinaryBuilderBase, Symbol("#clang_flags#80")),Linux})
    isdefined(BinaryBuilder.BinaryBuilderBase, Symbol("#clang_link_flags#82")) && precompile(Tuple{getfield(BinaryBuilder.BinaryBuilderBase, Symbol("#clang_link_flags#82")),Linux})
    isdefined(BinaryBuilder.BinaryBuilderBase, Symbol("#clangxx#89")) && precompile(Tuple{getfield(BinaryBuilder.BinaryBuilderBase, Symbol("#clangxx#89")),IOStream,Linux})
    isdefined(BinaryBuilder.BinaryBuilderBase, Symbol("#cxx#92")) && precompile(Tuple{getfield(BinaryBuilder.BinaryBuilderBase, Symbol("#cxx#92")),IOStream,Linux})
    isdefined(BinaryBuilder.BinaryBuilderBase, Symbol("#cxxfilt#104")) && precompile(Tuple{getfield(BinaryBuilder.BinaryBuilderBase, Symbol("#cxxfilt#104")),IOStream,Linux})
    isdefined(BinaryBuilder.BinaryBuilderBase, Symbol("#fortran_flags#78")) && precompile(Tuple{getfield(BinaryBuilder.BinaryBuilderBase, Symbol("#fortran_flags#78")),Linux})
    isdefined(BinaryBuilder.BinaryBuilderBase, Symbol("#gcc_flags#77")) && precompile(Tuple{getfield(BinaryBuilder.BinaryBuilderBase, Symbol("#gcc_flags#77")),Linux})
    isdefined(BinaryBuilder.BinaryBuilderBase, Symbol("#gcc_link_flags#83")) && precompile(Tuple{getfield(BinaryBuilder.BinaryBuilderBase, Symbol("#gcc_link_flags#83")),Linux})
    isdefined(BinaryBuilder.BinaryBuilderBase, Symbol("#gfortran#87")) && precompile(Tuple{getfield(BinaryBuilder.BinaryBuilderBase, Symbol("#gfortran#87")),IOStream,Linux})
    isdefined(BinaryBuilder.BinaryBuilderBase, Symbol("#go#96")) && precompile(Tuple{getfield(BinaryBuilder.BinaryBuilderBase, Symbol("#go#96")),IOStream,Linux})
    isdefined(BinaryBuilder.BinaryBuilderBase, Symbol("#meson#102")) && precompile(Tuple{getfield(BinaryBuilder.BinaryBuilderBase, Symbol("#meson#102")),IOStream,Linux})
    isdefined(BinaryBuilder.BinaryBuilderBase, Symbol("#objc#90")) && precompile(Tuple{getfield(BinaryBuilder.BinaryBuilderBase, Symbol("#objc#90")),IOStream,Linux})
    isdefined(BinaryBuilder.BinaryBuilderBase, Symbol("#patchelf#103")) && precompile(Tuple{getfield(BinaryBuilder.BinaryBuilderBase, Symbol("#patchelf#103")),IOStream,Linux})
    isdefined(BinaryBuilder.BinaryBuilderBase, Symbol("#rust_flags#98")) && precompile(Tuple{getfield(BinaryBuilder.BinaryBuilderBase, Symbol("#rust_flags#98")),Linux})
    isdefined(BinaryBuilder.BinaryBuilderBase, Symbol("#rustc#99")) && precompile(Tuple{getfield(BinaryBuilder.BinaryBuilderBase, Symbol("#rustc#99")),IOStream,Linux})
    isdefined(BinaryBuilder.BinaryBuilderBase, Symbol("#rustup#100")) && precompile(Tuple{getfield(BinaryBuilder.BinaryBuilderBase, Symbol("#rustup#100")),IOStream,Linux})
    isdefined(BinaryBuilder.OutputCollectors, Symbol("#1#3")) && precompile(Tuple{getfield(BinaryBuilder.OutputCollectors, Symbol("#1#3"))})
    isdefined(BinaryBuilder.OutputCollectors, Symbol("#14#15")) && precompile(Tuple{getfield(BinaryBuilder.OutputCollectors, Symbol("#14#15"))})
    isdefined(BinaryBuilder.OutputCollectors, Symbol("#2#4")) && precompile(Tuple{getfield(BinaryBuilder.OutputCollectors, Symbol("#2#4"))})
    let fbody = try __lookup_kwbody__(which(BinaryBuilder.Auditor.audit, (Prefix,String,))) catch missing end
        if !ismissing(fbody)
            precompile(fbody, (Base.TTY,Linux,Bool,Bool,Bool,Bool,Bool,typeof(audit),Prefix,String,))
        end
    end
    let fbody = try __lookup_kwbody__(which(BinaryBuilder.Auditor.check_cxxstring_abi, (ObjectFile.ELF.ELFHandle{IOStream},Linux,))) catch missing end
        if !ismissing(fbody)
            precompile(fbody, (Base.TTY,Bool,typeof(BinaryBuilder.Auditor.check_cxxstring_abi),ObjectFile.ELF.ELFHandle{IOStream},Linux,))
        end
    end
    let fbody = try __lookup_kwbody__(which(BinaryBuilder.BinaryBuilderBase.compress_dir, (String,))) catch missing end
        if !ismissing(fbody)
            precompile(fbody, (Type{T} where T,Int64,String,Bool,typeof(BinaryBuilder.BinaryBuilderBase.compress_dir),String,))
        end
    end
    let fbody = try __lookup_kwbody__(which(BinaryBuilder.BinaryBuilderBase.resolve_jlls, (Array{Dependency,1},))) catch missing end
        if !ismissing(fbody)
            precompile(fbody, (Pkg.Types.Context,Base.TTY,typeof(BinaryBuilder.BinaryBuilderBase.resolve_jlls),Array{Dependency,1},))
        end
    end
    let fbody = try __lookup_kwbody__(which(BinaryBuilder.build_tarballs, (Any,Any,Any,Any,Any,Any,Any,Any,))) catch missing end
        if !ismissing(fbody)
            precompile(fbody, (Any,typeof(build_tarballs),Any,Any,Any,Any,Any,Any,Any,Any,))
        end
    end
    let fbody = try __lookup_kwbody__(which(run, (BinaryBuilder.BinaryBuilderBase.UserNSRunner,Cmd,IOStream,))) catch missing end
        if !ismissing(fbody)
            precompile(fbody, (Bool,Base.TTY,typeof(run),BinaryBuilder.BinaryBuilderBase.UserNSRunner,Cmd,IOStream,))
        end
    end
    precompile(Tuple{Core.kwftype(typeof(Base.Broadcast.broadcasted_kwsyntax)),NamedTuple{(:verbose,),Tuple{Bool}},typeof(Base.Broadcast.broadcasted_kwsyntax),Function,Array{BinaryBuilder.BinaryBuilderBase.AbstractSource,1}})
    precompile(Tuple{Core.kwftype(typeof(BinaryBuilder.Auditor.check_cxxstring_abi)),NamedTuple{(:verbose,),Tuple{Bool}},typeof(BinaryBuilder.Auditor.check_cxxstring_abi),ObjectFile.ELF.ELFHandle{IOStream},Linux})
    precompile(Tuple{Core.kwftype(typeof(BinaryBuilder.Auditor.check_dynamic_linkage)),NamedTuple{(:platform, :silent, :verbose, :autofix),Tuple{Linux,Bool,Bool,Bool}},typeof(BinaryBuilder.Auditor.check_dynamic_linkage),ObjectFile.ELF.ELFHandle{IOStream},Prefix,Array{String,1}})
    precompile(Tuple{Core.kwftype(typeof(BinaryBuilder.Auditor.check_isa)),NamedTuple{(:verbose, :silent),Tuple{Bool,Bool}},typeof(BinaryBuilder.Auditor.check_isa),ObjectFile.ELF.ELFHandle{IOStream},Linux,Prefix})
    precompile(Tuple{Core.kwftype(typeof(BinaryBuilder.Auditor.check_libgfortran_version)),NamedTuple{(:verbose, :has_csl),Tuple{Bool,Bool}},typeof(BinaryBuilder.Auditor.check_libgfortran_version),ObjectFile.ELF.ELFHandle{IOStream},Linux})
    precompile(Tuple{Core.kwftype(typeof(BinaryBuilder.Auditor.check_libgomp)),NamedTuple{(:verbose, :has_csl),Tuple{Bool,Bool}},typeof(BinaryBuilder.Auditor.check_libgomp),ObjectFile.ELF.ELFHandle{IOStream},Linux})
    precompile(Tuple{Core.kwftype(typeof(BinaryBuilder.Auditor.check_license)),NamedTuple{(:verbose, :silent),Tuple{Bool,Bool}},typeof(BinaryBuilder.Auditor.check_license),Prefix,String})
    precompile(Tuple{Core.kwftype(typeof(BinaryBuilder.Auditor.check_os_abi)),NamedTuple{(:verbose,),Tuple{Bool}},typeof(BinaryBuilder.Auditor.check_os_abi),ObjectFile.ELF.ELFHandle{IOStream},Linux})
    precompile(Tuple{Core.kwftype(typeof(BinaryBuilder.Auditor.ensure_soname)),NamedTuple{(:verbose, :autofix),Tuple{Bool,Bool}},typeof(BinaryBuilder.Auditor.ensure_soname),Prefix,String,Linux})
    precompile(Tuple{Core.kwftype(typeof(BinaryBuilder.Auditor.relink_to_rpath)),NamedTuple{(:verbose,),Tuple{Bool}},typeof(BinaryBuilder.Auditor.relink_to_rpath),Prefix,Linux,String,SubString{String}})
    precompile(Tuple{Core.kwftype(typeof(BinaryBuilder.Auditor.symlink_soname_lib)),NamedTuple{(:verbose, :autofix),Tuple{Bool,Bool}},typeof(BinaryBuilder.Auditor.symlink_soname_lib),String})
    precompile(Tuple{Core.kwftype(typeof(BinaryBuilder.Auditor.translate_symlinks)),NamedTuple{(:verbose,),Tuple{Bool}},typeof(BinaryBuilder.Auditor.translate_symlinks),String})
    precompile(Tuple{Core.kwftype(typeof(BinaryBuilder.Auditor.update_linkage)),NamedTuple{(:verbose,),Tuple{Bool}},typeof(BinaryBuilder.Auditor.update_linkage),Prefix,Linux,String,SubString{String},String})
    precompile(Tuple{Core.kwftype(typeof(BinaryBuilder.BinaryBuilderBase.Type)),NamedTuple{(:cwd, :platform),Tuple{String,Linux}},Type{BinaryBuilder.BinaryBuilderBase.UserNSRunner},String})
    precompile(Tuple{Core.kwftype(typeof(BinaryBuilder.BinaryBuilderBase.Type)),NamedTuple{(:cwd, :platform, :verbose),Tuple{String,Linux,Bool}},Type{BinaryBuilder.BinaryBuilderBase.UserNSRunner},String})
    precompile(Tuple{Core.kwftype(typeof(BinaryBuilder.BinaryBuilderBase.Type)),NamedTuple{(:cwd, :platform, :verbose, :workspaces, :compiler_wrapper_dir, :src_name, :shards, :compilers),Tuple{String,Linux,Bool,Array{Pair{String,String},1},String,String,Array{BinaryBuilder.BinaryBuilderBase.CompilerShard,1},Array{Symbol,1}}},Type{BinaryBuilder.BinaryBuilderBase.UserNSRunner},String})
    precompile(Tuple{Core.kwftype(typeof(BinaryBuilder.BinaryBuilderBase.choose_shards)),NamedTuple{(:compilers,),Tuple{Array{Symbol,1}}},typeof(BinaryBuilder.BinaryBuilderBase.choose_shards),Linux})
    precompile(Tuple{Core.kwftype(typeof(BinaryBuilder.BinaryBuilderBase.generate_compiler_wrappers!)),NamedTuple{(:bin_path, :compilers),Tuple{String,Array{Symbol,1}}},typeof(BinaryBuilder.BinaryBuilderBase.generate_compiler_wrappers!),Linux})
    precompile(Tuple{Core.kwftype(typeof(BinaryBuilder.BinaryBuilderBase.generate_compiler_wrappers!)),NamedTuple{(:bin_path,),Tuple{String}},typeof(BinaryBuilder.BinaryBuilderBase.generate_compiler_wrappers!),Linux})
    precompile(Tuple{Core.kwftype(typeof(BinaryBuilder.BinaryBuilderBase.satisfied)),NamedTuple{(:verbose, :platform),Tuple{Bool,Linux}},typeof(satisfied),ExecutableProduct,Prefix})
    precompile(Tuple{Core.kwftype(typeof(BinaryBuilder.BinaryBuilderBase.setup_workspace)),NamedTuple{(:verbose,),Tuple{Bool}},typeof(BinaryBuilder.BinaryBuilderBase.setup_workspace),String,Array{BinaryBuilder.BinaryBuilderBase.SetupSource,1}})
    precompile(Tuple{Core.kwftype(typeof(BinaryBuilder.autobuild)),NamedTuple{(:verbose, :debug, :meta_json_stream),Tuple{Bool,Bool,IOStream}},typeof(autobuild),String,String,VersionNumber,Array{DirectorySource,1},String,Array{Linux,1},Array{Product,1},Array{Dependency,1}})
    precompile(Tuple{Core.kwftype(typeof(BinaryBuilder.autobuild)),NamedTuple{(:verbose, :debug, :meta_json_stream, :compilers),Tuple{Bool,Bool,Nothing,Array{Symbol,1}}},typeof(autobuild),String,String,VersionNumber,Array{BinaryBuilder.BinaryBuilderBase.AbstractSource,1},String,Array{Linux,1},Array{Product,1},Array{Dependency,1}})
    precompile(Tuple{Core.kwftype(typeof(BinaryBuilder.build_tarballs)),NamedTuple{(:compilers,),Tuple{Array{Symbol,1}}},typeof(build_tarballs),Array{String,1},String,VersionNumber,Array{BinaryBuilder.BinaryBuilderBase.AbstractSource,1},String,Array{Linux,1},Array{Product,1},Array{Dependency,1}})
    precompile(Tuple{Core.kwftype(typeof(BinaryBuilder.get_compilers_versions)),NamedTuple{(:compilers,),Tuple{Array{Symbol,1}}},typeof(BinaryBuilder.get_compilers_versions)})
    precompile(Tuple{Type{ArchiveSource},String,String})
    precompile(Tuple{Type{Array{BinaryBuilder.BinaryBuilderBase.AbstractSource,1}},UndefInitializer,Int64})
    precompile(Tuple{Type{Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1},Axes,F,Args} where Args<:Tuple where F where Axes},typeof(BinaryBuilder.BinaryBuilderBase.coerce_dependency),Tuple{Array{Dependency,1}}})
    precompile(Tuple{Type{Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1},Axes,F,Args} where Args<:Tuple where F where Axes},typeof(BinaryBuilder.BinaryBuilderBase.coerce_source),Tuple{Array{DirectorySource,1}}})
    precompile(Tuple{Type{BinaryBuilder.OutputCollectors.OutputCollector},Cmd,Base.Process,BinaryBuilder.OutputCollectors.LineStream,BinaryBuilder.OutputCollectors.LineStream,Base.GenericCondition{Base.AlwaysLockedST},Base.DevNull,Bool,Bool})
    precompile(Tuple{Type{BinaryBuilder.OutputCollectors.OutputCollector},Cmd,Base.Process,BinaryBuilder.OutputCollectors.LineStream,BinaryBuilder.OutputCollectors.LineStream,Base.GenericCondition{Base.AlwaysLockedST},Base.TTY,Bool,Bool})
    precompile(Tuple{Type{Dependency},String})
    precompile(Tuple{Type{DirectorySource},String})
    precompile(Tuple{Type{ExecutableProduct},String,Symbol})
    precompile(Tuple{typeof(Base.Broadcast.broadcasted),Function,Array{BinaryBuilder.BinaryBuilderBase.AbstractSource,1}})
    precompile(Tuple{typeof(Base.Broadcast.broadcasted),Function,Array{Dependency,1}})
    precompile(Tuple{typeof(Base.Broadcast.broadcasted),Function,Array{DirectorySource,1}})
    precompile(Tuple{typeof(Base.Broadcast.broadcasted),Function,Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1},Nothing,typeof(BinaryBuilder.BinaryBuilderBase.getname),Tuple{Array{Dependency,1}}},String})
    precompile(Tuple{typeof(Base.Broadcast.combine_eltypes),Function,Tuple{Array{BinaryBuilder.BinaryBuilderBase.AbstractSource,1}}})
    precompile(Tuple{typeof(Base.Broadcast.copyto_nonleaf!),Array{ArchiveSource,1},Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1},Tuple{Base.OneTo{Int64}},typeof(BinaryBuilder.BinaryBuilderBase.coerce_source),Tuple{Base.Broadcast.Extruded{Array{BinaryBuilder.BinaryBuilderBase.AbstractSource,1},Tuple{Bool},Tuple{Int64}}}},Base.OneTo{Int64},Int64,Int64})
    precompile(Tuple{typeof(Base.Broadcast.instantiate),Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1},Nothing,typeof(BinaryBuilder.BinaryBuilderBase.coerce_dependency),Tuple{Array{Dependency,1}}}})
    precompile(Tuple{typeof(Base.Broadcast.instantiate),Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1},Nothing,typeof(BinaryBuilder.BinaryBuilderBase.coerce_source),Tuple{Array{DirectorySource,1}}}})
    precompile(Tuple{typeof(Base.Broadcast.restart_copyto_nonleaf!),Array{BinaryBuilder.BinaryBuilderBase.AbstractSource,1},Array{ArchiveSource,1},Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1},Tuple{Base.OneTo{Int64}},typeof(BinaryBuilder.BinaryBuilderBase.coerce_source),Tuple{Base.Broadcast.Extruded{Array{BinaryBuilder.BinaryBuilderBase.AbstractSource,1},Tuple{Bool},Tuple{Int64}}}},DirectorySource,Int64,Base.OneTo{Int64},Int64,Int64})
    precompile(Tuple{typeof(Base._compute_eltype),Type{Tuple{Pair{String,String},Pair{String,String},Pair{String,Array{DirectorySource,1}},Pair{String,String},Pair{String,Array{Product,1}},Pair{String,Array{Dependency,1}},Pair{String,Bool}}}})
    precompile(Tuple{typeof(Base.allocatedinline),Type{BinaryBuilder.BinaryBuilderBase.CompilerShard}})
    precompile(Tuple{typeof(Base.allocatedinline),Type{Dependency}})
    precompile(Tuple{typeof(Base.grow_to!),Dict{String,Any},Tuple{Pair{String,String},Pair{String,String},Pair{String,Array{DirectorySource,1}},Pair{String,String},Pair{String,Array{Product,1}},Pair{String,Array{Dependency,1}},Pair{String,Bool}},Int64})
    precompile(Tuple{typeof(Base.grow_to!),Dict{String,String},Tuple{Pair{String,String},Pair{String,String},Pair{String,Array{DirectorySource,1}},Pair{String,String},Pair{String,Array{Product,1}},Pair{String,Array{Dependency,1}},Pair{String,Bool}},Int64})
    precompile(Tuple{typeof(Base.merge_types),NTuple{8,Symbol},Type{NamedTuple{(:cwd, :platform, :verbose, :workspaces, :compiler_wrapper_dir, :src_name, :shards),Tuple{String,Linux,Bool,Array{Pair{String,String},1},String,String,Array{BinaryBuilder.BinaryBuilderBase.CompilerShard,1}}}},Type{NamedTuple{(:compilers,),Tuple{Array{Symbol,1}}}}})
    precompile(Tuple{typeof(Base.promote_typeof),ArchiveSource,DirectorySource})
    precompile(Tuple{typeof(Base.vect),ArchiveSource,Vararg{Any,N} where N})
    precompile(Tuple{typeof(BinaryBuilder.Auditor.cppfilt),Array{SubString{String},1},Linux})
    precompile(Tuple{typeof(BinaryBuilder.Auditor.is_default_lib),SubString{String},ObjectFile.ELF.ELFHandle{IOStream}})
    precompile(Tuple{typeof(BinaryBuilder.Auditor.is_troublesome_library_link),SubString{String},Linux})
    precompile(Tuple{typeof(BinaryBuilder.Auditor.minimum_instruction_set),Dict{String,Int64},Bool})
    precompile(Tuple{typeof(BinaryBuilder.Auditor.should_ignore_lib),SubString{String},ObjectFile.ELF.ELFHandle{IOStream}})
    precompile(Tuple{typeof(BinaryBuilder.BinaryBuilderBase.storage_dir),String,Vararg{String,N} where N})
    precompile(Tuple{typeof(BinaryBuilder.BinaryBuilderBase.with_logfile),Function,Prefix,String})
    precompile(Tuple{typeof(BinaryBuilder.check_flag!),Array{String,1},String})
    precompile(Tuple{typeof(BinaryBuilder.extract_flag!),Array{String,1},String,String})
    precompile(Tuple{typeof(BinaryBuilder.extract_flag!),Array{String,1},String})
    precompile(Tuple{typeof(JSON.Writer.show_element),JSON.Writer.CompactContext{Base.GenericIOBuffer{Array{UInt8,1}}},JSON.Serializations.StandardSerialization,ExecutableProduct})
    precompile(Tuple{typeof(JSON.Writer.show_pair),JSON.Writer.CompactContext{Base.GenericIOBuffer{Array{UInt8,1}}},JSON.Serializations.StandardSerialization,String,Array{Dependency,1}})
    precompile(Tuple{typeof(JSON.Writer.show_pair),JSON.Writer.CompactContext{Base.GenericIOBuffer{Array{UInt8,1}}},JSON.Serializations.StandardSerialization,String,Array{DirectorySource,1}})
    precompile(Tuple{typeof(JSON.Writer.show_pair),JSON.Writer.CompactContext{Base.GenericIOBuffer{Array{UInt8,1}}},JSON.Serializations.StandardSerialization,String,Array{Product,1}})
    precompile(Tuple{typeof(build_tarballs),Any,Any,Any,Any,Any,Any,Any,Any})
    precompile(Tuple{typeof(copy),Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1},Tuple{Base.OneTo{Int64}},typeof(BinaryBuilder.BinaryBuilderBase.coerce_dependency),Tuple{Array{Dependency,1}}}})
    precompile(Tuple{typeof(copy),Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1},Tuple{Base.OneTo{Int64}},typeof(BinaryBuilder.BinaryBuilderBase.coerce_source),Tuple{Array{DirectorySource,1}}}})
    precompile(Tuple{typeof(copyto!),Array{BinaryBuilder.BinaryBuilderBase.AbstractSource,1},Tuple{ArchiveSource,DirectorySource}})
    precompile(Tuple{typeof(getindex),Type{Dependency},Dependency})
    precompile(Tuple{typeof(getindex),Type{Dependency}})
    precompile(Tuple{typeof(getindex),Type{DirectorySource}})
    precompile(Tuple{typeof(getindex),Type{Product},ExecutableProduct,ExecutableProduct,ExecutableProduct,ExecutableProduct,Vararg{ExecutableProduct,N} where N})
    precompile(Tuple{typeof(isequal),ExecutableProduct,ExecutableProduct})
    precompile(Tuple{typeof(merge),NamedTuple{(:cwd, :platform, :verbose, :workspaces, :compiler_wrapper_dir, :src_name, :shards),Tuple{String,Linux,Bool,Array{Pair{String,String},1},String,String,Array{BinaryBuilder.BinaryBuilderBase.CompilerShard,1}}},NamedTuple{(:compilers,),Tuple{Array{Symbol,1}}}})
    precompile(Tuple{typeof(setindex!),Array{ArchiveSource,1},ArchiveSource,Int64})
    precompile(Tuple{typeof(setindex!),Array{BinaryBuilder.BinaryBuilderBase.SetupSource{ArchiveSource},1},BinaryBuilder.BinaryBuilderBase.SetupSource{ArchiveSource},Int64})
    precompile(Tuple{typeof(setindex!),Dict{Product,Any},Dict{String,String},ExecutableProduct})
    precompile(Tuple{typeof(setindex!),Dict{String,Any},Array{Dependency,1},String})
    precompile(Tuple{typeof(setindex!),Dict{String,Any},Array{DirectorySource,1},String})
    precompile(Tuple{typeof(setindex!),Dict{String,Any},Array{Product,1},String})
    precompile(Tuple{typeof(similar),Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1},Tuple{Base.OneTo{Int64}},typeof(BinaryBuilder.BinaryBuilderBase.coerce_source),Tuple{Base.Broadcast.Extruded{Array{BinaryBuilder.BinaryBuilderBase.AbstractSource,1},Tuple{Bool},Tuple{Int64}}}},Type{ArchiveSource}})
    precompile(Tuple{typeof(strip),IOStream,Linux})
end


text:
using CairoMakie
using DiffEqSensitivity
using Zygote
using ForwardDiff
using LinearAlgebra
using DifferentialEquations

"""
  vanderpol!(du,u,μ,t)

Van der Pol oscillator dynamics for DifferentialEquations.jl
inline
"""
function vanderpol!(du,u,μ,t)
    x,y = u
    du[1] = μ * (x - x^3/3 - y)
    du[2] = x/μ
end

"""
  vanderpol(u,μ,t)

Van der Pol oscillator dynamics for DifferentialEquations.jl
not inline
"""
function vanderpol(u,μ,t)
    du = zeros(2)
    vanderpol!(du,u,μ,t)
    return du
end

"""
  phasemap(u, positions)

Applies to Van der Pol oscillator.

Assigns a point `u` in the statespace to a position on the
limit cycle of the oscillator. Given a list of `positions`
spaced equally in time along the limit cycle, this function
simulates the system from initial condition `u` then determines
the closest element in `positions`

Returns: The index of the closest element.
"""
function phasemap(u, positions)
    μ = 1.0
    tspan = (0.0,40.0)

    prob = ODEProblem(vanderpol!,u,tspan,μ)

    sol = solve(prob)
    final = sol.u[end]
    ~, out = findmin(sum(abs2,positions .- final,dims=1))
    return out[2]
end

"""
  simystem(u)

Compute the value of the Van der Pol oscillator initialized at
`u` 20 units of time in the future.
"""
function simsystem(u)
    μ = 1.0
    tspan = (0.0,20.0)

    prob = ODEProblem(vanderpol!,u,tspan,μ)

    sol = solve(prob,rtol=1e-6,atol=1e-6,sensealg=QuadratureAdjoint())
    return sol.u[end]
end

"""
  informativeVec(u)

Find the non-zero eigenvalue of dφ₂₀ in the Van Der Pol
oscillator around state space point `u`.
"""
function informativeVec(u)
    dx = ForwardDiff.jacobian(simsystem,u)
    ~, Σ, Vt = svd(dx)
    # Force counterclockwise
    test_vec = [u[2], -u[1]]
    if dot(test_vec,Vt[1,:]) > 0
        return Vt[1,:]
    else
        return -Vt[1,:]
    end
end


text:
using BeastUtils.XMLConstructor
using BeastUtils.TreeUtils, PhyloNetworks

xc = XMLConstructor

using DataFrames, CSV, Test
### Making sure https://github.com/suchard-group/incomplete_measurements/blob/master/scripts/xml_setup.jl works

function make_mbd(data_path::String, newick_path::String, xml_path::String,
                filename::String; dates_path::String = "")

    df = DataFrame(CSV.File(data_path))

    use_dates = false
    if length(dates_path) > 0
        dates_df = DataFrame(CSV.File(dates_path))
        @assert dates_df[!, :taxon] == df[!, :taxon]
        use_dates = true
    end

    newick = read(newick_path, String)

    taxa, data = XMLConstructor.df_to_matrix(df)

    bx = XMLConstructor.make_residual_xml(data, taxa, newick, chain_length = 100_000)
    if use_dates
        XMLConstructor.use_dates!(bx)
        xc.set_data_dates(bx, dates_df[!, :date])
    end
    xc.set_screen_logEvery(bx, 100)
    xc.set_file_logEvery(bx, 10)
    xc.set_filename(bx, filename)
    XMLConstructor.add_MBD_loggables!(bx)

    XMLConstructor.save_xml(xml_path, bx)
end

n = 10
p = 4
taxa = ["taxon$i" for i = 1:n]
data = randn(n, p)
dates = rand(n)

df = DataFrame()
df.taxon = taxa
for i = 1:p
    df[!, Symbol("trait$i")] = data[:, i]
end

data_path = "data.csv"
dates_path = "dates.csv"
newick_path = "newick.txt"
xml_path = "XMLConstructor.xml"
dates_xml_path = "xml_dates.xml"
filename = "test"

CSV.write(data_path, df)
CSV.write(dates_path, DataFrame(taxon = taxa, date = dates))

newick = writeTopology(TreeUtils.rtree(taxa))
write(newick_path, newick)

make_mbd(data_path, newick_path, xml_path, filename)

@test isfile(xml_path)

make_mbd(data_path, newick_path, dates_xml_path, filename, dates_path = dates_path)

@test isfile(dates_xml_path)

### Testing Factor xml
k = 3

# HMC, no shrinkage
bx = XMLConstructor.make_pfa_xml(data, taxa, newick, k, useHMC = false)
XMLConstructor.save_xml("facHMC.xml", bx)
@test isfile("facHMC.xml")

# Gibbs, no shrinkage
bx = XMLConstructor.make_pfa_xml(data, taxa, newick, k, useHMC = true)
XMLConstructor.save_xml("facGibbs.xml", bx)
@test isfile("facGibbs.xml")

# # HMC, shrinkage
# bx = XMLConstructor.make_pfa_xml(data, taxa, newick, k, useHMC = true,
#             shrink_loadings = true)
# XMLConstructor.save_xml("facHMCShrink.xml", bx)
# @test isfile("facHMCShrink.xml")

# # Gibbs, shrinkage
# bx = XMLConstructor.make_pfa_xml(data, taxa, newick, k, useHMC = false,
#             shrink_loadings = true)
# XMLConstructor.save_xml("facGibbsShrink.xml", bx)
# @test isfile("facGibbsShrink.xml")

# # Rotate prior
# fn = "facGibbsRotate.xml"
# bx = XMLConstructor.make_pfa_xml(data, taxa, newick, k, useHMC = false,
#             shrink_loadings = true, rotate_prior = true)
# XMLConstructor.save_xml(fn, bx)
# @test isfile(fn)

################################################################################
## Orthogonal factor analysis
################################################################################

fn = "facOrthogonal.xml"
bx = XMLConstructor.make_orthogonal_pfa_xml(data, taxa, newick, k)
XMLConstructor.save_xml(fn, bx)
@test isfile(fn)

################################################################################
## Joint models
################################################################################

n = 10
k = 2
p_res = 3
p_fac = 5
p_diff = k + p_res

taxa = ["taxon_$i" for i = 1:n]

tree = rtree(taxa)
newick = writeTopology(tree)

rm = ResidualVarianceModel(p_res)
fm = IntegratedFactorModel(k, p_fac)

jm = JointProcessModel([rm, fm])

dm = DataModel(taxa,
               [randn(n, p_res), randn(n, p_fac)],
               ["trait.res", "trait.fac"])

bx = make_joint_xml(newick, dm, jm)
joint_path = "joint.xml"
# XMLConstructor.save_xml(joint_path, bx)
# @test isfila(joint_path)
# rm(joint_path)

################################################################################
## Integrate two xml
################################################################################

test_path = joinpath(@__DIR__, "data", "sequence.xml")

bx_seq = BEASTXMLElement(test_path)

taxa = XMLConstructor.find_element(bx_seq, XMLConstructor.EmptyDataXMLElement).taxa

k = 2
p = 10
n = length(taxa)
data = randn(n, p)
newick = writeTopology(rtree(taxa))

bx = XMLConstructor.make_pfa_xml(data, taxa, newick, k, useHMC = false)
XMLConstructor.save_xml("facGibbs.xml", bx)
@test isfile("facGibbs.xml")

XMLConstructor.merge_xml!(bx, bx_seq)
# xml = XMLConstructor.make_xml(bx)
XMLConstructor.save_xml("merge.xml", bx);
# print(xml)

################################################################################
## Sampled factor model with integrated factors sampler
################################################################################
k = 2
p = 10
n = length(taxa)
data = randn(n, p)
newick = writeTopology(rtree(taxa))

bx = XMLConstructor.make_sampled_pfa_xml(data, taxa, newick, k)
save_xml("facSampledHMC.xml", bx, change_filename=true)


text:
import StatsBase

"""
Create a nonlinear state space model.

#### Parameters
- f : Function. Process function, maps (x_{t}, t+1) -> x_{t+1}.
- b : Funciton. (Optional) input function, mapping (x_{t}, u_{t}, t) -> x_{t}.
- obs_loglik : Function. Should take a state, observation, and time, and
	return a log-likelihood value, i.e. log(P(y | x, t))

#### Returns
- NonlinearSSM object.
"""
type NonlinearSSM
	f::Function
	b::Function
	obs_loglik::Function
	nx::Int
	nu::Int
	ny::Int
end

function NonlinearSSM(f::Function, obs_loglik::Function,
		nx::Int, ny::Int, nu::Int)
	b(x::Vector, u::Vector, t::Int) = x
	NonlinearSSM(f, b, obs_loglik, nx, ny, ny)
end

"""
Type representing a particle filter.

#### Parameters
- nx : Int. Size of the state vector.
- nparticles : Int.  Number of particles in the ensemble.
- jitter : Function.  Returns a vector of nx random numbers.
"""
type ParticleFilter
	nx::Int
	nparticles::Int
	jitter::Function
end

function jitter(filter::ParticleFilter)
	noise = zeros(filter.nx, filter.nparticles)
	for i in 1:filter.nparticles
		noise[:, i] = filter.jitter()
	end
	return noise
end

function predict(m::NonlinearSSM, x::Vector; u::Vector=zeros(m.nu), t::Int=1)
	x = m.f(x, t)
	x = m.b(x, u, t)
	return x
end

function predict(m::NonlinearSSM, ensemble::Matrix; u::Vector=zeros(m.nu),
		t::Int=1)
	output = similar(ensemble)
	for i in 1:size(ensemble, 2)
		output[:, i] = m.b(m.f(ensemble[:, i], t), u, t)
	end
	return output
end

# do prediction in place, modifying ensemble
function predict!(m::NonlinearSSM, ensemble::Matrix; u::Vector=zeros(m.nu),
		t::Int=1)
	for i in 1:size(ensemble, 2)
		ensemble[:, i] = m.f(ensemble[:, i], u, t)
	end
end

function simulate(m::NonlinearSSM, nt::Int, x0::Vector)
	x = zeros(m.nx, nt)
	x[:, 1] = predict(m, x0)
	for i in 2:nt
		x[:, i] = predict(m, x[:, i-1])
	end
	return x
end

function observe(m::NonlinearSSM, x::Vector, y::Vector, t::Int=1)
	return m.obs_loglik(x, y, t)
end

function observe(m::NonlinearSSM, ensemble::Matrix, y::Vector, t::Int=1)
	return Float64[observe(m, ensemble[:, i], y, t) for i in 1:size(ensemble, 2)]
end

function update(m::NonlinearSSM, pred::Matrix, y::Vector, filter::ParticleFilter,
		t::Int=1)
	n = filter.nparticles
	likelihoods = exp(observe(m, pred, y, t))
	weights = StatsBase.WeightVec(likelihoods / sum(likelihoods))
	indices = sample(1:n, weights, n)
	return pred[:, indices] + jitter(filter)
end

# update the the ensemble in place
function update!(m::NonlinearSSM, pred::Matrix, y::Vector, filter::ParticleFilter,
		t::Int=1)
	n = size(pred, 2)
	likelihoods = exp(observe(m, pred, y, t))
	weights = StatsBase.WeightVec(likelihoods / sum(likelihoods))
	indices = sample(1:n, weights, n)
	pred = pred[:, indices] + jitter(filter)
end

function filter{T}(m::NonlinearSSM, y::Matrix{T}, x0::Matrix{T},
		f::ParticleFilter, u::Matrix{T}=zeros(m.nu, size(y, 2)))
	@assert size(x0, 1) == m.nx == f.nx
	@assert size(x0, 2) == f.nparticles
	@assert size(u, 1) == m.nu
	@assert size(y, 1) == m.ny
	@assert size(u, 2) == size(y, 2)

	nt = size(y, 2)
	fs = zeros(f.nx, f.nparticles, nt)
	x_pred = predict(m, x0, u=u[:, 1], t=1)
	fs[:, :, 1] = update(m, x0, y[:, 1], f, 1)
	for i in 2:nt
		x_pred = predict(m, fs[:, :, i-1], u=u[:, i], t=i)
		fs[:, :, i] = update(m, x_pred, y[:, i], f, i)
	end
	return fs
end


text:
mutable struct LinearSurrogate{X,Y,C,L,U} <: AbstractSurrogate
    x::X
    y::Y
    coeff::C
    lb::L
    ub::U
end

function LinearSurrogate(x,y,lb::Number,ub::Number)
    ols = lm(reshape(x,length(x),1),y)
    LinearSurrogate(x,y,coef(ols),lb,ub)
end

function add_point!(my_linear::LinearSurrogate,new_x,new_y)
    if length(my_linear.lb) == 1
        #1D
        my_linear.x = vcat(my_linear.x,new_x)
        my_linear.y = vcat(my_linear.y,new_y)
        md = lm(reshape(my_linear.x,length(my_linear.x),1),my_linear.y)
        my_linear.coeff = coef(md)
    else
        #ND
        n_previous = length(my_linear.x)
        a = vcat(my_linear.x,new_x)
        n_after = length(a)
        dim_new = n_after - n_previous
        n = length(my_linear.x)
        d = length(my_linear.x[1])
        tot_dim = n + dim_new
        X = Array{Float64,2}(undef,tot_dim,d)
        for j = 1:n
            X[j,:] = vec(collect(my_linear.x[j]))
        end
        if dim_new == 1
            X[n+1,:] = vec(collect(new_x))
        else
            i = 1
            for j = n+1:tot_dim
                X[j,:] = vec(collect(new_x[i]))
                i = i + 1
            end
        end
        my_linear.x = vcat(my_linear.x,new_x)
        my_linear.y = vcat(my_linear.y,new_y)
        md = lm(X,my_linear.y)
        my_linear.coeff = coef(md)
    end
    nothing
end

function (lin::LinearSurrogate)(val)
    return lin.coeff'*[val...]
end

"""
LinearSurrogate(x,y,lb,ub)

Builds a linear surrogate using GLM.jl

"""
function LinearSurrogate(x,y,lb,ub)
    #X = Array{eltype(x[1]),2}(undef,length(x),length(x[1]))
    #=
    X = Array{eltype(x[1]),2}(undef,length(x),length(x[1]))
    for j = 1:length(x)
        X[j,:] = vec(collect(x[j]))
    end
    =#
    T = collect(reshape(collect(Base.Iterators.flatten(x)), (length(x[1]),length(x)))')
    #T = transpose(reshape(reinterpret(eltype(x[1]), x), length(x[1]), length(x)))
    X = Array{eltype(x[1]),2}(undef,length(x),length(x[1]))
    X = copy(T)
    ols = lm(X,y)
    return LinearSurrogate(x,y,coef(ols),lb,ub)
end


text:
abstract type InferenceType end
abstract type Results{I<:InferenceType} end
abstract type SamplingResults{I} <: Results{I} end

struct BareResults{I} <: Results{I}
    interpretation :: I
    values :: DefaultDict{Any, Vector{Any}}
end
export BareResults

struct StructuredResults{I} <: SamplingResults{I}
    interpretation :: I
    values :: Dict{String, Vector{Any}}
    log_weights :: Vector{Float64}
end
Base.keys(sr::StructuredResults) = keys(sr.values)
Base.getindex(sr::StructuredResults, k) = Base.getindex(sr.values, k)
export StructuredResults

### base types ###
@doc raw"""
    struct NonparametricSamplingResults{I} <: SamplingResults{I}
        interpretation :: I
        log_weights :: Array{Float64, 1}
        return_values :: Array{Any, 1}
        traces :: Array{Trace, 1}
    end

Wrapper for results of sampling. Implements the following methods
from `Base`: `getindex`, `length`, `keys`.
Intepretation of log weights is dependent on `I`.
"""
struct NonparametricSamplingResults{I} <: SamplingResults{I}
    interpretation :: I
    log_weights :: Array{Float64, 1}
    return_values :: Array{Any, 1}
    traces :: Array{Trace, 1}
end
function Base.keys(r :: NonparametricSamplingResults{I}) where I <: InferenceType
    k = Set{Any}()
    for t in r.traces
        for a in keys(t)
            union!(k, (a,))
        end
    end
    k
end
function Base.getindex(r :: NonparametricSamplingResults{I}, k) where I <: InferenceType
    results = []
    for t in r.traces
        if k in keys(t)
            push!(results, t.trace[k].value)
        end
    end
    convert(Vector{typeof(results[1])}, results)  # assuming address type stability
end
Base.getindex(r :: SamplingResults{I}) where I <: InferenceType = r.return_values
Base.length(r :: SamplingResults{I}) where I <: InferenceType = length(r.log_weights)

@doc raw"""
    function addresses(r::NonparametricSamplingResults{I}) where I <: InferenceType

Get all addresses associated with the `SamplingResults` object,
``A = \bigcup_{t\in \text{traces}}\mathcal A_t``
"""
function addresses(r::NonparametricSamplingResults{I}) where I <: InferenceType
    a = Set()
    for t in r.traces
        union!(a, collect(keys(t)))
    end
    a
end
function get_first_node(r::NonparametricSamplingResults{I}, a) where I<:InferenceType
    for t in r.traces
        if a in keys(t.trace)
            return t[a]
        end
    end
end

StatsBase.mean(s :: NonparametricSamplingResults{I}, k, n) where I<:InferenceType = StatsBase.mean(sample(s, k, n))
StatsBase.mean(s :: NonparametricSamplingResults{I}, k) where I<:InferenceType = StatsBase.mean(s, k, length(s.log_weights))
StatsBase.std(s :: NonparametricSamplingResults{I}, k, n) where I<:InferenceType = StatsBase.std(sample(s, k, n))
StatsBase.std(s :: NonparametricSamplingResults{I}, k) where I<:InferenceType = StatsBase.std(s, k, length(s.log_weights))

@doc raw"""
    struct ParametricSamplingResults{I} <: SamplingResults{I}
        interpretation :: I
        log_weights :: Array{Float64, 1}
        return_values :: Array{Any, 1}
        traces :: Array{Trace, 1}
        distributions :: Dict
    end

`distributions` maps from addresses to distributions, ``a \mapsto \pi^{(a)}_\psi(z)``, where
``\pi^{(a)}_\psi(z)`` solves

```math
\max_\psi E_{z \sim p(z|x)}[\log \pi^{(a)}_\psi(z)].
```

The distributions are not used to generate values but only to score sampled values; values are
still sampled from the posterior traces.
Right now, the parametric approximation is very simple: values with support over the
negative orthant of  ``\mathbb R^D`` are approximated by (multivariate) normal distributions, while
values with support over only the positive orthant of ``\mathbb{R}^D`` are approximated by
(multivariate) lognormal distributions. This behavior is expected to change in the future.
"""
struct ParametricSamplingResults{I} <: SamplingResults{I}
    interpretation :: I
    log_weights :: Array{Float64, 1}
    return_values :: Array{Any, 1}
    traces :: Array{Trace, 1}
    distributions :: Dict
end

@doc raw"""
    function to_parametric(r::NonparametricSamplingResults{I}) where I<:InferenceType

Converts a nonparametric sampling results object into one that additionally contains
a mapping from addresses to distributions.
"""
function to_parametric(r::NonparametricSamplingResults{I}) where I<:InferenceType
    a_set = addresses(r)
    distributions = Dict()
    for a in a_set
        rep_node = get_first_node(r, a)
        if !rep_node.observed
            parametric_dist = parametric_posterior(rep_node.address, rep_node.dist, r)
            distributions[a] = parametric_dist
        end
    end
    ParametricSamplingResults{I}(r.interpretation, r.log_weights, r.return_values, r.traces, distributions)
end
Base.getindex(r :: ParametricSamplingResults{I}, k) where I <: InferenceType = r.distributions[k]
getsampled(r :: ParametricSamplingResults{I}, k) where I <: InferenceType = [t.trace[k].value for t in r.traces]

@doc raw"""
    function addresses(r::ParametricSamplingResults{I}) where I <: InferenceType

Get all addresses associated with the `SamplingResults` object,
``A = \bigcup_{t\in \text{traces}}\mathcal A_t``
"""
addresses(r::ParametricSamplingResults{I}) where I <: InferenceType = Set(collect(keys(r.distributions)))

@doc raw"""
    function Distributions.logpdf(r :: A, v) where A <: AbstractArray

Interprets an array of objects as a delta distribution over those objects. If `v` is in the support set,
returns ``-log |r|``. Otherwise, returns ``-\infty``.
"""
function Distributions.logpdf(r :: A, v) where A <: AbstractArray
    if v in r
        -log(length(r))
    else
        -Inf
    end
end

@doc raw"""
    function sample(t :: Trace, a, r :: NonparametricSamplingResults{I}; pa = ()) where I <: InferenceType

Treat a marginal site of a `SamplingResults` as a distribution, sampling from it into a trace.
"""
function sample(t :: Trace, a, r :: NonparametricSamplingResults{I}) where I <: InferenceType
    s = sample(r, a, 1)  # relies on inference-specific implementations
    n = node(s, a, r[a], false, EMPIRICAL)
    t[a] = n
    s
end

@doc raw"""
    function sample(t :: Trace, a, r::ParametricSamplingResults{I}; pa = ()) where I <: InferenceType

Treat a marginal site of a `SamplingResults` as a distribution, sampling from it into a trace.
"""
function sample(t::Trace, a, r::ParametricSamplingResults{I}) where I <:InferenceType
    s = sample(r, a, 1)  # relies on inference-specific implementations
    n = node(s, a, r.distributions[a], false, EMPIRICAL)
    t[a] = n
    s
end

function sample(t :: Trace, a, d, i :: Empirical; pa = ())
    connect_pa_ch!(t, pa, a)
    t[a].value
end

function parametric_posterior(address, dist, result :: SamplingResults)
    sd = size(dist)
    lsd = length(sd)
    lsd > 1 && error("parametric posterior supported for only rank <= 1 tensors")
    etype = eltype(dist)
    if lsd == 0 && etype === Float64
        univariate_continuous_parametric_posterior(address, dist, result)
    elseif lsd == 0 && etype === Int64
        univariate_discrete_parametric_posterior(address, dist, result)
    elseif lsd == 1 && etype === Float64
        multivariate_continuous_parametric_posterior(address, dist, result)
    else
        error("Parametric posterior supported for only univariate and continuous multivariate distributions")
    end
end

parametric_posterior(address, dist::D, result::SamplingResults) where D <: Categorical = fit_mle(Categorical, length(dist.p), results[address])
parametric_posterior(address, dist::D, results::SamplingResults) where D <: Beta = fit(Beta, results[address])
parametric_posterior(address, dist::D, results::SamplingResults) where D <: Dirichlet = fit_mle(Dirichlet, hcat(results[address]...); maxiter=5000, tol=1.0e-12)
parametric_posterior(address, dist::D, results::SamplingResults) where D <: Binomial = fit_mle(Binomial, dist.n, results[address])

function univariate_continuous_parametric_posterior(address, dist, result)
    if insupport(dist, -1.0)
        fit(Normal, result[address])
    else
        if !insupport(dist, 1.1)
            return fit(Beta, result[address])
        end
        fit(LogNormal, result[address])
    end
end
univariate_discrete_parametric_posterior(address, dist, result) = insupport(dist, -1) ? fit(DiscreteUniform, result[address]) : fit(Poisson, result[address])
function multivariate_continuous_parametric_posterior(address, dist, result)
    if insupport(dist, -1.0 .* ones(size(dist)))
        mat = hcat(result[address]...)
        m = vec(mean(mat, dims=2))
        s = vec(std(mat, dims=2))
        MvNormal(m, s)
    else
        mat = log.(hcat(result[address]...))
        m = vec(mean(mat, dims=2))
        s = vec(std(mat, dims=2))
        MvLogNormal(MvNormal(m, s))
    end
end

function logprob(r::Union{NonparametricSamplingResults{I},ParametricSamplingResults{I}}, address::T) where {I<:InferenceType, T}
    count = 0
    denom = 0
    lps = Float64[]
    for trace in r.traces
        for (a, node) in trace.trace
            if a == address
                denom += 1
                push!(lps, node.logprob_sum)
            end
        end
        count += 1
    end
    p_in_trace = denom / count
    logsumexp(lps) - log(denom) + log(p_in_trace)
end
function logprob(r::Union{NonparametricSamplingResults{I},ParametricSamplingResults{I}}, addresses::Vector{T}) where {I<:InferenceType, T}
    counts = Dict(a => 0 for a in addresses)
    any_in_tr = false
    denom = 0
    lps = Float64[]
    for trace in r.traces
        for (a, node) in trace.trace
            if a in addresses
                any_in_tr = true
                counts[a] += 1
                push!(lps, node.logprob_sum)
            end
        end
        if any_in_tr
            denom += 1
            any_in_tr = false
        end
    end
    p_in_trace = collect(values(counts)) ./ length(r.traces)
    logsumexp(lps) - log(denom) + sum(log.(p_in_trace))
end
export logprob

export NonparametricSamplingResults, ParametricSamplingResults, sample, aic
export to_parametric, getsampled, addresses, get_first_node