text:
atomcoords(atoms) = coordarray(atoms) |> transpose |> collect
atomcolors(atoms, color = elecolors) = [color[BioStructures.element(x)] for x in atoms]
atomradii(atoms) = [vdwrad[BioStructures.element(x)] for x in atoms]
resids(residues) = resid.(residues)
function resatoms(res; typ = :OrderedDict)
	resvec1 = [k for k in res]
	if typ in (:OrderedDict, :ordereddict, :odict, OrderedDict)
		resvec2 = [(resvec1[i].name, coords(resvec1[i])) for i in 1:size(resvec1,1)]
		return OrderedDict(resvec2)
	elseif typ in (:Dict, :dict, Dict)
		resvec2 = [(resvec1[i].name, coords(resvec1[i])) for i in 1:size(resvec1,1)]
		return Dict(resvec2)
	elseif typ in (:ComponentArray, :comp, ComponentArray)
		resvec2 = [(resvec1[i].name, coords(resvec1[i])) for i in 1:size(resvec1,1)]
		resvec3 = dict2ntuple(Dict(resvec2))
		return ComponentArray(resvec3)
	else
		return [[resvec1[i].name, coords(resvec1[i])] for i in 1:size(resvec1,1)] |> combinedims |> _t
	end
	return res
end
"""
    viewstruc(args)

Create and return a Makie Figure for a PDB structure.
# Examples
```julia
using BioStructures
struc = retrievepdb("2vb1", dir = "data\\") |> Node
sv = viewstruc(struc)

struc = read("data\\2vb1_mutant1.pdb", BioStructures.PDB) |> Node
sv = viewstruc(struc)
```
Keyword arguments:
resolution ---- Default - (800,800)
atmcolors ----- Default - "element", define your own dict for atoms like: "N" => :blue
atmscale ------ Default - 1/3, size adjustment of atom radii
figdims ------- Default - [13,10], GridLayout dimensions for structure view: 'fig = Figure(); fig[ 1:(figdims[1]), 1:(figdims[2]) ]'
"""
function viewstruc( struc::T,
					selectors = [standardselector];
					resolution = (800,800),
					atmcolors = "element",
					atmscale = 1/3,
					figdims = [13,10]
					) where {T<:Node}
    atms = @lift BioStructures.collectatoms($struc,selectors...)
    atmcords = @lift atomcoords($atms)
    colr = @lift atomcolors($atms; color = atmcolors)
    marksize = @lift (atmscale).*atomradii($atms)
    fig = Figure(resolution = resolution)
    ly = fig[ 1:(figdims[1]), 1:(figdims[2]) ]
    plt = meshscatter(ly, atmcords; show_axis = false, color = colr, markersize = marksize)
	resshps = @lift bondshape(SplitApplyCombine.flatten(bonds(collectresidues($struc,selectors...))))
	bbshps = @lift bondshape(SplitApplyCombine.flatten(backbonebonds.(collectchains($struc))))
	resbnds = @lift normal_mesh.($resshps)
	bckbnds = @lift normal_mesh.($bbshps)
	mesh!(ly, resbnds, color = RGBAf0(0.5,0.5,0.5,0.8))
	mesh!(ly, bckbnds, color = RGBAf0(0.5,0.5,0.5,0.8))
    fig
end


text:
# Install RMS v0.3.2
import Pkg
Pkg.add(Pkg.PackageSpec(url="https://github.com/ReactionMechanismGenerator/ReactionMechanismSimulator.jl", rev="v0.3.2"))

ext: jl


text:
# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
# Since these FMAs can increase the performance of many numerical algorithms,
# we need to opt-in explicitly.
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
@muladd begin

# Composite type that represents a NDIMS-dimensional tree (parallel version).
#
# Implements everything required for AbstractContainer.
#
# Note: The way the data structures are set up and the way most algorithms
# work, it is *always* assumed that
#   a) we have a balanced tree (= at most one level difference between
#                                 neighboring cells, or 2:1 rule)
#   b) we may not have all children (= some children may not exist)
#   c) the tree is stored depth-first
#
# However, the way the refinement/coarsening algorithms are currently
# implemented, we only have fully refined cells. That is, a cell either has 2^NDIMS children or
# no children at all (= leaf cell). This restriction is also assumed at
# multiple positions in the refinement/coarsening algorithms.
#
# An exception to the 2:1 rule exists for the low-level `refine_unbalanced!`
# function, which is required for implementing level-wise refinement in a sane
# way. Also, depth-first ordering *might* not by guaranteed during
# refinement/coarsening operations.
mutable struct ParallelTree{NDIMS} <: AbstractTree{NDIMS}
  parent_ids::Vector{Int}
  child_ids::Matrix{Int}
  neighbor_ids::Matrix{Int}
  levels::Vector{Int}
  coordinates::Matrix{Float64}
  original_cell_ids::Vector{Int}
  mpi_ranks::Vector{Int}

  capacity::Int
  length::Int
  dummy::Int

  center_level_0::SVector{NDIMS, Float64}
  length_level_0::Float64
  periodicity::NTuple{NDIMS, Bool}

  function ParallelTree{NDIMS}(capacity::Integer) where NDIMS
    # Verify that NDIMS is an integer
    @assert NDIMS isa Integer

    # Create instance
    t = new()

    # Initialize fields with defaults
    # Note: length as capacity + 1 is to use `capacity + 1` as temporary storage for swap operations
    t.parent_ids = fill(typemin(Int), capacity + 1)
    t.child_ids = fill(typemin(Int), 2^NDIMS, capacity + 1)
    t.neighbor_ids = fill(typemin(Int), 2*NDIMS, capacity + 1)
    t.levels = fill(typemin(Int), capacity + 1)
    t.coordinates = fill(NaN, NDIMS, capacity + 1)
    t.original_cell_ids = fill(typemin(Int), capacity + 1)
    t.mpi_ranks = fill(typemin(Int), capacity + 1)

    t.capacity = capacity
    t.length = 0
    t.dummy = capacity + 1

    t.center_level_0 = SVector(ntuple(_ -> NaN, NDIMS))
    t.length_level_0 = NaN

    return t
  end
end

# Constructor for passing the dimension as an argument
ParallelTree(::Val{NDIMS}, args...) where NDIMS = ParallelTree{NDIMS}(args...)

# Create and initialize tree
function ParallelTree{NDIMS}(capacity::Int, center::AbstractArray{Float64},
                 length::Real, periodicity=true) where NDIMS
  # Create instance
  t = ParallelTree{NDIMS}(capacity)

  # Initialize root cell
  init!(t, center, length, periodicity)

  return t
end

# Constructor accepting a single number as center (as opposed to an array) for 1D
ParallelTree{1}(cap::Int, center::Real, len::Real, periodicity=true) = ParallelTree{1}(cap, [convert(Float64, center)], len, periodicity)

# Clear tree with deleting data structures, store center and length, and create root cell
function init!(t::ParallelTree, center::AbstractArray{Float64}, length::Real, periodicity=true)
  clear!(t)

  # Set domain information
  t.center_level_0 = center
  t.length_level_0 = length

  # Create root cell
  t.length += 1
  t.parent_ids[1] = 0
  t.child_ids[:, 1] .= 0
  t.levels[1] = 0
  set_cell_coordinates!(t, t.center_level_0, 1)
  t.original_cell_ids[1] = 0
  t.mpi_ranks[1] = typemin(Int)

  # Set neighbor ids: for each periodic direction, the level-0 cell is its own neighbor
  if all(periodicity)
    # Also catches case where periodicity = true
    t.neighbor_ids[:, 1] .= 1
    t.periodicity = ntuple(x->true, ndims(t))
  elseif !any(periodicity)
    # Also catches case where periodicity = false
    t.neighbor_ids[:, 1] .= 0
    t.periodicity = ntuple(x->false, ndims(t))
  else
    # Default case if periodicity is an iterable
    for dimension in 1:ndims(t)
      if periodicity[dimension]
        t.neighbor_ids[2 * dimension - 1, 1] = 1
        t.neighbor_ids[2 * dimension - 0, 1] = 1
      else
        t.neighbor_ids[2 * dimension - 1, 1] = 0
        t.neighbor_ids[2 * dimension - 0, 1] = 0
      end
    end

    t.periodicity = Tuple(periodicity)
  end
end

# Convenience output for debugging
function Base.show(io::IO, ::MIME"text/plain", t::ParallelTree)
  @nospecialize t # reduce precompilation time

  l = t.length
  println(io, '*'^20)
  println(io, "t.parent_ids[1:l] = $(t.parent_ids[1:l])")
  println(io, "transpose(t.child_ids[:, 1:l]) = $(transpose(t.child_ids[:, 1:l]))")
  println(io, "transpose(t.neighbor_ids[:, 1:l]) = $(transpose(t.neighbor_ids[:, 1:l]))")
  println(io, "t.levels[1:l] = $(t.levels[1:l])")
  println(io, "transpose(t.coordinates[:, 1:l]) = $(transpose(t.coordinates[:, 1:l]))")
  println(io, "t.original_cell_ids[1:l] = $(t.original_cell_ids[1:l])")
  println(io, "t.mpi_ranks[1:l] = $(t.mpi_ranks[1:l])")
  println(io, "t.capacity = $(t.capacity)")
  println(io, "t.length = $(t.length)")
  println(io, "t.dummy = $(t.dummy)")
  println(io, "t.center_level_0 = $(t.center_level_0)")
  println(io, "t.length_level_0 = $(t.length_level_0)")
  println(io, '*'^20)
end

# Check if cell is own cell, i.e., belongs to this MPI rank
is_own_cell(t::ParallelTree, cell_id) = t.mpi_ranks[cell_id] == mpi_rank()

# Return an array with the ids of all leaf cells for a given rank
leaf_cells_by_rank(t::ParallelTree, rank) = filter_leaf_cells(t) do cell_id
                                              t.mpi_ranks[cell_id] == rank
                                            end

# Return an array with the ids of all local leaf cells
local_leaf_cells(t::ParallelTree) = leaf_cells_by_rank(t, mpi_rank())

# Set information for child cell `child_id` based on parent cell `cell_id` (except neighbors)
function init_child!(t::ParallelTree, cell_id, child, child_id)
  t.parent_ids[child_id] = cell_id
  t.child_ids[child, cell_id] = child_id
  t.child_ids[:, child_id] .= 0
  t.levels[child_id] = t.levels[cell_id] + 1
  set_cell_coordinates!(t,
    child_coordinates(t, cell_coordinates(t, cell_id), length_at_cell(t, cell_id), child), child_id)
  t.original_cell_ids[child_id] = 0
  t.mpi_ranks[child_id] = t.mpi_ranks[cell_id]

  return nothing
end

# Reset range of cells to values that are prone to cause errors as soon as they are used.
#
# Rationale: If an invalid cell is accidentally used, we want to know it as soon as possible.
function invalidate!(t::ParallelTree, first::Int, last::Int)
  @assert first > 0
  @assert last <= t.capacity + 1

  # Integer values are set to smallest negative value, floating point values to NaN
  t.parent_ids[first:last] .= typemin(Int)
  t.child_ids[:, first:last] .= typemin(Int)
  t.neighbor_ids[:, first:last] .= typemin(Int)
  t.levels[first:last] .= typemin(Int)
  t.coordinates[:, first:last] .= NaN
  t.original_cell_ids[first:last] .= typemin(Int)
  t.mpi_ranks[first:last] .= typemin(Int)

  return nothing
end

# Raw copy operation for ranges of cells.
#
# This method is used by the higher-level copy operations for AbstractContainer
function raw_copy!(target::ParallelTree, source::ParallelTree, first::Int, last::Int, destination::Int)
  copy_data!(target.parent_ids, source.parent_ids, first, last, destination)
  copy_data!(target.child_ids, source.child_ids, first, last, destination,
             n_children_per_cell(target))
  copy_data!(target.neighbor_ids, source.neighbor_ids, first, last,
             destination, n_directions(target))
  copy_data!(target.levels, source.levels, first, last, destination)
  copy_data!(target.coordinates, source.coordinates, first, last, destination, ndims(target))
  copy_data!(target.original_cell_ids, source.original_cell_ids, first, last, destination)
  copy_data!(target.mpi_ranks, source.mpi_ranks, first, last, destination)
end

# Reset data structures by recreating all internal storage containers and invalidating all elements
function reset_data_structures!(t::ParallelTree{NDIMS}) where NDIMS
  t.parent_ids = Vector{Int}(undef, t.capacity + 1)
  t.child_ids = Matrix{Int}(undef, 2^NDIMS, t.capacity + 1)
  t.neighbor_ids = Matrix{Int}(undef, 2*NDIMS, t.capacity + 1)
  t.levels = Vector{Int}(undef, t.capacity + 1)
  t.coordinates = Matrix{Float64}(undef, NDIMS, t.capacity + 1)
  t.original_cell_ids = Vector{Int}(undef, t.capacity + 1)
  t.mpi_ranks = Vector{Int}(undef, t.capacity + 1)

  invalidate!(t, 1, capacity(t) + 1)
end

end # @muladd


text:
function find_passage(r)
    side=Ost
    while isborder(r,Nord)==true # прохода сверху нет
        putmarker!(r)
        move_by_markers!(r,side)
        side=inverse(side)
    end
end

move_by_markers(r::Robot,side::HorizonSide) =
    while ismarker(r)==true
        move!(r,side)
    end

inverse(side::HorizonSide) = HorizonSide(mod(Int(side)+2,4))


text:
module XLStats

using OrderedCollections
using Statistics
using Parameters
using Printf

export Link, getscore, read_all, filter, point_biserial, indomain

export name, resnames, indexes, ismatch, getlink,
       consistency, deuc, dtop, dmax, nscans,
       maxscore1, avgscore1,
       maxscore2, avgscore2,
       hasxic, maxxic, avgxic,
       nspecies, count_nspecies

#
# Convert one letter and three-letter amino acid codes
#
const oneletter = Dict("CYS"=>"C", "ASP"=>"D", "SER"=>"S", "GLN"=>"Q", "LYS"=>"K",
                       "ILE"=>"I", "PRO"=>"P", "THR"=>"T", "PHE"=>"F", "ASN"=>"N",
                       "GLY"=>"G", "HIS"=>"H", "LEU"=>"L", "ARG"=>"R", "TRP"=>"W",
                       "ALA"=>"A", "VAL"=>"V", "GLU"=>"E", "TYR"=>"Y", "MET"=>"M")
const threeletter = Dict( v => k for (k,v) in oneletter )

#
# List of possible data of the scans and the associated symbols in the Link struct
# needs to be ordered because it will define the substitution sequence when parsing
# the Scan line, and score1 and score2 conflict.
#
const data_list = OrderedDict( "Scan:"                => :index,
                               "Peptide:"             => :peptide,
                               "Secondary Score:"     => :score2,
                               "Score:"               => :score1,
                               "XL Position 1:"       => :pos1,
                               "XL Position 2:"       => :pos2,
                               "Source File:"         => :source_file,
                               "Experimental M+H:"    => :mplush,
                               "Precursor Charge:"    => :prec_charge,
                               "Peaks Matched alpha:" => :matched_alpha,
                               "Peaks Matched beta:"  => :matched_beta,
                               "Retention Time:"      => :retention_time,
                               "Assessment:"          => :assessment )

#
# Structure to contain scan data
#
@with_kw mutable struct Scan
  index :: Int = 0
  pep1 :: String = ""
  pep2 :: String = ""
  score1 :: Float64 = 0.
  score2 :: Float64 = 0.
  pos1 :: Int = 0
  pos2 :: Int = 0
  source_file :: String = ""
  mplush :: Float64 = 0
  prec_charge :: Int = 0
  matched_alpha :: Int = 0
  matched_beta :: Int = 0
  retention_time :: Float64 = 0
  xic :: Float64 = -1
end

#
# Structure that contain each link data
#
@with_kw mutable struct Link

  name::AbstractString = "None"

  consistency::Bool = false
  deuc::Float64 = 0.
  dtop::Float64 = 0.
  dmax::Float64 = 0.

  nspecies :: Int = 0
  hasxic::Bool = false
  scans :: Vector{Scan}

end

function Link(name::String,nscans::Int)
  str = split(name,'-')
  i = parse(Int,str[1][2:end])
  j = parse(Int,str[2][2:end])
  scans = [ Scan() for i in 1:nscans ]
  Link(name=name, scans=scans)
end

#
# syntax sugars for getfield on array of links
#
name(link::Link) = link.name
name(links::Vector{Link}) = name.(links)

function resnames(linkname::AbstractString;type=3)
  name = split(linkname,"-")
  if type == 3
    return threeletter[name[1][1:1]], threeletter[name[2][1:1]]
  elseif type == 1
    return name[1][1:1], name[2][1:1]
  else
    error(" type must be 1 or 3 in resnames. ")
  end
end
resnames(link::Link;type=3) = resnames(link.name,type=type)

function indexes(linkname::AbstractString)
  name = split(linkname,"-")
  return parse(Int,name[1][2:end]), parse(Int,name[2][2:end])
end
indexes(link::Link) = indexes(link.name)

function ismatch(name1::AbstractString,name2::AbstractString)
   ( ismatch(resnames(name1),resnames(name2)) &&
     ismatch(indexes(name1),indexes(name2)) )
end

function ismatch(x::Tuple{T,T},y::Tuple{T,T}) where T
  ( (x[1] == y[1] && x[2] == y[2]) ||
    (x[2] == y[1] && x[1] == y[2]) )
end

#
# Allow retrieving links by name (much slower, of course)
#
import Base.getindex
getindex(links::Vector{Link},name::AbstractString) = links[findfirst( link -> ismatch(link.name,name), links )]

consistency(link::Link;tol=0.) = (0. <= link.dtop <= (link.dmax + tol))
consistency(links::Vector{Link};tol=0) = consistency.(links,tol=tol)

deuc(link::Link) = link.deuc
deuc(links::Vector{Link}) = deuc.(links)

dtop(link::Link) = link.dtop
dtop(links::Vector{Link}) = dtop.(links)

dmax(link::Link) = link.dmax
dmax(links::Vector{Link}) = dmax.(links)

nscans(link::Link) = length(link.scans)
nscans(links::Vector{Link}) = nscans.(links)

#
# Get maximum value of a field of the scans of a link
#
function getmax(link::Link,field::Symbol)
  nscans(link) == 0 && return missing
  valmax = getfield(link.scans[1],field)
  for iscan in 2:nscans(link)
     valmax = max(valmax,getfield(link.scans[iscan],field))
  end
  return valmax
end

#
# Get the mean value of a field of the scans of a link
#
function getmean(link::Link,field::Symbol)
  nscans(link) == 0 && return missing
  meanval = zero(typeof(getfield(link.scans[1],field)))
  for iscan in 2:nscans(link)
     meanval += getfield(link.scans[iscan],field)
  end
  return meanval / nscans(link)
end

maxscore1(link::Link) = getmax(link,:score1)
maxscore1(links::Vector{Link}) = maxscore1.(links)

avgscore1(link::Link) = getmean(link,:score1)
avgscore1(links::Vector{Link}) = avgscore1.(links)

maxscore2(link::Link) = getmax(link,:score2)
maxscore2(links::Vector{Link}) = maxscore2.(links)

avgscore2(link::Link) = getmean(link,:score2)
avgscore2(links::Vector{Link}) = avgscore2.(links)

nspecies(link::Link) = link.nspecies
nspecies(links::Vector{Link}) = nspecies.(links)

# Not all data has xic values, so these function have to be special
hasxic(link::Link) = link.hasxic
hasxic(links::Vector{Link}) = hasxic.(links)
maxxic(link::Link) = getmax(link,:xic)
maxxic(links::Vector{Link}) = maxxic.(links)
function avgxic(links::Vector{Link})
  nxic = count(link -> link.hasxic, links)
  meanxic = zeros(nxic)
  ixic = 0
  for link in links
    if link.hasxic
      ixic += 1
      n = 0
      a  = 0.
      for scan in link.scans
        if scan.xic > 0
          n += 1
          a += scan.xic
        end
      end
      meanxic[ixic] = a/n
    end
  end
  return meanxic
end

"""

```
count_nspecies(link::Link;mtol=1.0)
```

Count the number of independent species identified in `link`, with a mass difference tolerance
`mtol`, which defaults to `1.0`. This default value is stored in the `link.nspecies` field.

## Example:

```julia-repl
julia> links
 Vector of Links with: 102 links.

julia> count_nspecies(links[1],mtol=1)
7

```

"""
function count_nspecies(link::Link;mtol=1.0)
  nscans(link) == 0 && return 0
  unique = Float64[link.scans[1].mplush]
  for i in 2:nscans(link)
    if findfirst( m -> abs(link.scans[i].mplush-m) < mtol, unique) == nothing
      push!(unique,link.scans[i].mplush)
    end
  end
  return length(unique)
end

"""

```
read_all(xml_file::String=nothing,
         topolink_log = nothing,
         topolink_input = nothing,
         xic_file_name = nothing,
         domain::UnitRange = -1000:1000)
```

Reads the `xml_file` SimXL file, `topolink` log file, `topolink` input file and,
optionally, a file containing `xic` data for each scan, and returns a vector of `Link`
structures containing all data. Optionally, one can set the `domain` to be considered,
consisting of a range of residue numbers.

## Example:

```julia_repl

julia> using XLStats

julia> links = read_all(xml_file="./data/salbiii_hitsDetail.dat",
                        topolink_log="./data/salbiii_topolink.log",
                        topolink_input="./data/topolink.inp",
                        #xic_file_name="./data/salbiii_xic.dat",
                        domain=2:134)
 Reading all data files...
 Reading XML file ...
 Reading XIC data...
 Vector of Links with: 102 links.

```

"""
function read_all(;xml_file = nothing,
                   topolink_log = nothing,
                   topolink_input = nothing,
                   xic_file_name = nothing,
                   domain::UnitRange = -1000:1000)

  println(" Reading all data files... ")

  error = false
  if xml_file == nothing
    println("ERROR: You need to provide the xml_file.")
    error = true
  end
  if topolink_log == nothing
    println("ERROR: You need to provide the topolink log file.")
    error = true
  end
  if topolink_input == nothing
    println("ERROR: You need to provide the topolink input file containing linktype information.")
    error = true
  end

  if error
    return nothing
  end

  # Read xml file from SIM-XL
  links = read_xml(xml_file,xic_file_name,domain)

  # Order links by residue index
  sort!(links,lt=isless_indexes)

  # Read topolink input to get the length of the linkers
  for link in links
    link.dmax = getdmax(link,topolink_input)
  end

  # Read topolink log file to get the euclidean and topological distances
  for link in links
    link.deuc, link.dtop = readlog(link,topolink_log)
  end

  # Set consistency with 0. tolerance:
  for link in links
    link.consistency = consistency(link,tol=0.)
  end

  return links
end
isless_indexes(link1,link2) = ( indexes(link1) < indexes(link2) )

#
# Function that checks if the link belongs to the chosen domain
#
function indomain(name,domain)
  name == "None" && return false
  r = split(name,'-')
  # Remove wrong assignments for which the two residues are the same
  r[1] == r[2] && return false
  # Check if pair is in domain
  i = parse(Int,r[1][2:end])
  j = parse(Int,r[2][2:end])
  if ( i >= domain[begin] &&
       i <= domain[end] ) &&
     ( j >= domain[begin] &&
       j <= domain[end] )
    return true
  else
    return false
  end
end
indomain(link::Link,domain) = indomain(link.name,domain)

#
# Function that reads the XIC file and fills the xic fields
#
function readxic!(links,xic_file_name)

  # Not every xic is provided, -1 means it was not
  for link in links
     link.hasxic = false
     for i in 1:nscans(link)
       link.scans[i].xic = -1
     end
  end

  xic_file_name == nothing && return
  xic_file = open(xic_file_name,"r")
  iread = 0
  xic_in_next_line = false
  for line in eachline(xic_file)
    data = split(line)
    if xic_in_next_line
      for link in links, i in 1:nscans(link)
        if iread == link.scans[i].index
          xic = parse(Float64,data[end])
          if xic > 0
            link.hasxic = true
            link.scans[i].xic = xic
          end
          xic_in_next_line = false
          break
        end
      end
    end
    if can_be_parsed_as(Int,data[1])
      iread = parse(Int,data[1])
      xic_in_next_line = true
    end
  end
  close(xic_file)

  nothing
end

function can_be_parsed_as(T::DataType,s::AbstractString)
  try
    parse(T,s)
    return true
  catch
    return false
  end
end

#
# Function that reads a topolink link log and sets link data
#
function readlog(link,logfile_name)

  file = open(logfile_name,"r")
  for line in eachline(file)
    comment(line) && continue
    if split(strip(line))[1] == "LINK:"
      residue1_name = oneletter[strip(line[8:12])]
      residue1_index = strip(line[15:19])
      residue2_name = oneletter[strip(line[25:29])]
      residue2_index = strip(line[32:36])
      name = residue1_name*residue1_index*'-'*residue2_name*residue2_index
      if ismatch(name,link.name)
        deuc = parse(Float64,line[43:51])
        dtop = parse(Float64,replace(line[52:61],">"=>""))
        close(file)
        return deuc, dtop
      end
    end
  end
  close(file)
  println(" Warning: Could not find LINK data in TopoLink log for link: $(link.name)")

  return -1., -1.
end

#
# Function that reads the linker length from a topolink input file with linktype
# definitions
#
function getdmax(link,linktype_file)
  file = open(linktype_file,"r")
  for line in eachline(file)
    comment(line) && continue
    if occursin("linktype",line)
      data = split(line)
      residue1_name = oneletter[data[2]]
      residue2_name = oneletter[data[6]]
      dmax = parse(Float64,data[10])
      name = split(link.name,'-')
      if ( name[1][1:1] == residue1_name && name[2][1:1] == residue2_name ) ||
         ( name[1][1:1] == residue2_name && name[2][1:1] == residue1_name )
        close(file)
        return dmax
      end
    end
  end
  close(file)
  return -1.
end

#
# Print link in the topolink format
#
function write(link;tol=0.)
  data = split(link.name,"-")
  res1_name = threeletter[data[1][1:1]]
  res2_name = threeletter[data[2][1:1]]
  res1_index = parse(Int,data[1][2:end])
  res2_index = parse(Int,data[2][2:end])
  cons = consistency(link,tol=tol)
  println(@sprintf("%3s %3i %3s %3i %3i %3.2f %3.2f",
                   res1_name,res1_index,res2_name,res2_index,cons,link.dtop,link.dmax))
end

"""

```
point_biserial(x::Vector{Bool},y::Vector{<:Real})
```

Compute point-biserial correlation (x assumes 0 or 1 values, y is continuous)
x is boolean (True or False for each group)

"""
function point_biserial(x::AbstractVector{Bool},y::AbstractVector{<:Real})

  ndata = length(x)
  g1 = zeros(Bool,ndata)
  g2 = zeros(Bool,ndata)
  for i in eachindex(x)
    if x[i]
      g1[i] = true
    else
      g2[i] = true
    end
  end
  data1 = y[g1]
  data2 = y[g2]

  n1 = length(data1)
  n2 = length(data2)
  if n1 == 0 || n2 == 0
    return 0.
  end

  avg1 = mean(data1)
  avg2 = mean(data2)
  sd = std(y)

  pbs = ( ( avg1 - avg2 ) / sd ) * sqrt( n1*n2 / ndata^2 )

  return pbs
end

#
# Function that reads the xml file generated by sim-xl
#
function read_xml(data_file_name,xic_file_name,domain)

  println(" Reading XML file ... ")

  # Count the number of links listed
  data_file = open(data_file_name,"r")
  nlinks = 0
  for line in eachline(data_file)
    comment(line) && continue
    if newlink(line)
      name = setname(line)
      if indomain(name,domain)
        nlinks = nlinks + 1
      end
    end
  end
  seek(data_file,0)

  # Number of scans per link
  names = String[ "None" for i in 1:nlinks ]
  nscans = zeros(Int,nlinks)

  # Count the number of scans per link
  ilink = 0
  consider_link = false
  for line in eachline(data_file)
    comment(line) && continue
    if newlink(line)
      name = setname(line)
      consider_link = indomain(name,domain)
      if consider_link
        ilink += 1
        names[ilink] = name
      end
    end
    if consider_link && occursin("Scan",line)
      nscans[ilink] += 1
    end
  end
  seek(data_file,0)

  # Initialize the scans per link
  links = [ Link(names[i],nscans[i]) for i in 1:nlinks ]

  # Read the score parameters for each scan
  ilink = 0
  iscan = 0
  for line in eachline(data_file)
    comment(line) && continue
    if newlink(line)
      name = setname(line)
      consider_link = indomain(name,domain)
      if consider_link
        ilink += 1
        iscan = 0
      end
    end
    if consider_link && occursin("Scan",line)
      iscan = iscan + 1

      for pair in data_list
        line = replace(line,pair)
      end
      data = split(line)

      for field in values(data_list)
        ival = findfirst(isequal("$field"),data)
        if ! isnothing(ival) && length(data) > ival
          if field == :peptide
            setfield!(links[ilink].scans[iscan],:pep1,String(data[ival+1]))
            setfield!(links[ilink].scans[iscan],:pep2,String(data[ival+3]))
          else
            if fieldtype(Scan,field) == String
              setfield!(links[ilink].scans[iscan],field,String(data[ival+1]))
            else
              fieldvalue = parse(fieldtype(Scan,field),data[ival+1])
              setfield!(links[ilink].scans[iscan],field,fieldvalue)
            end
          end
        end
      end
    end
  end
  close(data_file)

  # Read XIC data
  println(" Reading XIC data... ")
  readxic!(links,xic_file_name)

  # Compute the number of species of each link
  for link in links
    link.nspecies = count_nspecies(link)
  end

  return links
end

#
# Remove a specific link from the list
#
function remove(links,name)
  i=-1
  for link in links
    i=i+1
    if link.name == name
      deleteat!(links,i)
      return links
    end
  end
  error(" Error in remove function: link $(link.name) not found in link list. ")
end

#
# Check if a line is commented
#
comment(line) = (length(strip(line)) == 0) || (strip(line)[1] == '#')

#
# check if line defines a new link
#
newlink(line) = occursin("(",line)

#
# Check if a line in a xml file is a new link line, and return
# its name if so
#
function setname(line)
  ! occursin("-",line) && return "None"
  line = replace(line," ("=>"",count=2)
  line = replace(line,") - "=>"-")
  line = replace(line,")"=>"")
  name = strip(line)
  i1, i2 = indexes(name)
  if i2 < i1
    newname = split(name,"-")
    name = newname[2]*"-"*newname[1]
  end
  return name
end

function Base.show( io :: IO,::MIME"text/plain", links::Vector{Link} )
  println(" Vector of Links with: ", length(links)," links.")
end

end # module


text:
using SubClonalSelection
using Test
using Plots
using Statistics
using Random

#useful to check what the backend is, seems to cause some problems
println(backend())
tests = ["neutral", "1clone"]

println("Running tests ...")

for t in tests
    fn = "test_$t.jl"
    println("* $fn ...")
    include(fn)
end


text:
Tables.istable(::Type{TrialPerf}) = true
Tables.columnaccess(::Type{TrialPerf}) = true
Tables.columns(perf::TrialPerf) = Tables.CopiedColumns(fulltable(perf))


text:
"""
    LaplaceLikelihood(β::T=1.0)  #  Laplace likelihood with scale β

Laplace likelihood for regression:
```math
1/(2β) exp(-|y-f|/β)
```
see [wiki page](https://en.wikipedia.org/wiki/Laplace_distribution)
---
For the analytical solution, it is augmented via:
```math
p(y|f,ω) = N(y|f,ω⁻¹)
```
where ``ω ~ Exp(ω | 1/(2 β^2))``, and `Exp` is the [Exponential distribution](https://en.wikipedia.org/wiki/Exponential_distribution)
We use the variational distribution ``q(ω) = GIG(ω | a,b,p)``
"""
mutable struct LaplaceLikelihood{T<:Real,A<:AbstractVector{T}} <: RegressionLikelihood{T}
    β::T
    a::T
    p::T
    b::A #Variational parameter b of GIG
    θ::A #Expected value of ω
    function LaplaceLikelihood{T}(β::T) where {T<:Real}
        return new{T,Vector{T}}(β, β^-2, 0.5)
    end
    function LaplaceLikelihood{T}(
        β::T, b::A, θ::AbstractVector{T}
    ) where {T<:Real,A<:AbstractVector{T}}
        return new{T,A}(β, β^(-2), 0.5, b, θ)
    end
end

function LaplaceLikelihood(β::T=1.0) where {T<:Real}
    return LaplaceLikelihood{T}(β)
end

function implemented(
    ::LaplaceLikelihood, ::Union{<:AnalyticVI,<:QuadratureVI,<:GibbsSampling}
)
    return true
end

function init_likelihood(
    likelihood::LaplaceLikelihood{T},
    inference::AbstractInference{T},
    ::Int,
    nSamplesUsed::Int,
) where {T}
    if inference isa AnalyticVI || inference isa GibbsSampling
        LaplaceLikelihood{T}(likelihood.β, rand(T, nSamplesUsed), zeros(T, nSamplesUsed))
    else
        LaplaceLikelihood{T}(likelihood.β)
    end
end

function (l::LaplaceLikelihood)(y::Real, f::Real)
    return pdf(Laplace(f, l.β), y)
end

function Distributions.loglikelihood(l::LaplaceLikelihood, y::Real, f::Real)
    return logpdf(Laplace(f, l.β), y)
end

function Base.show(io::IO, l::LaplaceLikelihood{T}) where {T}
    return print(io, "Laplace likelihood (β=$(l.β))")
end

function compute_proba(
    l::LaplaceLikelihood{T}, μ::AbstractVector{<:Real}, σ²::AbstractVector{<:Real}
) where {T<:Real}
    return μ, max.(σ², 0.0) .+ 2 * l.β^2
end

## Local Updates ##

function local_updates!(
    l::LaplaceLikelihood{T}, y::AbstractVector, μ::AbstractVector, diagΣ::AbstractVector
) where {T}
    @. l.b = diagΣ + abs2(μ - y)
    @. l.θ = sqrt(l.a) / sqrt.(l.b)
end

function sample_local!(l::LaplaceLikelihood, y::AbstractVector, f::AbstractVector)
    @. l.b = rand(GeneralizedInverseGaussian(1 / l.β^2, abs2(f - y), 0.5))
    set_ω!(l, inv.(l.b))
    return nothing
end

@inline function ∇E_μ(l::LaplaceLikelihood{T}, ::AOptimizer, y::AbstractVector) where {T}
    return (l.θ .* y,)
end
@inline function ∇E_Σ(l::LaplaceLikelihood{T}, ::AOptimizer, ::AbstractVector) where {T}
    return (0.5 * l.θ,)
end

## ELBO ##
function expec_loglikelihood(
    l::LaplaceLikelihood{T},
    ::AnalyticVI,
    y::AbstractVector,
    μ::AbstractVector,
    diag_cov::AbstractVector,
) where {T}
    tot = -0.5 * length(y) * log(twoπ)
    tot += 0.5 * Zygote.@ignore(sum(log, l.θ))
    tot +=
        -0.5 * (
            dot(l.θ, diag_cov) + dot(l.θ, abs2.(μ)) - 2.0 * dot(l.θ, μ .* y) +
            dot(l.θ, abs2.(y))
        )
    return tot
end

AugmentedKL(l::LaplaceLikelihood, ::AbstractVector) = GIGEntropy(l) - expecExponentialGIG(l)

GIGEntropy(l::LaplaceLikelihood{T}) where {T} = GIGEntropy(l.a, l.b, l.p)

function expecExponentialGIG(l::LaplaceLikelihood{T}) where {T}
    return sum(
        -log(2 * l.β^2) .-
        0.5 * (l.a .* sqrt.(l.b) + l.b .* sqrt(l.a)) ./ (l.a .* l.b * l.β^2),
    )
end

## PDF and Log PDF Gradients ##

function grad_quad(
    likelihood::LaplaceLikelihood{T},
    y::Real,
    μ::Real,
    σ²::Real,
    inference::AbstractInference,
) where {T<:Real}
    nodes = inference.nodes * sqrt(σ²) .+ μ
    Edloglike = dot(inference.weights, ∇loglikehood.(likelihood, y, nodes))
    Ed²loglike = (1 / sqrt(twoπ * σ²)) / (likelihood.β^2)
    return -Edloglike::T, Ed²loglike::T
end

@inline function ∇loglikehood(l::LaplaceLikelihood{T}, y::Real, f::Real) where {T<:Real}
    return sign(y - f) ./ l.β
end

@inline hessloglikelihood(::LaplaceLikelihood{T}, ::Real, ::Real) where {T<:Real} = zero(T)


text:
function rgammanm(z, n, m=0)
    if n==m
        return one(z)
    elseif n < m
        return one(z) / rgammanm(z, m, n)
    end

    rg = one(z)
    for i = (n-1):-1:m
        rg = rg * (z + i)
    end

    return rg
end

function rgammaxy(x, y, n)
    rg = x/y * gamma(x)/gamma(y)

    for i = 1:(n-1)
        rg = rg * (x+i)/(y+i)
    end
    return rg
end

function rgammaxnym(x, n, y, m)

    if x == y
        return rgammanm(x, n, m)
    elseif n == m
        return rgammaxy(x, y, n)
    elseif n < m
        return one(x)/rgammaxnym(y, m, x, n)
    end

    rg = x/y * gamma(x)/gamma(y)

    for i = (m-1):-1:1
        rg = rg * (x+i)/(y+i)
    end

    for i = (n-1):-1:m
        rg = rg * (x+i)
    end

    return rg
end


text:
export Clock, RT_SEC, RT_MSEC, RT_USEC, RT_NSEC, TICK, SLEEP_TIME
export nsleep, usleep, msleep, ssleep, tick, job!
export clk

mutable struct Clock <: System
  started::Base.Event
  stopped::Bool
  periodic_jobs::Vector{Tuple{float,Function,String}}
  freq::AbstractFloat
end

Clock() = Clock(Base.Event(), true, [], 0.01667) # default frequency of approximately 60 Hz

const clk = Clock()

# RT == "real time"
# Δ carries the "actual" number of given time units elapsed
struct RT_SCALE{scale, signature}
  Δ::AbstractFloat
end
const SIG_TIME_COUNT = 0
const SIG_TIME_TICK = 1
RT_SEC = RT_SCALE{1.0, SIG_TIME_COUNT}
RT_MSEC = RT_SCALE{1e-3, SIG_TIME_COUNT}
RT_USEC = RT_SCALE{1e-6, SIG_TIME_COUNT}
RT_NSEC = RT_SCALE{1e-9, SIG_TIME_COUNT}
TICK = RT_SCALE{1.0, SIG_TIME_TICK}

struct SLEEP_TIME
  Δ::UInt # time in nanoseconds to sleep for
end

function Base.sleep(s::SLEEP_TIME)
  t1 = time_ns()
  while true
    if time_ns() - t1 >= s.Δ break end
    yield()
  end
  return time_ns() - t1
end

function sleep_with_message(sleep_time,timescale, message, debug)
  δ = sleep(SLEEP_TIME(sleep_time*timescale))
  sendMessage(message(δ/timescale))
  @debug debug
end

nsleep(Δ) = sleep_with_message(Δ,1.0, RT_NSEC, "nanosecond")
usleep(Δ) = sleep_with_message(Δ,1e3, RT_USEC, "microsecond")
msleep(Δ) = sleep_with_message(Δ,1e6, RT_MSEC, "millisecond")
ssleep(Δ) = sleep_with_message(Δ,1e9, RT_SEC, "second")
tick(Δ) = sleep_with_message(Δ,1e9, TICK, "tick")

function job!(c::Clock, f, arg=1)
  function job()
    Base.wait(c.started)
    while !c.stopped
      f(arg)
    end
  end
  schedule(Task(job))
end

function awake!(c::Clock)
  for i in c.periodic_jobs
    job!(c, Δ-> sleep_with_message(Δ,i[1], i[2], i[3]))
  end

  job!(c, tick, c.freq)

  c.stopped = false

  Base.notify(c.started)

  return true
end

function shutdown!(c::Clock)
  c.stopped = true
  c.started = Base.Event() # old one remains signaled no matter what, replace
  return false
end


text:
function dp(M, i, p, memo)
    # Base case
    if i == 1
        return M[i][p]
    end

    #Memorized
    if memo[i, p] > 0
        return memo[i, p]
    end

    # Recursion
    if p == 1
        return memo[i, p] = M[i][p] + max(dp(M, i-1, 2, memo), dp(M, i-1, 3, memo))
    elseif p == 2
        return memo[i, p] = M[i][p] + max(dp(M, i-1, 1, memo), dp(M, i-1, 3, memo))
    elseif p == 3
        return memo[i, p] = M[i][p] + max(dp(M, i-1, 1, memo), dp(M, i-1, 2, memo))
    end
end

function main()
    N = parse(Int, readline())
    M = [parse.(Int, split(readline())) for _ in 1:N]
    memo = zeros(Int, N, 3)

    println(max(dp(M, N, 1, memo), dp(M, N, 2, memo), dp(M, N, 3, memo)))
end

main()


text:
import ONNX: rename_keys

@testset "Utils" begin
    attrs = Dict(:a => 0, :b => 1)
    @test rename_keys(attrs, Dict(:a => :c)) == Dict(:c => 0, :b => 1)
    @test rename_keys(attrs, Dict(:f => :g)) == attrs
    @test rename_keys(attrs, Dict()) == attrs
end


text:
# ------------------------------------------------------------------------------------------
# # Iterators
# (c) Deniz Yuret, 2019
#
# * Objective: Learning how to construct and use Julia iterators.
# * Reading: [Interfaces](https://docs.julialang.org/en/v1/manual/interfaces/#man-interface-
# iteration-1),  [Collections](https://docs.julialang.org/en/v1/base/collections/#lib-
# collections-iteration-1), [Iteration
# Utilities](https://docs.julialang.org/en/v1/base/iterators) and [Generator
# expressions](https://docs.julialang.org/en/v1/manual/arrays/#Generator-Expressions-1) in
# the Julia manual.
# * Prerequisites: [minibatch,
# Data](https://github.com/denizyuret/Knet.jl/blob/master/src/data.jl) from the [MNIST
# notebook](20.mnist.ipynb)
# * New functions:
# [first](https://docs.julialang.org/en/v1/base/collections/#Base.first),
# [collect](https://docs.julialang.org/en/v1/base/collections/#Base.collect-Tuple{Any}),
# [repeat](https://docs.julialang.org/en/v1/base/arrays/#Base.repeat),
# [take](https://docs.julialang.org/en/v1/base/iterators/#Base.Iterators.take),
# [drop](https://docs.julialang.org/en/v1/base/iterators/#Base.Iterators.drop),
# [cycle](https://docs.julialang.org/en/v1/base/iterators/#Base.Iterators.cycle),
# [Stateful](https://docs.julialang.org/en/v1/base/iterators/#Base.Iterators.Stateful),
# [iterate](https://docs.julialang.org/en/v1/base/collections/#lib-collections-iteration-1)
#
# The `minibatch` function returns a `Knet.Data` object implemented as a Julia iterator that
# generates (x,y) minibatches. Iterators are lazy objects that only generate their next
# element when asked. This has the advantage of not wasting time and memory trying to create
# and store all the elements at once. We can even have infinite iterators! The training
# algorithms in Knet are also implemented as iterators so that:
# 1. We can monitor and report the training loss
# 2. We can take snapshots of the model during training
# 3. We can pause/terminate training when necessary
#
# Here are some things Julia can do with iterators:
# ------------------------------------------------------------------------------------------

# Set display width, load packages, import symbols
ENV["COLUMNS"]=72
using Pkg; haskey(Pkg.installed(),"Knet") || Pkg.add("Knet")
using Base.Iterators: take, drop, cycle, Stateful
using Knet

# Load data
include(Knet.dir("data","mnist.jl"))
xtrn,ytrn,xtst,ytst = mnist()
dtst = minibatch(xtst,ytst,100)

# We can peek at the first element using first()
summary.(first(dtst))

# Iterators can be used in for loops
# Let's count the elements in dtst:
n = 0
for (x,y) in dtst; n += 1; end
@show n;

# Iterators can be converted to arrays using `collect`
# (don't do this unless necessary, it just wastes memory. Use a for loop instead)
collect(dtst) |> summary

# We can generate an iterator for multiple epochs using `repeat`
# (an epoch is a single pass over the dataset)
n = 0
for (x,y) in repeat(dtst,5); n += 1; end
@show n;

# We can generate partial epochs using `take` which takes the first n elements
n = 0
for (x,y) in take(dtst,20); n += 1; end
@show n;

# We can also generate partial epochs using `drop` which drops the first n elements
n = 0
for (x,y) in drop(dtst,20); n += 1; end
@show n;

# We can repeat forever using `cycle` (this is useful to train until convergence)
# You do not want to collect a cycle or run a for loop without break!
n = 0; t = time_ns()
for (x,y) in cycle(dtst)
    n += 1
    time_ns() - t > 2e9 && break # Break after 2 seconds
end
@show n;

# We can make an iterator `Stateful` so it remembers where it left off.
# (by default iterators start from the beginning)
dtst1 = dtst            # dtst1 will start from beginning every time
dtst2 = Stateful(dtst)  # dtst2 will remember where we left off
for (x,y) in dtst1; println(Int.(y[1:5])); break; end
for (x,y) in dtst1; println(Int.(y[1:5])); break; end
for (x,y) in dtst2; println(Int.(y[1:5])); break; end
for (x,y) in dtst2; println(Int.(y[1:5])); break; end

# We can shuffle instances at every epoch using the keyword argument `shuffle=true`
# (by default elements are generated in the same order)
dtst1 = minibatch(xtst,ytst,100)              # dtst1 iterates in the same order
dtst2 = minibatch(xtst,ytst,100,shuffle=true) # dtst2 shuffles each time
for (x,y) in dtst1; println(Int.(y[1:5])); break; end
for (x,y) in dtst1; println(Int.(y[1:5])); break; end
for (x,y) in dtst2; println(Int.(y[1:5])); break; end
for (x,y) in dtst2; println(Int.(y[1:5])); break; end

# We can construct new iterators using [Generator expressions](https://docs.julialang.org/en/v1/manual/arrays/#Generator-Expressions-1)
# The following example constructs an iterator over the x norms in a dataset:
xnorm(data) = (sum(abs2,x) for (x,y) in data)
collect(xnorm(dtst))'

# Every iterator implements the `iterate` function which returns
# the next element and state (or nothing if no elements left).
# Here is how the for loop for dtst is implemented:
n = 0; next = iterate(dtst)
while next != nothing
    ((x,y), state) = next
    n += 1
    next = iterate(dtst,state)
end
@show n;

# You can define your own iterator by declaring a new type and overriding the `iterate` method.
# Here is another way to define an iterator over the x norms in a dataset:
struct Xnorm; itr; end

function Base.iterate(f::Xnorm, s...)
    next = iterate(f.itr, s...)
    next === nothing && return nothing
    ((x,y),state) = next
    return sum(abs2,x), state
end

Base.length(f::Xnorm) = length(f.itr) # collect needs this

collect(Xnorm(dtst))'

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
"""
Used to represent arguments to non-terminal nodes.
"""
type KozaNonTerminalArgument
  """
  The label used by this non-terminal.
  """
  label::AbstractString

  """
  The (name of the) type used by this non-terminal.
  """
  ty::AbstractString

  KozaNonTerminalArgument(def::AbstractString) = new(def[1:Base.search(def, ':')-1],
    def[Base.last(Base.search(def, "::"))+1:end])
end

abstract KozaNonTerminal <: KozaNode

function compose_non_terminal(inputs::Vector{KozaInput}, def::AbstractString)
  label = def[1:Base.search(def, '(')-1]
  returns = Base.strip(def[last(Base.search(def, ")::"))+1:Base.search(def, '=')-1])
  body = Base.strip(def[Base.search(def, '=')+1:end])
  accepts = [KozaNonTerminalArgument(Base.strip(ss))
    for ss in Base.split(def[Base.search(def, '(')+1:Base.search(def, ')')-1], ',')]

  # Compose the unique type for this non-terminal.
  src = "type <: KozaNonTerminal
    parent::KozaParent
    children::Vector{KozaNode}

    constructor() = new()
    constructor(p::KozaParent) =
      new(p, Array(KozaNode, $(length(accepts))))
    constructor($(
      join(vcat(["p::KozaParent"], ["c$(i)::KozaNode" for i in 1:length(inputs)]), ",")
    )) =
      new($(
      join(vcat(["p"], ["c$(i)" for i in 1:length(inputs)]), ",")
    ))
  end"
  t = anonymous_type(koza, src)

  # Compose the "fresh" function.
  eval(koza, Base.parse("fresh(::$(t), p::KozaParent) = $(t)(p)"))

  # Compose the label function.
  eval(koza, Base.parse("label(::$(t)) = $(label)"))

  # Compose the arity function.
  eval(koza, Base.parse("arity(::$(t)) = $(length(accepts))"))

  # Compose the execution function.
  input_list = [i.label for i in inputs]
  src = join(vcat(["n::$(t)"], ["$(i.label)::$(i.ty)" for i in inputs]), ",")
  src = "function execute($(src))\n"

  # Find each argument in the function body and replace with an execution call.
  for (i, a) in enumerate(accepts)
    body = replace(body, Regex("\\b$(a.label)\\b"),
      "execute($(join(vcat(["n.children[$(i)]"], input_list), ", ")))")
  end
  src = src * body * "\nend"
  eval(koza, Base.parse(src))

  # Compose the clone function.
  src = "clone(n::$(t), p::KozaParent) = $(t)($(
    join(vcat(["p"], ["clone(n.children[$(i)], n)" for i in 1:length(inputs)]), ",")
  ))"
  eval(koza, Base.parse(src))
  return t
end

isterminal(::KozaNonTerminal) = false
execute(::KozaNonTerminal) = error("No execution method for this non-terminal.")
arity(::KozaNonTerminal) = error("No arity method for this non-terminal.")
#num_nodes(n::KozaNonTerminal) = 1 + sum(Int[num_nodes(c) for c in n.children])
#height(n::KozaNonTerminal) = 1 + maximum(Int[height(c) for c in n.children])
describe(n::KozaNonTerminal) = "$(label(n))($(join([describe(n) for c in n.children], ",")))"

function num_nodes(n::KozaNonTerminal)
  i = 0
  for c in n.children
    i += num_nodes(c)
  end
  return i + 1
end

function height(n::KozaNonTerminal)
  i = 0
  for c in n.children
    i = max(i, height(c))
  end
  return i + 1
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

module TestBunchKaufman

using Test, LinearAlgebra, Random
using LinearAlgebra: BlasComplex, BlasFloat, BlasReal, QRPivoted
using Base: getproperty

n = 10

# Split n into 2 parts for tests needing two matrices
n1 = div(n, 2)
n2 = 2*n1

Random.seed!(12343210)

areal = randn(n,n)/2
aimg  = randn(n,n)/2
a2real = randn(n,n)/2
a2img  = randn(n,n)/2
breal = randn(n,2)/2
bimg  = randn(n,2)/2

@testset "$eltya argument A" for eltya in (Float32, Float64, ComplexF32, ComplexF64, Int)
    a = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex.(areal, aimg) : areal)
    a2 = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex.(a2real, a2img) : a2real)
    asym = transpose(a) + a                  # symmetric indefinite
    aher = a' + a                  # Hermitian indefinite
    apd  = a' * a                  # Positive-definite
    for (a, a2, aher, apd) in ((a, a2, aher, apd),
                               (view(a, 1:n, 1:n),
                                view(a2, 1:n, 1:n),
                                view(aher, 1:n, 1:n),
                                view(apd , 1:n, 1:n)))
        ε = εa = eps(abs(float(one(eltya))))

        # check that factorize gives a Bunch-Kaufman
        @test isa(factorize(asym), LinearAlgebra.BunchKaufman)
        @test isa(factorize(aher), LinearAlgebra.BunchKaufman)
        @testset "$uplo Bunch-Kaufman factor of indefinite matrix" for uplo in (:L, :U)
            bc1 = bunchkaufman(Hermitian(aher, uplo))
            @test LinearAlgebra.issuccess(bc1)
            @test logabsdet(bc1)[1] ≈ log(abs(det(bc1)))
            if eltya <: Real
                @test logabsdet(bc1)[2] == sign(det(bc1))
            else
                @test logabsdet(bc1)[2] ≈ sign(det(bc1))
            end
            @test inv(bc1)*aher ≈ Matrix(I, n, n)
            @testset for rook in (false, true)
                @test inv(bunchkaufman(Symmetric(transpose(a) + a, uplo), rook))*(transpose(a) + a) ≈ Matrix(I, n, n)
                if eltya <: BlasFloat
                    # test also bunchkaufman! without explicit type tag
                    # no bunchkaufman! method for Int ... yet
                    @test inv(bunchkaufman!(transpose(a) + a, rook))*(transpose(a) + a) ≈ Matrix(I, n, n)
                end
                @test size(bc1) == size(bc1.LD)
                @test size(bc1, 1) == size(bc1.LD, 1)
                @test size(bc1, 2) == size(bc1.LD, 2)
                if eltya <: BlasReal
                    @test_throws ArgumentError bunchkaufman(a)
                end
                # Test extraction of factors
                if eltya <: Real
                    @test getproperty(bc1, uplo)*bc1.D*getproperty(bc1, uplo)' ≈ aher[bc1.p, bc1.p]
                    @test getproperty(bc1, uplo)*bc1.D*getproperty(bc1, uplo)' ≈ bc1.P*aher*bc1.P'
                end

                bc1 = bunchkaufman(Symmetric(asym, uplo))
                @test getproperty(bc1, uplo)*bc1.D*transpose(getproperty(bc1, uplo)) ≈ asym[bc1.p, bc1.p]
                @test getproperty(bc1, uplo)*bc1.D*transpose(getproperty(bc1, uplo)) ≈ bc1.P*asym*transpose(bc1.P)
                @test_throws ErrorException bc1.Z
                @test_throws ArgumentError uplo == :L ? bc1.U : bc1.L
            end
            # test Base.iterate
            ref_objs = (bc1.D, uplo == :L ? bc1.L : bc1.U, bc1.p)
            for (bki, bkobj) in enumerate(bc1)
                @test bkobj == ref_objs[bki]
            end
            if eltya <: BlasFloat
                @test convert(LinearAlgebra.BunchKaufman{eltya}, bc1) === bc1
                @test convert(LinearAlgebra.Factorization{eltya}, bc1) === bc1
                if eltya <: BlasReal
                    @test convert(LinearAlgebra.Factorization{Float16}, bc1) == convert(LinearAlgebra.BunchKaufman{Float16}, bc1)
                elseif eltya <: BlasComplex
                    @test convert(LinearAlgebra.Factorization{ComplexF16}, bc1) == convert(LinearAlgebra.BunchKaufman{ComplexF16}, bc1)
                end
            end
            @test Base.propertynames(bc1) == (:p, :P, :L, :U, :D)
        end

        @testset "$eltyb argument B" for eltyb in (Float32, Float64, ComplexF32, ComplexF64, Int)
            b = eltyb == Int ? rand(1:5, n, 2) : convert(Matrix{eltyb}, eltyb <: Complex ? complex.(breal, bimg) : breal)
            for b in (b, view(b, 1:n, 1:2))
                εb = eps(abs(float(one(eltyb))))
                ε = max(εa,εb)

                @testset "$uplo Bunch-Kaufman factor of indefinite matrix" for uplo in (:L, :U)
                    bc1 = bunchkaufman(Hermitian(aher, uplo))
                    @test aher*(bc1\b) ≈ b atol=1000ε
                end

                @testset "$uplo Bunch-Kaufman factors of a pos-def matrix" for uplo in (:U, :L)
                    @testset "rook pivoting: $rook" for rook in (false, true)
                        bc2 = bunchkaufman(Hermitian(apd, uplo), rook)
                        @test LinearAlgebra.issuccess(bc2)
                        bks = split(sprint(show, "text/plain", bc2), "\n")
                        @test bks[1] == summary(bc2)
                        @test bks[2] == "D factor:"
                        @test bks[4+n] == "$uplo factor:"
                        @test bks[6+2n] == "permutation:"
                        @test logdet(bc2) ≈ log(det(bc2))
                        @test logabsdet(bc2)[1] ≈ log(abs(det(bc2)))
                        @test logabsdet(bc2)[2] == sign(det(bc2))
                        @test inv(bc2)*apd ≈ Matrix(I, n, n)
                        @test apd*(bc2\b) ≈ b rtol=eps(cond(apd))
                        @test ishermitian(bc2) == !issymmetric(bc2)
                    end
                end
            end
        end
    end
end

@testset "Singular matrices" begin
    R = Float64[1 0; 0 0]
    C = ComplexF64[1 0; 0 0]
    for A in (R, Symmetric(R), C, Hermitian(C))
        @test_throws SingularException bunchkaufman(A)
        @test_throws SingularException bunchkaufman!(copy(A))
        @test_throws SingularException bunchkaufman(A; check = true)
        @test_throws SingularException bunchkaufman!(copy(A); check = true)
        @test !issuccess(bunchkaufman(A; check = false))
        @test !issuccess(bunchkaufman!(copy(A); check = false))
    end
    F = bunchkaufman(R; check = false)
    @test sprint(show, "text/plain", F) == "Failed factorization of type $(typeof(F))"
end

@testset "test example due to @timholy in PR 15354" begin
    A = rand(6,5); A = complex(A'*A) # to avoid calling the real-lhs-complex-rhs method
    F = cholesky(A);
    v6 = rand(ComplexF64, 6)
    v5 = view(v6, 1:5)
    @test F\v5 == F\v6[1:5]
end

@testset "issue #32080" begin
    A = Symmetric([-5 -9 9; -9 4 1; 9 1 2])
    B = bunchkaufman(A, true)
    @test B.U * B.D * B.U' ≈ A[B.p, B.p]
end

@test_throws DomainError logdet(bunchkaufman([-1 -1; -1 1]))
@test logabsdet(bunchkaufman([8 4; 4 2]; check = false))[1] == -Inf

@testset "0x0 matrix" begin
    for ul in (:U, :L)
        B = bunchkaufman(Symmetric(ones(0, 0), ul))
        @test isa(B, BunchKaufman)
        @test B.D == Tridiagonal([], [], [])
        @test B.P == ones(0, 0)
        @test B.p == []
        if ul == :U
            @test B.U == UnitUpperTriangular(ones(0, 0))
            @test_throws ArgumentError B.L
        else
            @test B.L == UnitLowerTriangular(ones(0, 0))
            @test_throws ArgumentError B.U
        end
    end
end

end # module TestBunchKaufman

max_stars_repo_name: TimoLarson/julia

max_issues_repo_name: TimoLarson/julia

max_forks_repo_name: TimoLarson/julia


text:
toml_str(a; kwargs...) = sprint(io -> TOML.print(io, a; kwargs...))
toml_str(f, a; kwargs...) = sprint(io -> TOML.print(f, io, a; kwargs...))

@test toml_str(Dict("b" => 1, "c" => 2, "a" => 3); sorted=true) ==
    """
    a = 3
    b = 1
    c = 2
    """

@test toml_str(Dict("b" => 1, "ac" => 2, "abc" => 3); sorted=true, by=length) ==
    """
    b = 1
    ac = 2
    abc = 3
    """

struct MyStruct
    a::Int
end
@test_throws ErrorException toml_str(Dict("foo" => MyStruct(1)))
# simple value
@test toml_str(Dict("foo" => MyStruct(1))) do x
        x isa MyStruct && return x.a
    end == """
        foo = 1
        """

# tabular values
@test toml_str(Dict("foo" => MyStruct(1)); sorted=true) do x
         x isa MyStruct && return [x.a]
     end == """
         foo = [1]
         """
@test toml_str(Dict("foo" => MyStruct(1)); sorted=true) do x
        x isa MyStruct && return Dict(:bar => x.a)
    end == """
        [foo]
        bar = 1
        """

# validation against the usual case
@test toml_str(Dict("foo" => MyStruct(1)); sorted=true) do x
         x isa MyStruct && return [x.a]
     end == toml_str(Dict("foo" => [1]); sorted=true)
@test toml_str(Dict("foo" => MyStruct(1)); sorted=true) do x
        x isa MyStruct && return Dict(:bar => x.a)
    end == toml_str(Dict("foo" => Dict(:bar => 1)); sorted=true)

@test toml_str(Dict("b" => SubString("foo"))) == "b = \"foo\"\n"

@testset "empty dict print" begin
    s = """
    user = "me"
    [julia]
    [option]
    """
    d = TOML.parse(s)
    @test toml_str(d) == "user = \"me\"\n\n[julia]\n\n[option]\n"
end

@testset "special characters" begin
    s = """
    "\U1f355 \0 \x0 \x1 \t \b" = "\U1f355 \0 \x0 \x1 \t \b"
    "\x7f" = "\x7f"
    """
    @test roundtrip(s)

    d = Dict("str" => string(Char(0xd800)))
    @test_throws ErrorException TOML.print(d)
end


text:
# Contingency Tables
# ==================

"""
A `ContingencyTable` is a multidimensional array. It stores the contingency matrix, its
marginal values and total. The type also has an internal and private temporal array and an
alphabet object. It's a parametric type, taking three ordered parameters:

- `T` : The element type of the multidimensional array.
- `N` : It's the dimension of the array and should be an `Int`.
- `A` : This should be a type, subtype of `ResidueAlphabet`, i.e.: `UngappedAlphabet`,
`GappedAlphabet` or `ReducedAlphabet`.

A `ContingencyTable` can be created from an alphabet if all the parameters are given.
Otherwise, you need to give a type, a number (`Val`) and an alphabet. You can also create a
`ContingencyTable` using a matrix and a alphabet. For example:

```julia
ContingencyTable{Float64, 2, UngappedAlphabet}(UngappedAlphabet())
ContingencyTable(Float64, Val{2}, UngappedAlphabet())
ContingencyTable(zeros(Float64,20,20), UngappedAlphabet())
```
"""
mutable struct ContingencyTable{T,N,A} <: AbstractArray{T,N}
    alphabet::A
    temporal::Array{T,N}
    table::NamedArray{T,N,Array{T,N},NTuple{N,OrderedDict{String,Int}}}
    marginals::NamedArray{T,2,Array{T,2},NTuple{2,OrderedDict{String,Int}}}
    total::T
end

# Probability and Counts
# ----------------------

"""
A `Probabilities` object wraps a `ContingencyTable` storing probabilities. It doesn't
perform any check. If the total isn't one, you must use `normalize` or `normalize!`on the
`ContingencyTable` before wrapping it to make the sum of the probabilities equal to one.
"""
mutable struct Probabilities{T,N,A} <: AbstractArray{T,N}
    table::ContingencyTable{T,N,A}
end

"""
A `Counts` object wraps a `ContingencyTable` storing counts/frequencies.
"""
mutable struct Counts{T,N,A} <: AbstractArray{T,N}
    table::ContingencyTable{T,N,A}
end

# Getters

"""
`getcontingencytable` allows to access the wrapped `ContingencyTable` in a `Probabilities`
or `Counts` object.
"""
@inline getcontingencytable(p::Probabilities{T,N,A}) where {T,N,A} = p.table
@inline getcontingencytable(n::Counts{T,N,A}) where {T,N,A} = n.table

for f in (:getalphabet, :gettable, :getmarginals, :gettotal,
          :gettablearray, :getmarginalsarray)
    @eval $(f)(p::Probabilities{T,N,A}) where {T,N,A} = $(f)(getcontingencytable(p))
    @eval $(f)(n::Counts{T,N,A}) where {T,N,A} = $(f)(getcontingencytable(n))
end

# AbstractArray

for f in (:size, :getindex, :setindex!)
    @eval Base.$(f)(p::Probabilities{T,N,A},args...) where {T,N,A} = $(f)(getcontingencytable(p),args...)
    @eval Base.$(f)(n::Counts{T,N,A},args...) where {T,N,A} = $(f)(getcontingencytable(n),args...)
end

# Getters
# -------

"`getalphabet` allows to access the stored alphabet object."
@inline getalphabet(table::ContingencyTable) = table.alphabet

"`gettable` allows to access the table (`NamedArray`)."
@inline gettable(table::ContingencyTable) = table.table

"`getmarginals` allows to access the array with the marginal values (`NamedArray`)."
@inline getmarginals(table::ContingencyTable) = table.marginals

"`gettotal` allows to access the stored total value."
@inline gettotal(table::ContingencyTable) = table.total

Base.sum(table::ContingencyTable) = gettotal(table)

"`gettablearray` allows to access the table (`Array` without names)."
@inline gettablearray(table::ContingencyTable) = getarray(table.table)

"""
`getmarginalsarray` allows to access the array with the marginal values
(`Array` without names).
"""
@inline getmarginalsarray(table::ContingencyTable) = getarray(table.marginals)

# Cartesian (helper functions)
# ----------------------------

"""
`_marginal(1,:A,:i,:value)` generates the expression: `A[i_1, 1] += value`
"""
function _marginal(N::Int, marginal::Symbol, index::Symbol, value::Symbol)
    aexprs = [Expr(:escape, Expr(:(+=), Expr(:ref, marginal, Symbol(index,'_',i), i), value)) for i = 1:N]
    Expr(:block, aexprs...)
end

macro _marginal(N, marginal, index, value)
    _marginal(N, marginal, index, value)
end

"""
`_test_index(1, i, continue)` generates the expression: `i_1 >= 22 && continue`
"""
function _test_index(N::Int, index::Symbol, expr::Expr)
    aexprs = Expr[]
    for i = 1:N
        push!(aexprs,
            Expr(:escape,
                Expr(:&&,
                    Expr(:call, :>=, Symbol(index, "_", i), 22),
                    expr
                )
            )
        )
    end
    Expr(:block, aexprs...)
end

macro _test_index(N, index, expr)
    _test_index(N, index, expr)
end

# AbstractArray Interface
# -----------------------

Base.size(table::ContingencyTable) = size(table.table)

Base.getindex(table::ContingencyTable, i...) = getindex(table.table, i...)

@generated function Base.getindex(table::ContingencyTable, I::Residue...)
    N = length(I)
    quote
        alphabet = getalphabet(table)
        matrix = getarray(gettable(table))
        # index_1 = alphabet[I[1]]
        # index_2 ...
        @nextract $N index d->alphabet[I[d]]
        # index_1 >= 22 error("There is a Residue outside the alphabet")
        # index_2 ...
        @_test_index $N index error("There is a Residue outside the alphabet")
        # getindex(matrix, index_1, index_2...
        @nref $N matrix index
    end
end

function Base.setindex!(table::ContingencyTable, value, i...)
    setindex!(table.table, value, i...)
    update_marginals!(table)
end

@generated function Base.setindex!(table::ContingencyTable, value, I::Residue...)
    N = length(I)
    quote
        alphabet = getalphabet(table)
        matrix = getarray(gettable(table))
        @nextract $N index d->alphabet[I[d]]
        @_test_index $N index error("There is a Residue outside the alphabet")
        @nref($N, matrix, index) = value
        update_marginals!(table)
    end
end

# Similar
# -------

Base.similar(table::ContingencyTable{T,N,A}) where {T,N,A} = ContingencyTable(T, Val{N}, table.alphabet)

function Base.similar(table::ContingencyTable{T,N,A}, ::Type{S}) where {T,S,N,A}
    ContingencyTable(S, Val{N}, table.alphabet)
end

# Show
# ====

function Base.show(io::IO, ::MIME"text/plain",
                   table::Union{ContingencyTable{T,N,A},Probabilities{T,N,A},Counts{T,N,A}}) where {T,N,A}
    show(io, table)
end

function Base.show(io::IO, table::ContingencyTable{T,N,A}) where {T,N,A}
    println(io, typeof(table), " : ")
    print(io, "\ntable : ")
    show(io, gettable(table))
    if length(size(table)) != 1
        print(io, "\n\nmarginals : ")
        show(io, getmarginals(table))
    end
    print(io, "\n\ntotal : $(gettotal(table))")
end

function Base.show(io::IO, table::Union{Probabilities{T,N,A},Counts{T,N,A}}) where {T,N,A}
    print(io, typeof(table), " wrapping a ")
    show(io, getcontingencytable(table))
end

# Creation
# --------

@generated function (::Type{ContingencyTable{T,N,A}})(alphabet::A) where {T,N,A}
    @assert N > 0 "The dimension should be a natural number"
    quote
        n = length(alphabet)
        # residue_names = names(alphabet)
        # namedict = OrderedDict{String,Int}(residue_names[i] => i for i in 1:n)
        namedict = getnamedict(alphabet)
        dimtable = @ntuple $N k -> n # (n, n, ...)
        dimtemporal = @ntuple $N k -> 22 # (22, 22, ...)
        table = NamedArray(zeros($T, dimtable),
                           @ntuple($N, k -> namedict), # (namedict, namedict, ...)
                           @ntuple($N, k -> "Dim_k")) # ("Dim_1", "Dim_2", ...)
        marginals = NamedArray(zeros($T, n, $N),
                               # OrderedDict{String,Int}("Dim_$i" => i for i in 1:N)
                               (namedict, OrderedDict{String,Int}(@ntuple $N k -> "Dim_k"=>k)),
                               ("Residue","Dim"))
        ContingencyTable{T,N,A}(alphabet,
                                zeros(T, dimtemporal),
                                table,
                                marginals,
                                zero(T))
        end
end

function ContingencyTable(::Type{T}, ::Type{Val{N}}, alphabet::A) where {T,A,N}
    ContingencyTable{T,N,A}(alphabet)
end

function ContingencyTable(matrix::AbstractArray{T,N}, alphabet::A) where {T,N,A}
    n = length(alphabet)
    @assert size(matrix) == ((n for i in 1:N)...,) "Matrix size doesn't match alphabet length"
    table = ContingencyTable(T, Val{N}, alphabet)
    getarray(table.table)[:] = matrix
    _update_marginals!(table)
    _update_total!(table)
end

# Update!
# -------

# Update the table, marginal and total using temporal

@generated function _update_table!(table::ContingencyTable{T,N,A}) where {T,N,A}
    if A <: ReducedAlphabet
        quote
            temporal = table.temporal::Array{T,N}
            alphabet = getalphabet(table)::A
            freqtable = getarray(table.table)::Array{T,N}
            @inbounds @nloops $N i temporal begin
                @_test_index $N i continue
                @nextract $N a d->alphabet[i_d]
                @_test_index $N a continue
                @nref($N, freqtable, a) += @nref($N, temporal, i)
            end
            table
        end
    elseif (A === UngappedAlphabet) || (A === GappedAlphabet)
        quote
            temporal = table.temporal::Array{T,N}
            freqtable = getarray(table.table)::Array{T,N}
            @inbounds @nloops $N i freqtable begin
                @nref($N, freqtable, i) += @nref($N, temporal, i)
            end
            table
        end
    end
end

@generated function _update_marginals!(table::ContingencyTable{T,N,A}) where {T,N,A}
        quote
            freqtable = getarray(table.table)::Array{T,N}
            marginal  = getarray(table.marginals)::Matrix{T}
            @inbounds @nloops $N i freqtable begin
                value = @nref $N freqtable i
                @_marginal($N, marginal, i, value)
            end
            table
        end
end

function _update_total!(table::ContingencyTable{T,N,A}) where {T,N,A}
    marginals = getarray(table.marginals)::Matrix{T}
    n = size(marginals,1)
    total = zero(T)
    @inbounds @simd for i in 1:n
        total += marginals[i] # marginals[i,1]
    end
    table.total = total
    table
end

"`update_marginals!` updates the marginal and total values using the table."
function update_marginals!(table::ContingencyTable{T,N,A}) where {T,N,A}
    fill!(getarray(table.marginals)::Matrix{T}, zero(T))
    _update_marginals!(table)
    _update_total!(table)
    table
end

function _update!(table::ContingencyTable{T,N,A}) where {T,N,A}
    _update_table!(table)
    update_marginals!(table)
end

function _cleanup_table!(table::ContingencyTable{T,N,A}) where {T,N,A}
    fill!(getarray(table.table)::Array{T,N}, zero(T))
end

function _cleanup_temporal!(table::ContingencyTable{T,N,A}) where {T,N,A}
    fill!(table.temporal, zero(T))
end

"""
`cleanup!` fills the temporal, table and marginals arrays with zeros.
It also sets total to zero.
"""
function cleanup!(table::ContingencyTable{T,N,A}) where {T,N,A}
    _cleanup_temporal!(table)
    _cleanup_table!(table)
    fill!(getarray(table.marginals)::Matrix{T}, zero(T))
    table.total = zero(T)
    table
end

# Fill
# ====

function Base.fill!(table::ContingencyTable{T,N,A}, value::T) where {T,N,A}
    fill!(getarray(table.table)::Array{T,N}, value)
    update_marginals!(table)
end

Base.fill!(table::ContingencyTable{T,N,A}, p::AdditiveSmoothing{T}) where {T,N,A} = fill!(table, p.λ)

@inline Base.fill!(table::ContingencyTable{T,N,A}, p::NoPseudocount) where {T,N,A} = table

# Apply pseudocount
# =================

# This is faster than array[:] += value
function _sum!(matrix::NamedArray, value)
    matrix_array = getarray(matrix)
    @inbounds for i in eachindex(matrix_array)
        matrix_array[i] += value
    end
    matrix
end

"It adds the `pseudocount` value to the table cells."
function apply_pseudocount!(table::ContingencyTable{T,N,A}, pseudocount::T) where {T,N,A}
    _sum!(table.table, pseudocount)
    update_marginals!(table)
end

function apply_pseudocount!(table::ContingencyTable{T,N,A}, p::AdditiveSmoothing{T}) where {T,N,A}
    apply_pseudocount!(table, p.λ)
end

@inline apply_pseudocount!(table::ContingencyTable, p::NoPseudocount) = table

# Normalize
# =========

# This is faster than array[:] /= value
function _div!(matrix::NamedArray, value)
    matrix_array = getarray(matrix)
    @inbounds for i in eachindex(matrix_array)
        matrix_array[i] /= value
    end
    matrix
end

"`normalize!` makes the sum of the frequencies to be one, in place."
function LinearAlgebra.normalize!(table::ContingencyTable{T,N,A}) where {T,N,A}
    if table.total != one(T)
        _div!(table.table, table.total)
        update_marginals!(table)
    end
    table
end

"`normalize` returns another table where the sum of the frequencies is one."
LinearAlgebra.normalize(table::ContingencyTable{T,N,A}) where {T,N,A} = normalize!(deepcopy(table))

# Delete Dimensions
# =================

function _list_without_dimensions(len::Int, output_len::Int, dimensions::Int...)
  ndim = length(dimensions)
  @assert (len-ndim) == output_len "$output_len should be = $(len-ndim)"
  index_list = Array{Int}(undef, output_len)
  j = 1
  @inbounds for i in 1:len
    if ! (i in dimensions)
      index_list[j] = i
      j += 1
    end
  end
  index_list
end

"""
`delete_dimensions!(out::ContingencyTable, in::ContingencyTable, dimensions::Int...)`

This function fills a ContingencyTable with the counts/probabilities on `in` after the
deletion of `dimensions`. i.e. This is useful for getting Pxy from Pxyz.
"""
function delete_dimensions!(output::ContingencyTable{T,S,A},
                            input::ContingencyTable{T,N,A},
                            dimensions::Int...) where {T,N,S,A}
  output_marginals = getarray(output.marginals)
  output_table = getarray(output.table)
  input_marginals = getarray(input.marginals)
  input_table = getarray(input.table)
  output_marginals[:] = input_marginals[:,_list_without_dimensions(N, S, dimensions...)]
  output_table[:] = sum(input_table, dims=dimensions)
  output.total = input.total
  output
end

"""
`delete_dimensions(in::ContingencyTable, dimensions::Int...)`

This function creates a ContingencyTable with the counts/probabilities on `in` after the
deletion of `dimensions`. i.e. This is useful for getting Pxy from Pxyz.
"""
function delete_dimensions(input::ContingencyTable{T,N,A}, dims::Vararg{Int,I}) where {T,N,A,I}
    delete_dimensions!(ContingencyTable(T, Val{N-I}, input.alphabet), input, dims...)
end

for tp in (:Probabilities, :Counts)
    @eval begin
        function delete_dimensions!(output::$(tp){T,S,A}, input::$(tp){T,N,A},
                                    dimensions::Int...) where {T,N,S,A}
            delete_dimensions!(getcontingencytable(output), getcontingencytable(input), dimensions...)
            output
        end

        function delete_dimensions(input::$(tp){T,N,A}, dims::Vararg{Int,I}) where {T,N,A,I}
            output = delete_dimensions(getcontingencytable(input), dims...)
            $(tp)(output)
        end
    end
end


text:
using LinearAlgebra
using LinearManifoldCluster
using Distributions
using Distances
using LMCLUS
using Test
using Random

@testset "Random direction" begin
    n = 1000
    N = 10
    pinf, psup = extrema(LinearManifoldCluster.on_ball(N, n))
    @test pinf > -1.0
    @test psup <  1.0
    pinf, psup = extrema(LinearManifoldCluster.in_ball(N, n))
    @test pinf > -1.0
    @test psup <  1.0
    for d in colwise(Euclidean(), LinearManifoldCluster.on_ball(N,n), zeros(N))
        @test d ≈ 1.
    end
end

@testset "Translation vectors" begin
    n = 1000
    N = 10
    μ_bounds = 0.25 => 0.75
    μ = LinearManifoldCluster.translations(N, n, μ_bounds)
    # distance from center of unit hypercube to
    # any translation vector should be within bounds
    for d in colwise(Euclidean(), μ, ones(N)/2.)
        @test μ_bounds[1]/2. <= d <= μ_bounds[2]/2.
    end
end

@testset "Basis formation" begin
    N = 10
    M = 3
    B, B′ = LinearManifoldCluster.lm_basis(N, M)
    @test size(B) == (N,M)
    @test size(B′) == (N,N-M)
    P = B*B'
    @test P*P ≈ P
    @test P' ≈ P
    P′ = B′*(B′)'
    @test P′*P′ ≈ P′
    @test (P′)' ≈ P′
end

@testset "Manifold cluster generator" begin
    n = 1000
    N = 3
    M = 2
    μ_bounds = 0.25 => 0.75
    t = LinearManifoldCluster.translations(N, 1, μ_bounds)[:]
    D = Matrix{Float64}(I, N, N)
    B = D[:,1:2]
    B′ = D[:, 3:3]
    Φ = 0.5
    E = 0.05

    Random.seed!(230898573857)
    mfld_s = LinearManifoldCluster.generate_cluster(n, t, B, B′, Φ, E)

    Random.seed!(230898573857)
    Φ = fill(0.5, M)
    E = fill(0.05, N-M)
    mfld = LinearManifoldCluster.generate_cluster(n, t, B, B′, Φ, E)
    # one scalar bounds and array bounds return values
    @test mfld_s == mfld

    # manifold size & space dimension
    @test size(mfld) == (N,n)
    # if Φ is scalar, then manifold forms a hypercube
    # with side equal to Φ (± E)?
    mbound_max = norm(fill(Φ[1]+E[1],M),2)
    for d in colwise(Euclidean(), mfld.-t, zeros(N))
        @test 0.0 <= d <= mbound_max
    end
    # if E is scalar, then a distance from any point to
    # the manifold should be within ||[ɛ, ...,ɛ]||
    mdist_max = norm(E,2)
    for d in LMCLUS.distance_to_manifold(mfld, t, B)
        @test 0.0 <= d <= mdist_max
    end

    # manifold dimension generator
    @test LinearManifoldCluster.mdims(10,13,7) == [1, 2, 3, 4, 5, 6, 7, 7, 7, 7, 7, 7, 7]
end

@testset "Multi-manifold cluster generator" begin
    n = 1000
    N = 3
    M = [1,2,2]
    τ = 0.25 => 0.75
    Φ = Vector{Float64}[   # bounds on manifold points (per manifold wrt dims)
        [.6],
        [.2,  .73],
        [.4, .65]
    ]
    E = Vector{Float64}[   # bound of a point extent (per manifold wrt dims)
        [.05,  .03],
        [.01],
        [.04]
    ]
    lmcs, Bs, μs = generate(n,N,M,τ,Φ,E)
    @test length(lmcs) == length(M) # number of manifold
    @test size(lmcs[1]) == (N,n)      # number of manifold points
    for i in length(Φ)
        @test size(Bs[i]) == (N,length(Φ[i]))  # bases dimensions
    end
    @test size(μs) == (N,length(M))

    # multi-manifold generator (scalar)
    n = 1000
    N = 3
    m = 12
    τ = 0.25 => 0.75
    κ = 1.5
    lmcs, Bs, μs = generate(n,N,m,τ,κ)
    @test length(lmcs) == m
end


text:
#=
This file is part of the replication code for: Hasenzagl, T., Pellegrino, F., Reichlin, L., & Ricco, G. (2020). A Model of the Fed's View on Inflation.
Please cite the paper if you are using any part of the code for academic work (including, but not limited to, conference and peer-reviewed papers).
=#

function ex_blkdiag(x, y, z=[])

# ----------------------------------------------------------------------------------------------------------------------
# Custom blkdiag function
# ----------------------------------------------------------------------------------------------------------------------

    # Initialise
    rx, cx = size(x);
    ry, cy = size(y);

    if size(z)[1] == 0
        rz = 0;
        cz = 0;
    else
        rz, cz = size(z);
    end

    # Blkdiag
    out                                     = zeros(rx+ry+rz, cx+cy+cz);
    out[1:rx, 1:cx]                         = x;
    out[rx+1:rx+ry, cx+1:cx+cy]             = y;
    out[rx+ry+1:rx+ry+rz, cx+cy+1:cx+cy+cz] = z;

    return out;
end


text:
module FranklinTheorems

import Franklin
export lx_fakebiblabel, lx_newcounter, lx_stepcounter, lx_arabic

module ExportConfigPath

	export css_path, config_path

	"""
	Returns a path to FranklinTheorem's default CSS file. Copy this file to the Franklin `\\_css\\` directory.
	"""
	function css_path()
		return normpath(joinpath(@__FILE__, "..", "FranklinTheorems.css"))
	end

	"""
	Returns a path to FranklinTheorem's markdown file, full of LaTeX-style definitions.
	By default, the theorem, lemma, and definition blocks are defined.
	If a list of strings is given, then the associated counters and environments are created.
	Since this is dynamic, because of the optional parameters, this returns a path to a tempfile.
	"""
	function config_path()
		config_path([])
	end
	function config_path(theorem_class_names)
		all_class_names = [["theorem", "lemma", "definition"]; theorem_class_names]
		baseMarkdownText = read(normpath(joinpath(@__FILE__, "..", "FranklinTheorems.md")), String)

		f = tempname()
		write(f, baseMarkdownText * "\n" * make_custom_environments(all_class_names))
		return f
	end

	"""
	Returns a path to a temporary file that includes Franklin-Markdown code for a custom amstheorem-style environment.
	"""
	function make_custom_environments(theorem_class_names)

		env_string = ""
		enable_string = "\\newcommand{\\enabletheorems}{"

		for i=1:length(theorem_class_names)
			defn_name = theorem_class_names[i]
			defn_titlecase = titlecase(defn_name)
			defn_body = """
\\newcommand{\\$(defn_name)ref}[1]{\\cite{!#1}}

\\newenvironment{$(defn_name)}[2]{
	\\stepcounter{Num$(defn_name)}
	\\begin{thmBlock}{$(defn_titlecase) \\arabic{Num$(defn_name)}}{#1}{$(defn_titlecase)}{ \\fakebiblabel{!#2 @ $(defn_titlecase) \\arabic{Num$(defn_name)}} }
}{
	\\end{thmBlock}
}
			"""
			env_string = env_string * "\n" * defn_body
			enable_string = enable_string * "\n" * "\\newcounter{Num$(defn_name)}"
		end

		enable_string = enable_string * "\n}"

		return env_string * "\n" * enable_string
	end

end

using .ExportConfigPath

"""
Specialized code for automatic numbering by abuse of the biblabel function.

- Given an input of form `\\fakebiblable{@ Legible Name}`, this returns the string `"Legible Name"`
- Given an input of form `\\fakebiblabel{goob @ Legible Name}`, this runs `\\biblabel{goob}{Legible Name}`, so that Franklin's referencing system keeps track of `goob` as a name that points to `"Legible Name"`, and the function returns the string `"goob"`.

This `@` symbol approach is a hack to allow a single optional input variable to the function.
The `@` symbol was chosen because of its relevance to LaTeX.

Further, the string `"Legible Name"` is always parsed by Franklin, which has two notable consequences
1. Other lx_functions like `\\arabic{counter}` can be used in the full name
2. Markdown syntaxing like `_italics_` will be parsed into `<em>italics<em>`, which breaks the references

So, while `Legible Name` can reference counters, it should not have any markdown styling.
"""
function lx_fakebiblabel(lxc, _)
	args = Franklin.content(lxc.braces[1]) |> strip |> split

	nvars = length(args)

	if nvars < 2
		error("Bad fakebiblabel input: Less than 2 inputs were given. Input was:\n\\fakebiblablel{" * join(args, " ") * "}")
	end

	if args[1] == "@"
		return Franklin.fd2html(join(args[2:end], " "), internal=true, nop=true)
	elseif args[2] == "@"
		ref_name = args[1]
		legible_name = Franklin.fd2html(join(args[3:end], " "), internal=true, nop=true)
		Franklin.fd2html("\\biblabel{" * ref_name * "}{" * legible_name * "}", internal=true)
		return ref_name
	else
		error("Bad fakebiblabel input: @ was not the second input. Input was:\n\\fakebiblablel{" * join(args, " ") * "}")
	end
end

# Module for the counters. By making this into a module, this hides
# access to the counters dictionary. If you need a new way to interact
# with the system of counters, impliment a new method in this module.
module LxCounters

	# Used for latex-style parsing
	using FranklinUtils

	export lx_newcounter, lx_stepcounter, lx_arabic

	# The core dictionary that stores all the things we want to count.
	counters = Dict{String,Int64}()

	"""
	Call `\\newcounter{name}` to create a new counter. Returns an empty string.
	"""
	function lx_newcounter(lxc, _)
		args, kwargs = lxargs(lxc)

		name = string(args[1])

		delete!(counters, name)
		get!(counters, name, 0)
		return ""
	end

	"""
	Call `\\stepcounter{name}` to incriment the value of a counter by 1. Returns an empty string.
	"""
	function lx_stepcounter(lxc, _)
		args, kwargs = lxargs(lxc)

		name = string(args[1])

		counters[name] = counters[name] + 1
		return ""
	end

	"""
	Call `\\arabic{name}` to return the value of a counter, as an arabic number.
	"""
	function lx_arabic(lxc, _)
		args, kwargs = lxargs(lxc)

		name = string(args[1])

		return string(counters[name])
	end

end

using .LxCounters

end # Module


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

using Markdown
using InteractiveUtils

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

# ╔═╡ 0e77075e-3b2b-11eb-3ccc-b5566e9738bf
begin
	using Markdown
	using Luxor
	using Colors
	using PlutoUI
end

# ╔═╡ 473d669e-3b2b-11eb-0cf8-6525810c72ea
f((x,y)) = (cos(y)^5*3^x, x^3+y) .|> p -> mod(p,1);

# ╔═╡ f2cdbc9e-3bc2-11eb-03da-bf1ad0245bb2
g((x,y)) = let a=2.1, b=5.9
	(sin(x*y/b)*y + cos(a*x - y), x + sin(y)/b) .|> p -> mod(p,1)
end;

# ╔═╡ 9b32c812-3b57-11eb-242b-0bb52e31a065
arnold((x,y)) = (2x + y, x + y) .% 1;

# ╔═╡ d7240ae8-3b57-11eb-188b-3f0ee0efecad
function compose(f, n)
    function (x)
        val = x
        for _ in 1:n
            val = f(val)
        end
        return val
    end
end;

# ╔═╡ 1173cd5a-3b5a-11eb-376d-99a5b00641bf
Base.:^(f::Function, n::Integer) = compose(f,n);

# ╔═╡ 4663058e-3b59-11eb-2929-dbaea3ec1c78
dist((x1,x2),(y1,y2)) = sqrt((x1-y1)^2 + (x2-y2)^2);

# ╔═╡ 168b884a-3b59-11eb-1e46-8d67f5ead8f8
function dyndist(x,y,f,n)
	d = dist(x,y)
	for _ in 1:n
		x = f(x)
		y = f(y)
		d = max(d, dist(x,y))
	end
	d
end;

# ╔═╡ 8284ce6c-3b31-11eb-2796-835ee6d006ca
Np = 800

# ╔═╡ 552c91cc-3b5c-11eb-3301-6dd3e136df1d
@bind cut Slider(0:0.01:1)

# ╔═╡ b03b0738-3b5c-11eb-26c3-173861b22912
@bind n Slider(0:1:20)

# ╔═╡ 3f0b2320-3b5e-11eb-24ea-7d0a59c99f3e
p = (.671,.81);

# ╔═╡ e68f1386-3b5e-11eb-0838-d95297ccef28
fun = f;

# ╔═╡ 66ac124e-3b52-11eb-1bea-339a709ddf25
let	s = 800
	Drawing(s, s, "my-drawing.png")
	background("antiquewhite")
	for x in 0:(1/Np):1, y in 0:(1/Np):1
		x = rand() # * .01 + p[1] - .005
		y = rand() # * .01 + p[2] - .005
		let (tx,ty) = (fun^n)((x,y))
			setcolor(HSV(360tx,.5ty,dyndist((x,y),p,fun,n)))
		end
		if dyndist((x,y),p,fun,n) < cut
			setcolor("black")
		end
		circle(Point(s .* (x, y)), sqrt(s)*0.026, :fill) # .+ (.005,.005) .- p
	end
	setcolor("red")
	circle(s .* Point(p), 5, :fill) # Point(.5,.5)
	finish()
	preview()
end

# ╔═╡ Cell order:
# ╠═0e77075e-3b2b-11eb-3ccc-b5566e9738bf
# ╠═473d669e-3b2b-11eb-0cf8-6525810c72ea
# ╠═f2cdbc9e-3bc2-11eb-03da-bf1ad0245bb2
# ╠═9b32c812-3b57-11eb-242b-0bb52e31a065
# ╠═d7240ae8-3b57-11eb-188b-3f0ee0efecad
# ╠═1173cd5a-3b5a-11eb-376d-99a5b00641bf
# ╠═4663058e-3b59-11eb-2929-dbaea3ec1c78
# ╠═168b884a-3b59-11eb-1e46-8d67f5ead8f8
# ╠═8284ce6c-3b31-11eb-2796-835ee6d006ca
# ╠═552c91cc-3b5c-11eb-3301-6dd3e136df1d
# ╠═b03b0738-3b5c-11eb-26c3-173861b22912
# ╠═3f0b2320-3b5e-11eb-24ea-7d0a59c99f3e
# ╠═e68f1386-3b5e-11eb-0838-d95297ccef28
# ╠═66ac124e-3b52-11eb-1bea-339a709ddf25


text:
"""
```
augment_states(m::AnSchorfheide, TTT::Matrix{T}, RRR::Matrix{T}, CCC::Matrix{T}) where {T<:AbstractFloat}
```

### Arguments

-`m`: the m object
-`TTT`, `RRR`, and `CCC`: matrices of the state transition equation

### Return values

- `TTT_aug`, `RRR_aug`, and `CCC_aug`: extend the corresponding input matrices to include
  jobservables which are growth rates.

### Description

Some observables in the model are growth rates, which are calculated as a linear combination
of a present and lagged state (which is not yet accounted for in the `TTT`, `RRR`,and `CCC`
matrices). To improve the performance of `gensys`, these additional states are added after
the model is solved. `augment_states` assigns an index to each lagged state, and extends the
input `TTT`, `RRR`, and `CCC` matrices to accommodate the additional states and capture the
lagged state value in the current state vector. `RRR` and `CCC` are mostly augmented with
zeros.

The diagram below shows how `TTT` is extended to `TTT_aug`.

                TTT_aug
     (m.endogenous_states_additional
                  x
     m.endogenous_states_additional)
     _________________________________
    |                     |           |
    |          TTT        | endog_    |
    | (endogenous_states  | states_   |
    |          x          | augmented |
    |  endogenous_states) |           |
    |_____________________|           |
    |                                 |
    |    endogenous_states_augmented  |
    |_________________________________|

"""
function augment_states(m::AnSchorfheide, TTT::Matrix{T}, RRR::Matrix{T},
                        CCC::Vector{T}) where {T<:AbstractFloat}
    return TTT, RRR, CCC
end


text:
using AutoGrad

# === building blocks for NN ===

function init_weights(;n_in=1, n_hidden=10, n_out=1)
    W1 = randn(n_hidden, n_in) # for left multiply W1*x
    b1 = zeros(n_hidden)
    W2 = randn(n_out, n_hidden)
    b2 = zeros(n_out)
    params = [W1, b1, W2, b2]
    return params
end

function sizes_and_length(params)
    sizes = map(size, params)
    lengths = map(length, params)
    total_l = sum(lengths)
    return sizes, lengths, total_l
end

function predict(params, t, y0, t0; act=tanh)
    W1, b1, W2, b2 = params

    # normal NN calculation
    a = act.(W1*t .+ b1)
    out = W2*a .+ b2

    # force intial condition
    y = y0 .+ (t .- t0) .* out
    return y
end

# gradient w.r.t t
# get diagonal of jacobian to vectorize over t
predict_sum(params, t, y0, t0) = sum(predict(params, t, y0, t0))
predict_dt = grad(predict_sum, 2)

# === flatten & unflatten ===

flat_opt = p->collect(Iterators.flatten(p)) # flatten operation for a single NN params

function unflatten(params_flat, sizes, lengths)
    params = []
    i1 = 1
    for j in 1:length(sizes)
        s = sizes[j]
        l = lengths[j]
        i2 = i1+l
        #p = reshape(params_flat[i1:i2-1], s)
        p = reshape(view(params_flat,i1:i2-1), s)
        push!(params, p)
        i1 = i2
    end
    return params
end

function unflatten_all(all_flattened_params, Nvar, total_l, sizes, lengths)
    params_list = []
    all_params_reshape = reshape(all_flattened_params, total_l, Nvar)
    for i in 1:Nvar
        params = unflatten(all_params_reshape[:,i], sizes, lengths)
        push!(params_list, params)
    end
    return params_list
end


text:
"YAML configuration"
type YamlConfiguration <: Configuration
    name::AbstractString
    source::AbstractString
    state::LifeCycle.State
    root::LoggerConfig

    properties::PROPERTIES
    appenders::APPENDERS
    loggers::LOGCONFIGS
    filter::FILTER

    data::Dict # Configuration data
end

# Register configuration type
LOG4JL_CONFIG_TYPES[:YAML] = YamlConfiguration
LOG4JL_CONFIG_EXTS[:YAML]  = [".yaml", ".yml"]

function YamlConfiguration(cfgloc::AbstractString, cfgname::AbstractString="YAML")
    eval(:(import YAML)) # Package lazy eval
    conf = YAML.load(open(cfgloc))
    fltr = FILTER()

    # Set status logger parameters
    if haskey(conf, "configuration")
        stat = conf["configuration"]
        cfgname = get(stat, "name", cfgname)
        # custom levels
        haskey(conf["configuration"], "customlevels") && for (l,v) in conf["configuration"]["customlevels"]
            Level.add(symbol(l), Int32(v))
        end
        # status
        haskey(stat, "status") && level!(LOGGER, evaltype((stat["status"] |> uppercase), "Level"))
        # filters
        fltrs = filterconfig(stat, "filter")
        if length(fltrs) > 0
            fltr = parsefilters(fltrs)
        end
    else
        error(LOGGER, "Malformed configuration: `configuration` node does not exist.")
    end
    YamlConfiguration(cfgname, cfgloc, LifeCycle.INITIALIZED, LoggerConfig(),
                      PROPERTIES(), APPENDERS(), LOGCONFIGS(), fltr, conf)
end
getconfig(::Type{YamlConfiguration}, cfgloc::AbstractString, cfgname::AbstractString="YAML") = YamlConfiguration(cfgloc, cfgname)

appender(cfg::YamlConfiguration, name::AbstractString) = get(cfg.appenders, name, nothing)
appenders(cfg::YamlConfiguration) = cfg.appenders
loggers(cfg::YamlConfiguration) = cfg.loggers

function configure(cfg::YamlConfiguration)
    # if configuration is malformed
    if !haskey(cfg.data, "configuration")
        warn(LOGGER, "Malformed configuration: creating default ERROR-level Root logger with Console appender")
        default!(cfg)
        return
    end
    conf = cfg.data["configuration"]

    # parse properties
    haskey(conf, "properties") && for p in values(conf["properties"])
        cfg.properties[p["name"]] = p["value"]
    end

    # update appenders and loggers with pattern values
    for (k,v) in cfg.properties
        haskey(conf, "appenders") && subst_value(conf["appenders"], k, v)
        haskey(conf, "loggers")   && subst_value(conf["loggers"], k, v)
    end

    # parse appenders
    haskey(conf, "appenders") && for (atype, aconf) in conf["appenders"]
        # get appender type
        apndType = contains(atype, ".") ?
                   evaltype(atype, "") :
                   evaltype(atype, "Appenders")
        apndType === Void && continue
        apndCfg = Dict{AbstractString,Any}(aconf)

        # find first layout description and create layout object
        lyts = filterconfig(aconf, "layout")
        if length(lyts) > 0 # take first layout description
            strLType = first(keys(lyts))
            # get type from internal
            lytType = contains(strLType, ".") ?
                      evaltype(strLType, "") :
                      evaltype(strLType, "Layouts")
            if lytType !== Void
                apndCfg["layout"] = lytType(aconf[strLType])
            end
            delete!(apndCfg, strLType)
        end

        # check filters
        fltrs = filterconfig(apndCfg, "filter")
        if length(fltrs) > 0
            fltr = parsefilters(fltrs)
            apndCfg["filter"] = fltr
            for strFType in keys(fltrs)
                delete!(apndCfg, strFType)
            end
        end

        # create appender object
        cfg.appenders[aconf["name"]] = apndType(apndCfg)
    end

    refs = Dict[]
    # add loggers
    if haskey(conf, "loggers")
        lconf = conf["loggers"]

        # if there is a root logger configuration
        if haskey(lconf, "root")
            rconf = lconf["root"]
            cfg.root.level = configlevel(rconf)

            if haskey(rconf, "appenderref")
                push!(refs, Dict("name"=>"root", "appenderref" => rconf["appenderref"]))
            end
        else
            warn(LOGGER, "No Root logger was configured, creating default ERROR-level Root logger with Console appender")
            default!(cfg)
        end

        # if there are other logger configuration
        if haskey(lconf, "logger")
            for lcconf in lconf["logger"]
                lcname = lcconf["name"]
                lcadd = get(lcconf, "additivity", false)
                lclvl = configlevel(lcconf)
                #TODO: add filters
                logger!(cfg, lcname, LoggerConfig(lcname, level=lclvl, additive=lcadd))
            end
            append!(refs, lconf["logger"])
        end
    else
        warn(LOGGER, "No Loggers were configured, using default. Is the Loggers element missing?")
        default!(cfg)
    end

    # configure references
    for lcconf in refs
        lcname = lcconf["name"]
        lc = lcname == "root" ? cfg.root : cfg.loggers[lcname]
        lcconf_arefs = lcconf["appenderref"]
        lcconf_arefs = isa(lcconf_arefs, Dict) ? [lcconf_arefs] : lcconf_arefs
        for lcconf_aref in lcconf_arefs
            if haskey(lcconf_aref, "ref")
                aref = lcconf_aref["ref"]
                if haskey(cfg.appenders, aref)
                    apnd = cfg.appenders[aref]
                    #TODO: add filters
                    reference!(lc, apnd, configlevel(lcconf_aref))
                else
                    error(LOGGER, "Unable to locate appender '$aref' for logger '$lcname'")
                end
            else
                error(LOGGER, "No reference provided for logger '$lcname'")
            end
        end
    end

    # setup parents
    parents!(cfg)
end

function subst_value(conf::Dict, k::AbstractString, v::AbstractString)
    for (dk,dv) in conf
        if isa(dv, AbstractString)
            pat = "\${$k}"
            if contains(dv, pat)
                conf[dk] = replace(dv, pat, v)
            end
        elseif isa(dv, Dict)
            subst_value(dv, k, v)
        end
    end
    return conf
end

function filterconfig(conf::Dict, stype::AbstractString)
    ckeys = collect(keys(conf))
    lidxs = find(s->contains(s,stype), map(lowercase, ckeys))
    return if length(lidxs) > 0
        filter((k,v)->k in ckeys[lidxs], conf)
    else
        Dict()
    end
end

function parsefilters(conf::Dict)
    fltrs = Filter[]
    for (fltrstrtype, fltrcfg) in conf
        fltrtype = evaltype(fltrstrtype, "")
        fltr = (fltrcfg === nothing) ? fltrtype() : fltrtype(Dict{AbstractString,Any}(fltrcfg))
        fltr !== nothing && push!(fltrs, fltr)
    end
    fltrcount = length(fltrs)
    return if fltrcount == 0
        FILTER()
    elseif fltrcount == 1
        FILTER(fltrs[1])
    else
        FILTER(CompositeFilter(fltrs))
    end
end


text:
""" model function for an isotope """
function alphamodel(x,var)
    μ = var[1]
    σ = var[2]
    τ = var[3]
    A = var[4]
    amp = A / 2τ
    oneOverSqrt2 = 1.0 / √2
    sigScaledByTau  = σ / τ
    return @.  amp * exp((x - μ ) / τ + (sigScaledByTau^2) ) * erfc( oneOverSqrt2 * ((x - μ ) / σ  +  sigScaledByTau ))
end