text:
""" Data structure for (directed) wiring diagrams, aka string diagrams.

A (directed) wiring diagram consists of a collection of boxes with input and
output ports connected by wires. A box can be atomic (possessing no internal
structure) or can itself be a wiring diagram. Thus, wiring diagrams can be
nested recursively. Wiring diagrams are closely related to what the CS
literature calls "directed graphs with ports" or more simply "port graphs". The
main difference is that a wiring diagram has an "outer box": a wiring diagram
has its own ports that can be connected to the ports of its boxes.

This module provides a generic data structure for wiring diagrams. Arbitrary
data can be attached to the boxes, ports, and wires of a wiring diagram. The
diagrams are "abstract" in the sense that they cannot be directly rendered as
raster or vector graphics. However, they form a useful intermediate
representation that can be serialized to and from GraphML or translated into
Graphviz or other declarative diagram languages.
"""
module DirectedWiringDiagrams
export AbstractBox, Box, WiringDiagram, Wire, Port, PortKind,
  InputPort, OutputPort, input_ports, output_ports, input_id, output_id,
  outer_ids, boxes, box_ids, nboxes, nwires, box, wires, has_wire, graph,
  add_box!, add_boxes!, add_wire!, add_wires!, rem_box!, rem_boxes!, rem_wire!,
  rem_wires!, port_value, validate_ports, is_permuted_equal,
  all_neighbors, neighbors, outneighbors, inneighbors, in_wires, out_wires,
  singleton_diagram, induced_subdiagram, encapsulated_subdiagram,
  ocompose, substitute, encapsulate

using Compat
using AutoHashEquals

using ...Present, ...CSetDataStructures, ...Graphs.BasicGraphs
using ...Graphs.BasicGraphs: TheoryGraph
import ...Graphs: all_neighbors, neighbors, outneighbors, inneighbors

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

""" Kind of port: input or output.
"""
@enum PortKind InputPort OutputPort

""" A port on a box to which wires can be connected.
"""
@auto_hash_equals struct Port
  box::Int
  kind::PortKind
  port::Int
end
set_box(port::Port, box::Int) = Port(box, port.kind, port.port)

function Base.isless(p1::Port, p2::Port)::Bool
  # Lexicographic order.
  p1.box < p2.box ||
    (p1.box == p2.box &&
      (p1.kind < p2.kind || (p1.kind == p2.kind && p1.port < p2.port)))
end

""" A wire connecting one port to another.
"""
@auto_hash_equals struct Wire{Value}
  value::Value
  source::Port
  target::Port
end

Wire(value, src::Tuple{Int,PortKind,Int}, tgt::Tuple{Int,PortKind,Int}) =
  Wire(value, Port(src[1],src[2],src[3]), Port(tgt[1],tgt[2],tgt[3]))
Wire(value, src::Tuple{Int,Int}, tgt::Tuple{Int,Int}) =
  Wire(value, Port(src[1],OutputPort,src[2]), Port(tgt[1],InputPort,tgt[2]))
Wire(value, pair::Pair) = Wire(value, first(pair), last(pair))

Wire(src::Port, tgt::Port) = Wire(nothing, src, tgt)
Wire(src::Tuple, tgt::Tuple) = Wire(nothing, src, tgt)
Wire(pair::Pair) = Wire(nothing, first(pair), last(pair))

function Base.show(io::IO, wire::Wire)
  skip_kind = wire.source.kind == OutputPort && wire.target.kind == InputPort
  show_port = (io::IO, port::Port) -> begin
    if skip_kind
      print(io, "($(port.box),$(port.port))")
    else
      print(io, "($(port.box),$(string(port.kind)),$(port.port)")
    end
  end
  print(io, "Wire(")
  if !isnothing(wire.value)
    show(io, wire.value)
    print(io, ", ")
  end
  show_port(io, wire.source)
  print(io, " => ")
  show_port(io, wire.target)
  print(io, ")")
end

function Base.isless(w1::Wire, w2::Wire)::Bool
  # Lexicographic order.
  isless(w1.source, w2.source) ||
    (w1.source == w2.source && isless(w1.target, w2.target))
end

""" Internal wiring diagram data corresponding to `Port`. Do not use directly.
"""
@auto_hash_equals struct PortData
  kind::PortKind
  port::Int
end
PortData(port::Port) = PortData(port.kind, port.port)
Port(port::PortData, v::Int) = Port(v, port.kind, port.port)

""" Internal wiring diagram data corresponding to `Wire`. Do not use directly.
"""
@auto_hash_equals struct WireData{Value}
  value::Value
  source::PortData
  target::PortData
end
function WireData(wire::Wire)
  WireData(wire.value, PortData(wire.source), PortData(wire.target))
end
function Wire(wire::WireData, src::Int, tgt::Int)
  Wire(wire.value, Port(wire.source, src), Port(wire.target, tgt))
end

""" Base type for any box (node) in a wiring diagram.

This type represents an arbitrary black box with inputs and outputs.
"""
abstract type AbstractBox end

input_ports(box::AbstractBox)::Vector = box.input_ports
output_ports(box::AbstractBox)::Vector = box.output_ports

""" An atomic box in a wiring diagram.

These boxes have no internal structure.
"""
@auto_hash_equals struct Box{Value} <: AbstractBox
  value::Value
  input_ports::Vector
  output_ports::Vector
end

Box(inputs::Vector, outputs::Vector) = Box(nothing, inputs, outputs)

function Base.show(io::IO, box::Box)
  print(io, "Box(")
  if !isnothing(box.value)
    show(io, box.value)
    print(io, ", ")
  end
  print(io, "[")
  join(io, [sprint(show, port) for port in box.input_ports], ",")
  print(io, "], [")
  join(io, [sprint(show, port) for port in box.output_ports], ",")
  print(io, "])")
end

@present TheoryWiringDiagramGraph <: TheoryGraph begin
  Box::Data
  WireData::Data

  box::Attr(V,Box)
  wire::Attr(E,WireData)
end

const WiringDiagramGraphUnionAll =
  ACSetType(TheoryWiringDiagramGraph, index=[:src, :tgt])

""" Internal datatype for graph underlying a directed wiring diagram.

Boxes and wires are attached to vertices and edges, respectively.
"""
const WiringDiagramGraph = WiringDiagramGraphUnionAll{
  Union{AbstractBox,Nothing},WireData}

""" A directed wiring diagram, also known as a string diagram.

The wiring diagram is implemented using the following internal data structure.
The "skeleton" of the diagram is an instance of `Catlab.Graphs.AbstractGraph`: a
directed multigraph whose vertices correspond to boxes and whose edges
correspond to wires. There are two special vertices, accessible via `input_id`
and `output_id`, that represent the input and output ports of the outer box.
"""
mutable struct WiringDiagram{Theory} <: AbstractBox
  graph::WiringDiagramGraph
  value::Any
  input_ports::Vector
  output_ports::Vector

  function WiringDiagram{T}(value, inputs::Vector, outputs::Vector) where T
    graph = WiringDiagramGraph()
    diagram = new{T}(graph, value, inputs, outputs)
    add_vertices!(graph, 2, box=nothing)
    return diagram
  end
  function WiringDiagram(d::WiringDiagram{T}) where T
    # Copy constructor for shallow copy.
    new{T}(copy(d.graph), d.value, d.input_ports, d.output_ports)
  end
end

function WiringDiagram{T}(inputs::Vector, outputs::Vector) where T
  WiringDiagram{T}(nothing, inputs, outputs)
end
WiringDiagram(args...) = WiringDiagram{Any}(args...)

input_id(::WiringDiagram) = 1
output_id(::WiringDiagram) = 2
outer_ids(::WiringDiagram) = (1,2)

""" Check equality of wiring diagrams.

Warning: This method checks equality of the underlying graph representation, not
mathematical equality which involves graph isomorphism.

See also: `is_permuted_equal`
"""
function Base.:(==)(d1::WiringDiagram, d2::WiringDiagram)
  (input_ports(d1) == input_ports(d2) &&
   output_ports(d1) == output_ports(d2) && d1.value == d2.value &&
   boxes(d1) == boxes(d2) && sort!(wires(d1)) == sort!(wires(d2)))
end

""" Check equality of wiring diagram under permutation of boxes.

When the boxes in the first diagram `d1` are permuted according to `σ`,
does it become identical to the second diagram `d2`?
"""
function is_permuted_equal(d1::WiringDiagram, d2::WiringDiagram, σ::Vector{Int})
  @assert nboxes(d1) == length(σ) && nboxes(d2) == length(σ)
  d1_ids, d2_ids = box_ids(d1), box_ids(d2)
  box_map = Dict{Int,Int}(d1_ids[σ[i]] => d2_ids[i] for i in eachindex(σ))
  is_induced_equal(d1, d2, box_map)
end
function is_induced_equal(d1::WiringDiagram, d2::WiringDiagram, box_map::Dict{Int,Int})
  box_map[input_id(d1)] = input_id(d2)
  box_map[output_id(d1)] = output_id(d2)
  map_wire = wire::Wire -> Wire(wire.value,
    set_box(wire.source, box_map[wire.source.box]),
    set_box(wire.target, box_map[wire.target.box]))
  (input_ports(d1) == input_ports(d2) && output_ports(d1) == output_ports(d2) &&
   all(box(d1,v) == box(d2,box_map[v]) for v in box_ids(d1)) &&
   sort!(map(map_wire, wires(d1))) == sort!(wires(d2)))
end

Base.copy(diagram::WiringDiagram) = WiringDiagram(diagram)

function Base.show(io::IO, diagram::WiringDiagram{T}) where T
  sshowcompact = x -> sprint(show, x, context=:compact => true)
  print(io, "WiringDiagram")
  if T != Any
    print(io, "{$T}")
  end
  print(io, "(")
  if !isnothing(diagram.value)
    show(io, diagram.value)
    print(io, ", ")
  end
  print(io, "[")
  join(io, map(sshowcompact, input_ports(diagram)), ",")
  print(io, "], [")
  join(io, map(sshowcompact, output_ports(diagram)), ",")
  print(io, "], ")
  if get(io, :compact, false)
    print(io, "{$(nboxes(diagram)) boxes}, {$(nwires(diagram)) wires}")
  else
    print(io, "\n[ $(input_id(diagram)) => {inputs},\n  ")
    print(io, "$(output_id(diagram)) => {outputs},\n  ")
    join(io, [ "$v => $(sshowcompact(box(diagram, v)))"
               for v in box_ids(diagram) ], ",\n  ")
    print(io, " ],\n[ ")
    join(io, map(sshowcompact, wires(diagram)), ",\n  ")
    print(io, " ]")
  end
  print(io, ")")
end

# Imperative interface
######################

# Basic accessors.

box(f::WiringDiagram, v::Int) = subpart(f.graph, v, :box)

function boxes(f::WiringDiagram)
  collect(AbstractBox, subpart(f.graph, box_ids(f), :box))
end

nboxes(f::WiringDiagram) = nv(f.graph) - 2

function box_ids(f::WiringDiagram)
  Int[ v for v in 1:nv(f.graph) if !(v in outer_ids(f)) ]
end

function wires(f::WiringDiagram, src::Int, tgt::Int)
  [ Wire(subpart(f.graph, e, :wire), src, tgt)
    for e in edges(f.graph, src, tgt) ]
end
function wires(f::WiringDiagram)
  g = f.graph
  [ Wire(subpart(g, e, :wire), src(g, e), tgt(g, e)) for e in edges(g) ]
end
nwires(f::WiringDiagram) = ne(f.graph)

has_wire(f::WiringDiagram, src::Int, tgt::Int) =
  !isempty(edges(f.graph, src, tgt))
has_wire(f::WiringDiagram, wire::Wire) =
  wire in wires(f, wire.source.box, wire.target.box)
has_wire(f::WiringDiagram, pair::Pair) = has_wire(f, Wire(pair))

function input_ports(f::WiringDiagram, v::Int)
  if v == input_id(f)
    error("Input vertex does not have input ports within wiring diagram")
  elseif v == output_id(f)
    output_ports(f)
  else
    input_ports(box(f, v))
  end
end

function output_ports(f::WiringDiagram, v::Int)
  if v == input_id(f)
    input_ports(f)
  elseif v == output_id(f)
    error("Output vertex does not have output ports within wiring diagram")
  else
    output_ports(box(f, v))
  end
end

function port_value(f::WiringDiagram, port::Port)
  get_ports = port.kind == InputPort ? input_ports : output_ports
  get_ports(f, port.box)[port.port]
end

# Graph mutation.

function add_box!(f::WiringDiagram, box::AbstractBox)
  add_vertex!(f.graph, box=box)
  return nv(f.graph)
end

function add_boxes!(f::WiringDiagram, boxes)
  boxes = collect(boxes)
  add_vertices!(f.graph, length(boxes), box=boxes)
end

function rem_box!(f::WiringDiagram, v::Int)
  @assert v ∉ outer_ids(f)
  rem_vertex!(f.graph, v)
end

function rem_boxes!(f::WiringDiagram, vs)
  @assert all(v ∉ outer_ids(f) for v in vs)
  rem_vertices!(f.graph, sort!(collect(vs)))
end

function add_wire!(f::WiringDiagram, wire::Wire)
  validate_ports(port_value(f, wire.source), port_value(f, wire.target))
  add_edge!(f.graph, wire.source.box, wire.target.box, wire=WireData(wire))
end
add_wire!(f::WiringDiagram, pair::Pair) = add_wire!(f, Wire(pair))

function add_wires!(f::WiringDiagram, wires)
  for wire in wires
    add_wire!(f, wire)
  end
end

function rem_wire!(f::WiringDiagram, wire::Wire)
  g, wire_data = f.graph, WireData(wire)
  for e in edges(g, wire.source.box, wire.target.box)
    subpart(g, e, :wire) == wire_data && return rem_edge!(g, e)
  end
  error("Wire $wire does not exist, so cannot be removed")
end
rem_wire!(f::WiringDiagram, pair::Pair) = rem_wire!(f, Wire(pair))

function rem_wires!(f::WiringDiagram, wires)
  for wire in wires
    rem_wire!(f, wire)
  end
end

function rem_wires!(f::WiringDiagram, src::Int, tgt::Int)
  rem_edges!(f.graph, sort!(collect(edges(f.graph, src, tgt))))
end

""" Check compatibility of source and target ports.

The default implementation is a no-op.
"""
function validate_ports(source_port, target_port) end

# Graph properties.

""" Retrieve the graph underlying the wiring diagram.

The graph is an instance of `Catlab.Graphs.AbstractGraph`. Do not mutate it! All
mutations should use the wiring diagrams API: `add_box!`, `rem_box!`, and so on.
"""
graph(diagram::WiringDiagram) = diagram.graph

# Convenience methods delegated to underlying graph.
all_neighbors(d::WiringDiagram, v::Int) = all_neighbors(graph(d), v)
neighbors(d::WiringDiagram, v::Int) = neighbors(graph(d), v)
outneighbors(d::WiringDiagram, v::Int) = outneighbors(graph(d), v)
inneighbors(d::WiringDiagram, v::Int) = inneighbors(graph(d), v)

""" Get all wires coming into or out of the box.
"""
function wires(d::WiringDiagram, v::Int)
  g = graph(d)
  [ Wire(subpart(g, e, :wire), src(g, e), tgt(g, e))
    for e in unique!(sort!([incident(g, v, :src); incident(g, v, :tgt)])) ]
end

""" Get all wires coming into the box.
"""
function in_wires(d::WiringDiagram, v::Int)
  g = graph(d)
  [ Wire(subpart(g, e, :wire), src(g, e), v) for e in incident(g, v, :tgt) ]
end

""" Get all wires coming into the port.
"""
function in_wires(d::WiringDiagram, port::Port)
  filter(wire -> wire.target == port, in_wires(d, port.box))
end
function in_wires(d::WiringDiagram, v::Int, port::Int)
  in_wires(d, Port(v, InputPort, port))
end

""" Get all wires coming out of the box.
"""
function out_wires(d::WiringDiagram, v::Int)
  g = graph(d)
  [ Wire(subpart(g, e, :wire), v, tgt(g, e)) for e in incident(g, v, :src) ]
end

""" Get all wires coming out of the port.
"""
function out_wires(d::WiringDiagram, port::Port)
  filter(wire -> wire.source == port, out_wires(d, port.box))
end
function out_wires(d::WiringDiagram, v::Int, port::Int)
  out_wires(d, Port(v, OutputPort, port))
end

# Other constructors
#-------------------

""" Wiring diagram with a single box connected to the outer ports.
"""
function singleton_diagram(T::Type, box::AbstractBox)
  inputs, outputs = input_ports(box), output_ports(box)
  d = WiringDiagram{T}(inputs, outputs)
  v = add_box!(d, box)
  add_wires!(d, ((input_id(d),i) => (v,i) for i in eachindex(inputs)))
  add_wires!(d, ((v,i) => (output_id(d),i) for i in eachindex(outputs)))
  return d
end
singleton_diagram(box::AbstractBox) = singleton_diagram(Any, box)

""" The wiring diagram induced by a subset of its boxes.

See also `encapsulated_subdiagram`.
"""
function induced_subdiagram(d::WiringDiagram{T}, vs::Vector{Int}) where T
  sub = WiringDiagram{T}(input_ports(d), output_ports(d))
  vmap = Dict(input_id(d) => input_id(sub), output_id(d) => output_id(sub))
  for v in vs
    vmap[v] = add_box!(sub, box(d, v))
  end
  for wire in wires(d)
    src, tgt = wire.source, wire.target
    if haskey(vmap, src.box) && haskey(vmap, tgt.box)
      add_wire!(sub,
        Wire(set_box(src, vmap[src.box]), set_box(tgt, vmap[tgt.box])))
    end
  end
  return sub
end

# Operadic interface
####################

""" Operadic composition of wiring diagrams.

This generic function has two different signatures, corresponding to the "full"
and "partial" notions of operadic composition (Yau, 2018, *Operads of Wiring
Diagrams*, Definitions 2.3 and 2.10).

This operation is a simple wrapper around [`substitute`](@ref).
"""
function ocompose(f::WiringDiagram, gs::Vector{<:WiringDiagram})
  @assert length(gs) == nboxes(f)
  substitute(f, box_ids(f), gs)
end
function ocompose(f::WiringDiagram, i::Int, g::WiringDiagram)
  @assert 1 <= i <= nboxes(f)
  substitute(f, box_ids(f)[i], g)
end

# Substitution
##############

""" Substitute wiring diagrams for boxes.

Performs one or more substitutions. When performing multiple substitutions, the
substitutions are simultaneous.

This operation implements the operadic composition of wiring diagrams, see also
[`ocompose`](@ref).
"""
function substitute(d::WiringDiagram; kw...)
  substitute(d, filter(v -> box(d,v) isa WiringDiagram, box_ids(d)); kw...)
end
function substitute(d::WiringDiagram, v::Int; kw...)
  substitute(d, v, box(d,v)::WiringDiagram; kw...)
end
function substitute(d::WiringDiagram, vs::AbstractVector{Int}; kw...)
  substitute(d, vs, WiringDiagram[ box(d,v) for v in vs ]; kw...)
end
function substitute(d::WiringDiagram, v::Int, sub::WiringDiagram; kw...)
  substitute(d, [v], [sub]; kw...)
end

function substitute(d::WiringDiagram{T}, vs::AbstractVector{Int},
                    subs::Vector{<:WiringDiagram};
                    merge_wire_values=default_merge_wire_values) where T
  # In outline, the algorithm is:
  #
  # 1. Create an empty wiring diagram.
  # 2. Add *all* boxes of original diagram and the diagrams to be substituted
  #    (in the appropriate order).
  # 3. Add *all* wires of original diagram and the diagrams to be substituted.
  # 4. Remove the boxes that were substituted (hence also removing extraneous
  #    wires from step 3).
  #
  # This may seem convoluted, but it is the simplest way I know to handle the
  # problem of *instantaneous wires*. Some authors ban instantaneous wires, but
  # we need them to represent identities, braidings, etc.
  @assert length(vs) == length(subs)
  result = WiringDiagram{T}(d.value, input_ports(d), output_ports(d))

  # Add boxes by interleaving, in the correct order, the non-substituted boxes
  # of the original diagram and the internal boxes of the substituted diagrams.
  # At the very end, add the substituted boxes too.
  vmap = Dict(input_id(d) => input_id(result), output_id(d) => output_id(result))
  sub_maps = Dict(zip(vs, ((sub, Dict{Int,Int}()) for sub in subs)))
  for v in box_ids(d)
    if haskey(sub_maps, v)
      sub, sub_map = sub_maps[v]
      for u in box_ids(sub)
        sub_map[u] = add_box!(result, box(sub, u))
      end
    else
      vmap[v] = add_box!(result, box(d, v))
    end
  end
  for v in vs
    vmap[v] = add_box!(result, box(d, v))
  end

  # Add the wires of the original diagram, then add the internal wires of the
  # substituted diagrams.
  for wire in wires(d)
    add_wire!(result, Wire(wire.value,
      set_box(wire.source, vmap[wire.source.box]),
      set_box(wire.target, vmap[wire.target.box])))
  end
  for v in vs
    substitute_wires!(result, vmap[v], sub_maps[v]...;
      merge_wire_values=merge_wire_values)
  end

  # Finally, remove the substituted boxes. Because they were added last, this
  # will not change the IDs of the other boxes.
  rem_boxes!(result, (vmap[v] for v in vs))
  result
end

""" Substitute wires of sub-diagram into containing wiring diagram.
"""
function substitute_wires!(d::WiringDiagram, v::Int,
                           sub::WiringDiagram, sub_map::Dict{Int,Int};
                           merge_wire_values=default_merge_wire_values)
  for wire in wires(sub)
    src = get(sub_map, wire.source.box, 0)
    tgt = get(sub_map, wire.target.box, 0)
    # Special case: wire from input port to output port.
    if wire.source.box == input_id(sub) && wire.target.box == output_id(sub)
      for in_wire in in_wires(d, v, wire.source.port)
        for out_wire in out_wires(d, v, wire.target.port)
          add_wire!(d, Wire(
            merge_wire_values(in_wire.value, wire.value, out_wire.value),
            in_wire.source, out_wire.target))
        end
      end
    # Special case: wire from input port to internal box.
    elseif wire.source.box == input_id(sub)
      for in_wire in in_wires(d, v, wire.source.port)
        add_wire!(d, Wire(
          merge_wire_values(in_wire.value, wire.value, nothing),
          in_wire.source, set_box(wire.target, tgt)))
      end
    # Special case: wire from internal box to output port.
    elseif wire.target.box == output_id(sub)
      for out_wire in out_wires(d, v, wire.target.port)
        add_wire!(d, Wire(
          merge_wire_values(nothing, wire.value, out_wire.value),
          set_box(wire.source, src), out_wire.target))
      end
    # Default case: wire between two internal boxes.
    else
      add_wire!(d, Wire(
        merge_wire_values(nothing, wire.value, nothing),
        set_box(wire.source, src), set_box(wire.target, tgt)))
    end
  end
end

default_merge_wire_values(::Any, middle::Any, ::Any) = middle

# Encapsulation
###############

""" Encapsulate multiple boxes within a single sub-diagram.

This operation is a (one-sided) inverse to subsitution, see
[`substitute`](@ref).
"""
function encapsulate(d::WiringDiagram, vs::Vector{Int}; value=nothing, kw...)
  encapsulate(d, [vs]; values=[value], kw...)
end

function encapsulate(d::WiringDiagram{T}, vss::Vector{Vector{Int}};
    discard_boxes::Bool=false, make_box=Box, values=nothing) where T
  if isempty(vss); return d end
  if any(isempty(vs) for vs in vss)
    error("Cannot encapsulate an empty set of boxes")
  end
  if !allunique(reduce(vcat, vss))
    error("Cannot encapsulate overlapping sets of boxes")
  end
  if isnothing(values)
    values = repeat([nothing], length(vss))
  end
  result = WiringDiagram{T}(d.value, input_ports(d), output_ports(d))

  # Add boxes, both encapsulated and non-encapsulated, to new wiring diagram.
  encapsulated_representatives = Dict(
    minimum(vs) => (vs, value) for (vs, value) in zip(vss, values))
  all_encapsulated = Set(v for vs in vss for v in vs)
  vmap = Dict(input_id(d) => input_id(result), output_id(d) => output_id(result))
  port_map = Dict{Port,Port}()
  for v in box_ids(d)
    if haskey(encapsulated_representatives, v)
      vs, value = encapsulated_representatives[v]
      sub, sub_map = encapsulated_subdiagram(d, vs;
        discard_boxes=discard_boxes, make_box=make_box, value=value)
      subv = add_box!(result, sub)
      merge!(port_map, Dict(port => Port(data, subv)
                            for (port, data) in sub_map))
    elseif v ∉ all_encapsulated
      vmap[v] = add_box!(result, box(d, v))
    end
  end

  # Add wires to new wiring diagram.
  for wire in wires(d)
    src, tgt = wire.source, wire.target
    new_src = if haskey(vmap, src.box); set_box(src, vmap[src.box])
      elseif haskey(port_map, src); port_map[src]
      else; continue end
    new_tgt = if haskey(vmap, tgt.box); set_box(tgt, vmap[tgt.box])
      elseif haskey(port_map, tgt); port_map[tgt]
      else; continue end
    add_wire!(result, Wire(new_src, new_tgt))
  end
  result
end

""" Create an encapsulating box for a set of boxes in a wiring diagram.

To a first approximation, the union of input ports of the given boxes will
become the inputs ports of the encapsulating box and likewise for the output
ports. However, when copies or merges occur, as in a cartesian or cocartesian
category, a simplification procedure may reduce the number of ports on the
encapsulating box.

Specifically:

1. Each input port of an encapsulated box will have at most one incoming wire
from the encapsulating outer box, and each output port of an encapsulated box
will have at most one outgoing wire to the encapsulating outer box.

2. A set of ports connected to the same outside (non-encapsulated) ports will be
simplified into a single port of the encapsulating box.

See also `induced_subdiagram`.
"""
function encapsulated_subdiagram(d::WiringDiagram{T}, vs::Vector{Int};
    discard_boxes::Bool=false, make_box=Box, value=nothing) where T
  # Add encapsulated box to new wiring diagram.
  inputs, outputs = [], []
  result = discard_boxes ? nothing : WiringDiagram{T}(value, inputs, outputs)
  vmap = if discard_boxes
    Dict(v => nothing for v in vs)
  else
    Dict(v => add_box!(result, box(d, v)) for v in vs)
  end

  # Process wires into, or out of, encapsulated boxes.
  port_map = Dict{Port,PortData}()
  inner_port_map = Dict{Tuple{Vector{Port},Any},Int}()
  for v in vs
    # Add input ports to encapsulating box and corresponding wire.
    for (port, value) in enumerate(input_ports(d, v))
      tgt = Port(v, InputPort, port)
      srcs = sort!([ wire.source for wire in in_wires(d, tgt)
                     if !haskey(vmap, wire.source.box) ])
      if isempty(srcs) continue end
      src = get!(inner_port_map, (srcs, value)) do
        push!(inputs, value)
        port_data = port_map[tgt] = PortData(InputPort, length(inputs))
        port_data.port
      end
      if discard_boxes; continue end
      add_wire!(result,
        Wire(Port(input_id(result), OutputPort, src), set_box(tgt, vmap[v])))
    end

    # Add output ports to encapsulating box and corresponding wire.
    for (port, value) in enumerate(output_ports(d, v))
      src = Port(v, OutputPort, port)
      tgts = sort([ wire.target for wire in out_wires(d, src)
                    if !haskey(vmap, wire.target.box) ])
      if isempty(tgts) continue end
      tgt = get!(inner_port_map, (tgts, value)) do
        push!(outputs, value)
        port_data = port_map[src] = PortData(OutputPort, length(outputs))
        port_data.port
      end
      if discard_boxes; continue end
      add_wire!(result,
        Wire(set_box(src, vmap[v]), Port(output_id(result), InputPort, tgt)))
    end

    # Add wires between two encapsulated boxes.
    if discard_boxes; continue end
    for wire in out_wires(d, v)
      src, tgt = wire.source, wire.target
      if haskey(vmap, src.box) && haskey(vmap, tgt.box) # Clause #1 always true.
        add_wire!(result,
          Wire(set_box(src, vmap[src.box]), set_box(tgt, vmap[tgt.box])))
      end
    end
  end

  # Yield input and output port lists with the tightest possible types.
  inputs, outputs = [ x for x in inputs ], [ x for x in outputs ]
  if discard_boxes
    result = make_box(value, inputs, outputs)
  else
    result.input_ports, result.output_ports = inputs, outputs
  end
  (result, port_map)
end

end


text:
#!/usr/bin/env julia

const dir = abspath( splitdir(@__FILE__)[1] )

tic()
println( STDERR, "EvoDuplexes.jl loading and compiling... " )

using ArgParse
using BufferedStreams
using Libz
using GenomicFeatures
using BioSequences
using DataFrames
using Gadfly
using Distances
using Clustering
using MultivariateStats
using StatsBase

unshift!( LOAD_PATH, dir * "/../src" )
using EvoDuplexes
using EvoDuplexes: score, score!

function parse_cmd()
   s = ArgParseSettings()
   @add_arg_table s begin
      "--tree"
       help     = "Phylogenetic tree with neutral branch lengths (in newick format)"
       arg_type = String
       required = true
      "--cons-regions"
       help     = "BED file containing conserved regions"
       arg_type = String
       required = false
      "--gene-regions"
       help     = "BED file containing gene units to allow long-range folding within"
       arg_type = String
       required = true
      "--maf"
       help     = "Directory with MAF files named by chromosome (chr1.maf.gz...)"
       arg_type = String
       required = true
       default  = "../maf"
      "--model-load"
       help     = "Load pre-trained IsolationForest models, .evt.jls"
       arg_type = String
       required = false
      "--model-data"
       help     = "Load training data from `.jlt` file, output `.evt.jls` file"
       arg_type = String
       required = false
      "--model-train"
       help     = "Train the IsolationForest models from scratch using input regions, output `.evt.jls` file"
       action   = :store_true
      "--output","-o"
       help     = "Output prefix"
       arg_type = String
       default  = "evoduplex"
      "--cons-branch"
       help     = "When conserved regions are given, use this branch length multiplier"
       arg_type = Float64
       default  = 0.33
      "--sig-ratio"
       help     = "Output the top fraction of outliers"
       arg_type = Float64
       default  = 0.05
      "--output-all"
       help     = "Output all folds regardless of prediction status"
       action   = :store_true
      "--max-distance"
       help     = "Set a limit on the maximum distance between a left/right arm of a duplex"
       arg_type = Int64
       default  = 2000
      "--max-duplex-length"
       help     = "Set the maximum duplex length that will be accessible from the suffix array"
       arg_type = Int64
       default  = 50
      "--max-bulges"
       help     = "Set the maximum number of bulges to allow for a duplex"
       arg_type = Int64
       default  = 3
      "--max-mismatches"
       help     = "Set the maximum number of mismatches to allow for a duplex"
       arg_type = Int64
       default  = 3
      "--max-deltag"
       help     = "Set the maximum allowable deltaG value"
       arg_type = Float64
       default  = -0.8
    end
   return parse_args(s)
end

function main()

   istrue( arg ) = (arg != nothing)

   args = parse_cmd()

   println(STDERR, " $( round( toq(), 6 ) ) seconds" )

   neutraltree = parsenewick(String(chomp(String(read(open(args["tree"]))))))
   pairtree = deepcopy(neutraltree)
   constree = deepcopy(neutraltree)

   extend_branches!( neutraltree, 1.0 )
   extend_branches!( constree,    args["cons-branch"] )
   extend_branches!( pairtree,    args["cons-branch"] )

   set_prob_mat!( neutraltree, GTR_SINGLE_Q )
   set_prob_mat!( constree,    GTR_SINGLE_Q )
   set_prob_mat!( pairtree,    GTR_PAIRED_Q )

   secbool   = false
   conserved = false
   if istrue( args["cons-regions"] )
      regs  = loadbed( args["cons-regions"], expandfirst=25, expandlast=25 )
      conserved = true
      if istrue( args["gene-regions"] )
         genes   = loadbed( args["gene-regions"] )#, expandfirst=2500, expandlast=2500 )
         secbool = true
      else
         error("No specific `--gene-regions` given! Can't fold entire chromosomes..")
      end
   elseif istrue( args["gene-regions"] )
      genes = loadbed( args["gene-regions"] )
      regs  = genes
   else
      error("No specific regions given! Can't fold entire chromosomes..")
   end

   train = false
   if istrue( args["model-load"] )
      println("Loading model from $(args["model-load"]) file...")
      model = open(deserialize, args["model-load"])
   elseif istrue( args["model-data"] )
      println("Training model from .jlt datafile $(args["model-data"])...")
      tabdata = open(readdlm, args["model-data"])
      model = DistanceForest( args["sig-ratio"] )
      train!( model, tabdata )
   elseif args["model-train"]
      train = true
   else
      error("No --model-file or --model-data given!, if training is desired use --model-train!")
   end

   datout = BufferedOutputStream(open("$(args["output"])-raw.jlt", "w"))
   bedout = open("$(args["output"]).bed", "w")
   duplexes = Vector{DuplexInterval}()

   for chr in keys(genes.trees)
      println("Reading $(args["maf"])/$chr.maf.gz")
      reader = ZlibInflateInputStream(open("$(args["maf"])/$chr.maf.gz"), bufsize=129000) |> MAFReader
      @time maf = readmaf!( reader, neutraltree.index, minlength=18, minspecies=2, minscore=-Inf, regionbool=true, regions=genes, secbool=secbool, secregions=regs )

      println(maf)

      for r in collect(genes)
         haskey( maf.trees, r.seqname ) || continue
         col = collect(intersect( maf.trees[r.seqname], r ))
         length(col) >= 1 || continue
         if r.strand == STRAND_NEG
            for i in col
               EvoDuplexes.reverse_complement!( i.metadata )
            end
         end
         println("Building DuplexArray for $(r.metadata)")
         rda = RNADuplexArray{DNAAlphabet{2},UInt16,UInt16}( col[1].metadata, neutraltree, args["max-duplex-length"] )
         for i in 2:length(col)
            push!( rda, col[i].metadata, neutraltree )
         end

         @time trav = collect(traverse( rda, 1:args["max-distance"],
                                         single=neutraltree,
                                         bulge_max=args["max-bulges"],
                                         mismatch_max=args["max-mismatches"],
                                         minfold=args["max-deltag"] ))

         println("Calculating structural conservation across the phylogeny...")
         for i in trav
            #const conscore = (score( i.duplex, constree, gapdenom=1.25 ) - score( i.duplex, neutraltree )) / npairs( i.duplex )
            if !conserved
               score!( i.duplex, neutraltree, pairtree )
               const cons = score( i.duplex ) / npairs( i.duplex )
            else
               score!( i.duplex, constree, pairtree )
               const cons = score( i.duplex ) / npairs( i.duplex )
            end
            const ent     = dinucleotide_entropy(i.duplex.duplex) / (npairs( i.duplex ) + 5)
            const meancov = covariance( i.duplex ) * 10
            const len     = npairs( i.duplex )
            const dist    = distance( i )
            const deltag  = energy( i.duplex )

            write( datout, string(dist) ); write( datout, '\t' )
            write( datout, string(cons) ); write( datout, '\t' )
            write( datout, string(deltag) ); write( datout, '\t' )
            write( datout, string(len) ); write( datout, '\t' )
            write( datout, string(ent) ); write( datout, '\t' )
            write( datout, string(meancov) ); write( datout, '\n' )
         end

         duplexes = vcat(duplexes, trav)

      end
   end

   full_char = [collect('0':'9'); collect('A':'Z')]

   close(datout)
   tabdata = open(readdlm, "$(args["output"])-raw.jlt")

   if train
      println("Training IsolationForest models...")
      model = DistanceForest( args["sig-ratio"] )
      train!( model, tabdata )
   end

   println("Identifying statistical anomalies...")
   if !args["output-all"]
      sig = EvoDuplexes.predict( model, tabdata )
   else
      sig = collect(1:length(duplexes))
   end

   println("Found $(length(sig)) significant duplexes..")

   for i in sig
      name = randcode(full_char)
      writebed( bedout, duplexes[i], name )
      s,b = brackets( duplexes[i].duplex )
      println("$name: $b")
      println("$name: $s")
   end

   modelhndl = open( args["output"], "w" )
   println("Serializing model to $( args["output"] ).evt.jls")
   serialize(modelhndl, model)
   close(modelhndl)

   close(bedout)
end

@time main()

ext: jl


text:
######################################################################
# dipole.jl
#
# allows user to easily make dipole antennas
######################################################################

struct Dipole <: AbstractAntenna
    wire::Wire
end

function Dipole(fmhz::Real, rad::Real=0.001; n_lambdas::Real=0.5)
    lambda = get_lambda(fmhz)
    L = n_lambdas * lambda       # length of

    a = (0,0,-L/2)
    b = (0,0,L/2)
    n_segments = divide(L, fmhz)

    return Dipole(Wire(a, b, rad, n_segments))
end

function VoltageSource(dipole::Dipole; v::Float64=1.0)
    return VoltageSource(dipole.wire, middle(dipole.wire), v)
end

function RadiationPattern(d::Dipole, fmhz::Real)
    return RadiationPattern([d.wire], VoltageSource(d), fmhz)
end


text:
import Printf: ini_dec, fix_dec, ini_hex, ini_HEX

if VERSION < v"1.1"
using Base.Grisu: DIGITSs
fix_dec(out, d::Float128, flags::String, width::Int, precision::Int, c::Char) = fp128_printf(out, d, flags, width, precision, c, DIGITSs[Threads.threadid()])
ini_dec(out, d::Float128, ndigits::Int, flags::String, width::Int, precision::Int, c::Char) = fp128_printf(out, d, flags, width, precision, c, DIGITSs[Threads.threadid()])
ini_hex(out, d::Float128, ndigits::Int, flags::String, width::Int, precision::Int, c::Char) = fp128_printf(out, d, flags, width, precision, c, DIGITSs[Threads.threadid()])
ini_HEX(out, d::Float128, ndigits::Int, flags::String, width::Int, precision::Int, c::Char) = fp128_printf(out, d, flags, width, precision, c, DIGITSs[Threads.threadid()])
ini_hex(out, d::Float128, flags::String, width::Int, precision::Int, c::Char) = fp128_printf(out, d, flags, width, precision, c, DIGITSs[Threads.threadid()])
ini_HEX(out, d::Float128, flags::String, width::Int, precision::Int, c::Char) = fp128_printf(out, d, flags, width, precision, c, DIGITSs[Threads.threadid()])
else
fix_dec(out, d::Float128, flags::String, width::Int, precision::Int, c::Char, digits) = fp128_printf(out, d, flags, width, precision, c, digits)
ini_dec(out, d::Float128, ndigits::Int, flags::String, width::Int, precision::Int, c::Char, digits) = fp128_printf(out, d, flags, width, precision, c, digits)
ini_hex(out, d::Float128, ndigits::Int, flags::String, width::Int, precision::Int, c::Char, digits) = fp128_printf(out, d, flags, width, precision, c, digits)
ini_HEX(out, d::Float128, ndigits::Int, flags::String, width::Int, precision::Int, c::Char, digits) = fp128_printf(out, d, flags, width, precision, c, digits)
ini_hex(out, d::Float128, flags::String, width::Int, precision::Int, c::Char, digits) = fp128_printf(out, d, flags, width, precision, c, digits)
ini_HEX(out, d::Float128, flags::String, width::Int, precision::Int, c::Char, digits) = fp128_printf(out, d, flags, width, precision, c, digits)
end

function fp128_printf(out, d::Float128, flags::String, width::Int, precision::Int, c::Char, digits)
    fmt_len = sizeof(flags)+4
    if width > 0
        fmt_len += ndigits(width)
    end
    if precision >= 0
        fmt_len += ndigits(precision)+1
    end
    fmt = IOBuffer(maxsize=fmt_len)
    print(fmt, '%')
    print(fmt, flags)
    if width > 0
        print(fmt, width)
    end
    if precision == 0
        print(fmt, '.')
        print(fmt, '0')
    elseif precision > 0
        print(fmt, '.')
        print(fmt, precision)
    end
    print(fmt, 'Q')
    print(fmt, c)
    write(fmt, UInt8(0))
    printf_fmt = take!(fmt)
    @assert length(printf_fmt) == fmt_len
    bufsiz = length(digits)
    lng = @ccall(libquadmath.quadmath_snprintf(digits::Ptr{UInt8}, bufsiz::Csize_t, printf_fmt::Ptr{UInt8}, d::(Cfloat128...))::Cint)
    lng > 0 || error("invalid printf formatting for Float128")
    unsafe_write(out, pointer(digits), min(lng, bufsiz-1))
    return (false, ())
end


text:
# Using Mosek as the SDP solver
model = SOSModel(Mosek.Optimizer)

max_stars_repo_name: sums-of-squares/tools

max_issues_repo_name: sums-of-squares/tools

max_forks_repo_name: sums-of-squares/tools


text:
let mode, n, A, p, k, A
# test two properties of oscillating matrix
# For an n x n oscillating matrix A, the following holds:
# 1. A has n distinct and postive eigenvalues λ_1 > λ_2 > ... > λ_n > 0
# 2. The ith eigenvector, corresponding to λ_i in the above ordering has
#    exactly i-1 sign changes.

n = 8 # rand(2:10)
@test_throws ArgumentError matrixdepot("oscillate", n, 4)
for mode = 1:2
A = matrixdepot("oscillate", n, mode)

eva, evc = eigen(A)
@test all([ei > 0 for ei in eva])

p = sortperm(eva, rev = true)
evc = evc[:,p]

# compute the number of sign change
function num_sign_change(v)
    signv = sign.(v)
    change = signv[1]
    num = 0
    for i in signv
        if i != change
            num += 1
            change = i
        end
    end
    return num
end

k = rand(1:n)
@test num_sign_change(evc[:,k]) == k - 1
end

@test matrixdepot("oscillate", n) !== nothing
end
println("'oscillate' passed test...")


text:
mutable struct Cuboid
  x :: Tuple{Int,Int}
  y :: Tuple{Int,Int}
  z :: Tuple{Int,Int}
end

mutable struct Instruction
  state  :: Bool
  cuboid :: Cuboid
end

function apply(I::Instruction, M::Array{Bool,3})
  M[I.x[1]:I.x[2],I.y[1]:I.y[2],I.z[1]:I.z[2]] = I.state
end

function parseCuboid(s::String) :: Cuboid
  return Cuboid(map(y -> (parse(Int,y[1]),parse(Int,y[2])), map(x->split(x[3:end],".."),split(s,',')))...)
end

function parseLine(line::String) :: Instruction
  parts = split(line, ' ')
  state  :: Bool    = (parts[1] == "on" ? true : false)
  cuboid :: Cuboid = parseCuboid(String(parts[2]))
  return Instruction(state,cuboid)
end

function readInput() :: Vector{Instruction}
  f = open("inputs/day22.txt")
  data = map(line -> parseLine(line), eachline(f))
  return data
end

data = readInput()


text:
using Glob

function add_icc(filename,
                 icc_prof = joinpath(dirname(@__FILE__), "sRGB_v4_ICC_preference.icc"))
    run(`convert $(filename) -profile $(icc_prof) $(filename)`)
end

function compile_pgf(filebase)
    w,h = [gcf()[s]() for s in
           [:get_figwidth,:get_figheight]]
    latex_source = """\\documentclass[12pt]{minimal}
\\usepackage[paperwidth=$(w)in, paperheight=$(h)in, margin=0in]{geometry}
\\usepackage{pdfx}
$(pgf_backend[:get_preamble]())
$(pgf_backend[:get_fontspec]())
\\usepackage{pgf}
\\begin{document}
\\centering
\\input{$(filebase).pgf}
\\end{document}"""
    open("$(filebase).tex", "w") do file
        write(file, latex_source)
    end
    tex_cmd = pgf_backend.get_texcommand()
    args = ["-interaction=nonstopmode", "-halt-on-error"]
    if tex_cmd in ["xetex", "xelatex"]
        println("WARNING: XeTeX is problematic")
        push!(args, "-shell-escape")
        push!(args, "-output-driver=xdvipdfmx -z 0")
    end
    run(`$(tex_cmd) $(args) $(filebase).tex`)
end

function save_pgf_with_icc(filename, args...;
                           icc_prof = joinpath(dirname(@__FILE__), "sRGB_v4_ICC_preference.icc"),
                           compile_pdf = true, kwargs...)
    savefig(filename, args...; kwargs...)
    dir,file,ext=dirname(filename),splitext(filename)...
    imgs = glob("$(file)-img*.png", dir)
    length(imgs) > 0 && println("Applying ICC profile to $(length(imgs)) image(s)")
    for img in imgs
        println(img)
        add_icc(img, icc_prof)
    end
    compile_pdf && compile_pgf(joinpath(dir, file))
end

export add_icc, save_pgf_with_icc


text:
using LoopVectorization
using BenchmarkTools

function set_ghost_cells_bot_kernel!(
    a, a_bot, a_gradbot, dzh, ks)

    # fac = 0.5f0 # This works.
    fac = 1//2 # This does not work.

    @tturbo for j in 1:size(a, 2)
        for i in 1:size(a, 1)
            a[i, j, ks-1] = -a_gradbot[i, j]*dzh[ks] + a[i, j, ks]
            a_bot[i, j] = fac * (a[i, j, ks-1] + a[i, j, ks])
        end
    end
end

TF = Float32
# TF = Float64

a = rand(TF, 50, 50, 50)
a_bot = rand(TF, 50, 50)
a_gradbot = rand(TF, 50, 50)
dzh = rand(TF, 50)
ks = 2

@btime set_ghost_cells_bot_kernel!($a, $a_bot, $a_gradbot, $dzh, $ks)


text:
"""
From sandreza/Learning/sandbox/oceananigans_converter.jl
  https://github.com/sandreza/Learning/blob/master/sandbox/oceananigans_converter.jl

custom_avg(Φ, n)
# Description
- Average a field down to n.
- Requires field to have evenly spaced points. Size of N leq length(Φ).
- Furthermore requires
# Arguments
- `Φ` :(vector) The field, an array
- `n` :(Int) number of grid points to average down to.
# Return
- `Φ2` :(vector) The field with values averaged, an array
"""
function custom_avg(Φ, n)
  m = length(Φ)
  scale = Int(floor(m/n))
  if ( abs(Int(floor(m/n)) - m/n) > eps(1.0))
      return error
  end
  Φ2 = zeros(n)
  for i in 1:n
      Φ2[i] = 0
          for j in 1:scale
              Φ2[i] += Φ[scale*(i-1) + j] / scale
          end
  end
  return Φ2
end

max_stars_repo_name: adelinehillier/LearnConvection

max_issues_repo_name: adelinehillier/LearnConvection

max_forks_repo_name: adelinehillier/LearnConvection


text:
Dict(
    "hello" => "c'a va latn"
   )


text:
# ROUGE score implementation
#Lin C.Y. , 2004
#Rouge: A package for automatic evaluation of summaries
#Proceedings of the workshop on text summarization branches out (WAS 2004) (2004), pp. 25-26
#Link to paper:
#http://www.aclweb.org/anthology/W04-1013

#It is a n-gram recall between a candidate summary
#and a set of reference summaries.

# param references : list of reference strings
#type references : Array{String,1}

# param candidate :  the candidate string
#type (candidate) : Array{String,1}

# param n : length of ngram
#type (n) : int

#ngram_cand : list of ngrams in candidate
#ngram_ref : list of ngrams in reference
#r_lcs : recall factor
#p_lcs : precision factor
#rouge_recall : list containing all the rouge-n scores for
#               every reference against the candidate

function rouge_n(references, candidate, n, averaging = true)

    ngram_cand = listify_ngrams(ngrams(StringDocument(candidate), n))
    rouge_recall = []

    for ref in references
        matches = 0  #variable counting the no.of matching ngrams
        ngram_ref = listify_ngrams(ngrams(StringDocument(ref), n))
        print(ngram_ref)

        for ngr in ngram_cand
            if ngr in ngram_ref
                matches += 1
            end

        end

        push!(rouge_recall, matches/length(ngram_ref))

    end

    if averaging == true
        rouge_recall = jackknife_avg(rouge_recall)
    end

    return(rouge_recall)
end

# It calculates the rouge-l score between the candidate
#and the reference at the sentence level.

# param references : list of reference strings
#type references : Array{String,1}

# param candidate :  the candidate string
#type (candidate) : Array{String,1}

# param beta : user-defined parameter. Default value = 8
#type (beta) : float

#rouge_l_list : list containing all the rouge scores for
#                every reference against the candidate
#r_lcs : recall factor
#p_lcs : precision factor
#score : rouge-l score between the reference sentence and
#        the candidate sentence

function rouge_l_sentence(references, candidate, beta=8, averaging = true)

    ngram_cand = tokenize(Languages.English(), candidate)
    rouge_l_list = []

    for ref in references
        ngram_ref = tokenize(Languages.English(), ref)
        r_lcs = weighted_lcs(ngram_ref, ngram_cand,true, false, sqrt)/length(ngram_ref)
        p_lcs = weighted_lcs(ngram_ref, ngram_cand,true, false, sqrt)/length(ngram_cand)
        score = fmeasure_lcs(r_lcs, p_lcs, beta)
        push!(rouge_l_list,score)

    end
    if averaging == true
        rouge_l_list = jackknife_avg(rouge_l_list)
    end
    return rouge_l_list
end

#It calculates the rouge-l score between the candidate
#and the reference at the summary level.
# param references : list of reference summaries. Each of the summaries
#                    must be tokenized list of words
#type (references) : list

# param candidate : candidate summary tokenized into list of words
#type (candidate) : list
# param beta : user-defined parameter
#type (beta) : float

#rouge_l_list : list containing all the rouge scores for
#                every reference against the candidate

#r_lcs : recall factor
#p_lcs : precision factor
#score : rouge-l score between a reference and the candidate

function rouge_l_summary(references, candidate, beta, averaging=true)

    rouge_l_list = []
    cand_sent_list = split_sentences(candidate)

    for ref in references
        ref_sent_list = split_sentences(ref)
        sum_value = 0

        for ref_sent in ref_sent_list
            l_ = []
            arg1 = tokenize(Languages.English(), ref)

            for cand_sent in cand_sent_list
                arg2 = tokenize(Languages.English(), cand_sent)
                d = tokenize(Languages.English(), weighted_lcs(arg1, arg2, false, true, sqrt))
                append!(l_,d)
            end

            print(l_)
            sum_value = sum_value+length(unique(l_))

        end

        r_lcs = sum_value/length(tokenize(Languages.English(), ref))
        p_lcs = sum_value/length(tokenize(Languages.English(), candidate))
        score = fmeasure_lcs(r_lcs, p_lcs, beta)
        push!(rouge_l_list,score)

    end

    if averaging == true
        rouge_l_list = jackknife_avg(rouge_l_list)
    end

    return rouge_l_list
end


text:
"""
#### Get the simulation overview from RAMSES info, descriptor and output header files
```julia
getinfo(; output::Real=1, path::String="", namelist::String="", verbose::Bool=verbose_mode)
return InfoType
```

#### Keyword Arguments
- **`output`:** timestep number (default=1)
- **`path`:** the path to the output folder relative to the current folder or absolute path
- **`namelist`:** give the path to a namelist file (by default the namelist.txt-file in the output-folder is read)
- **`verbose:`:** informations are printed on the screen by default: gloval variable `verbose_mode=true`

#### Examples
```julia
# read simulation information from output `1` in current folder
julia> info = getinfo()

# read simulation information from output `420` in given folder (relative path to the current working folder)
julia> info = getinfo(output=420, path="../MySimFolder/")

# or simply use
julia> info = getinfo(420, "../MySimFolder/")

# get an overview of the returned field-names
julia> propertynames(info)

# a more detailed overview
julia> viewfields(info)
julia> viewallfields(info)
julia> namelist(info)
```

"""
function getinfo(output::Real; path::String="", namelist::String="", verbose::Bool=verbose_mode)
    return getinfo(output=output, path=path, namelist=namelist, verbose=verbose)
end

function getinfo(output::Real, path::String; namelist::String="", verbose::Bool=verbose_mode)
    return getinfo(output=output, path=path, namelist=namelist, verbose=verbose)
end

function getinfo(path::String; output::Real=1, namelist::String="", verbose::Bool=verbose_mode)
    return getinfo(output=output, path=path, namelist=namelist, verbose=verbose)
end

function getinfo(; output::Real=1, path::String="", namelist::String="", verbose::Bool=verbose_mode)

    printtime("",verbose)

    info = InfoType()   # predeclare InfoType
    info.descriptor = DescriptorType()
    info.files_content = FilesContentType()

    info.output     = output
    info.path       = joinpath(pwd(), path) # store full path
    info.simcode    = "RAMSES"
    info.Narraysize = 0

    createpath!(info, namelist=namelist)    # create path to the files
    readinfofile!(info)     # infofile overview
    createconstants!(info)  # predefined constants (create before scales)
    createscales!(info)     # get predefined scales and the corresponding cgs units

    readnamelistfile!(info) # read namelist, check for ASCII file

    readamrfile1!(info)     # amr overview
    readhydrofile1!(info)   # hydro overview
    isgravityfile1!(info)   # gravity overview
    readparticlesfile1!(info)   # particles overview
    readrtfile1!(info)      # rt overview
    readclumpfile1!(info)   # clumps overview

    # todo: check for sinks
    info.sinks = false
    info.sinks_variable_list = Symbol[]
    info.descriptor.sinks = Symbol[]
    info.descriptor.usesinks = false
    info.descriptor.sinksfile = false

    readtimerfile!(info)      # check for timer-file
    readcompilationfile!(info)  # compilation overview
    readmakefile!(info)       # check for makefile
    readpatchfile!(info)      # check for patchfile

    printsimoverview(info, verbose) # print overview on screen
    return info
end

# ============================================
# 'getinfo' related functions
# ============================================

function readinfofile!(dataobject::InfoType)
    if  !isfile(dataobject.fnames.info)
      println()
      error("[Mera]:  File or folder does not exist: $(dataobject.fnames.info) !")
    end
    mtime = Dates.unix2datetime(stat(dataobject.fnames.info).mtime)
    ctime = Dates.unix2datetime(stat(dataobject.fnames.info).ctime)

    f = open(dataobject.fnames.info)
    lines = readlines(f)

    ncpu      = parse(Int32, rsplit(lines[1],"=")[2] )
    ndim      = parse(Int32, rsplit(lines[2],"=")[2] )
    levelmin  = parse(Int32, rsplit(lines[3],"=")[2] )
    levelmax  = parse(Int32, rsplit(lines[4],"=")[2] )

    dataobject.grid_info = GridInfoType() #0,0,0,0,0,0,0,0,zeros(Float64, ncpu+1) ,falses(ncpu) )
    dataobject.grid_info.bound_key = zeros(Float64, ncpu+1)
    dataobject.grid_info.cpu_read = falses(ncpu)
    dataobject.grid_info.ngridmax  = parse(Int32, rsplit(lines[5],"=")[2] )
    dataobject.grid_info.nstep_coarse = parse(Int32, rsplit(lines[6],"=")[2] )

    boxlen    = parse(Float64, rsplit(lines[8],"=")[2] )
    time      = parse(Float64, rsplit(lines[9],"=")[2] )
    aexp      = parse(Float64, rsplit(lines[10],"=")[2] )
    H0        = parse(Float64, rsplit(lines[11],"=")[2] )
    omega_m   = parse(Float64, rsplit(lines[12],"=")[2] )
    omega_l   = parse(Float64, rsplit(lines[13],"=")[2] )
    omega_k   = parse(Float64, rsplit(lines[14],"=")[2] )
    omega_b   = parse(Float64, rsplit(lines[15],"=")[2] )

    unit_l    = parse(Float64, rsplit(lines[16],"=")[2] )
    unit_d    = parse(Float64, rsplit(lines[17],"=")[2] )
    unit_t    = parse(Float64, rsplit(lines[18],"=")[2] )
    unit_v    = unit_l / unit_t
    unit_m    = unit_d * unit_l^3

    hilbert_ordering = occursin(r"(?i)hilbert", lines[20])

    if  ndim != 3
      error("[Mera]: Program only works with 3D data!")
    end

    dataobject.grid_info.bound_key[1]   = parse( Float64, split(lines[22])[2] )
    for i = 1:ncpu
      #println( parse( Float64, split(lines[21+i])[1] ), " ", parse( Float64, split(lines[21+i])[2] ), " ",  parse( Float64, split(lines[21+i])[3] ) )
      dataobject.grid_info.bound_key[i+1] = parse( Float64, split(lines[21+i])[3] )
    end
    close(f)

    dataobject.ndim       = ndim
    #dataobject.grid_info  = grid_info
    dataobject.mtime      = mtime
    dataobject.ctime      = ctime
    dataobject.ncpu       = ncpu
    dataobject.levelmin   = levelmin
    dataobject.levelmax   = levelmax

    dataobject.boxlen     = boxlen
    dataobject.time       = time
    dataobject.aexp       = aexp
    dataobject.H0         = H0
    dataobject.omega_m    = omega_m
    dataobject.omega_l    = omega_l
    dataobject.omega_k    = omega_k
    dataobject.omega_b    = omega_b
    dataobject.unit_l     = unit_l
    dataobject.unit_d     = unit_d
    dataobject.unit_m     = unit_m
    dataobject.unit_v     = unit_v
    dataobject.unit_t     = unit_t

    return dataobject
end

function readamrfile1!(dataobject::InfoType)
    dataobject.amr = true

    if  !isfile(dataobject.fnames.amr * "out00001")
      dataobject.amr = false
      #println()
      #error("""[Mera]:  File or folder does not exist: $(dataobject.fnames.amr * "out00001")!""")
    else

        f = FortranFile(dataobject.fnames.amr * "out00001")
        #read(f) # ncpu    = read(f, Int32)
        #read(f) # ndim    = read(f, Int32)
        skiplines(f, 2)
        dataobject.grid_info.nx, dataobject.grid_info.ny, dataobject.grid_info.nz = read(f, Int32,Int32,Int32 )
        dataobject.grid_info.nlevelmax = read(f, Int32)
        dataobject.grid_info.ngridmax  = read(f, Int32)
        dataobject.grid_info.nboundary = read(f, Int32)
        dataobject.grid_info.ngrid_current = read(f, Int32)
        #info.boxlen = read(f, Float64)
        close(f)
    end
    return dataobject

end

function readhydrofile1!(dataobject::InfoType)

    # read descriptor file
    descriptor_file = false
    variable_descriptor_list = Symbol[]
    variable_types = String[]
    version = 0
    if  isfile(dataobject.fnames.hydro_descriptor)
        descriptor_file = true
        f = open(dataobject.fnames.hydro_descriptor)
        lines = readlines(f)

        # check for descriptor version
        if occursin("nvar", String(lines[1]))   # =< stable_17_09
            version = 0
        elseif occursin("version", String(lines[1])) # > stable_17_09
            version = parse(Int, rsplit(lines[1], ":" )[2])
        end

        # read descriptor variables
        if version == 0
            dnvar = parse(Int, rsplit(lines[1],"=")[2] )
            for i = 1:dnvar
                ivar = String(rsplit(lines[i+1], ":" )[2])
                ivar = strip(ivar)
                append!(variable_descriptor_list, [Symbol(ivar)])
            end
        elseif version == 1
            dnvar = length(lines)
            for i = 3:dnvar
                ivar = String(rsplit(lines[i], "," )[2])
                itype = String(rsplit(lines[i], "," )[3])

                ivar = strip(ivar)
                itype = strip(itype)
                append!(variable_descriptor_list, [Symbol(ivar)])
                append!(variable_types, [itype])
            end
        else
        # version not supported,
        # descriptor variables not read
        end
        close(f)
    end

    #read header from first cpu file
    hydro_files = false
    nvarh = 0
    variable_list = Symbol[]
    if  isfile(dataobject.fnames.hydro * "out00001")
        f_hydro = FortranFile(dataobject.fnames.hydro * "out00001")
        skiplines(f_hydro, 1)
        nvarh = read(f_hydro, Int32)
        skiplines(f_hydro, 3)
        ###skiplines(f_hydro, 3)
        dataobject.gamma = read(f_hydro, Float64)
        close(f_hydro)

        variable_list = [:rho,:vx,:vy,:vz,:p]
        if nvarh > 5
            for ivar=6:nvarh
                append!(variable_list, [Symbol("var$ivar")])
            end
        end
       hydro_files = true
    end

    if !isfile(dataobject.fnames.hydro_descriptor)
        variable_descriptor_list = variable_list
    end

    dataobject.hydro            = hydro_files
    dataobject.nvarh            = nvarh
    dataobject.variable_list    = variable_list

    # descriptor
    dataobject.descriptor.hversion      = version
    dataobject.descriptor.hydro         = variable_descriptor_list
    dataobject.descriptor.htypes        = variable_types
    dataobject.descriptor.usehydro      = false
    dataobject.descriptor.hydrofile     = descriptor_file
    return dataobject
end

# in work
function readrtfile1!(dataobject::InfoType)

    # read descriptor file
    descriptor_file = false
    rtPhotonGroups = Dict()
    descriptor_list= Dict()
    version = 0

    if  isfile(dataobject.fnames.rt_descriptor)
        version = 1
        descriptor_file = true
    elseif isfile(dataobject.fnames.rt_descriptor_v0)
        descriptor_file = true
    end

    if descriptor_file
        if version == 0
            f = open(dataobject.fnames.rt_descriptor_v0)
        elseif version == 1
            f = open(dataobject.fnames.rt_descriptor)
        end
        lines = readlines(f)

        if length(lines) != 0 # check for empty file

            # read descriptor variables
            descriptor_list[Symbol("nRTvar")]  = parse(Int, rsplit(lines[1],"=")[2] )
            descriptor_list[Symbol("nIons")]   = parse(Int, rsplit(lines[2],"=")[2] )
            descriptor_list[Symbol("nGroups")] = parse(Int, rsplit(lines[3],"=")[2] )
            descriptor_list[Symbol("iIons")]   = parse(Int, rsplit(lines[4],"=")[2] )

            descriptor_list[Symbol("X_fraction")] = parse(Float64, rsplit(lines[6],"=")[2] )
            descriptor_list[Symbol("Y_fraction")] = parse(Float64, rsplit(lines[7],"=")[2] )

            descriptor_list[Symbol("unit_np")] = parse(Float64, rsplit(lines[9],"=")[2] )
            descriptor_list[Symbol("unit_pf")] = parse(Float64, rsplit(lines[10],"=")[2] )
            descriptor_list[Symbol("rt_c_frac")] = parse(Float64, rsplit(lines[11],"=")[2] )

            descriptor_list[Symbol("n_star")] = parse(Float64, rsplit(lines[13],"=")[2] )
            descriptor_list[Symbol("T2_star")] = parse(Float64, rsplit(lines[14],"=")[2] )
            descriptor_list[Symbol("g_star")] = parse(Float64, rsplit(lines[15],"=")[2] )

            #todo read photon groups
            #rtPhotonGroups

        else # if file is empty
            descriptor_file = false

        end

    end

    #read header from first cpu file
    rt_files = false
    nvarrt = 0
    if  isfile(dataobject.fnames.rt * "out00001")
        rt_files = true
        f_rt = FortranFile(dataobject.fnames.rt * "out00001")
        skiplines(f_rt, 1)
        nvarrt = read(f_rt, Int32)
        skiplines(f_rt, 3)
        #println(read(f_rt, Float64)) # gamma

        close(f_rt)
    end

    dataobject.rt                   = rt_files
    dataobject.nvarrt               = nvarrt
    dataobject.rt_variable_list     = Symbol[]

    # descriptor
    dataobject.descriptor.rtversion     = version
    dataobject.descriptor.rt            = descriptor_list
    dataobject.descriptor.rtPhotonGroups  = rtPhotonGroups
    dataobject.descriptor.usert         = false
    dataobject.descriptor.rtfile        = descriptor_file

    return dataobject
end

function isgravityfile1!(dataobject::InfoType)
    grav_files = false
    # grav file of first cpu
    if  isfile(dataobject.fnames.gravity * "out00001")
        grav_files = true
    end

    dataobject.gravity                  = grav_files
    dataobject.gravity_variable_list    = [:epot,:ax,:ay,:az]

    # descriptor
    dataobject.descriptor.gravity       = [:epot,:ax,:ay,:az]
    dataobject.descriptor.usegravity    = false
    dataobject.descriptor.gravityfile   = false

    return dataobject
end

# todo: introduce new RAMSES version
function readparticlesfile1!(dataobject::InfoType)
    Npart = 0
    Ndm = 0
    Nstars = 0
    Nsinks = 0
    other_tracer1   = 0
    debris_tracer   = 0
    cloud_tracer    = 0
    star_tracer     = 0
    other_tracer2   = 0
    gas_tracer      = 0

    Ncloud  = 0
    Ndebris = 0
    Nother  = 0
    Nundefined = 0

    part_files = false
    part_header = false
    version=0
    if  isfile(dataobject.fnames.particles * "out00001")
        part_files = true
        if isfile(dataobject.fnames.header)
            part_header = true
            f = open(dataobject.fnames.header)
            line = readline(f)

            # check for header version
            if occursin("Total", String(line))   # =< stable_18_09
                version = 0
            elseif occursin("Family", String(line)) # > stable_18_09
                version = 1
            else
                version =-1
            end

            if version == 0
                Npart = parse(Int, readline(f) )
                line = readline(f)
                Ndm = parse(Int, readline(f) )
                line = readline(f)
                Nstars = parse(Int, readline(f) )
                line = readline(f)
                Nsinks = parse(Int, readline(f) )

            elseif version == 1
                other_tracer1 = parse(Int, rsplit(readline(f))[2] )
                debris_tracer = parse(Int, rsplit(readline(f))[2]  )
                cloud_tracer = parse(Int, rsplit(readline(f))[2]  )
                star_tracer = parse(Int, rsplit(readline(f))[2] )
                other_tracer2 = parse(Int, rsplit(readline(f))[2]  )
                gas_tracer = parse(Int, rsplit(readline(f))[2]  )

                Ndm     = parse(Int, rsplit(readline(f))[2]  )
                Nstars  = parse(Int, rsplit(readline(f))[2]  )
                Ncloud  = parse(Int, rsplit(readline(f))[2]  )
                Ndebris = parse(Int, rsplit(readline(f))[2]  )
                Nother  = parse(Int, rsplit(readline(f))[2]  )
                Nundefined = parse(Int, rsplit(readline(f))[2]  )
            end
            close(f)

            #if Npart != 0
            #    part_files = true
            #end
        end
    end

    dataobject.part_info = PartInfoType()
    dataobject.part_info.eta_sn = 0.
    dataobject.part_info.age_sn = 10. / dataobject.scale.Myr
    dataobject.part_info.f_w = 0.

    # overwrite some default parameters from namelist file
    if dataobject.namelist
        keylist_header = keys(dataobject.namelist_content)
        for i in keylist_header
            icontent = dataobject.namelist_content[i]
            keylist_parameters = keys(icontent)
            for j in keylist_parameters
                if j == "eta_sn"
                    dataobject.part_info.eta_sn = parse(Float64,icontent[j])
                elseif j == "age_sn"
                    dataobject.part_info.age_sn = parse(Float64,icontent[j])
                elseif j == "f_w"
                    dataobject.part_info.f_w = parse(Float64,icontent[j])
                end
            end

        end
    end

    # read descriptor file
    descriptor_file = false
    variable_descriptor_list = Symbol[]
    variable_types = String[]
    dnvar = 0
    version = 0
    if  isfile(dataobject.fnames.part_descriptor)
        descriptor_file = true
        f = open(dataobject.fnames.part_descriptor)
        lines = readlines(f)

        # read descriptor version # > stable_18_09
        version = parse(Int, rsplit(lines[1], ":" )[2])

        # read descriptor variables
        if version == 1
            dnvar = length(lines)
            for i = 3:dnvar
                ivar = String(rsplit(lines[i], "," )[2])
                itype = String(rsplit(lines[i], "," )[3])

                ivar = strip(ivar)
                itype = strip(itype)
                append!(variable_descriptor_list, [Symbol(ivar)])
                append!(variable_types, [itype])
            end
        else
        # version not supported,
        # descriptor variables not read
        end
        close(f)
    end

    dataobject.nvarp = 5
    if version <= 0
        dataobject.particles_variable_list=[:vx, :vy, :vz, :mass, :birth]
    else

        if in(:metallicity, variable_descriptor_list)
            pre_variable_list = [:position_x, :position_y, :position_z, :velocity_x, :velocity_y, :velocity_z,
                                 :mass, :identity, :levelp, :birth_time, :metallicity, :family, :tag]
            addvar_index = Int[]
            for (i,ival) in enumerate(variable_descriptor_list)
                if !in(ival, pre_variable_list)
                    append!(addvar_index, [i])
                end
            end

            particles_variable_list=[:vx, :vy, :vz, :mass, :family, :tag, :birth, :metals]
            if length(addvar_index) != 0
                addvar = variable_descriptor_list[addvar_index]
                append!(particles_variable_list, addvar)
            end
            dataobject.particles_variable_list=particles_variable_list
            dataobject.nvarp = length(particles_variable_list)
        else
            dataobject.particles_variable_list=[:vx, :vy, :vz, :mass, :family, :tag, :birth]
            dataobject.nvarp = 7
        end

        # todo: automatic detection of more variables
        #if (dnvar-3) > dataobject.nvarp + 7  # variables + (id,x,y,z,level,family,tag)
        #    for invar = (dataobject.nvarp+7):(dnvar-3)
        #        icount = invar - 3
        #        append!(dataobject.particles_variable_list, [Symbol("var$icount")])
        #        dataobject.nvarp = dataobject.nvarp + 1
        #    end
        #end
    end

    dataobject.part_info.Npart      = Npart
    dataobject.part_info.Ndm        = Ndm
    dataobject.part_info.Nstars     = Nstars
    dataobject.part_info.Nsinks     = Nsinks
    dataobject.part_info.Ncloud     = Ncloud
    dataobject.part_info.Ndebris    = Ndebris
    dataobject.part_info.Nother     = Nother
    dataobject.part_info.Nundefined = Nundefined

    dataobject.part_info.other_tracer1  = other_tracer1
    dataobject.part_info.debris_tracer  = debris_tracer
    dataobject.part_info.cloud_tracer   = cloud_tracer
    dataobject.part_info.star_tracer    = star_tracer
    dataobject.part_info.other_tracer2  = other_tracer2
    dataobject.part_info.gas_tracer     = gas_tracer

    dataobject.particles  = part_files
    dataobject.headerfile = part_header

    dataobject.descriptor.pversion = version
    if version == 0
        dataobject.descriptor.particles = dataobject.particles_variable_list
        dataobject.descriptor.ptypes = String[]
    elseif version > 0
        dataobject.descriptor.particles = variable_descriptor_list
        dataobject.descriptor.ptypes = variable_types
    end

    dataobject.descriptor.useparticles=false
    dataobject.descriptor.particlesfile=descriptor_file
    return dataobject
end

function readnamelistfile!(dataobject::InfoType)
    namelist_file = false
    asciifile = false
    namelist = Dict()
    variables = Dict()
    if  isfile(dataobject.fnames.namelist)

        f = open(dataobject.fnames.namelist )
        lines = readlines(f)
        #check for ASCII content

        for i in lines
            if occursin("&RUN_PARAMS", i)
                asciifile = true
                namelist_file = true
            end
        end

        if asciifile
            iheader = "false"
            Nlines = length(lines)
            for (j,i) in enumerate(lines)
                if iheader != "false"
                    if occursin("=", i)
                        variable = String(rsplit(i, "=" )[1])
                        content = String(rsplit(i, "=" )[2])
                        variables[variable] = content
                    end
                end

                if occursin("&", i) || j == Nlines
                    if iheader != "false"
                        namelist[iheader]= variables
                        variables = Dict()
                    end
                    iheader = i
                end

            end
        end
    end

    dataobject.namelist = namelist_file
    dataobject.namelist_content = namelist
    return dataobject
end

function readclumpfile1!(dataobject::InfoType)
    # clump file of first cpu
    clump_files = false
    header= Symbol[]
    if  isfile(dataobject.fnames.clumps * "txt00001")
        clump_files = true
        header, NColumns = getclumpvariables(dataobject, [:all], dataobject.fnames)
    end

    dataobject.clumps     = clump_files
    dataobject.clumps_variable_list = header

    # descriptor
    dataobject.descriptor.clumps = header
    dataobject.descriptor.useclumps = false
    dataobject.descriptor.clumpsfile = false

    return dataobject
end

function readtimerfile!(dataobject::InfoType)
    timer_file = false
    dataobject.files_content.timerfile = [""]
    if isfile(dataobject.fnames.timer )
        timer_file = true
        f = open(dataobject.fnames.timer);
        dataobject.files_content.timerfile = readlines(f)
    end
    dataobject.timerfile = timer_file

    return dataobject
end

function readcompilationfile!(dataobject::InfoType)
    compilation_file = false
    compile_date = ""
    patch_dir = ""
    remote_repo = ""
    local_branch= ""
    last_commit= ""
    if  isfile(dataobject.fnames.compilation )
        compilation_file = true
        f = open(dataobject.fnames.compilation )

        lines = readlines(f)
        if !occursin("\0", String(rsplit(lines[1], "=" )[2])) # skip embedded NULs
            compile_date = String(rsplit(lines[1], "=" )[2])
            patch_dir  =  String(rsplit(lines[2], "=" )[2])
            remote_repo =  String(rsplit(lines[3], "=" )[2])
            local_branch =  String(rsplit(lines[4], "=" )[2])
            last_commit =  String(rsplit(lines[5], "=" )[2])
        end
        close(f)
    end

    dataobject.compilation = CompilationInfoType() #"", "" ,"" , "", "" )
    dataobject.compilation.compile_date = compile_date
    dataobject.compilation.patch_dir = patch_dir
    dataobject.compilation.remote_repo = remote_repo
    dataobject.compilation.local_branch= local_branch
    dataobject.compilation.last_commit= last_commit

    dataobject.compilationfile = compilation_file
    return dataobject
end

function readmakefile!(dataobject::InfoType)
    make_file = false
    dataobject.files_content.makefile = [""]
    if  isfile(dataobject.fnames.makefile)
        make_file = true
        f = open(dataobject.fnames.makefile);
        dataobject.files_content.makefile = readlines(f)
    end

    dataobject.makefile = make_file
    return dataobject
end

function readpatchfile!(dataobject::InfoType)
    patch_file = false
    dataobject.files_content.patchfile = [""]
    if  isfile(dataobject.fnames.patchfile)
        patch_file = true
        f = open(dataobject.fnames.patchfile);
        dataobject.files_content.patchfile = readlines(f)
    end
    dataobject.patchfile = patch_file
    return dataobject
end

function printsimoverview(info::InfoType, verbose::Bool)
    if verbose
        println("Code: ", info.simcode)
        println("output [$(info.output)] summary:")
        println("mtime: ", info.mtime)
        println("ctime: ", info.ctime)
        printstyled("=======================================================\n", bold=true, color=:normal)

        time_val, time_unit  = humanize(info.time, info.scale, 2, "time")
        println("simulation time: ",   time_val, " [$time_unit]")

        boxlen_val, boxlen_unit  = humanize(info.boxlen, info.scale, 2, "length")
        println("boxlen: ",   boxlen_val, " [$boxlen_unit]")

        println("ncpu: ",     info.ncpu)
        println("ndim: ",     info.ndim)

        println("-------------------------------------------------------")
        min_cellsize, min_unit  = humanize(info.boxlen / 2^info.levelmin, info.scale, 2, "length")
        max_cellsize, max_unit  = humanize(info.boxlen / 2^info.levelmax, info.scale, 2, "length")
        println("amr:           ",     info.amr)
        if info.levelmin != info.levelmax # if AMR
            println("level(s): ", info.levelmin, " - ", info.levelmax, " --> cellsize(s): ", min_cellsize ," [$min_unit]", " - ", max_cellsize," [$max_unit]")
        else
            println("level of uniform grid: ", info.levelmax, " --> cellsize(s): ", max_cellsize," [$max_unit]")
        end

        if info.hydro
            println("-------------------------------------------------------")
        else
            println()
        end
        println("hydro:         ",     info.hydro)

        if info.hydro
            println("hydro-variables:  ",     info.nvarh, "  --> ", tuple(info.variable_list...) )
            if info.descriptor.hydrofile
                println("hydro-descriptor: ", tuple(info.descriptor.hydro...) )
            end
            println("γ: ", info.gamma)

        end

        if info.gravity
            println("-------------------------------------------------------")
        end

        println("gravity:       ",     info.gravity)
        if info.gravity
            println("gravity-variables: ", tuple(info.gravity_variable_list...))
        end

        if info.particles
            println("-------------------------------------------------------")
        end
        print("particles:     ",     info.particles)
        if info.particles
            if info.headerfile == false
                print("  (no particle header file) \n")
            else
                print("\n")

                if info.part_info.other_tracer1  != 0 @printf "- other_tracer1:    %e \n" info.part_info.other_tracer1 end
                if info.part_info.debris_tracer  != 0 @printf "- debris_tracer:    %e \n" info.part_info.debris_tracer end
                if info.part_info.cloud_tracer   != 0 @printf "- cloud_tracer:    %e \n" info.part_info.cloud_tracer end
                if info.part_info.star_tracer    != 0 @printf "- star_tracer:    %e \n" info.part_info.star_tracer end
                if info.part_info.other_tracer2  != 0 @printf "- other_tracer2:    %e \n" info.part_info.other_tracer2 end
                if info.part_info.gas_tracer     != 0 @printf "- gas_tracer:    %e \n" info.part_info.gas_tracer end

                if info.part_info.Npart != 0 @printf "- Npart:    %e \n" info.part_info.Npart end
                if info.part_info.Nstars != 0 @printf "- Nstars:   %e \n" info.part_info.Nstars end
                if info.part_info.Ndm != 0 @printf "- Ndm:      %e \n" info.part_info.Ndm end
                if info.part_info.Nsinks != 0 @printf "- Nsinks:   %e \n" info.part_info.Nsinks end
                if info.part_info.Ncloud != 0 @printf "- Ncloud:   %e \n" info.part_info.Ncloud end
                if info.part_info.Ndebris != 0 @printf "- Ndebris:   %e \n" info.part_info.Ndebris end
                if info.part_info.Nother != 0 @printf "- Nother:   %e \n" info.part_info.Nother end
                if info.part_info.Nundefined != 0 @printf "- Nundefined:   %e \n" info.part_info.Nundefined end
            end
        else
            print("\n")
        end

        if info.particles

            println("particle-variables: ", info.nvarp, "  --> ", tuple(info.particles_variable_list...) )
            if info.descriptor.particlesfile
                println("particle-descriptor: ", tuple(info.descriptor.particles...) )
            end
            if !info.rt
                println("-------------------------------------------------------")
            end
        end

        if info.rt
            println("-------------------------------------------------------")
        end
        println("rt:            ",     info.rt)
        if info.rt
            println("rt-variables: ", info.nvarrt)
            if info.descriptor.rtfile
                #println("nRTvar: ", info.descriptor.rt[:nRTvar] )
                println("nIons: ", info.descriptor.rt[:nIons] )
                println("nGroups: ", info.descriptor.rt[:nGroups] )
                println("iIons: ", info.descriptor.rt[:iIons] )
            end
            if !info.clumps
                println("-------------------------------------------------------")
            end
        end

        if info.clumps
            println("-------------------------------------------------------")
        end
        println("clumps:           ", info.clumps)
        if info.clumps
            println("clump-variables: ", tuple(info.clumps_variable_list...) )
            println("-------------------------------------------------------")
        end

        if info.namelist
            if !info.clumps
                println("-------------------------------------------------------")
            end
            println("namelist-file: ", tuple(keys(info.namelist_content)...) )
            println("-------------------------------------------------------")
        else
            println("namelist-file:    ", info.namelist)
        end

        println("timer-file:       ", info.timerfile)
        println("compilation-file: ", info.compilationfile)
        println("makefile:         ", info.makefile)
        println("patchfile:        ", info.patchfile)
        #println("nx, ny, nz: ", overview.nx, ", ", overview.ny,  ", ",overview.nz)
        #println("infofile creation-time: ", DateTime(ctime) ) #todo: format
        #println("infofile modification-time: ", DateTime(mtime) )  #todo: format
        printstyled("=======================================================\n", bold=true, color=:normal)
        println()
    end
end


text:
# this constants come from and are inspired by the pyprocessing project
# (https://code.google.com/p/pyprocessing/source/browse/trunk/pyprocessing/constants.py)
# with some slight adjustments

#modes
const CORNER = 1
export CORNER
const CORNERS = 2
export CORNERS
const CENTER = 3
export CENTER
const RADIUS = 4
export RADIUS

# arc modes
const OPEN = 1
export OPEN
const PIE = 2
export PIE
const CHORD = 3
export CHORD

#text alignment
const TOP = 7
export TOP
const BOTTOM = 8
export BOTTOM
const BASELINE = 9
export BASELINE

# const textAlignConst = [LEFT=>"left", RIGHT=>"right",
const textAlignConst = [CENTER=>"center", TOP=>"top",
                  		BOTTOM=>"bottom", BASELINE=>"baseline"]
export textAlignConst

# colorMode
#const RGB=0
#export RGB
#const HSB=1
#export HSB

# math
const PI = pi
export PI
const TWO_PI = pi*2
export TWO_PI
const HALF_PI = pi/2
export HALF_PI
const QUARTER_PI = pi/4
export QUARTER_PI
const TAU = Tau.tau
export TAU

# hints
#const ENABLE_DEPTH_TEST = "ENABLE_DEPTH_TEST"
#export ENABLE_DEPTH_TEST
#const DISABLE_DEPTH_TEST = "DISABLE_DEPTH_TEST"
#export DISABLE_DEPTH_TEST
#const DISABLE_POLYGON_SMOOTH = "DISABLE_POLYGON_SMOOTH"
#export DISABLE_POLYGON_SMOOTH
#const DOUBLE_FLIP_POLICY = "DOUBLE_FLIP_POLICY"
#export DOUBLE_FLIP_POLICY
#const SINGLE_FLIP_POLICY = "SINGLE_FLIP_POLICY"
#export SINGLE_FLIP_POLICY
#const FBO_FLIP_POLICY = "FBO_FLIP_POLICY"
#export FBO_FLIP_POLICY
#const ACCUM_FLIP_POLICY = "ACCUM_FLIP_POLICY"
#export ACCUM_FLIP_POLICY
#const BACKUP_FLIP_POLICY = "BACKUP_FLIP_POLICY"
#export BACKUP_FLIP_POLICY

#filter modes
const THRESHOLD = "THRESOLD"
export THRESHOLD
const GRAY = "GRAY"
export GRAY
const INVERT = "INVERT"
export INVERT
const POSTERIZE = "POSTERIZE"
export POSTERIZE
const BLUR = "BLUR"
export BLUR
const OPAQUE = "OPAQUE"
export OPAQUE
const ERODE = "ERODE"
export ERODE
const DILATE = "DILATE"
export DILATE

# shapes
#const POINTS = OpenGLStd.POINTS
#export POINTS
#const LINES = OpenGLStd.LINES
#export LINES
#const TRIANGLES = OpenGLStd.TRIANGLES
#export TRIANGLES
#const TRIANGLE_FAN = OpenGLStd.TRIANGLE_FAN
#export TRIANGLE_FAN
#const TRIANGLE_STRIP = OpenGLStd.TRIANGLE_STRIP
#export TRIANGLE_STRIP
#const QUADS = OpenGLStd.QUADS
#export QUADS
#const QUAD_STRIP = OpenGLStd.QUAD_STRIP
#export QUAD_STRIP
#const CLOSE = 1
#export CLOSE

#blend modes
const BLEND = 1
export BLEND
const ADD = 2
export ADD
const SUBTRACT = 3
export SUBTRACT
const DARKEST = 4
export DARKEST
const LIGHTEST = 5
export LIGHTEST
const DIFFERENCE = 6
export DIFFERENCE
const EXCLUSION = 7
export EXCLUSION
const MULTIPLY = 8
export MULTIPLY
const SCREEN = 9
export SCREEN
const OVERLAY = 10
export OVERLAY
const HARD_LIGHT = 11
export HARD_LIGHT
const SOFT_LIGHT = 12
export SOFT_LIGHT
const DODGE = 13
export DODGE
const BURN = 14
export BURN

# texture modes
const NORMALIZED = "NORMALIZED"
export NORMALIZED
const IMAGE = "IMAGE"
export IMAGE

# color formats
#const RGB = "RGB"
#export RGB
#const ARGB = "RGBA"
#export ARGB
#const ALPHA = "A"
#export ALPHA

# stroke join modes
const MITER = "MITER"
export MITER
const BEVEL = "BEVEL"
export BEVEL

# stroke cap modes
const SQUARE = "SQUARE"
export SQUARE
const PROJECT = "PROJECT"
export PROJECT
const ROUND = "ROUND"
export ROUND

# cursor types
const ARROW = "arrow"
export ARROW
const CROSS = "crosshair"
export CROSS
const HAND = "hand"
export HAND
const MOVE = "size"
export MOVE
const TEXT = "text"
export TEXT
const WAIT = "wait"
export WAIT

#perlin noise table
const perlin = [
151,160,137,91,90,15,131,13,201,95,96,53,194,233,7,225,140,36,103,
30,69,142,8,99,37,240,21,10,23,190,6,148,247,120,234,75,0,26,197,
62,94,252,219,203,117,35,11,32,57,177,33,88,237,149,56,87,174,20,
125,136,171,168,68,175,74,165,71,134,139,48,27,166,77,146,158,231,
83,111,229,122,60,211,133,230,220,105,92,41,55,46,245,40,244,102,
143,54,65,25,63,161,1,216,80,73,209,76,132,187,208,89,18,169,200,
196,135,130,116,188,159,86,164,100,109,198,173,186,3,64,52,217,226,
250,124,123,5,202,38,147,118,126,255,82,85,212,207,206,59,227,47,16,
58,17,182,189,28,42,223,183,170,213,119,248,152,2,44,154,163,70,
221,153,101,155,167,43,172,9,129,22,39,253,19,98,108,110,79,113,
224,232,178,185,112,104,218,246,97,228,251,34,242,193,238,210,144,
12,191,179,162,241,81,51,145,235,249,14,239,107,49,192,214,31,181,
199,106,157,184,84,204,176,115,121,50,45,127,4,150,254,138,236,
205,93,222,114,67,29,24,72,243,141,128,195,78,66,215,61,156,180,
151,160,137,91,90,15,131,13,201,95,96,53,194,233,7,225,140,36,103,
30,69,142,8,99,37,240,21,10,23,190,6,148,247,120,234,75,0,26,197,
62,94,252,219,203,117,35,11,32,57,177,33,88,237,149,56,87,174,20,
125,136,171,168,68,175,74,165,71,134,139,48,27,166,77,146,158,231,
83,111,229,122,60,211,133,230,220,105,92,41,55,46,245,40,244,102,
143,54,65,25,63,161,1,216,80,73,209,76,132,187,208,89,18,169,200,
196,135,130,116,188,159,86,164,100,109,198,173,186,3,64,52,217,226,
250,124,123,5,202,38,147,118,126,255,82,85,212,207,206,59,227,47,16,
58,17,182,189,28,42,223,183,170,213,119,248,152,2,44,154,163,70,
221,153,101,155,167,43,172,9,129,22,39,253,19,98,108,110,79,113,
224,232,178,185,112,104,218,246,97,228,251,34,242,193,238,210,144,
12,191,179,162,241,81,51,145,235,249,14,239,107,49,192,214,31,181,
199,106,157,184,84,204,176,115,121,50,45,127,4,150,254,138,236,
205,93,222,114,67,29,24,72,243,141,128,195,78,66,215,61,156,180]
export perlin


text:
"""
    module Simulate_BouncingBall

Simulation of a Modia3D model of a ball, bouncing on ground using an
elastic response calculation. The result is compared
with the simulation of a direct implementation of this system in ModiaMath
using an impulsive response calculation.

There is a bug with the event handling (if interval=0.01 is selected, the simulation is wrong.)
"""
module Simulate_BouncingBall

using  Modia3D
import Modia3D.ModiaMath

mutable struct Model <: ModiaMath.AbstractSimulationModel
    simulationState::ModiaMath.SimulationState

    # Parameters
    h0::Float64   # initial height
    cor::Float64    # coefficient of restitution
    g::Float64    # gravity constant
    vsmall::Float64

    # Discrete variables
    flying::Bool

    function Model(;h0=1.0, cor=0.7, g=9.81, vsmall=0.1)
        @assert(h0 > 0.0)
        @assert(0.0 <= cor <= 1.0)
        @assert(vsmall > 0.0)
        simulationState = ModiaMath.SimulationState("BouncingBall1", getModelResidues!, [h0;0.0], getVariableName;
                                                    nz=1, nw=2)
        new(simulationState, h0, cor, g, vsmall, true)
    end
end

getVariableName(model, vcat, vindex) = ModiaMath.getVariableName(model, vcat, vindex;
                                                                 xNames=["h", "v"],
                                                                 wNames=["cor_res", "flying"])

function getModelResidues!(m::Model, t::Float64, _x::Vector{Float64}, _derx::Vector{Float64}, _r::Vector{Float64}, _w::Vector{Float64})
    sim = m.simulationState
    if ModiaMath.isInitial(sim)
        println("... h0 = ", _x[1])
        if ModiaMath.isLogInfos(sim)
            println("        flying = ", m.flying)
        end
        _w[1] = m.cor
    end
    if ModiaMath.isTerminal(sim)
        if ModiaMath.isLogInfos(sim)
            println("\n      BouncingBall model is terminated (flying = ", m.flying, ")")
        end
        return
    end

    h = _x[1]
    v = _x[2]

    # println("+++ time = ", t, ", h = ", h)
    if ModiaMath.edge!(sim, 1, -h, "-h")
        cor_res = Modia3D.resultantCoefficientOfRestitution(m.cor, abs(v), m.vsmall, cor_min=0.0)
        _w[1] =  cor_res
            v = -cor_res * v
        _x[2] =  v    # re-initialize state vector x
        println("... v = $v")
    end
#=
    if ModiaMath.edge!(sim, 2, v, "v")
        if h <= 0.0
            m.flying = false
            v = 0.0
            if ModiaMath.isLogInfos(sim)
                println("        flying = ", m.flying)
            end
            _x[2] = v    # re-initialize state vector x
        end
    end
=#

    derh = v
    derv = m.flying ? -m.g : 0.0

    _r[1] = derh - _derx[1]
    _r[2] = derv - _derx[2]
    _w[2] = m.flying

    return nothing
end

cor      = 0.7
vsmall   = 0.1
radius   = 0.01
h0       = 20*radius
stopTime = 1.0

model1  = Model(h0=h0, cor=cor, vsmall=vsmall )
#result1 = ModiaMath.simulate!(model1, stopTime=stopTime, log=true)
result1 = ModiaMath.simulate!(model1, stopTime=stopTime, interval=0.02, log=true)  # interval=0.01 gives wrong behaviour
ModiaMath.plot(result1, ["h", ("v", "flying"), "cor_res"], heading="Bouncing ball with impulsive response calculation", figure=1)

ballMaterial  = Modia3D.Material(color="Red"       , transparency=0.5)
tableMaterial = Modia3D.Material(color="LightBlue" , transparency=0.5)

@assembly BouncingBall2 begin
  world = Object3D()
  ball  = Object3D(world, Solid(SolidSphere(2*radius) , "Steel"  , ballMaterial ; contactMaterial="Steel"); r=[0.0, 0.0, h0], fixed=false)
  box   = Object3D(world, Solid(SolidBox(20*radius, 20*radius, 2*radius), "DryWood", tableMaterial; contactMaterial="Steel");
                   r=[0.0, 0.0,-2*radius], fixed=true)
end

gravField     = Modia3D.UniformGravityField(g=9.81, n=[0,0,-1.0])
bouncingBall2 = BouncingBall2(sceneOptions=Modia3D.SceneOptions(gravityField=gravField, visualizeFrames=false,
                                                                defaultFrameLength=0.1, enableContactDetection=true))
model2        = SimulationModel( bouncingBall2 )
result2       = ModiaMath.simulate!(model2; stopTime=stopTime, tolerance=1e-8, log=false)
ModiaMath.plot(result2, ["ball.r[3]", "ball.v[3]"], heading="Bouncing ball with compliant response calculation", figure=2)

# Need to be included in ModiaMath
import Base
function Base.getproperty(result::ModiaMath.ResultWithVariables, name::AbstractString)
    (sig, dummy1, dummy2) = ModiaMath.Result.getSignal(result.series, name)
    return sig
end

#=
# The following commands can only be used if PyPlot is available
using PyPlot
using PyCall

fig, ax = PyPlot.subplots(figsize=(3,9))

pyplot_rc = PyCall.PyDict(PyPlot.matplotlib."rcParams")
pyplot_rc["font.family"]      = "sans-serif"
pyplot_rc["font.sans-serif"]  = ["Calibri", "Arial", "Verdana", "Lucida Grande"]
pyplot_rc["font.size"]        = 12.0
pyplot_rc["lines.linewidth"]  = 1.5
pyplot_rc["grid.linewidth"]   = 0.5
pyplot_rc["axes.grid"]        = true
pyplot_rc["axes.titlesize"]   = "medium"
pyplot_rc["figure.titlesize"] = "medium"

t1      = result1["time"]
h1      = result1["h"]
v1      = result1["v"]
cor_res = result1["cor_res"]

t2 = result2."time"
h2 = result2."ball.r[3]"
v2 = result2."ball.v[3]"

figure(3)
clf()
plot(t1,h1,"b", t2,h2,"r", t1,cor_res,"g")
grid(true)
xlabel("time [s]")
ylabel("height [m]")
legend(["impulsive response",
        "compliant response",
        "\$cor_{res}\$"])

plot(t1,h1,"b", t2,h2,"r")
grid(true)
xlabel("time [s]")
ylabel("height [m]")
legend(["impulsive response",
        "compliant response"])
=#

println("... success of examples/collisions/Simulate_BouncingBall.jl")
end


text:
include("splitcell.jl")
include("lineage.jl");
include("segmentation3d.jl")
include("normalization3d.jl")
include("julia2ims.jl")
using Dates

data_dir = "/datahub/rawdata/tandeng/mRNA_imaging/mRNA_confocal_hamamatsu-60X-TIRF/20200318";
ret_dir = "/datahub/rawdata/tandeng/mRNA_imaging/mRNA_confocal_hamamatsu-60X-TIRF/20200318_result/output";

#for s in 3:19
#for s in 2:20
function playing(s)
    println("~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~")
    println("$(Dates.now()) Loading $s")
	@time img = load(File(format"TIFF","$data_dir/HE7-1NLS-22uw-const_1_s$s.ome.tiff"));
	GC.gc()
	println("Load done")

    t_len = size(img)[3]÷20
	@time markers = split_cell_LoG(img, t_len);
    GC.gc()
	@time time_line, longlived_labels, livingtime, time_line_whole =
		find_time_line(markers);
	if length(longlived_labels) == 0
		println("No longlived cells found, skipping")
		return nothing
	end
	@time split_contacted_cell!(time_line, longlived_labels, livingtime, time_line_whole);
	tracks = walking(time_line, longlived_labels, livingtime);
	@time longlived_maps, watershed_maps = grant_domain(img, time_line, longlived_labels, livingtime, time_line_whole);
    #save("$ret_dir/d21s$(s)-watershed.jld", "maps", watershed_maps)
	for index in 1:length(longlived_labels)
		println("tracking cell $index")
		@time local cell = pick_cell(img, longlived_maps, longlived_labels[index], tracks[:,:,index],
								livingtime[index]);
		@time cell_nu, nucleus_size, nucleus_th = extract3dnucleus(cell);
		@time cell_nu_nor, nor_para = normalize(cell_nu);
        save("$ret_dir/d18s$(s)-$index.jld", "threshold", nucleus_th, "normal", nor_para,
             "livingtime", livingtime[index], "track", tracks[:,:,index], "size", nucleus_size)
        save2ims(reinterpret.(reshape(cell_nu_nor, (512,512,20,t_len))),
                 "$ret_dir/d18s$(s)-$index.ims")
	end
	GC.gc()
end

println("Processing $data_dir, output $ret_dir")
for i in parse(Int, ARGS[1]):parse(Int, ARGS[2])
    playing(i)
end


text:
#: ==========================================================

function preprocess_raw_input(raw::String)
    s = "$raw"
    for r ∈ __preproc__
        s = replace(s, r)
    end
    return s
end

#: ==========================================================

function mask(s0::String)::String
    s = "$s0"
    for r ∈ MASK_RULES
        s = replace(s, r)
    end
    return s
end

#: ==========================================================

splt_ln(l) = split(l,r"[^\S\n\r]",keepempty=false)

encode_line(l,dic) = Int[dic[w] for w ∈ splt_ln(l)]

function tokenize(
    lines::Vector{S},
    code::Dict{String,Int}
    )::Vector{TPattern} where {S<:AbstractString}
    enc(l) = encode_line(l,code)
    patts = TPattern[]
    try
        #: line *1
        patts = TPattern[((unique(p)==[0,]) ? __PATT__all_number_line__ : p) for p ∈ enc.(lines)]
    catch _e_
        @warn "tokenize failed."
        @warn _e_
        return TPattern[]
    end
    return patts
end

function tokenize0(S0::String)::Tuple{Vector{TPattern},Dict{String,Int}}
    lines = split(S0,"\n",keepempty=false) ; #TODO
    unique_words = unique(vcat(splt_ln.(lines)...)) ;
    code = Dict{String,TCode}(w=>i for (i,w) ∈ enumerate(unique_words))
    code["£"] = 0  # reserved token 0 for number
    conflict_kw = [k for (k,v) ∈ code if "£"!=k && occursin("£",k)]
    if length(conflict_kw)>0
        @warn "conflict_kw = $(conflict_kw)"
    end
    patts = tokenize(lines, code)
    return patts, code
end

function tokenize(
    S0::String;
    REF_CODE = __REFDIC__
    )::Tuple{Vector{TPattern},Dict{String,TCode}}

    lines = split(S0,"\n",keepempty=false) ; #TODO
    unique_words = unique(vcat(splt_ln.(lines)...)) ;
    code = copy(REF_CODE)
    k_REF_CODE = keys(REF_CODE)
    i = length(k_REF_CODE)+2
    for w ∈ unique_words
        if w∉k_REF_CODE &&  w!="£"
            code[w] = i
            i += 1
        end
    end
    code["£"] = 0  # reserved token 0 for number
    conflict_kw = [k for (k,v) ∈ code if "£"!=k && occursin("£",k)]
    if length(conflict_kw)>0
        @warn "conflict_kw = $(conflict_kw)"
    end
    patts = tokenize(lines, code)
    return patts, code
end

#: ==========================================================

function revert(code::Dict{String,Int})
    dic = Dict{Int,String}(i=>k for (k,i) ∈ code)
    dic[__PATT__all_number_line__[1]] = "£++++++++++++++++"
    return dic
end

#: ==========================================================

function add_newline(lines)
    @inline headspace(l) = ((length(l)>0 && startswith(l," "))
                                ? (findfirst(x->x!=' ',l)!==nothing ? (findfirst(x->x!=' ',l)-1)
                                : findlast(x->x==' ',l)) : 0)
    hstack = Stack{Int}()
    lines1 = []
    for l ∈ lines
        h = headspace(l)
        if length(hstack)>0
            if h > top(hstack)
                push!(lines1,"")
            elseif h == top(hstack)
                nothing
            else
                ht = top(hstack)
                push!(lines1,"")
                while length(hstack)>0 && top(hstack)>h
                    x = pop!(hstack)
                    if x != ht
                        ht = x
                        push!(lines1,"")
                    end
                end
            end
        end
        push!(lines1, l)
        push!(hstack, h)
    end
    return lines1
end

function load_file(
    fn
    )::Tuple{Vector{String},Vector{String},Vector{TPattern},Dict{String,Int}}
    if isfile(fn)
        S0 = read(fn,String) |> preprocess_raw_input
        MS0 = mask(S0)
        patts, code = tokenize(MS0)
        lines = split(S0,"\n",keepempty=false) #TODO
        lines_masked = split(MS0,"\n",keepempty=false) #TODO
        return lines, lines_masked, patts, code
    else
        return String[], String[], TPattern[], Dict{String,Int}()
    end
end

max_stars_repo_name: algorithmx/SmartParser

max_issues_repo_name: algorithmx/SmartParser

max_forks_repo_name: algorithmx/SmartParser


text:
using PtFEM

data = Dict(
  # Beam(ndim, nst, nxe, nip, direction, finite_element(nod, nodof), axisymmet
  :struc_el => Beam(3, 1, 4, 1, :x, Line(2, 3), false),
  :properties => [
    4.0e6 1.0e6 0.3e6 0.3e6],
  :x_coords => [0.0, 1.0, 2.0, 3.0, 4.0],
  :y_coords => zeros(5),
  :z_coords => zeros(5),
  :gamma => [0.0, 0.0, 90.0],
  :g_num => [
    1 2 3 4;
    2 3 4 5],
  :support => [
    (1, [0 0 0])
    ],
  :loaded_nodes => [
    (5, [0.0 -100.0 0.0])]
)

data |> display
println()

@time fem, dis_df, fm_df = p43(data)
println()

display(dis_df)
println()
display(fm_df)


text:
using FractionalDiffEq
h=0.01
T=50
α=0.97
τ=2
q=0.5
f(t, ϕ, y) = 2*ϕ/(1+ϕ^9.65)-y

prob = FDDEProblem(f, q, α, τ, T)
x, y=solve(prob, h, DelayABM())

using Plots
plot(x[500:end], y[500:end])


text:
module TestParameters

# using Revise
using MLJ
using Test

import MLJ: Scale, scale, transform, inverse_transform

mutable struct DummyModel <: Deterministic
    K::Int
    metric::Float64
    kernel::Char
end

dummy_model = DummyModel(4, 9.5, 'k')
@test params(dummy_model) ==
    (K = 4, metric = 9.5, kernel = 'k')

mutable struct SuperModel <: Deterministic
    lambda::Float64
    model1::DummyModel
    model2::DummyModel
end

dummy1 = DummyModel(1, 9.5, 'k')
dummy2 = DummyModel(2, 9.5, 'k')
super_model = SuperModel(0.5, dummy1, dummy2)
params(super_model)

tree = params(super_model)

p1 = range(dummy_model, :K, lower=1, upper=10, scale=:log10)
p2 = range(dummy_model, :kernel, values=['c', 'd'])
p3 = range(super_model, :lambda, lower=0.1, upper=1, scale=:log2)
p4 = range(dummy_model, :K, lower=1, upper=3, scale=x->2x)
@test_throws ErrorException range(dummy_model, :K, lower=1, values=['c', 'd'])
@test_throws ErrorException range(dummy_model, :kernel, upper=10)

@test MLJ.scale(p1) == :log10
@test MLJ.scale(p2) == :none
@test MLJ.scale(p3) == :log2
@test MLJ.scale(p4) == :custom
@test scale(sin) === sin
@test transform(Scale, scale(:log), ℯ) == 1
@test inverse_transform(Scale, scale(:log), 1) == float(ℯ)

@test MLJ.iterator(p1, 5)  == [1, 2, 3, 6, 10]
@test MLJ.iterator(p2) == collect(p2.values)
u = 2^(log2(0.1)/2)
@test MLJ.iterator(p3, 3) ≈ [0.1, u, 1]
@test MLJ.iterator(p4, 3) == [2, 4, 6]

# test ranges constructed from neste parameters specified with dots:
q1 = range(super_model, :(model1.K) , lower=1, upper=10, scale=:log10)
@test iterator(q1, 5) == iterator(p1, 5)
q2 = range

# test unwinding of iterators
iterators = ([1, 2], ["a","b"], ["x", "y", "z"])
@test MLJ.unwind(iterators...) ==
[1  "a"  "x";
 2  "a"  "x";
 1  "b"  "x";
 2  "b"  "x";
 1  "a"  "y";
 2  "a"  "y";
 1  "b"  "y";
 2  "b"  "y";
 1  "a"  "z";
 2  "a"  "z";
 1  "b"  "z";
 2  "b"  "z"]

end
true


text:
using LinearAlgebra
using Test
using Random
using MyGEMM.TransposeNoPacking: multiply_transpose_no_packing!
using MyGEMM.TransposeMultithreading: multiply_transpose!
using MyGEMM.Types: TilingParams

Random.seed!(777)

params = TilingParams{96, 2048, 256, 8, 6, 4}()

@testset "TransposeNoPacking" begin
  # Square test
  n, m, k = 48, 48, 48
  A, B, C = rand(m, k), rand(n, k), rand(m, n)
  C0 = C + A * B'
  multiply_transpose_no_packing!(C, A, B, params)
  @test C0 ≈ C

  # Rectangle test
  n, m, k = 48, 2*48, 3*48
  A, B, C = rand(m, k), rand(n, k), rand(m, n)
  C0 = C + A * B'
  multiply_transpose_no_packing!(C, A, B, params)
  @test C0 ≈ C
end

@testset "TransposeMultithreading" begin
  # Square test
  n, m, k = 48, 48, 48
  A, B, C = rand(m, k), rand(n, k), rand(m, n)
  C0 = C + A * B'
  multiply_transpose!(C, A, B, params)
  @test C0 ≈ C

  # Rectangle test
  n, m, k = 48, 2*48, 3*48
  A, B, C = rand(m, k), rand(n, k), rand(m, n)
  C0 = C + A * B'
  multiply_transpose!(C, A, B, params)
  @test C0 ≈ C
end

max_stars_repo_name: bischtob/MyGEMM

max_issues_repo_name: bischtob/MyGEMM

max_forks_repo_name: bischtob/MyGEMM