text:
""" Generalized algebraic theories (GATs) in Julia.

GATs generalize (multisorted) algebraic theories by incorporating a fragment of
dependent type theory. They allow type and term constructors to be partially
defined. GATs provide a convenient formal syntax for categorical structures.

At present, this module only supports defining the *signature* of a GAT. In the
future we may allow *axioms* to be expressed as well. Regardless, we will
persist in calling this module "GAT". Signatures are defined using the
`@signature` macro.

There are several useful things to do with a GAT signature:

1. Define *instances* of the signature (models of the theory) in Julia code
   using the `@instance` macro
2. Define *syntax* systems using the `@syntax` macro (see `Syntax` module)
3. Define models of the theory by finite *presentation* using a syntax system
   and the `@present` macro (see `Present` module)

References:

- (Cartmell, 1986, "Generalized algebraic theories and contextual categoies")
- (Pitts, 1995, "Categorical logic", Sec 6: "Dependent types")
"""
module GAT
export @signature, @instance, invoke_term

using Base.Meta: ParseError
using Compat
using AutoHashEquals
using DataStructures: OrderedDict
using Match

using ..Meta

# Data types
############

""" Base type for method stubs in GAT signature.
"""
abstract type Stub end

const Context = OrderedDict{Symbol,Expr0}

""" Type constructor in a GAT.
"""
@auto_hash_equals struct TypeConstructor
  name::Symbol
  params::Vector{Symbol}
  context::Context
  doc::Union{String,Nothing}

  function TypeConstructor(name::Symbol, params::Vector,
                           context::Context, doc=nothing)
    new(name, params, context, doc)
  end
end

""" Term constructor in a GAT.
"""
@auto_hash_equals struct TermConstructor
  name::Symbol
  params::Vector{Symbol}
  typ::Expr0
  context::Context
  doc::Union{String,Nothing}

  function TermConstructor(name::Symbol, params::Vector, typ::Expr0,
                           context::Context, doc=nothing)
    new(name, params, typ, context, doc)
  end
end

""" Signature for a generalized algebraic theory (GAT).
"""
@auto_hash_equals struct Signature
  types::Vector{TypeConstructor}
  terms::Vector{TermConstructor}
end

""" Typeclass = GAT signature + Julia-specific content.
"""
struct Typeclass
  name::Symbol
  type_params::Vector{Symbol}
  signature::Signature
  functions::Vector{JuliaFunction}
end

struct SignatureBinding
  name::Symbol
  params::Vector{Symbol}
end
struct SignatureHead
  main::SignatureBinding
  base::Vector{SignatureBinding}
  SignatureHead(main, base=[]) = new(main, base)
end

# Signatures
############

""" Define a signature of a generalized algebraic theory (GAT).

Three kinds of things can go in the signature body:

1. Type constructors, indicated by the special type `TYPE`, e.g.
   `Hom(X::Ob,Y::Ob)::TYPE`
2. Term constructors, e.g.,
   `id(X::Ob)::Hom(X,X)`
3. Julia functions operating on the term constructors to provide additional
   functionality

A signature can extend existing signatures (at present only one).
"""
macro signature(head, body)
  # Parse signature header.
  head = parse_signature_head(head)
  @assert length(head.base) <= 1 "Multiple signature extension not supported"
  if length(head.base) == 1
    base_name, base_params = head.base[1].name, head.base[1].params
    @assert all(p in head.main.params for p in base_params)
  else
    base_name, base_params = nothing, []
  end

  # Parse signature body: GAT types/terms and extra Julia functions.
  types, terms, functions = parse_signature_body(body)
  signature = Signature(types, terms)
  class = Typeclass(head.main.name, head.main.params, signature, functions)

  # We must generate and evaluate the code at *run time* because the base
  # signature, if specified, is not available at *parse time*.
  expr = :(signature_code($class, $(esc(base_name)), $base_params))
  Expr(:block,
    Expr(:call, esc(:eval), expr),
    :(Core.@__doc__ $(esc(head.main.name))))
end
function signature_code(main_class, base_mod, base_params)
  # Add types/terms/functions from base class, if provided.
  if isnothing(base_mod)
    class = main_class
  else
    base_class = base_mod.class()
    bindings = Dict(zip(base_class.type_params, base_params))
    main_sig = main_class.signature
    base_sig = replace_types(bindings, base_class.signature)
    sig = Signature([base_sig.types; main_sig.types],
                    [base_sig.terms; main_sig.terms])
    functions = [ [ replace_symbols(bindings, f) for f in base_class.functions ];
                  main_class.functions ]
    class = Typeclass(main_class.name, main_class.type_params, sig, functions)
  end
  signature = class.signature

  # Generate module with stub types.
  mod = Expr(:module, true, class.name,
    Expr(:block, [
      # Prevents error about export not being at toplevel.
      # https://github.com/JuliaLang/julia/issues/28991
      LineNumberNode(0);
      Expr(:export, [cons.name for cons in signature.types]...);
      map(gen_abstract_type, signature.types);
      :(class() = $(class));
    ]...))

  # Generate method stubs.
  # (We put them outside the module, so the stub type names must be qualified.)
  bindings = Dict(cons.name => Expr(:(.), class.name, QuoteNode(cons.name))
                  for cons in signature.types)
  fns = interface(class)
  toplevel = [ generate_function(replace_symbols(bindings, f)) for f in fns ]

  # Modules must be at top level:
  # https://github.com/JuliaLang/julia/issues/21009
  Expr(:toplevel, mod, toplevel...)
end

function parse_signature_head(expr::Expr)::SignatureHead
  parse = parse_signature_binding
  @match expr begin
    (Expr(:call, [:(=>), Expr(:tuple, bases), main])
      => SignatureHead(parse(main), map(parse, bases)))
    Expr(:call, [:(=>), base, main]) => SignatureHead(parse(main), [parse(base)])
    _ => SignatureHead(parse(expr))
  end
end

function parse_signature_binding(expr::Expr)::SignatureBinding
  @match expr begin
    Expr(:call, [name::Symbol, params...]) => SignatureBinding(name, params)
    _ => throw(ParseError("Ill-formed signature binding $expr"))
  end
end

""" Parse the body of a GAT signature declaration.
"""
function parse_signature_body(expr::Expr)
  @assert expr.head == :block
  types = OrderedDict{Symbol,TypeConstructor}()
  terms = TermConstructor[]
  funs = JuliaFunction[]
  for elem in strip_lines(expr).args
    head = last(parse_docstring(elem)).head
    if head in (:(::), :call)
      cons = parse_constructor(elem)
      if isa(cons, TypeConstructor)
        if haskey(types, cons.name)
          throw(ParseError("Duplicate type constructor $elem"))
        else
          types[cons.name] = cons
        end
      else
        push!(terms, cons)
      end
    elseif head in (:(=), :function)
      push!(funs, parse_function(elem))
    else
      throw(ParseError("Ill-formed signature element $elem"))
    end
  end
  return (collect(values(types)), terms, funs)
end

""" Julia functions for type parameter accessors.
"""
function accessors(sig::Signature)::Vector{JuliaFunction}
  vcat(map(accessors, sig.types)...)
end
function accessors(cons::TypeConstructor)::Vector{JuliaFunction}
  [ JuliaFunction(Expr(:call, param, Expr(:(::), cons.name)),
                  strip_type(cons.context[param]))
    for param in cons.params ]
end

""" Julia functions for term constructors of GAT.
"""
function constructors(sig::Signature)::Vector{JuliaFunction}
  [ constructor(cons, sig) for cons in sig.terms ]
end
function constructor(cons::TermConstructor, sig::Signature)::JuliaFunction
  arg_names = cons.params
  arg_types = [ strip_type(cons.context[name]) for name in arg_names ]
  args = [ Expr(:(::), name, typ) for (name,typ) in zip(arg_names, arg_types) ]
  return_type = strip_type(cons.typ)

  call_expr = Expr(:call, cons.name, args...)
  if !any(has_type(sig, typ) for typ in arg_types)
    call_expr = add_type_dispatch(call_expr, return_type)
  end
  JuliaFunction(call_expr, return_type)
end

""" Complete set of Julia functions for a type class.
"""
function interface(class::Typeclass)::Vector{JuliaFunction}
  [ accessors(class.signature);
    constructors(class.signature);
    class.functions ]
end

""" Get type constructor by name.

Unlike term constructors, type constructors cannot be overloaded, so there is at
most one type constructor with a given name.
"""
function get_type(sig::Signature, name::Symbol)::TypeConstructor
  indices = findall(cons -> cons.name == name, sig.types)
  @assert length(indices) == 1
  sig.types[indices[1]]
end
function has_type(sig::Signature, name::Symbol)::Bool
  findfirst(cons -> cons.name == name, sig.types) != nothing
end

""" Add a type-valued first argument to a Julia function signature.

We need this to avoid ambiguity in method dispatch when a term constructor has
no arguments (e.g., `munit()`) or more generally has no arguments that are types
in the signature (e.g., object generators in a category).

The fundamental reason why these cases must be treated specially is that Julia
does not (yet?) support
[dispatching on return type](https://github.com/JuliaLang/julia/issues/19206).
"""
function add_type_dispatch(call_expr::Expr, type_expr::Expr0)::Expr
  @match call_expr begin
    (Expr(:call, [name, args...]) =>
      Expr(:call, name, :(::Type{$type_expr}), args...))
    _ => throw(ParseError("Ill-formed call expression $call_expr"))
  end
end

# GAT expressions
#################

""" Parse a raw expression in a GAT.

A "raw expression" is a just composition of function and constant symbols.
"""
function parse_raw_expr(expr)
  @match expr begin
    Expr(:call, args) => map(parse_raw_expr, args)
    head::Symbol => nothing
    _ => throw(ParseError("Ill-formed raw expression $expr"))
  end
  expr # Return the expression unmodified. This function just checks syntax.
end

""" Parse context for term or type in a GAT.
"""
function parse_context(expr::Expr)::Context
  context = Context()
  args = expr.head == :tuple ? expr.args : [ expr ]
  for arg in args
    name, typ = @match arg begin
      Expr(:(::), [name::Symbol, typ]) => (name, parse_raw_expr(typ))
      _ => throw(ParseError("Ill-formed context expression $expr"))
    end
    push_context!(context, name, typ)
  end
  return context
end
function push_context!(context, name, expr)
  if haskey(context, name)
    throw(ParseError("Name $name already defined"))
  end
  context[name] = expr
end

""" Parse type or term constructor in a GAT.
"""
function parse_constructor(expr::Expr)::Union{TypeConstructor,TermConstructor}
  # Context is optional.
  doc, expr = parse_docstring(expr)
  cons_expr, context = @match expr begin
    Expr(:call, [:<=, inner, context]) => (inner, parse_context(context))
    _ => (expr, Context())
  end

  # Allow abbreviated syntax where tail of context is included in parameters.
  function parse_param(param::Expr0)::Symbol
    @match param begin
      Expr(:(::), [name::Symbol, typ]) => begin
        push_context!(context, name, parse_raw_expr(typ))
        name
      end
      name::Symbol => name
      _ => throw(ParseError("Ill-formed type/term parameter $param"))
    end
  end

  @match cons_expr begin
    (Expr(:(::), [name::Symbol, :TYPE])
      => TypeConstructor(name, [], context, doc))
    (Expr(:(::), [Expr(:call, [name::Symbol, params...]), :TYPE])
      => TypeConstructor(name, map(parse_param, params), context, doc))
    (Expr(:(::), [Expr(:call, [name::Symbol, params...]), typ])
      => TermConstructor(name, map(parse_param, params), parse_raw_expr(typ),
                         context, doc))
    _ => throw(ParseError("Ill-formed type/term constructor $cons_expr"))
  end
end

""" Generate abstract type definition from a GAT type constructor.
"""
function gen_abstract_type(cons::TypeConstructor)::Expr
  stub_name = GlobalRef(GAT, :Stub)
  expr = :(abstract type $(cons.name) <: $stub_name end)
  generate_docstring(expr, cons.doc)
end

""" Replace names of type constructors in a GAT.
"""
function replace_types(bindings::Dict, sig::Signature)::Signature
  Signature([ replace_types(bindings, t) for t in sig.types ],
            [ replace_types(bindings, t) for t in sig.terms ])
end
function replace_types(bindings::Dict, cons::TypeConstructor)::TypeConstructor
  TypeConstructor(replace_symbols(bindings, cons.name), cons.params,
                  replace_types(bindings, cons.context), cons.doc)
end
function replace_types(bindings::Dict, cons::TermConstructor)::TermConstructor
  TermConstructor(cons.name, cons.params,
                  replace_symbols(bindings, cons.typ),
                  replace_types(bindings, cons.context), cons.doc)
end
function replace_types(bindings::Dict, context::Context)::Context
  GAT.Context(((name => @match expr begin
    (Expr(:call, [sym::Symbol, args...]) =>
      Expr(:call, replace_symbols(bindings, sym), args...))
    sym::Symbol => replace_symbols(bindings, sym)
  end) for (name, expr) in context))
end

""" Remove type parameters from dependent type.
"""
function strip_type(expr)::Symbol
  @match expr begin
    Expr(:call, [head::Symbol, args...]) => head
    sym::Symbol => sym
  end
end

# GAT expressions in a signature
################################

""" Expand context variables that occur implicitly in an expression.

Reference: (Cartmell, 1986, Sec 10: 'Informal syntax').
"""
function expand_in_context(expr, params::Vector{Symbol},
                           context::Context, sig::Signature)
  @match expr begin
    Expr(:call, [name::Symbol, args...]) => begin
      expanded = [expand_in_context(e, params, context, sig) for e in args]
      Expr(:call, name, expanded...)
    end
    name::Symbol => begin
      if name in params
        name
      elseif haskey(context, name)
        expand_symbol_in_context(name, params, context, sig)
      else
        error("Name $name missing from context $context")
      end
    end
    _ => throw(ParseError("Ill-formed raw expression $expr"))
  end
end
function expand_symbol_in_context(sym::Symbol, params::Vector{Symbol},
                                  context::Context, sig::Signature)
  # This code expands symbols that occur as direct arguments to type
  # constructors. If there are term constructors in between, it does not work:
  # indeed, it cannot work in general because the term constructors are not
  # necessarily injective. For example, we can expand :X in
  #   (:X => :Ob, :f => :(Hom(X)))
  # but not in
  #   (:X => :Ob, :Y => :Ob, :f => :(Hom(otimes(X,Y))))
  names = collect(keys(context))
  start = findfirst(names .== sym)
  for name in names[start+1:end]
    expr = context[name]
    if isa(expr, Expr) && expr.head == :call && sym in expr.args[2:end]
      cons = get_type(sig, expr.args[1])
      accessor = cons.params[findfirst(expr.args[2:end] .== sym)]
      expanded = Expr(:call, accessor, name)
      return expand_in_context(expanded, params, context, sig)
    end
  end
  error("Name $sym does not occur explicitly among $params in context $context")
end

""" Expand context variables that occur implicitly in the type expression
of a term constructor.
"""
function expand_term_type(cons::TermConstructor, sig::Signature)
  isa(cons.typ, Symbol) ? cons.typ :
    expand_in_context(cons.typ, cons.params, cons.context, sig)
end

""" Implicit equations defined by a context.

This function allows a generalized algebraic theory (GAT) to be expressed as
an essentially algebraic theory, i.e., as partial functions whose domains are
defined by equations.

References:
 - (Cartmell, 1986, Sec 6: "Essentially algebraic theories and categories with
    finite limits")
 - (Freyd, 1972, "Aspects of topoi")
"""
function equations(context::Context, sig::Signature)::Vector{Pair}
  # The same restrictions as `expand_symbol_in_context` apply here.
  eqs = Pair[]
  names = collect(keys(context))
  for (start, var) in enumerate(names)
    for name in names[start+1:end]
      expr = context[name]
      if isa(expr, Symbol) && !has_type(sig, expr)
        # If the constructor is a symbol and there isn't a matching type in
        # the signature, assume it's a Julia type. For now, these are
        # completely ignored by the syntax system.
        continue
      end
      expr = isa(expr, Symbol) ? Expr(:call, expr) : expr
      cons = get_type(sig, expr.args[1])
      accessors = cons.params[findall(expr.args[2:end] .== var)]
      append!(eqs, (Expr(:call, a, name) => var for a in accessors))
    end
  end
  eqs
end

""" Implicit equations defined by context, allowing for implicit variables.
"""
function equations(params::Vector{Symbol}, context::Context,
                   sig::Signature)::Vector{Pair}
  eqs = [ (expand_in_context(lhs, params, context, sig) =>
           expand_in_context(rhs, params, context, sig))
          for (lhs, rhs) in equations(context, sig) ]
  # Remove tautologies (expr == expr) resulting from expansions.
  # FIXME: Should we worry about redundancies from the symmetry of equality,
  # i.e., (expr1 == expr2) && (expr2 == expr1)?
  filter(eq -> eq.first != eq.second, eqs)
end

""" Implicit equations for term constructor.
"""
function equations(cons::TermConstructor, sig::Signature)::Vector{Pair}
  equations(cons.params, cons.context, sig)
end

# Instances
###########

""" Define an *instance* of a generalized algebraic theory (GAT).

These are perfectly analagous to instances of a type class in Haskell. See also
the Typeclass.jl library for Julia.
"""
macro instance(head, body)
  # Parse the instance definition.
  head = parse_signature_binding(head)
  functions = parse_instance_body(body)

  # We must generate and evaluate the code at *run time* because the signature
  # module is not defined at *parse time*.
  # Also, we "throw away" any docstring.
  # FIXME: Is there a better place to put the docstring?
  expr = :(instance_code($(esc(head.name)), $(esc(head.params)), $functions))
  Expr(:block,
    Expr(:call, esc(:eval), expr),
    :(Core.@__doc__ abstract type $(esc(gensym(:instance_doc))) end)) # /dev/null
end
function instance_code(mod, instance_types, instance_fns)
  code = Expr(:block)
  class = mod.class()
  bindings = Dict(zip(class.type_params, instance_types))
  bound_fns = [ replace_symbols(bindings, f) for f in interface(class) ]
  bound_fns = OrderedDict(parse_function_sig(f) => f for f in bound_fns)
  instance_fns = Dict(parse_function_sig(f) => f for f in instance_fns)
  for (sig, f) in bound_fns
    if haskey(instance_fns, sig)
      f_impl = instance_fns[sig]
    elseif !isnothing(f.impl)
      f_impl = f
    else
      error("Method $(f.call_expr) not implemented in $(class.name) instance")
    end
    push!(code.args, generate_function(f_impl))
  end
  return code
end

""" Parse the body of a GAT instance definition.
"""
function parse_instance_body(expr::Expr)::Vector{JuliaFunction}
  @match strip_lines(expr) begin
    Expr(:block, args) => map(parse_function, args)
    _ => throw(ParseEror("Ill-formed instance definition"))
  end
end

""" Invoke a term constructor by name on an instance.

This method provides reflection for GAT signatures. In everyday use the generic
method for the constructor should be called directly, not through this function.

Cf. Julia's builtin `invoke()` function.
"""
function invoke_term(signature_module::Module, instance_types::Tuple,
                     constructor_name::Symbol, args...)
  # Get the corresponding Julia method from the parent module.
  method = getfield(parentmodule(signature_module), constructor_name)
  args = collect(Any, args)

  # Add dispatch on return type, if necessary.
  if !any(typeof(arg) <: typ for typ in instance_types for arg in args)
    # Case 1: Name refers to type constructor, e.g., generator constructor
    # in syntax system.
    signature = signature_module.class().signature
    index = findfirst(cons -> cons.name == constructor_name, signature.types)
    if isnothing(index)
      # Case 2: Name refers to term constructor.
      # FIXME: Terms constructors can be overloaded, so there may be multiple
      # term constructors with the same name. Distinguishing them requires type
      # inference. I am punting on that right now.
      constructor = signature.terms[
        findfirst(cons -> cons.name == constructor_name, signature.terms)]
      return_name = strip_type(constructor.typ)
      index = findfirst(cons -> cons.name == return_name, signature.types)
    end
    insert!(args, 1, instance_types[index])
  end

  # Invoke the method!
  method(args...)
end

end


text:
# Autogenerated wrapper script for mimalloc_jll for x86_64-linux-gnu
export libmimalloc

JLLWrappers.@generate_wrapper_header("mimalloc")
JLLWrappers.@declare_library_product(libmimalloc, "libmimalloc.so.2.0")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libmimalloc,
        "lib/libmimalloc.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()

ext: jl


text:
# call with julia --track-allocation=user test/runprofiler.jl

using VortexCollisions

include("test_functions.jl")

function run_profiler()
    M = 64
    N = 64

    ℳcut = M
    𝒩cut = N

    grid = Grid2d(M,N)
    ft   = FourierTransform(grid; ℳcut=ℳcut, 𝒩cut=𝒩cut)
    op   = FokkerPlanckOperator(grid, ft)

    u    = get_field(grid)
    divJ = zeros(u)

    evaluate_function_on_grid(grid, u_test, u)

    collision_operator!(op, u, divJ)

    Profile.clear_malloc_data()

    collision_operator!(op, u, divJ)

end

run_profiler()


text:
function init_face_data(elem::Tri, N; quad_rule_face = gauss_quad(0,0,N))
    # Nodes on faces, and face node coordinate
    r1D, w1D = quad_rule_face
    e = ones(size(r1D)) # vector of all ones
    z = zeros(size(r1D)) # vector of all zeros
    rf,sf = map_face_nodes(elem,r1D)
    wf = vec(repeat(w1D,3,1));
    nrJ = [z; e; -e]
    nsJ = [-e; e; z]
    return rf,sf,wf,nrJ,nsJ
end

function init_face_data(elem::Quad, N; quad_rule_face = gauss_quad(0,0,N))
    r1D,w1D = quad_rule_face
    e = ones(size(r1D))
    z = zeros(size(r1D))
    rf,sf = map_face_nodes(elem,r1D)
    wf = vec(repeat(w1D,4,1)); # 4 faces
    nrJ = [z; e; z; -e]
    nsJ = [-e; z; e; z]
    return rf,sf,wf,nrJ,nsJ
end

function init_face_data(elem::Hex, N; quad_rule_face = quad_nodes(Quad(), N))
    rquad, squad, wquad = quad_rule_face
    e = ones(size(rquad))
    zz = zeros(size(rquad))
    rf,sf,tf = map_face_nodes(elem,rquad,squad)
    Nfaces = 6
    wf = vec(repeat(wquad,Nfaces,1));
    nrJ = [-e; e; zz;zz; zz;zz]
    nsJ = [zz;zz; -e; e; zz;zz]
    ntJ = [zz;zz; zz;zz; -e; e]
    return rf,sf,tf,wf,nrJ,nsJ,ntJ
end

function init_face_data(elem::Tet, N; quad_rule_face = quad_nodes(Tri(), N))
    rquad,squad,wquad = quad_rule_face
    e = ones(size(rquad))
    zz = zeros(size(rquad))
    rf,sf,tf = map_face_nodes(elem,rquad,squad)
    Nfaces = 4
    wf = vec(repeat(wquad,Nfaces,1));
    nrJ = [zz; e; -e; zz]
    nsJ = [-e; e; zz; zz]
    ntJ = [zz; e; zz; -e]
    return rf,sf,tf,wf,nrJ,nsJ,ntJ
end

"""
    RefElemData(elem::Line, N;
                quad_rule_vol = quad_nodes(elem, N+1))
    RefElemData(elem::Union{Tri, Quad}, N;
                 quad_rule_vol = quad_nodes(elem, N),
                 quad_rule_face = gauss_quad(0, 0, N))
    RefElemData(elem::Union{Hex, Tet}, N;
                 quad_rule_vol = quad_nodes(elem, N),
                 quad_rule_face = quad_nodes(Quad(), N))
    RefElemData(elem; N, kwargs...) # version with keyword args

Constructor for `RefElemData` for different element types.
"""
function RefElemData(elem::Line, approxType::Polynomial, N;
                     quad_rule_vol = quad_nodes(elem, N+1),
                     Nplot = 10)

    fv = face_vertices(elem)

    # Construct matrices on reference elements
    r = nodes(elem,N)
    Fmask = [1 N+1]
    VDM = vandermonde(elem, N, r)
    Dr = grad_vandermonde(elem, N, r)/VDM

    V1 = vandermonde(elem,1,r)/vandermonde(elem,1,[-1;1])

    rq,wq = quad_rule_vol
    Vq = vandermonde(elem,N,rq)/VDM
    M = Vq'*diagm(wq)*Vq
    Pq = M\(Vq'*diagm(wq))

    rf = [-1.0;1.0]
    nrJ = [-1.0;1.0]
    wf = [1.0;1.0]
    Vf = vandermonde(elem,N,rf)/VDM
    LIFT = M\(Vf') # lift matrix

    # plotting nodes
    rp = equi_nodes(elem,Nplot)
    Vp = vandermonde(elem,N,rp)/VDM

    return RefElemData(elem,approxType,N,fv,V1,
                       tuple(r),VDM,vec(Fmask),
                       Nplot,tuple(rp),Vp,
                       tuple(rq),wq,Vq,
                       tuple(rf),wf,Vf,tuple(nrJ),
                       M,Pq,tuple(Dr),LIFT)
end

function RefElemData(elem::Union{Tri, Quad},  approxType::Polynomial, N;
                     quad_rule_vol = quad_nodes(elem, N),
                     quad_rule_face = quad_nodes(face_type(elem), N),
                     Nplot = 10)

    fv = face_vertices(elem) # set faces for triangle

    # Construct matrices on reference elements
    r,s = nodes(elem,N)
    Fmask = hcat(find_face_nodes(elem,r,s)...)

    VDM,Vr,Vs = basis(elem,N,r,s)
    Dr = Vr/VDM
    Ds = Vs/VDM

    # low order interpolation nodes
    r1,s1 = nodes(elem,1)
    V1 = vandermonde(elem,1,r,s)/vandermonde(elem,1,r1,s1)

    rf,sf,wf,nrJ,nsJ = init_face_data(elem, N, quad_rule_face = quad_rule_face)

    rq,sq,wq = quad_rule_vol
    Vq = vandermonde(elem,N,rq,sq)/VDM
    M = Vq'*diagm(wq)*Vq
    Pq = M\(Vq'*diagm(wq))

    Vf = vandermonde(elem,N,rf,sf)/VDM # interpolates from nodes to face nodes
    LIFT = M\(Vf'*diagm(wf)) # lift matrix used in rhs evaluation

    # plotting nodes
    rp, sp = equi_nodes(elem,Nplot)
    Vp = vandermonde(elem,N,rp,sp)/VDM

    Drs = (Dr,Ds)

    return RefElemData(elem, approxType, N, fv, V1,
                       tuple(r, s), VDM, vec(Fmask),
                       Nplot, tuple(rp, sp), Vp,
                       tuple(rq, sq), wq, Vq,
                       tuple(rf, sf), wf, Vf, tuple(nrJ, nsJ),
                       M, Pq, Drs, LIFT)
end

function RefElemData(elem::Tet, approxType::Polynomial, N;
                     quad_rule_vol = quad_nodes(elem, N),
                     quad_rule_face = quad_nodes(face_type(elem), N),
                     Nplot = 10)

    fv = face_vertices(elem)

    # Construct matrices on reference elements
    r,s,t = nodes(elem,N)
    Fmask = hcat(find_face_nodes(elem,r,s,t)...)
    VDM,Vr,Vs,Vt = basis(elem,N,r,s,t)
    Dr,Ds,Dt = (A->A/VDM).((Vr,Vs,Vt))

    # low order interpolation nodes
    r1,s1,t1 = nodes(elem,1)
    V1 = vandermonde(elem,1,r,s,t)/vandermonde(elem,1,r1,s1,t1)

    # Nodes on faces, and face node coordinate
    rf,sf,tf,wf,nrJ,nsJ,ntJ = init_face_data(elem, N, quad_rule_face = quad_rule_face)

    # quadrature nodes - build from 1D nodes.
    rq,sq,tq,wq = quad_rule_vol
    Vq = vandermonde(elem,N,rq,sq,tq)/VDM
    M = Vq'*diagm(wq)*Vq
    Pq = M\(Vq'*diagm(wq))

    Vf = vandermonde(elem,N,rf,sf,tf)/VDM
    LIFT = M\(Vf'*diagm(wf))

    # plotting nodes
    rp,sp,tp = equi_nodes(elem,Nplot)
    Vp = vandermonde(elem,N,rp,sp,tp)/VDM

    Drst = (Dr,Ds,Dt)

    return RefElemData(elem, approxType, N, fv, V1,
                       tuple(r, s, t), VDM, vec(Fmask),
                       Nplot, tuple(rp, sp, tp), Vp,
                       tuple(rq, sq, tq), wq, Vq,
                       tuple(rf, sf, tf), wf, Vf, tuple(nrJ, nsJ, ntJ),
                       M, Pq, Drst, LIFT)
end

# specialize constructor for `Hex` to allow for higher polynomial degrees `N`
function RefElemData(elem::Hex, approxType::Polynomial, N;
                     quad_rule_vol = quad_nodes(elem, N),
                     quad_rule_face = quad_nodes(face_type(elem), N),
                     Nplot = 10)

    fv = face_vertices(elem)

    # Construct matrices on reference elements
    r, s, t = nodes(elem, N)
    Fmask = hcat(find_face_nodes(elem, r, s, t)...)

    # construct 1D operator for faster Kronecker solves
    r1D = nodes(Line(), N)
    rq1D, wq1D = quad_nodes(Line(), N)
    VDM_1D = vandermonde(Line(), N, r1D)
    Vq1D = vandermonde(Line(), N, rq1D) / VDM_1D
    invVDM_1D = inv(VDM_1D)
    invM_1D = VDM_1D * VDM_1D'
    M1D = Vq1D' * diagm(wq1D) * Vq1D

    # form kronecker products of multidimensional matrices to invert/multiply
    M = kronecker(M1D, M1D, M1D)
    VDM = kronecker(VDM_1D, VDM_1D, VDM_1D)
    invVDM = kronecker(invVDM_1D, invVDM_1D, invVDM_1D)
    invM = kronecker(invM_1D, invM_1D, invM_1D)

    _, Vr, Vs, Vt = basis(elem, N, r, s, t)
    Dr, Ds, Dt = (A -> A * invVDM).((Vr, Vs, Vt))

    # low order interpolation nodes
    r1, s1, t1 = nodes(elem, 1)
    V1 = vandermonde(elem, 1, r, s, t) / vandermonde(elem, 1, r1, s1, t1)

    # Nodes on faces, and face node coordinate
    rf, sf, tf, wf, nrJ, nsJ, ntJ = init_face_data(elem, N, quad_rule_face=quad_rule_face)

    # quadrature nodes - build from 1D nodes.
    rq, sq, tq, wq = quad_rule_vol
    Vq = vandermonde(elem, N, rq, sq, tq) * invVDM
    # M = Vq' * diagm(wq) * Vq
    Pq = invM * (Vq' * diagm(wq))

    Vf = vandermonde(elem, N, rf, sf, tf) * invVDM
    LIFT = invM * (Vf' * diagm(wf))

    # plotting nodes
    rp, sp, tp = equi_nodes(elem, Nplot)
    # Vp = vandermonde(elem, N, rp, sp, tp) * invVDM
    rp1D = LinRange(-1, 1, Nplot + 1)
    Vp1D = vandermonde(Line(), N, rp1D) / VDM_1D
    Vp = kronecker(Vp1D, Vp1D, Vp1D)

    Drst = (Dr, Ds, Dt)

    return RefElemData(elem, approxType, N, fv, V1,
                       tuple(r, s, t), VDM, vec(Fmask),
                       Nplot, tuple(rp, sp, tp), Vp,
                       tuple(rq, sq, tq), wq, Vq,
                       tuple(rf, sf, tf), wf, Vf, tuple(nrJ, nsJ, ntJ),
                       M, Pq, Drst, LIFT)
end


text:
module GitHub

import Base.show

using Requests
import JSON
using HttpCommon

abstract GitHubType
abstract Owner <: GitHubType

# types
export User,
       Organization,
       Repo,
       Issue,
       HttpError,
       AuthError,
       StatsError

# methods
export authenticate,
       set_api_endpoint,
       set_web_endpoint,
       user,
       star,
       unstar,
       stargazers,
       starred,
       forks,
       fork,
       contributors,
       commit_activity,
       code_frequency,
       participation,
       punch_card,
       collaborators,
       iscollaborator,
       add_collaborator,
       remove_collaborator,
       watchers,
       watched,
       watching,
       watch,
       unwatch,
       followers,
       following,
       org,
       orgs,
       repo,
       issue,
       create_issue,
       edit_issue,
       issues

include("utils.jl")
include("endpoint.jl")
include("error.jl")
include("auth.jl")
include("users.jl")
include("organizations.jl")
include("repos.jl")
include("issues.jl")
include("starring.jl")
include("forks.jl")
include("statistics.jl")
include("collaborators.jl")
include("watching.jl")

end


text:
using Documenter
using SymmetricSparseMatrices

makedocs(
    sitename = "SymmetricSparseMatrices",
    format = Documenter.HTML(),
    modules = [SymmetricSparseMatrices]
)

deploydocs(repo = "github.com/victorsndvg/SymmetricSparseMatrices.jl.git")


text:
for abstract type AbstractVector{T}
where T
  length(using) -> Int
  index(using, Int) -> T
end


text:
using Flux, Gym, Printf, Zygote
using Zygote: @adjoint
using Flux.Optimise: update!
using Statistics: mean
#using CuArrays

import Base.sign

#Load game environment

env = make("CartPole-v0")
reset!(env)

#ctx = Ctx(env)

#display(ctx.s)
#using Blink# when not on Juno
#body!(Blink.Window(), ctx.s)

# ----------------------------- Parameters -------------------------------------

STATE_SIZE = length(env._env.state)
ACTION_SIZE = length(env._env.action_space)
MAX_TRAIN_REWARD = env._env.x_threshold * env._env.θ_threshold_radians
SEQ_LEN = 8

# Optimiser params
η = 3f-2
# ------------------------------ Model Architecture ----------------------------
sign(x) = Base.sign.(x)
@adjoint sign(x) = sign(x), x̄ -> (x̄,)

model = Chain(Dense(STATE_SIZE, 24, relu),
              Dense(24, 48, relu),
              Dense(48, 1, tanh), x->sign(x)) |> gpu

opt = ADAM(η)

action(state) = state |> model |> (model_output) -> (3 .+ model_output) / 2

loss(rewards) = Flux.mse(rewards, MAX_TRAIN_REWARD)

# ----------------------------- Helper Functions -------------------------------

function train_reward(env::EnvWrapper)
    s = env._env.state
    x, ẋ, θ, θ̇  = s
    # Custom reward for training
    # Product of Triangular function over x-axis and θ-axis
    # Min reward = 0, Max reward = env.x_threshold * env.θ_threshold_radians
    x_upper = env._env.x_threshold - x
    x_lower = env._env.x_threshold + x

    r_x     = max(0f0, min(x_upper, x_lower))

    θ_upper = env._env.θ_threshold_radians - θ
    θ_lower = env._env.θ_threshold_radians + θ

    r_θ     = max(0f0, min(θ_upper, θ_lower))

    return r_x * r_θ
end

function μEpisode(env::EnvWrapper)
    l = 0
    for frames ∈ 1:SEQ_LEN
        #render(env, ctx)
        #sleep(0.01)
        a = action(env._env.state)
        s′, r, done, _ = step!(env, a)

        if trainable(env)
            l += loss(train_reward(env))
        end
        game_over(env) && break
    end
    return l
end

function episode!(env::EnvWrapper)
    reset!(env)
    while !game_over(env)
        if trainable(env)
            grads = gradient(()->μEpisode(env), params(model))
            update!(opt, params(model), grads)
        else
            μEpisode(env)
        end
    end
    env.total_reward
end

# -------------------------------- Testing -------------------------------------

function test(env::EnvWrapper)
    score_mean = 0f0
    testmode!(env)
    for _=1:100
        total_reward = episode!(env)
        score_mean += total_reward / 100
    end
    testmode!(env, false)
    return score_mean
end

# ------------------------------ Training --------------------------------------

e = 1
while true
    global e
    total_reward = @sprintf "%6.2f" episode!(env)
    print("Episode: $e | Score: $total_reward | ")

    score_mean = test(env)
    score_mean_str = @sprintf "%6.2f" score_mean
    print("Mean score over 100 test episodes: " * score_mean_str)

    println()

    if score_mean > env.reward_threshold
        println("CartPole-v0 solved!")
    break
    end
    e += 1
end

max_stars_repo_name: jstrube/model-zoo

max_issues_repo_name: jstrube/model-zoo

max_forks_repo_name: jstrube/model-zoo


text:
import Base.show

function show(io::IO, x::Select)
    write(io, "SELECT ",
        join(UTF8String[sprint(io->show(io, x)) for x in x.relation.attrs], ", "))
    show(io, x.from)
end

show(io::IO, x::DataStreamRelation) = fieldset(x).name

function show(io::IO, x::Attr)
    write(io, string(name(x)))
end

function show(io::IO, x::QualifiedAttr)
    write(io, string(x.relation.name), ".", string(name(x.attr)))
end

function show(io::IO, x::Where)
    show(io, x.source)
    write(io, " WHERE ")
    show(io, x.predicate)
end

function show(io::IO, x::GroupBy)
    show(io, x.source)
    write(io, " GRPOUPBY ")
    show(io, x.grouping)
end

function show(io::IO, x::Aggregate)
    show(io, x.source)
    write(io, ' ')
    show(io, x.aggregator)
end

function show(io::IO, agg::Count)
    write(io, "count()")
end

function show(io::IO, agg::Avg)
    write(io, "avg(")
    show(io, agg.attr)
    write(io, ")")
end

function show(io::IO, agg::Sum)
    write(io, "sum(")
    show(io, agg.attr)
    write(io, ")")
end

function show(io::IO, a::QualifiedAttr)
    write(io, string(a.relation.name))
    write(io, '.')
    write(io, string(name(a)))
end

function show{op}(io::IO, a::BinOp{op})
    write(io, '(')
    show(io, a.l)
    write(io, ' ')
    write(io, op)
    write(io, ' ')
    show(io, a.r)
    write(io, ')')
end


text:
function bisection(f::Function, a::Number, b::Number;
    tol::AbstractFloat=1e-5, maxiter::Integer=100)
    fa = f(a)
    fa*f(b) <= 0 || error("No real root in [a,b]")
    i = 0
    local c
    while b-a > tol
        i += 1
        i != maxiter || error("Max iteration exceeded")
        c = (a+b)/2
        fc = f(c)
        if fc == 0
            break
        elseif fa*fc > 0
            a = c  # Root is in the right half of [a,b].
            fa = fc
        else
            b = c  # Root is in the left half of [a,b].
        end
    end
    return c
end


text:
using LinearTransforms, BenchmarkTools

m, n = 100, 100
A = randn(m, n)
prod = v -> A*v
tprod = w -> A'*w
op = LinearTransform(m, n, prod, tprod)
x = randn(n)
y = randn(m)

println("Forward native multiply")
@benchmark A*x
println("Forward operator multiply")
@benchmark op*x

println("Ajoint native multiply")
@benchmark A'*y
println("Adjoint operator multiply")
@benchmark op'*y


text:
# Autogenerated wrapper script for Libgpg_error_jll for aarch64-apple-darwin
export libgpg_error

JLLWrappers.@generate_wrapper_header("Libgpg_error")
JLLWrappers.@declare_library_product(libgpg_error, "@rpath/libgpg-error.0.dylib")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libgpg_error,
        "lib/libgpg-error.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# ## Basis function definitions

mutable struct PGBF
    expn::Float64
    x::Float64
    y::Float64
    z::Float64
    I::Int64
    J::Int64
    K::Int64
    norm::Float64
end

function pgbf(expn,x=0,y=0,z=0,I=0,J=0,K=0,norm=1)
    p = PGBF(expn,x,y,z,I,J,K,norm)
    normalize!(p)
    return p
end

function amplitude(bf::PGBF,x,y,z)
    dx,dy,dz = x-bf.x,y-bf.y,z-bf.z
    r2 = dist2(dx,dy,dz)
    return bf.norm*(dx^bf.I)*(dy^bf.J)*(dz^bf.K)*exp(-bf.expn*r2)
end

function normalize!(pbf::PGBF)
    pbf.norm /= sqrt(overlap(pbf,pbf))
end

mutable struct CGBF
    x::Float64
    y::Float64
    z::Float64
    I::Int64
    J::Int64
    K::Int64
    norm::Float64
    pgbfs::Array{PGBF,1}
    coefs::Array{Float64,1}
end

cgbf(x=0,y=0,z=0,I=0,J=0,K=0) = CGBF(x,y,z,I,J,K,1.0,PGBF[],Float64[])

function amplitude(bf::CGBF,x,y,z)
    s = 0
    for (c,pbf) in primitives(bf)
        s += c*amplitude(pbf,x,y,z)
    end
    return bf.norm*s
end

function normalize!(bf::CGBF)
    bf.norm /= sqrt(overlap(bf,bf))
end

primitives(a::CGBF) = zip(a.coefs,a.pgbfs)

function push!(cbf::CGBF,expn,coef)
    Base.push!(cbf.pgbfs,pgbf(expn,cbf.x,cbf.y,cbf.z,cbf.I,cbf.J,cbf.K))
    Base.push!(cbf.coefs,coef)
    normalize!(cbf)
end

function contract(f,a::CGBF,b::CGBF)
    s = 0
    for (ca,abf) in primitives(a)
        for (cb,bbf) in primitives(b)
            s += ca*cb*f(abf,bbf)
        end
    end
    return a.norm*b.norm*s
end

function contract(f,a::CGBF,b::CGBF,c::CGBF,d::CGBF)
    s = 0
    for (ca,abf) in primitives(a)
        for (cb,bbf) in primitives(b)
            for (cc,cbf) in primitives(c)
                for (cd,dbf) in primitives(d)
                    s += ca*cb*cc*cd*f(abf,bbf,cbf,dbf)
                end
            end
        end
    end
    return a.norm*b.norm*c.norm*d.norm*s
end

function test_pgbf()
    @testset "test_pgbf" begin
    s = pgbf(1.0)
    px = pgbf(1.0,0,0,0,1,0,0)
    @test amplitude(s,0,0,0) ≈ 0.71270547
    @test amplitude(px,0,0,0) ≈ 0
    end # testset
end
function test_cgbf()
    @testset "test_cgbf" begin
    c = cgbf(0.0,0.0,0.0)
    push!(c,1,1)
    @test amplitude(c,0,0,0) ≈ 0.71270547
    c2 = cgbf(0,0,0)
    push!(c2,1,0.2)
    push!(c2,0.5,0.2)
    @test overlap(c2,c2) ≈ 1
    end # testset
end


text:
module errors

	function _next_perm(v::Int64)
	    v = unsigned(v);
	    t = (v | unsigned(v - 1)) + 1;
	    t = unsigned(t);
	    w = t | ((div((t & -t), (v & -v)) >> 1) - 1);
	    return convert(Int64, w)
	end

	function _all_bitstrings_of_weight(n, w)
		K = binomial(n, w);
		perm = 2^w-1;
		perms = [];
		for k=1:K
		    perms = [perms; perm];
		    perm = _next_perm(perm);
		end

		E = zeros(n, K);
		for k=1:K
		    E[:, k] = digits(perms[k], 2, n)
		end
		return convert(Array{Bool, 2}, E)
	end

	function _all_bitstrings_up_to_weight(n, w)
		w_0 = convert(Array{Bool, 2}, zeros(n, 1))
		w_1 = convert(Array{Bool, 2}, eye(n))
		# All errors up to weight w_x
		if w > 1
			# All weight 1 errors
			out = [w_0 w_1]
			# All other weights
			for j=2:w
				out = [out _all_bitstrings_of_weight(n, j)]
			end
		elseif w == 1
			out = [w_0 w_1]
		elseif w == 0
			out = w_0
		end
		return out
	end

	# Generate all errors up to and including weights w_x and w_z
	function _errors(n, w_x, w_z)
		# All X errors up to weight w_x
		out_x = _all_bitstrings_up_to_weight(n, w_x)'

		# All Z errors up to weight w_z
		out_z = _all_bitstrings_up_to_weight(n, w_z)'

		# Take the outer product of the errors
		n_x = size(out_z)[1];
		n_z = size(out_x)[1];
		rep_x = repeat(out_x, inner=[n_x; 1]);
		rep_z = repeat(out_z, outer=[n_z; 1]);
		return (rep_x, rep_z)
	end
	export _errors, _all_bitstrings_up_to_weight, _all_bitstrings_of_weight

end


text:
using JSServe, FileIO, Colors, Observables
using JSServe.DOM
using JSServe: App
JSServe.browser_display()
dir = mktempdir()
for i in 1:5
    save(joinpath(dir, "img$i.jpg"), rand(RGB{Colors.N0f8}, 200, 200))
end

app = App() do
    files = JSServe.Asset.(joinpath.(dir, filter(x-> endswith(x, ".jpg"), readdir(dir))))
    slider = JSServe.Slider(1:length(files))
    image_obs = Observable(first(files))
    img = DOM.img(src=image_obs, style="height: 200px;")
    on(slider) do idx
        image_obs[] = files[idx]
    end
    return DOM.div(slider, img)
end

display(app)


text:
# Include this file, that is,
# julia> include("Setup.jl")
# in order to use the package without installing the package...
function addpath!(this)
  new_path = joinpath(@__DIR__, this)
  if !any(LOAD_PATH .== new_path)
    push!(LOAD_PATH, new_path)
  end
end

addpath!("src")
nothing


text:
using LowRankModels

m,n,k = 10,20,3
Y = rand(3,20)
A = rand(10,20)

ry = Regularizer[fixed_latent_features(Y[:,i]) for i=1:n]
glrm = GLRM(A, QuadLoss(), SimplexConstraint(), ry, k+1)
X, Yp, ch = fit!(glrm)

@assert(Yp[1:k,end] == Y)

m,n,k = 10,20,3
Y = rand(3,20)
A = rand(10,20)

ry = Regularizer[fixed_last_latent_features(Y[:,i]) for i=1:n]
glrm = GLRM(A, QuadLoss(), SimplexConstraint(), ry, k+1)
X, Yp, ch = fit!(glrm)

@assert(Yp[2:end,:] == Y)


text:
# Copyright 2019 Tobias Frilling
#
# 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.

"""
    plot_matrix(m; <keyword arguments>)

Visualize the given U-, P- or U*-matrix.

# Arguments
- `projection`: A projection of data points to ESOM neuron to show on the map.
- `classes`: The class of each data point. Used to color the dots for the projection.
- `normalize`: Whether to normalize the given matrix. Default: `true`.
- `colorStyle`: Either `:umatrix` or `:pmatrix`. Default: `:umatrix`.
"""
@userplot Plot_Matrix

@recipe function f(p::Plot_Matrix;
                   settings = defaultSettings,
                   projection = Dict(),
                   key = 1:length(projection),
                   classes = Dict(),
                   normalize = true,
                   colorStyle = :umatrix)

    m = p.args[1]
    if normalize
        (min_m, max_m) = quantile(m[:], [.01, .99])
        m = (m .- min_m) ./ (max_m - min_m)
    end
    (rows, columns) = size(m)

    grid := false
    framestyle := :none
    yflip --> true
    colorbar --> true
    aspect_ratio --> 1

    @series begin
        seriestype := :contourf
        (min_m, max_m) = quantile(m[:], [.01, .99])
        levels --> round(Int, max_m / max(min_m, .05))
        seriescolor --> if colorStyle == :pmatrix colormap_pmatrix
                        else colormap_umatrix end
        linewidth --> 0.1
        if settings.toroid
            [m m; m m]
        else
            m
        end
    end

    if !isempty(projection)
        seriestype := :scatter
        widen := false
        legend := false
        colorbar_entry := false
        markersize --> 2.5
        markerstrokewidth --> 0.4

        if isempty(classes)
            classes[1] = keys(projection) |> collect
        end
        colors = distinguishable_colors(length(classes), RGB(1,0,1))

        for (class, indices) in classes
            @series begin
                markercolor := colors[class]
                v = [projection[key[index]] for index in indices]
                if settings.toroid
                    v = vcat(v, map(i -> i + CartesianIndex(rows, 0), v))
                    v = vcat(v, map(i -> i + CartesianIndex(0, columns), v))
                end
                getindex.(v,2), getindex.(v,1)
            end
        end
    end

end

colormap_umatrix = begin
    local c = ["#3C6DF0", "#3C6DF0", "#3C6DF0", "#006602", "#006A02", "#006D01",
               "#007101", "#007501", "#007901", "#007C00", "#008000", "#068103",
               "#118408", "#0B8305", "#17860A", "#1D870D", "#228810", "#288A12",
               "#2E8B15", "#348D18", "#398E1A", "#3F8F1D", "#45911F", "#4A9222",
               "#509325", "#569527", "#5C962A", "#61982C", "#67992F", "#6D9A32",
               "#729C34", "#789D37", "#7E9F39", "#84A03C", "#89A13F", "#8FA341",
               "#95A444", "#9AA547", "#A0A749", "#A6A84C", "#ACAA4E", "#B1AB51",
               "#B7AC54", "#BDAE56", "#C3AF59", "#C8B15B", "#CEB25E", "#CBAF5C",
               "#C8AC59", "#C5A957", "#C3A654", "#C0A352", "#BDA050", "#BA9D4D",
               "#B7994B", "#B49648", "#B29346", "#AF9044", "#AC8D41", "#A98A3F",
               "#A6873C", "#A3843A", "#A08138", "#9E7E35", "#9B7B33", "#987830",
               "#95752E", "#92722B", "#8F6E29", "#8C6B27", "#8A6824", "#876522",
               "#84621F", "#815F1D", "#7E5C1B", "#7B5918", "#795616", "#765313",
               "#714E0F", "#6C480B", "#674307", "#6F4D15", "#785822", "#806230",
               "#896D3E", "#91774C", "#998159", "#A28C67", "#AA9675", "#B3A183",
               "#BBAB90", "#C3B59E", "#CCC0AC", "#D4CABA", "#DDD5C7", "#E5DFD5",
               "#E7E1D8", "#E9E4DB", "#EBE6DE", "#ECE8E1", "#EEEAE4", "#F0EDE7",
               "#F2EFEA", "#F4F1ED", "#F6F4F0", "#F8F6F3", "#F9F8F6", "#FBFAF9",
               "#FDFDFC", "#FFFFFF", "#FFFFFF", "#FEFEFE", "#FEFEFE", "#FEFEFE",
               "#FDFDFD", "#FDFDFD", "#FDFDFD", "#FCFCFC", "#FCFCFC", "#FCFCFC",
               "#FBFBFB", "#FBFBFB", "#FBFBFB", "#FAFAFA", "#FAFAFA", "#FAFAFA",
               "#F9F9F9", "#F9F9F9", "#FFFFFF", "#FFFFFF"]
    (s -> parse(RGB, s)).(c) |> ColorGradient
end

colormap_pmatrix = begin
    local c = ["#FFFFFF", "#FFFFF7", "#FFFFEF", "#FFFFE7", "#FFFFDF", "#FFFFD7",
               "#FFFFCF", "#FFFFC7", "#FFFFBF", "#FFFFB7", "#FFFFAF", "#FFFFA7",
               "#FFFF9F", "#FFFF97", "#FFFF8F", "#FFFF87", "#FFFF80", "#FFFF78",
               "#FFFF70", "#FFFF68", "#FFFF60", "#FFFF58", "#FFFF50", "#FFFF48",
               "#FFFF40", "#FFFF38", "#FFFF30", "#FFFF28", "#FFFF20", "#FFFF18",
               "#FFFF10", "#FFFF08", "#FFFF00", "#FFFA00", "#FFF400", "#FFEF00",
               "#FFEA00", "#FFE400", "#FFDF00", "#FFDA00", "#FFD400", "#FFCF00",
               "#FFCA00", "#FFC500", "#FFBF00", "#FFBA00", "#FFB500", "#FFAF00",
               "#FFAA00", "#FFA500", "#FF9F00", "#FF9A00", "#FF9500", "#FF8F00",
               "#FF8A00", "#FF8500", "#FF8000", "#FF7A00", "#FF7500", "#FF7000",
               "#FF6A00", "#FF6500", "#FF6000", "#FF5A00", "#FF5500", "#FF5000",
               "#FF4A00", "#FF4500", "#FF4000", "#FF3A00", "#FF3500", "#FF3000",
               "#FF2B00", "#FF2500", "#FF2000", "#FF1B00", "#FF1500", "#FF1000",
               "#FF0B00", "#FF0500", "#FF0000", "#FA0000", "#F40000", "#EF0000",
               "#EA0000", "#E40000", "#DF0000", "#DA0000", "#D40000", "#CF0000",
               "#CA0000", "#C50000", "#BF0000", "#BA0000", "#B50000", "#AF0000",
               "#AA0000", "#A50000", "#9F0000", "#9A0000", "#950000", "#8F0000",
               "#8A0000", "#850000", "#800000", "#7A0000", "#750000", "#700000",
               "#6A0000", "#650000", "#600000", "#5A0000", "#550000", "#500000",
               "#4A0000", "#450000", "#400000", "#3A0000", "#350000", "#300000",
               "#2B0000", "#250000", "#200000", "#1B0000", "#150000", "#100000",
               "#0B0000", "#050000"]
    (s -> parse(RGB, s)).(c) |> ColorGradient
end


text:
function add_compat_section!(project::Dict)
    if !haskey(project, "compat")
        project["compat"] = Dict{Any, Any}()
    end
    return project
end

function get_random_string()
    return string(utc_to_string(now_localzone()), "-", rand(UInt32))::String
end

function generate_pr_title_parenthetical(keep_or_drop::Symbol,
                                         parenthetical_in_pr_title::Bool)::String
    if parenthetical_in_pr_title
        if keep_or_drop == :keep
            return " (keep existing compat)"
        elseif keep_or_drop == :drop
            return " (drop existing compat)"
        else
            return ""
        end
    else
        return ""
    end
end


text:
function make_model(grid, ν; timestepper = :RungeKutta3, advection = WENO5())

    v_bcs, w_bcs = no_slip_bc(grid)

    return IncompressibleModel(timestepper = timestepper,
                                 advection =  advection,
                                      grid = grid,
                                  buoyancy = BuoyancyTracer(),
                                   tracers = :b,
                       boundary_conditions = (v=v_bcs, w=w_bcs),
                                   closure = IsotropicDiffusivity(ν=ν, κ=0))
end


text:
# Some help functions to generate interesting test cases
circle_segment(N) = sin.(acos.(range(-1, stop=1, length=N)))
linear_trend(N) = (0:N-1)

""" brute_force_multi(g, M, t=1:length(g); tol=1e-3)
    Finds the solution to the cardinality constrained problem for
    m=1:M segments using brute force with possible breakpoints on grid t.
    `t` has to be supplied if `g` is not `AbstractArray`
"""
function brute_force_multi(g, M, t=1:length(g); tol=1e-3)

    if length(g) <= M
        @warn "Specified M too large. Setting M=length(g)-1."
        M = length(g)-1
    end

    ℓ, V_N = EllZeroTrendFiltering.get_transition_costs(g, t, true, tol=tol)

    I_vec = Vector{Vector{Int64}}(undef, M)
    Y_vec = Vector{Vector{Float64}}(undef, M)
    f_vec = Vector{Float64}(undef, M)

    for m=1:M
        (I_bf, Y_bf, f_bf) = brute_force_search(ℓ, V_N, m)
        I_vec[m] = t[I_bf]
        Y_vec[m] = Y_bf
        f_vec[m] = f_bf
    end

    return I_vec, Y_vec, f_vec
end


text:
using JuMP, Clp

rpomdp = Baby3RPOMDP()
dim = n_observations(rpomdp)
uncset = observation(rpomdp,true,true)
cdfl = cumsum(pdf.(uncset.lower))
cdfu = cumsum(pdf.(uncset.upper))

coeff = [5,2,3]
m = Model(solver = ClpSolver())
@variable(m, 0 <= p[i = 1:dim] <= 1)
@objective(m, Min, dot(coeff,p))
@constraint(m, sum(p) == 1)
@constraint(m, cdfl[1] <= p[1] <= cdfu[1])
for i = 2:dim
    pl = cdfl[i] - cdfu[i-1]
    pu = cdfu[i] - cdfl[i-1]
    @constraint(m, pl <= p[i] <= pu)
end
JuMP.solve(m)
argmin = getvalue(p)
obj = getobjectivevalue(m)

rpomdp = Baby3RPOMDP()
dist_t = transition(rpomdp,true,true)
dist = observation(rpomdp,true,true)
dim = n_observations(rpomdp)
cdfl = cumsum(pdf.(uncset.lower))
cdfu = cumsum(pdf.(uncset.upper))
pt = [0.5,0.6]
α = [1,2]
m = Model(solver = ClpSolver())
@variable(m, 0 <= p[i = 1:dim] <= 1)
@objective(m, Min, dot(α, p .* pt))
@constraint(m, cdfl[zind] <= p <= cdfu[zind])
JuMP.solve(m)
argmin = getvalue(p)
obj = getobjectivevalue(m)

dist = observation(prob,a,sp)
dim = n_observations(rpomdp)
function robustdist(uncset::PBox, dim::Int, α::Vector{Float64})
    cdfl = cumsum(pdf.(uncset.lower))
    cdfu = cumsum(pdf.(uncset.upper))
    m = Model(solver = ClpSolver())
    @variable(m, 0 <= p[i = 1:dim] <= 1)
    @objective(m, Min, dot(α, p .* pt))
    @constraint(m, sum(p) == 1)
    @constraint(m, cdfl[1] <= p[1] <= cdfu[1])
    for i = 2:dim
        pl = cdfl[i] - cdfu[i-1]
        pu = cdfu[i] - cdfl[i-1]
        @constraint(m, pl <= p[i] <= pu)
    end
    JuMP.solve(m)
    argmin = getvalue(p)
    obj = getobjectivevalue(m)
end


text:
# 1
wavelet_matrices(ws, J, L; ξ0=maximum(ws) / 2, c=1) =
    [wavelet_matrix(ws, j, θ; ξ0, L, c) for j in 0:J - 1, θ in range(0., length=L, step=2π / L)]

# 2
e_θ(l, L) = (cos(π / 2 + 2π * l / L), sin(π / 2 + 2π * l / L))
τ_θ_j(l, j_, L) = round.(Int, 2^j_ .* e_θ(l, L))

# ! k > 1 doesn't work if not in Γ_H

function default_Γ_H(J, L, K)
    Γ_H = NamedTuple{
        (:j, :l, :θ, :k, :j_, :l_, :θ_, :k_, :τ_),
        Tuple{Int64, Int64, Float64, Int64, Int64,
        Int64, Float64, Int64, Tuple{Int64, Int64}}
    }[]
    p = Iterators.product(
        0:L - 1,
        0:L - 1,
        0:J - 1,
        0:J - 1,
        [false, true],
        0:K,
        0:K
    )

    for (l, l_, j, j_, isshift, k, k_) in p
        l, l_, j, j_, isshift, k, k_
        if l_ > l
            continue
        end
        if (j > j_) || (j < j_ - 2)
            continue
        end
        if k > min(k_, 1)
            continue
        end

        if isshift
            τ_ = τ_θ_j(l, j_, L)
        else
            τ_ = (0, 0)
        end

        θ = 2π * l / L
        θ_ = 2π * l_ / L
        dl = abs(l - l_)

        if j == j_
            if k == 0 && k_ > 1
                continue
            end
            if k == 1
                if k_ != 1
                    continue
                elseif dl > 2
                    continue
                end
            end
        else
            if k == 0 && k_ == 4
                continue
            end
            if k == 1
                if k_ != 2^(j_ - j)
                    continue
                elseif dl > 2
                    continue
                end
            end
        end

        push!(Γ_H, (;j, l, θ, k, j_, l_, θ_, k_, τ_))
    end
    Γ_H
end

# 3
"Fourier transform of δ(x - x0)"
# Δ(w, x0) = exp(-2π * im * dot(x0, w))
"Fourier transform of point process"
# M(Ws, μ::AbstractMatrix) = [sum(Δ(w, x) for x in eachcol(μ)) for w in Ws]

function Δ(x0, ws)
    a = x0[1] .* ws
    b = x0[2] .* ws
    -2pi * (a' .+ b) .|> cis
end

M(x, ws) = sum(Δ(xi, ws) for xi in eachcol(x))

# main
struct WaveletParams{Tw}
    s::Float64
    N::Int
    J::Int
    L::Int
    K::Int
    σ::Float64

    ws::Tw
    Ws::Matrix{Vector{Float64}}
    FN::Function

    Ψs::Matrix{Matrix{Float64}}

    Γ_H::Vector{NamedTuple{
        (:j, :l, :θ, :k, :j_, :l_, :θ_, :k_, :τ_),
        Tuple{Int64, Int64, Float64, Int64, Int64,
        Int64, Float64, Int64, Tuple{Int64, Int64}}}
        }
    jlk::Array{Tuple{Int64, Int64, Int64}, 3}
end

function WaveletParams(s, N, J, L, K, σ, Γ_H = default_Γ_H(J, L, K))
    dx = 2s / N
    xs = range(-s; length=N, step=dx)
    ws = fftfreq(N, 1 / dx)
    Ws = [[w1, w2] for w1 in ws, w2 in ws]

    g = MvNormal([0, 0], σ^2 * I)
    gaus = [pdf(g, [x, y]) for x in xs, y in xs]
    Gaus = gaus |> ifftshift |> fft |> real

    # FN(μ) = M(Ws, μ)
    FN(μ) = M(μ, ws)

    Ψs = [Ψ .* Gaus for Ψ in wavelet_matrices(ws, J, L)]

    jlk = collect(Iterators.product(1:J, 1:L, 1:K + 1))

    WaveletParams(s, N, J, L, K, σ, ws, Ws, FN, Ψs, Γ_H, jlk)
end

phase_harmonics(z::Complex, k::Int) = abs(z) * exp(im * angle(z) * k)

function wavelet_phase_harmonics(μ, wp, vs)
    M = wp.FN(μ)
    # cs = [ifft(M .* Ψ) for Ψ in wp.Ψs]
    cs = map(Ψ -> ifft(M .* Ψ), wp.Ψs)
    m = ifft(M) .- vs[1]
    # csk = [phase_harmonics.(cs[j, l], k - 1) .- vs[2][j, l, k] for (j, l, k) in wp.jlk]
    csk = map(wp.jlk) do jlk
        (j, l, k) = jlk
        phase_harmonics.(cs[j, l], k - 1) .- vs[2][j, l, k]
    end
    m, csk
end

function v_λ_k_all(μ, wp::WaveletParams)
    m, csk = wavelet_phase_harmonics(μ, wp, (0., zeros(wp.J, wp.L, wp.K + 1)))
    mean(m), mean.(csk)
end

function K_μ(csk; j, l, k, j_, l_, k_, τ_)
    A = csk[j + 1, l + 1, k + 1]
    B = circshift(csk[j_ + 1, l_ + 1, k_ + 1], τ_)
    mean(@. A * conj(B))
end

function K_μ_0(m, csk)
    map(c -> mean(m .* conj.(c)), csk)
end

K_μ(csk, p) = K_μ(csk; p.j, p.l, p.k, p.j_, p.l_, p.k_, p.τ_)

function K_all(μ, wp::WaveletParams, vs)
    m, csk = wavelet_phase_harmonics(μ, wp, vs)
    K_main = map(p -> K_μ(csk, p), wp.Γ_H)
    K_0 = K_μ_0(m, csk)[:]
    [K_main; K_0]
end

function lossW(x; wp, vs, K0, weights, index)
    k = (K_all(x, wp, vs) .* weights)[index]
    k0 = (K0 .* weights)[index]
    sqL2dist(k, k0)
end
lossW(x, p) = lossW(x; p...)

function lossW_params(x, wp, weights=ones(length(wp.Γ_H) + length(wp.jlk)), index=:)
    vs = v_λ_k_all(x, wp)
    K0 = K_all(x, wp, vs)
    (; wp, vs, K0, weights, index)
end


text:
function read_struct_tag!(io::IO, v::Int64=0)
  s = read(io, UInt8)
  s == 0x53 || (close(io); error("damaged or scrambled SUDS file; can't continue."))
  m = read(io, UInt8)
  SB.sid = read(io, Int16)
  SB.nbs = read(io, Int32)
  SB.nbx = read(io, Int32)

  (v > 0) && println("suds_structtag: machine code = ", Char(m),
                      ", ID = ",          SB.sid,
                      ", struct size = ", SB.nbs, " B",
                      ", data size = ",   SB.nbx, " B",
                      )

  return nothing
end

function staident!(io::IO, v::Int64=0)
  read!(io, SB.hdr)
  fill!(SB.id, 0x00)

  # Set ID string
  j = 1
  for i = 1:2
    if SB.hdr[i] != 0x00
      SB.id[j] = SB.hdr[i]
      j += 1
    end
  end
  SB.id[j] = 0x2e
  j += 1
  for i = 5:9
    if SB.hdr[i] != 0x00
      SB.id[j] = SB.hdr[i]
      j += 1
    end
  end
  SB.id[j] = 0x2e
  SB.id[j+1] = 0x2e
  SB.id[j+2] = SB.hdr[10]
  j += 2
  setfield!(SB, :id_str, unsafe_string(pointer(getfield(SB, :id)), j))

  # Print new ID to STDOUT if v > 1
  (v > 1) && println( "id  = ",         SB.id_str,
                      ", inst code = ", Int16(SB.hdr[12]) << 8 | Int16(SB.hdr[11]))
  return nothing
end

function read_chansetentry!(S::GphysData, io::IO, v::Int64)
  SB.C.inst    = read(io, Int32)
  SB.C.stream  = read(io, Int16)
  SB.C.chno    = read(io, Int16)
  staident!(io, 0)
  if v > 2
    println("id = ",          SB.id_str,
            ", inst code = ", Int16(SB.hdr[12]) << 8 | Int16(SB.hdr[11]),
            ", inst = ",      SB.C.inst,
            ", stream = ",    SB.C.stream,
            ", chno = ",      SB.C.chno)
  end
  return nothing
end

# #  STAT_IDENT:  Station identification
# read_1!(S::GphysData, io::IO, v::Int64, full::Bool) = staident!(io, v)
#
# #  STRUCTTAG:  Structure to identify structures when archived together
# read_2!(S::GphysData, io::IO, v::Int64, full::Bool) = read_struct_tag!(io)

# STATIONCOMP:  Generic station component information
function read_5!(S::GphysData, io::IO, v::Int64, full::Bool)
  staident!(io, 0)
  SB.S05.az             = read(io, Int16)    # azimuth, N°E
  SB.S05.inc            = read(io, Int16)    # incidence, from vertical
  SB.S05.lat            = read(io, Float64)  # latitude, N = +
  SB.S05.lon            = read(io, Float64)  # longitude, E = +
  SB.S05.ele            = read(io, Float32)  # elevation, meters
  read!(io, SB.S05.codes)
  read!(io, SB.S05.gain)
  read!(io, SB.S05.a2d)
  SB.t_i32              = read(io, Int32)    # date/time values became effective
  read!(io, SB.S05.t_corr)

  # Post-read processing
  SB.data_type  = SB.S05.codes[9]
  sensor_type   = SB.S05.codes[8]

  # Does this station exist in S?
  i = findid(SB.id_str, S)
  if i == 0
    loc = GeoLoc( "",
                  SB.S05.lat,
                  SB.S05.lon,
                  Float64(SB.S05.ele),
                  0.0,
                  Float64(SB.S05.az),
                  Float64(SB.S05.inc)
                  )
    misc = if full
            Dict{String, Any}(
              "ic"             => Int16(SB.hdr[12]) << 8 | Int16(SB.hdr[11]),
              "enclosure"      => Char(SB.S05.codes[1]),
              "annotation"     => Char(SB.S05.codes[2]),
              "recorder"       => Char(SB.S05.codes[3]),
              "rockclass"      => Char(SB.S05.codes[4]),
              "rocktype"       => Int16(SB.S05.codes[6]) << 8 | Int16(SB.S05.codes[5]),
              "sitecondition"  => Char(SB.S05.codes[7]),
              "sensor_type"    => Char(sensor_type),
              "data_type"      => Char(SB.data_type),
              "data_units"     => Char(SB.S05.codes[10]),
              "polarity"       => Char(SB.S05.codes[11]),
              "st_status"      => Char(SB.S05.codes[12]),
              "max_gain"       => SB.S05.gain[1],
              "clip_value"     => SB.S05.gain[2],
              "con_mvolts"     => SB.S05.gain[3],
              "a2d_ch"         => SB.S05.a2d[1],
              "a2d_gain"       => SB.S05.a2d[2],
              "dt_eff"         => SB.t_i32,
              "clock_corr"     => SB.S05.t_corr[1],
              "sta_delay"      => SB.S05.t_corr[2]
              )
          else
            Dict{String, Any}(
              "data_type"      => Char(SB.data_type),
              "sensor_type"    => Char(sensor_type),
              )
          end

    push!(S, SeisChannel( id = SB.id_str,
                          loc = loc,
                          misc = misc,
                          gain = Float64(SB.S05.gain[1]),
                          units = sensor_types[sensor_type]
                          )
         )
    i = S.n
  end
  return i
end

# MUXDATA:  Header for (possibly) multiplexed data
function read_6!(S::GphysData, io::IO, v::Int64, full::Bool)
  SB.T.net      = read(io, UInt16)
  skip(io, 2)
  SB.t_f64      = read(io, Float64)
  SB.t_i16      = read(io, Int16)
  SB.T.n_ch     = read(io, Int16)
  SB.fs         = read(io, Float32)
  SB.data_type  = read(io, UInt8)
  skip(io, 3)
  SB.T.ns       = read(io, Int32)
  SB.nz         = read(io, Int32)

  if v > 1
    println("net = ",         Char(SB.T.net & 0x00ff), Char(SB.T.net >> 8),
            ", btime = ",     u2d(SB.t_f64),
            ", loctime = ",   SB.t_i16,
            ", n = ",         SB.T.n_ch,
            ", fs = ",        SB.fs,
            ", data_type = ", Char(SB.data_type),
            ", n sweeps = ",  SB.T.ns,
            ", nx = ",        SB.nz
            )
  end

  # return (net, data_type, btime, loctime, numchans, fs, data_type, nx)
  return nothing
end

#=  DESCRIPTRACE:  Descriptive information about a seismic trace.
                   Normally followed by waveform data =#
function read_7!(S::GphysData, io::IO, v::Int64, full::Bool)
  staident!(io, 0)
  SB.t_f64       = read(io, Float64)
  SB.t_i16       = read(io, Int16)
  SB.data_type   = read(io, UInt8)
  SB.T.desc      = read(io, UInt8)
  skip(io, 4)
  SB.nz          = read(io, Int32)
  SB.fs          = read(io, Float32)
  skip(io, 16)
  SB.t_f64      += read(io, Float64)
  SB.rc          = read(io, Float32)

  if v > 2
    tl = SB.t_i16
    println("id = ",          SB.id_str,
            ", inst code = ", Int16(SB.hdr[12]) << 8 | Int16(SB.hdr[11]),
            ", begin = ",     u2d(SB.t_f64 + tl*60),
            " (GMT ", tl < 0 ? "-" : "+", div(tl, 60), ")",
            ", data_type = ", Char(SB.data_type),
            ", nx = ",        SB.nx,
            ", fs = ",        SB.fs + SB.rc)
  end

  # return (id, data_type, nx, Float64(fs+rc), btime+tc)
  return nothing
end

# FEATURE:  Observed phase arrival time, amplitude, and period
function read_10!(S::GphysData, io::IO, v::Int64, full::Bool)
  staident!(io, 0)
  SB.P.pc           = read(io, Int16)
  SB.P.onset        = read(io, UInt8)
  SB.P.fm           = read(io, UInt8)
  SB.P.snr          = read(io, Int16)
  skip(io, 4)
  SB.P.gr           = read(io, Int16)
  SB.t_f64          = read(io, Float64)
  SB.P.amp          = read(io, Float32)
  SB.t_f32          = read(io, Float32)
  SB.t_i32          = read(io, Int32)
  skip(io, 4)

  phase = pick_types[SB.P.pc]
  (v > 1) && println(
            "id  = ",           SB.id_str,
            ", inst code = ",   Int16(SB.hdr[12]) << 8 | Int16(SB.hdr[11]),
            ", phase: ",        phase, Char(SB.P.onset), Char(SB.P.fm),
            ", t = ",           u2d(SB.t_f64),
            ", A = ",           SB.P.amp,
            ", τ = ",           SB.t_f32,
            ", SNR = ",         SB.P.snr,
            ", gain range = ",  SB.P.gr)

  if isa(S, EventTraceData)
    i = findid(SB.id_str, S)
    if i == 0
      push!(S, SeisChannel(id = SB.id_str))
      i = S.n
    end
    S.pha[i][phase] = SeisPha(amp = Float64(SB.P.amp),
                              tt = SB.t_f64,
                              pol = Char(SB.P.fm),
                              )
  end

  return nothing
end

#  ORIGIN: Information about a specific solution for a given event
function read_14!(S::GphysData, io::IO, v::Int64, full::Bool)
  SB.H.evno         = read(io, Int32)
  SB.H.auth         = read(io, Int16)
  read!(io, SB.H.chars)
  SB.H.reg          = read(io, Int32)
  SB.H.ot           = read(io, Float64)
  SB.H.lat          = read(io, Float64)
  SB.H.lon          = read(io, Float64)
  read!(io, SB.H.floats)
  read!(io, SB.H.model)
  SB.H.gap          = read(io, Int16)
  SB.H.d_min        = read(io, Float32)
  read!(io, SB.H.shorts)
  read!(io, SB.H.mag)
  SB.t_i32          = read(io, Int32)
  if v > 1
    println("evno = ",      SB.H.evno,
            ", auth = ",    SB.H.auth,
            ", codes = ",   map(Char, SB.H.chars),
            ", reg = ",     SB.H.reg,
            ", ot = ",      u2d(SB.H.ot),
            ", origin = ",  SB.H.lat, "N, ", SB.H.lon, "E, z ", SB.H.floats[1], " km",
            ", δx = ",      SB.H.floats[2], " km",
            ", δz = ",      SB.H.floats[3], " km",
            ", rms = ",     SB.H.floats[4],
            ", model = ",   String(copy(SB.H.model)),
            ", Δ = ",       SB.H.gap, "∘",
            ", d_min = ",   SB.H.d_min, " km",
            ", shorts = ",  SB.H.shorts,
            ", mag = ",     SB.H.mag,
            ", t_eff = ",   u2d(Float64(SB.t_i32))
            )
  end
  return nothing
end

# COMMENT:  Comment tag to be followed by the bytes of comment
function read_20!(S::GphysData, io::IO, v::Int64, full::Bool)
  read!(io, SB.comm_i)
  L = read(io, Int16)
  skip(io, 2)
  if v > 1
    println("structure ref. ID = ",   SB.comm_i[1],
            ", item in struct ID = ", SB.comm_i[2],
            ", L = ", L)
  end
  SB.comm_s = String(copy(read(io, L)))
  for i in 1:S.n
    if haskey(S.misc[i], "comment")
      append!(S.misc[i]["comment"], SB.comm_s)
    else
      S.misc[i]["comment"] = [SB.comm_s]
    end
  end
  (v > 2) && (printstyled("comment: \n", color=:green); println(SB.comm_s); printstyled("--\n", color=:green))
  return nothing
end

# TRIGGERS:  Earthquake detector trigger statistics
function read_25!(S::GphysData, io::IO, v::Int64, full::Bool)
  staident!(io, 0)
  shorts    = read!(io, zeros(Int16, 6))
  trig_time = read(io, Float64)
  if v > 1
    println("id = ",          SB.id_str,
            ", inst code = ", Int16(SB.hdr[12]) << 8 | Int16(SB.hdr[11]),
            ", shorts = ", shorts,
            ", btime = ", u2d(trig_time))
  end
  return nothing
end

# TRIGSETTING:  Settings for earthquake trigger system
function read_26!(S::GphysData, io::IO, v::Int64, full::Bool)
  net       = read!(io, zeros(UInt8, 4))
  btime     = read(io, Float64)
  shorts    = read!(io, zeros(Int16, 6))
  t_sweep   = read(io, Float32)
  t_aper    = read(io, Float32)
  alg       = read(io, UInt8)
  skip(io, 3)
  if v > 1
    println("net = ", String(copy(net)),
            ", btime = ", u2d(btime),
            ", shorts = ", shorts,
            ", t_sweep = ", t_sweep,
            ", t_aperture = ", t_aper,
            ", alg = ", Char(alg))
  end

  return nothing
end

# EVENTSETTING:  Settings for earthquake trigger system
function read_27!(S::GphysData, io::IO, v::Int64, full::Bool)
  net       = read!(io, zeros(UInt8, 4))
  btime     = read(io, Float64)
  shorts    = read!(io, zeros(Int16, 4))
  dur_min   = read(io, Float32)
  dur_max   = read(io, Float32)
  alg       = read(io, UInt8)
  skip(io, 3)
  if v > 1
    println("net = ", String(copy(net)),
            ", btime = ", u2d(btime),
            ", shorts = ", shorts,
            ", dur_min = ", dur_min,
            ", dur_max = ", dur_max,
            ", alg = ", Char(alg))
  end

  return nothing
end

# DETECTOR
function read_28!(S::GphysData, io::IO, v::Int64, full::Bool)
  alg           = read(io, UInt8)
  event_type    = read(io, UInt8)
  net_node_id   = read!(io, Array{UInt8,1}(undef, 10))
  version       = read(io, Float32)
  event_num     = read(io, Int32)
  skip(io, 4)

  if v > 1
    println("algorithm = ", Char(alg),
            ", event_type = ", Char(event_type),
            ", net_node_id = ", String(net_node_id),
            ", version = ", version,
            ", event_num = ", event_num
            )

  end
  return nothing
end

# ATODINFO
function read_29!(S::GphysData, io::IO, v::Int64, full::Bool)
  io_addr       = read(io, Int16)
  dev_id        = read(io, Int16)
  dev_flags     = read(io, UInt16)
  ext_bufs      = read(io, Int16)
  ext_mux       = read(io, Int16)
  timing_src    = read(io, UInt8)
  trig_src      = read(io, UInt8)

  if v > 1
    println("io_addr = ", io_addr,
            ", dev_id = ", dev_id,
            ", dev_flags = ", bitstring(dev_flags),
            ", ext_bufs = ", ext_bufs,
            ", ext_mux = ", ext_mux,
            ", timing_src = ", timing_src,
            ", trig_src = ", trig_src
            )

  end
  return nothing
end

# TIMECORRECTION:  Time correction information
function read_30!(S::GphysData, io::IO, v::Int64, full::Bool)
  staident!(io, 0)
  SB.tc          = read(io, Float64)
  SB.rc          = read(io, Float32)
  SB.sync_code   = Char(read(io, UInt8))
  skip(io,1)
  SB.t_i32       = read(io, Int32)
  skip(io, 2)

  if v > 1
    println("tc = ",          SB.tc,
            ", rc = ",        SB.rc,
            ", sync_code = ", SB.sync_code,
            ", t_eff = ",     SB.t_i32)
  end

  # Check: is the second part of ID "IRIG" (i.e., a local time clock)?
  SB.irig = false
  id_b = SB.id
  i1 = findfirst(SB.id.==0x2e)+1
  i2 = findnext(SB.id.==0x2e, i1+1)-1
  if SB.id[i1:i2] == [0x49, 0x52, 0x49, 0x47] # "IRIG"
    SB.irig = true
  end

  # return (round(Int64, tc*1.0e6), rc, t_eff*1000000, SB.id_str, netflg)
  return nothing
end

#= INSTRUMENT: Instrument hardware settings, mainly PADS related
                    added by R. Banfill, Jan 1991 =#
function read_31!(S::GphysData, io::IO, v::Int64, full::Bool)
  staident!(io, 0)
  serial    = read(io, Int16)
  ncmp      = read(io, Int16)
  chno      = read(io, Int16)
  sens_type = read(io, UInt8)
  data_type = read(io, UInt8)
  void_samp = read(io, Int32)
  floats    = read!(io, zeros(Float32,10))
  t_eff     = read(io, Int32)
  pre_evt   = read(io, Float32)
  trig_num  = read(io, Int16)
  study     = read(io, 6)
  sn_serial = read(io, Int16)
  if v > 1
    println("id = ",            SB.id_str,
            ", inst code = ",   Int16(SB.hdr[12]) << 8 | Int16(SB.hdr[11]),
            ", serial = ",      serial,
            ", ncmp = ",        ncmp,
            ", chno = ",        chno,
            ", sens_type = ",   Char(sens_type),
            ", data_type = ",   Char(data_type),
            ", void_samp = ",   void_samp,
            ", floats = ",      floats,
            ", t_eff = ",       u2d(t_eff),
            ", pre_evt = ",     pre_evt,
            ", trig_num = ",    trig_num,
            ", study = ",       String(study),
            ", sn_serial = ",   sn_serial
            )
  end
  return nothing
end

# CHANSET:  Associate stations and components with sets...???
function read_32!(S::GphysData, io::IO, v::Int64, full::Bool)
  SB.C.typ        = read(io, Int16)
  SB.C.n          = read(io, Int16)
  read!(io, SB.C.sta)
  skip(io, 1)
  SB.C.tu         = read(io, Int32)
  SB.C.td         = read(io, Int32)
  if v > 1
    println("type = ",    SB.C.typ,
            ", n = ",     SB.C.n,
            ", sta = ",   String(copy(SB.C.sta)),
            ", t_up = ",  u2d(SB.C.tu),
            ", t_dn = ",  u2d(SB.C.td)
            )
    (v > 2) && printstyled("channel entries:\n", color=:green)
  end

  for i = 1:SB.C.n
    read_chansetentry!(S, io, v)
  end

  return nothing
end

# # CHANSETENTRY
# read_33!(S::GphysData, io::IO, v::Int64, full::Bool) = read_chansetentry!(S, io, v)
#
function suds_support()
  println("\nCurrent support for SUDS structures\n")
  printstyled("CODE  STRUCTURE       \n", color=:green, bold=true)
  printstyled("====  =========       \n", color=:green, bold=true)
  for i in 1:33
    str = lpad(i, 4) * "  " * suds_codes[i] * "\n"
    if i in unsupported
      printstyled(str, color=:red)
    elseif i in disp_only
      printstyled(str, color=11)
    else
      printstyled(str, color=:green, bold=true)
    end
  end
  printstyled("\n(supported = ")
  printstyled("GREEN", color=:green, bold=true)
  printstyled(", logging only = ")
  printstyled("YELLOW", color=11)
  printstyled(", unsupported = ")
  printstyled("RED", color=:red)
  printstyled(")\n")
  return nothing
end


text:
#!/usr/bin/env julia
using Plots, Random, Printf, Plots.Measures

function initCondGen() #get random initial conditions for mass/radius, position, and velocity
    function getMass(nBodies) #generate random masses that better reflect actual stellar populations
        mList=zeros(nBodies)
        N=(0.5^(-1.3)-150^(-1.3))/1.3 #crude approximation of IMF integral assuming alpha = 2.3, stellar mass range of 0.5:150 solar masses
        rescale=1e6
        max=floor(Int,N*rescale)
        for i=1:nBodies
            intTarget=rand(0:max,1)[1]/rescale
            m=(0.5^(-1.3)-intTarget*1.3)^(-1/1.3) #just algebra from above
            mList[i]=round(m,digits=2)
        end
        return mList
    end
    m=rand(1:1500,3)./10 #3 random masses between 0.1 and 150 solar masses, uniform distribution
    #m=getMass(3) #get mass from IMF -- this way is kind of boring...so not using it, but left here in case I change my mind?
    rad=m.^0.8 #3 radii based on masses in solar units
    m=m.*2e30 #convert to SI kg
    rad=rad.*7e8 #convert to SI m
    pos1=rand(-10:10,3) #random initial coordinates x & y for first body, AU
    function genPos2(pos1)
        accept2=false
        while accept2==false
            pos2=rand(-10:10,3) #random initial coordinates for second body, AU
            dist21=sqrt((pos1[1]-pos2[1])^2+(pos1[2]-pos2[2])^2+(pos1[3]-pos2[3])^2)
            if (dist21*1.5e11)>(rad[1]+rad[2]) #they aren't touching
                accept2=true
                return pos2
            end
        end
    end
    pos2=genPos2(pos1)
    function genPos3(pos1,pos2)
        accept3=false
        while accept3==false
            pos3=rand(-10:10,3) #random initial coordinates for third body, AU
            dist31=sqrt((pos1[1]-pos3[1])^2+(pos1[2]-pos3[2])^2+(pos1[3]-pos3[3])^2)
            dist32=sqrt((pos2[1]-pos3[1])^2+(pos2[2]-pos3[2])^2+(pos2[3]-pos3[3])^2)
            if (dist31*1.5e11)>(rad[1]+rad[3]) && (dist32*1.5e11)>(rad[2]+rad[3]) #3rd isn't touching either
                accept3=true
                return pos3
            end
        end
    end
    pos3=genPos3(pos1,pos2)
    pos=[pos1[1],pos1[2],pos1[3],pos2[1],pos2[2],pos2[3],pos3[1],pos3[2],pos3[3]].*1.5e11 #convert accepted positions to SI, m
    v=rand(-7e3:7e3,9) #random xyz velocities with mag between -10 & 10 km/s, totally arbitrary...
    #r=[x1,y1,x2,y2,x3,y3,v1x,v1y,v2x,v2y,v3x,v3y]
    r=[pos[1],pos[2],pos[3],pos[4],pos[5],pos[6],pos[7],pos[8],pos[9],v[1],v[2],v[3],v[4],v[5],v[6],v[7],v[8],v[9]]
    open("initCond.txt","w") do f #save initial conditions to file in folder where script is run
        write(f,"m1=$(@sprintf("%.1f",(m[1]/2e30))) m2=$(@sprintf("%.1f",(m[2]/2e30))) m3=$(@sprintf("%.1f",(m[3]/2e30))) (solar masses)\nv1x=$(v[1]/1e3) v1y=$(v[2]/1e3) v1z=$(v[3]/1e3) v2x=$(v[4]/1e3) v2y=$(v[5]/1e3) v2z=$(v[6]/1e3) v3x=$(v[7]/1e3) v3y=$(v[8]/1e3) v3z=$(v[9]/1e3) (km/s)\nx1=$(pos1[1]) y1=$(pos1[2]) z1=$(pos1[3]) x2=$(pos2[1]) y2=$(pos2[2]) z2=$(pos2[3]) x3=$(pos3[1]) y3=$(pos3[2]) z3=$(pos3[3]) (AU from center)")
    end
    return r, rad, m
end

function dR(r,m) #function we will use RK4 on to approximate solution
    G=6.67408313131313e-11# Nm^2/kg^2
    M1,M2,M3=m[1],m[2],m[3] #kg
    x1,x2,x3=r[1],r[4],r[7] #m
    y1,y2,y3=r[2],r[5],r[8] #m
    z1,z2,z3=r[3],r[6],r[9]
    c1,c2,c3=G*M1,G*M2,G*M3 #Nm^2/kg
    r1_2=sqrt((x1-x2)^2+(y1-y2)^2+(z1-z2)^2) #distance from 1->2, m
    r1_3=sqrt((x1-x3)^2+(y1-y3)^2+(z1-z3)^2) #distance from 1->3, m
    r2_3=sqrt((x2-x3)^2+(y2-y3)^2+(z2-z3)^2) #distance from 2->3, m

    v1X,v2X,v3X=r[10],r[13],r[16] #these are our change in position after dt (dr/dt*dt=dr)
    v1Y,v2Y,v3Y=r[11],r[14],r[17] #m after * dt
    v1Z,v2Z,v3Z=r[12],r[15],r[18]

    #get change in velocity from accelerations (d^2r/dt^2*dt=dv/dt*dt=dv)
    dx1=-(c2*(x1-x2)/(r1_2^3))-(c3*(x1-x3)/(r1_3^3)) #d^2x/dt^2 for 1, m/s after * dt
    dx2=-(c1*(x2-x1)/(r1_2^3))-(c3*(x2-x3)/(r2_3^3)) #d^2x/dt^2 for 2, m/s
    dx3=-(c1*(x3-x1)/(r1_3^3))-(c2*(x3-x2)/(r2_3^3)) #d^2x/dt^2 for 3, m/s
    dy1=-(c2*(y1-y2)/(r1_2^3))-(c3*(y1-y3)/(r1_3^3)) #d^2y/dt^2 for 1, m/s
    dy2=-(c1*(y2-y1)/(r1_2^3))-(c3*(y2-y3)/(r2_3^3)) #d^2y/dt^2 for 2, m/s
    dy3=-(c1*(y3-y1)/(r1_3^3))-(c2*(y3-y2)/(r2_3^3)) #d^2y/dt^2 for 3, m/s
    dz1=-(c2*(z1-z2)/(r1_2^3))-(c3*(z1-z3)/(r1_3^3)) #d^2y/dt^2 for 1, m/s
    dz2=-(c1*(z2-z1)/(r1_2^3))-(c3*(z2-z3)/(r2_3^3)) #d^2y/dt^2 for 2, m/s
    dz3=-(c1*(z3-z1)/(r1_3^3))-(c2*(z3-z2)/(r2_3^3)) #d^2y/dt^2 for 3, m/s

    return [v1X,v1Y,v1Z,v2X,v2Y,v2Z,v3X,v3Y,v3Z,dx1,dy1,dz1,dx2,dy2,dz2,dx3,dy3,dz3]
end

function gen3Body(stopCond=[10,100],numSteps=10000) #default stop conditions of 10 yrs and 100 AU sep
    tStop=stopCond[1]*365*24*3600 #convert to SI s
    sepStop=stopCond[2]*1.5e11 #convert to SI m
    stop=false
    currentT=0
    t=range(0,stop=tStop,length=(numSteps+1)) #+1 because I don't want 0 to count
    stepSize=tStop/numSteps
    x1=zeros(length(t))
    y1=zeros(length(t))
    z1=zeros(length(t))
    x2=zeros(length(t))
    y2=zeros(length(t))
    z2=zeros(length(t))
    x3=zeros(length(t))
    y3=zeros(length(t))
    z3=zeros(length(t))
    r,rad,m=initCondGen()
    min12=rad[1]+rad[2]
    min13=rad[1]+rad[3]
    min23=rad[2]+rad[3]
    i=1
    stopT=maximum(t)
    collisionBool=false
    #implement RK4 to model solutions to differential equations
    while stop==false
        if currentT==stopT || currentT>stopT #in case of rounding error or something
            stop=true
        elseif i>numSteps+1 #inf loop failsafe
            stop=true
            println("error: shouldn't have gotten here")
        else
            x1[i]=r[1]
            y1[i]=r[2]
            z1[i]=r[3]
            x2[i]=r[4]
            y2[i]=r[5]
            z2[i]=r[6]
            x3[i]=r[7]
            y3[i]=r[8]
            z3[i]=r[9]

            k1=stepSize*dR(r,m)
            k2=stepSize*dR(r.+0.5.*k1,m)
            k3=stepSize*dR(r.+0.5.*k2,m)
            k4=stepSize*dR(r.+k3,m)
            r+=(k1.+2.0*k2.+2.0.*k3.+k4)./6

            #check separation after each dt step
            sep12=sqrt((x1[i]-x2[i])^2+(y1[i]-y2[i])^2+(z1[i]-z2[i])^2)
            sep13=sqrt((x1[i]-x3[i])^2+(y1[i]-y3[i])^2+(z1[i]-z3[i])^2)
            sep23=sqrt((x3[i]-x2[i])^2+(y3[i]-y2[i])^2+(z3[i]-z2[i])^2)

            if sep12<min12 || sep13<min13 || sep23<min23 || sep12>sepStop || sep13>sepStop || sep23>sepStop
                if sep12<min12 || sep13<min13 || sep23<min23
                    collisionBool=true
                else
                    collisionBool=false
                end
                stop=true #stop if collision happens or body is ejected
                t=range(0,stop=currentT,length=i) #t should match pos vectors
                x1=x1[1:i] #don't want trailing zeros
                y1=y1[1:i]
                z1=z1[1:i]
                x2=x2[1:i]
                y2=y2[1:i]
                z2=z2[1:i]
                x3=x3[1:i]
                y3=y3[1:i]
                z3=z3[1:i]
            end
            i+=1
            currentT+=stepSize #next step
        end
    end
    return [x1,y1,z1,x2,y2,z2,x3,y3,z3], t, m, rad, collisionBool
end

function getInteresting3Body(minTime=0) #in years, defaults to 0
    #sometimes random conditions result in a really short animation where things
    #just crash into each other/fly away, so this function throws away those
    yearSec=365*24*3600
    interesting=false
    i=1
    while interesting==false
        plotData,t,m,rad,collisionBool=gen3Body([60,150],600000)
        if (maximum(t)/yearSec)>minTime #only return if simulation runs for longer than minTime
            println(maximum(t)/yearSec) #tell me how many years we are simulating
            open("cron_log.txt","a") do f #for cron logging, a flag = append
                write(f,"$(maximum(t)/yearSec)\n")
            end
            return plotData,t,m,rad,collisionBool
            interesting=true
        elseif i>5000 #computationally expensive so don't want to go forever
            interesting=true #render it anyways I guess because sometimes it's fun?
            println("did not find interesting solution in number of tries allotted, running anyways")
            println(maximum(t)/yearSec) #how many years simulation runs for
            open("cron_log.txt","a") do f #for cron logging
                write(f,"$(maximum(t)/yearSec)\n")
            end
            return plotData,t,m,rad,collisionBool
        end
        i+=1
    end
end

function getLims(pos,padding) #determines plot limits at each frame, padding in units of pos
    x=[pos[1],pos[4],pos[7]]
    xMin=minimum(x)
    xMax=maximum(x)
    dx=xMax-xMin
    y=[pos[2],pos[5],pos[8]]
    yMin=minimum(y)
    yMax=maximum(y)
    dy=yMax-yMin
    z=[pos[3],pos[6],pos[9]]
    zMin=minimum(z)
    zMax=maximum(z)
    dz=zMax-zMin
    dList=[dx,dy,dz]
    if maximum(dList)==dx
        #use x for square
        xlims=[xMin-padding,xMax+padding]
        ylims=[yMin-padding,yMin+dx+padding]
        zlims=[zMin-padding,zMin+dx+padding]
    elseif maximum(dList)==dy
        #use y for square
        xlims=[xMin-padding,xMin+dy+padding]
        ylims=[yMin-padding,yMax+padding]
        zlims=[zMin-padding,zMin+dy+padding]
    else
        #use z for cube
        xlims=[xMin-padding,xMin+dz+padding]
        ylims=[yMin-padding,yMin+dz+padding]
        zlims=[zMin-padding,zMax+padding]
    end
    return xlims,ylims,zlims
end

function getColors(m,c) #places colors of objects according to mass/size
    #c=[:biggest,:medium,:smallest] (order of input colors)
    maxM=maximum(m)
    minM=minimum(m)
    colors=[:blue,:blue,:blue] #testing
    if m[1]==maxM
        colors[1]=c[1]
        if m[2]==minM
            colors[2]=c[3]
            colors[3]=c[2]
        else
            colors[3]=c[3]
            colors[2]=c[2]
        end
    elseif m[2]==maxM
        colors[2]=c[1]
        if m[1]==minM
            colors[1]=c[3]
            colors[3]=c[2]
        else
            colors[3]=c[3]
            colors[1]=c[2]
        end
    else
        colors[3]=c[1]
        if m[1]==minM
            colors[1]=c[3]
            colors[2]=c[2]
        else
            colors[2]=c[3]
            colors[1]=c[2]
        end
    end
    return colors
end

function makeCircleVals(r,center=[0,0])
    xOffset=center[1]
    yOffset=center[2]
    xVals=[r*cos(i)+xOffset for i=0:(pi/64):(2*pi)]
    yVals=[r*sin(i)+yOffset for i=0:(pi/64):(2*pi)]
    return xVals,yVals
end

plotData,t,m,rad,collisionBool=getInteresting3Body(15)
c=[:DodgerBlue,:Gold,:Tomato] #most massive to least massive, also roughly corresponds to temp
colors=getColors(m,c)
#adding fake stars
numStars=2500
starsX=zeros(numStars)
starsY=zeros(numStars)
starsZ=zeros(numStars)
for i=1:numStars
    num=rand(-200:200,3) #box size is 70 AU but we need some extra padding for movement
    starsX[i]=num[1]
    starsY[i]=num[2]
    starsZ[i]=num[3]
end
stars = [starsX,starsY,starsZ]
global frameNum=1
stop=length(t)
if collisionBool==true
    stop=length(t)-600
end

function makePanel(i1,i2,i,rad,plotData,stars,x,y,lims)
    p=plot(plotData[i1][1:33:i]./1.5e11,plotData[i2][1:33:i]./1.5e11,label="",linewidth=2,linecolor=colors[1],linealpha=max.((1:33:i) .+ 10000 .- i,2500)/10000) #plot orbits up to i
    p=plot!(plotData[i1+3][1:33:i]./1.5e11,plotData[i2+3][1:33:i]./1.5e11,label="",linewidth=2,linecolor=colors[2],linealpha=max.((1:33:i) .+ 10000 .- i,2500)/10000) #linealpha argument causes lines to decay
    p=plot!(plotData[i1+6][1:33:i]./1.5e11,plotData[i2+6][1:33:i]./1.5e11,label="",linewidth=2,linecolor=colors[3],linealpha=max.((1:33:i) .+ 10000 .- i,2500)/10000) #example: alpha=max.((1:i) .+ 100 .- i,0) causes only last 100 to be visible
    p=scatter!(stars[i1],stars[i2],markercolor=:white,markersize=:1,label="") #fake background stars
    star1=makeCircleVals(rad[1],[plotData[i1][i],plotData[i2][i]]) #generate circles with appropriate sizes for each star
    star2=makeCircleVals(rad[2],[plotData[i1+3][i],plotData[i2+3][i]]) #at current positions
    star3=makeCircleVals(rad[3],[plotData[i1+6][i],plotData[i2+6][i]])
    p=plot!(star1[1]./1.5e11,star1[2]./1.5e11,label="$(@sprintf("%.1f", m[1]./2e30))",color=colors[1],fill=true) #plot star circles with labels
    p=plot!(star2[1]./1.5e11,star2[2]./1.5e11,label="$(@sprintf("%.1f", m[2]./2e30))",color=colors[2],fill=true)
    p=plot!(star3[1]./1.5e11,star3[2]./1.5e11,label="$(@sprintf("%.1f", m[3]./2e30))",color=colors[3],fill=true)
    p=plot!(background_color=:black,background_color_legend=:transparent,foreground_color_legend=:transparent,
        background_color_outside=:white,aspect_ratio=:equal,legendtitlefontcolor=:white,legendfontfamily="Courier") #formatting for plot frame
    p=plot!(xlabel="$x: AU",ylabel="$y: AU",title="$x$y plane",
        legend=:best,xaxis=("$x: AU",(lims[i1][1],lims[i1][2]),font(9,"Courier")),yaxis=("$y: AU",(lims[i2][1],lims[i2][2]),font(9,"Courier")),
        grid=false,titlefont=font(14,"Courier"),size=(720,721),legendfontsize=8,legendtitle="Mass (in solar masses)",legendtitlefontsize=8,legendtitlefont="Courier") #add in axes/title/legend with formatting
    return p
end

for i=1:333:stop #this makes animation scale ~1 sec/year with other conditions
    GR.inline("png") #added to eneable cron/jobber compatibility, also this makes frames generate WAY faster? Prior to adding this when run from cron/jobber frames would stop generating at 408 for some reason.
    gr(legendfontcolor = plot_color(:white)) #legendfontcolor=:white plot arg broken right now (at least in this backend)
    print("$(@sprintf("%.2f",i/length(t)*100)) % complete\r") #output percent tracker
    pos=[plotData[1][i],plotData[2][i],plotData[3][i],plotData[4][i],plotData[5][i],plotData[6][i],plotData[7][i],plotData[8][i],plotData[9][i]] #current pos
    lims=getLims(pos./1.5e11,15) #convert to AU, 10 AU padding
    pXY=makePanel(1,2,i,rad,plotData,stars,"x","y",lims)
    pXZ=makePanel(1,3,i,rad,plotData,stars,"x","z",lims)
    pYZ=makePanel(2,3,i,rad,plotData,stars,"y","z",lims)
    title = plot(title = "Random Three-Body Problem: t = $(@sprintf("%0.2f",t[i]/365/24/3600)) years after start", grid = false, showaxis = false,titlefontcolor=:black,titlefont=(font(18,"Courier")),
    background_color=:white,foreground_color=:white,background_color_outside=:white,ticks=false,background_color_legend=:white,foreground_color_legend=:white,background_color_subplot=:white,background_color_inside=:white,foreground_color_subplot=:white)
    P = plot(title,pXY,pXZ,pYZ,layout=@layout([A{0.05h}; [B C D]]),size=(720*3,720),left_margin=5mm,bottom_margin=5mm)
    png(P,@sprintf("tmpPlots/frame_%06d.png",frameNum))
    global frameNum+=1
    closeall() #close plots
end
if collisionBool==true #this condition makes 2 seconds of slo-mo right before the collision
    println("making collision cam")
    for i=1:10:600
        GR.inline("png") #added to eneable cron/jobber compatibility, also this makes frames generate WAY faster? Prior to adding this when run from cron/jobber frames would stop generating at 408 for some reason.
        gr(legendfontcolor = plot_color(:white)) #legendfontcolor=:white plot arg broken right now (at least in this backend)
        print("$(@sprintf("%.2f",i/600*100)) % complete\r") #output percent tracker
        I = length(plotData[1])-(600-i)
        pos=[plotData[1][end-(600-i)],plotData[2][end-(600-i)],plotData[3][end-(600-i)],plotData[4][end-(600-i)],plotData[5][end-(600-i)],plotData[6][end-(600-i)],plotData[7][end-(600-i)],plotData[8][end-(600-i)],plotData[9][end-(600-i)]] #current pos
        lims=getLims(pos./1.5e11,15) #convert to AU, 15 AU padding
        pXY=makePanel(1,2,I,rad,plotData,stars,"x","y",lims)
        pXZ=makePanel(1,3,I,rad,plotData,stars,"x","z",lims)
        pYZ=makePanel(2,3,I,rad,plotData,stars,"y","z",lims)#add in axes/title/legend with formatting
        title = plot(title = "Random Three-Body Problem: t = $(@sprintf("%0.2f",t[i]/365/24/3600)) years after start", grid = false, showaxis = false,titlefontcolor=:black,titlefont=(font(18,"Courier")),
        background_color=:white,foreground_color=:white,background_color_outside=:white,ticks=false,background_color_legend=:white,foreground_color_legend=:white,background_color_subplot=:white,background_color_inside=:white,foreground_color_subplot=:white)
        P = plot(title,pXY,pXZ,pYZ,layout=@layout([A{0.05h}; [B C D]]),size=(720*3,720),left_margin=5mm,bottom_margin=5mm)
        P=annotate!((lims[1][1]+(lims[1][2]-lims[1][1])/40,lims[2][2]-(lims[2][2]-lims[2][1])/40,Plots.text("COLLISION CAM (slo-mo x 33)",10,"Courier",:orange,:left,:bold)))
        png(P,@sprintf("tmpPlots/frame_%06d.png",frameNum))
        global frameNum+=1
        closeall() #close plots
    end
end

run(`ffmpeg -y -framerate 30 -i "tmpPlots/frame_%06d.png" -c:v libx264 -preset slow -coder 1 -movflags +faststart -g 15 -crf 18 -pix_fmt yuv420p -profile:v high -y -bf 2 "/home/kirk/Documents/3Body/3Body3Frame_fps30.mp4"`)


text:
struct UnionFind
    parent::Vector{Int} # Parent index of each index, when -1, it is root
    size::Vector{Int} # group size of each index
    UnionFind(N) = new(fill(-1, N), fill(1, N))
end

function root!(uf::UnionFind, x::Int)
    if uf.parent[x] == -1
        return x
    else
        # Path compuression
        # Set the parent of the values traced from x as the root of x.
        return uf.parent[x] = root!(uf, uf.parent[x])
    end
end

function unite!(uf::UnionFind, x::Int, y::Int)
    x_root = root!(uf, x)
    y_root = root!(uf, y)
    x_root == y_root && return false  # these has already same root
    if uf.size[x_root] < uf.size[y_root] # union by size
        uf.parent[x_root] = y_root
        uf.size[y_root] += uf.size[x_root]
        uf.size[x_root] = uf.size[y_root]
    else
        uf.parent[y_root] = x_root
        uf.size[x_root] += uf.size[y_root]
        uf.size[y_root] = uf.size[x_root]
    end
    return true
end

function members(uf::UnionFind, x::Int)
    x_root = root!(uf, x)
    return [x for x in 1:length(uf.parent) if root!(uf, x) == x_root]
end

function all_group_members(uf::UnionFind)
    group_members = Dict{Int, Vector{Int}}()
    for x in 1:length(uf.parent)
        x_root = root!(uf, x)
        haskey(group_members, x_root) || (group_members[x_root] = Int[])
        push!(group_members[x_root], x)
    end
    return group_members
end

roots(uf::UnionFind)=[x for x in 1:length(uf.parent) if uf.parent[x] == -1]
issame(uf::UnionFind, x::Int, y::Int)=root!(uf, x) == root!(uf, y)
size(uf::UnionFind, x::Int)=uf.size[x]
group_size(uf::UnionFind)=length(roots(uf))


text:
using Documenter, HigherPrecision

makedocs(
    format = :html,
    sitename = "HigherPrecision.jl",
    #modules = [HigherPrecision],
    pages = [
        "Introduction" => "index.md",
        "Compute mode" => "computemode.md",
        "Reference" => "reference.md"
        ]
)

deploydocs(
    repo   = "github.com/saschatimme/HigherPrecision.jl.git",
    target = "build",
    julia = "0.6",
    osname = "linux",
    deps   = nothing,
    make   = nothing
)


text:
@defcomp co2cycle begin
    e_globalCO2emissions=Parameter(index=[time],unit="Mtonne/year")
    e0_globalCO2emissions=Parameter(unit="Mtonne/year", default=38191.0315797948)
    c_CO2concentration=Variable(index=[time],unit="ppbv")
    pic_preindustconcCO2=Parameter(unit="ppbv", default=278000.)
    exc_excessconcCO2=Variable(unit="ppbv")
    c0_CO2concbaseyr=Parameter(unit="ppbv", default=395000.)
    re_remainCO2=Variable(index=[time],unit="Mtonne")
    re_remainCO2base=Variable(unit="Mtonne")
    renoccf_remainCO2wocc=Variable(index=[time],unit="Mtonne")
    air_CO2fractioninatm=Parameter(unit="%", default=62.00)
    stay_fractionCO2emissionsinatm=Parameter(default=0.3)
    tea_CO2emissionstoatm=Variable(index=[time],unit="Mtonne/year")
    teay_CO2emissionstoatm=Variable(index=[time],unit="Mtonne/t")
    ccf_CO2feedback=Parameter(unit="%/degreeC", default=9.66666666666667)
    ccfmax_maxCO2feedback=Parameter(unit="%", default=53.3333333333333)
    cea_cumCO2emissionsatm=Variable(index=[time],unit="Mtonne")
    ce_0_basecumCO2emissions=Parameter(unit="Mtonne", default=2050000.)
    y_year=Parameter(index=[time],unit="year")
    y_year_0=Parameter(unit="year")
    res_CO2atmlifetime=Parameter(unit="year", default=73.3333333333333)
    den_CO2density=Parameter(unit="Mtonne/ppbv", default=7.8)
    rt_g0_baseglobaltemp=Parameter(unit="degreeC", default=0.735309967925382)
    rt_g_globaltemperature=Parameter(index=[time],unit="degreeC")

    function run_timestep(p, v, d, t)

        if is_first(t)
            #CO2 emissions gain calculated based on PAGE 2009
            gain=p.ccf_CO2feedback*p.rt_g0_baseglobaltemp
            #eq.6 from Hope (2006) - emissions to atmosphere depend on the sum of natural and anthropogenic emissions
            tea0=p.e0_globalCO2emissions*p.air_CO2fractioninatm/100
            v.tea_CO2emissionstoatm[t]=(p.e_globalCO2emissions[t])*p.air_CO2fractioninatm/100
            v.teay_CO2emissionstoatm[t]=(v.tea_CO2emissionstoatm[t]+tea0)*(p.y_year[t]-p.y_year_0)/2
            #adapted from eq.1 in Hope(2006) - calculate excess concentration in base year
            v.exc_excessconcCO2=p.c0_CO2concbaseyr-p.pic_preindustconcCO2
            #Eq. 2 from Hope (2006) - base-year remaining emissions
            v.re_remainCO2base=v.exc_excessconcCO2*p.den_CO2density
            #PAGE 2009 initial remaining emissions without CO2 feedback
            renoccf0_remainCO2wocc=v.re_remainCO2base/(1+gain/100)
            #eq. 8 from Hope (2006) - baseline cumulative emissions to atmosphere
            ceabase=p.ce_0_basecumCO2emissions*p.air_CO2fractioninatm/100
            #eq.9 from Hope(2006) - cumulative emissions in atmosphere
            v.cea_cumCO2emissionsatm[t]=ceabase+v.teay_CO2emissionstoatm[t]
            #eq.11 from Hope (2006) - anthropogenic remaining emissions
            v.renoccf_remainCO2wocc[t]=p.stay_fractionCO2emissionsinatm*ceabase*
                (1-exp(-(p.y_year[t]-p.y_year_0)/
                p.res_CO2atmlifetime))+renoccf0_remainCO2wocc*
                exp(-(p.y_year[t]-p.y_year_0)/p.res_CO2atmlifetime)+
                v.teay_CO2emissionstoatm[t]*exp(-(p.y_year[t]-p.y_year_0)/
                (2*p.res_CO2atmlifetime))
            #Hope 2009 - remaining emissions with CO2 feedback
            v.re_remainCO2[t]=v.renoccf_remainCO2wocc[t]*(1+gain/100)
        else
            #CO2 emissions gain calculated based on PAGE 2009
            gain=min(p.ccf_CO2feedback*p.rt_g_globaltemperature[t-1],p.ccfmax_maxCO2feedback)
            #eq.6 from Hope (2006) - emissions to atmosphere depend on the sum of natural and anthropogenic emissions
            v.tea_CO2emissionstoatm[t]=(p.e_globalCO2emissions[t])*p.air_CO2fractioninatm/100
            #eq.7 from Hope (2006) - total emissions over time period
            v.teay_CO2emissionstoatm[t]=(v.tea_CO2emissionstoatm[t]+v.tea_CO2emissionstoatm[t-1])*
                (p.y_year[t]-p.y_year[t-1])/2
            #eq.9 from Hope(2006) - cumulative emissions in atmosphere
            v.cea_cumCO2emissionsatm[t]=v.cea_cumCO2emissionsatm[t-1]+v.teay_CO2emissionstoatm[t]
            #eq.11 from Hope (2006) - anthropogenic remaining emissions
            v.renoccf_remainCO2wocc[t]=p.stay_fractionCO2emissionsinatm*v.cea_cumCO2emissionsatm[t-1]*
                (1-exp(-(p.y_year[t]-p.y_year[t-1])/
                p.res_CO2atmlifetime))+v.renoccf_remainCO2wocc[t-1]*
                exp(-(p.y_year[t]-p.y_year[t-1])/p.res_CO2atmlifetime)+
                v.teay_CO2emissionstoatm[t]*exp(-(p.y_year[t]-p.y_year[t-1])/
                (2*p.res_CO2atmlifetime))
            #Hope 2009 - remaining emissions with CO2 feedback
            v.re_remainCO2[t]=v.renoccf_remainCO2wocc[t]*(1+gain/100)
        end
        #eq.11 from Hope(2006) - CO2 concentration
        v.c_CO2concentration[t]=p.pic_preindustconcCO2+v.exc_excessconcCO2 * v.re_remainCO2[t]/v.re_remainCO2base
    end
end


text:
include("../src/TeukEvolution.jl")

import .TeukEvolution as TE

@time TE.launch(ARGS[1])