text:
"""
    abstract type Grid{Dc,Dp}

Abstract type that represents mesh of a domain of parametric dimension `Dc` and
physical dimension `Dp`.

The interface of `Grid` is defined by overloading:

- [`get_node_coordinates(trian::Grid)`](@ref)
- [`get_cell_node_ids(trian::Grid)`](@ref)
- [`get_reffes(trian::Grid)`](@ref)
- [`get_cell_type(trian::Grid)`](@ref)

The `Grid`  interface has the following traits

- [`OrientationStyle(::Type{<:Grid})`](@ref)
- [`RegularityStyle(::Type{<:Grid})`](@ref)

The interface of `Grid` is tested with

- [`test_grid`](@ref)

"""
abstract type Grid{Dc,Dp} <: GridapType end

# Traits

abstract type OrientationStyle end
struct Oriented <: OrientationStyle end
struct NonOriented <: OrientationStyle end

"""
    OrientationStyle(::Type{<:Grid})
    OrientationStyle(::Grid)

`Oriented()` if has oriented faces, `NonOriented()` otherwise (default).
"""
OrientationStyle(a::Grid) = OrientationStyle(typeof(a))
OrientationStyle(::Type{<:Grid}) = NonOriented()

abstract type RegularityStyle end
struct Regular <: RegularityStyle end
struct Irregular <: RegularityStyle end

"""
    RegularityStyle(::Type{<:Grid})
    RegularityStyle(::Grid)

`Regular()` if no hanging-nodes default), `Irregular()` otherwise.
"""
RegularityStyle(a::Grid) = RegularityStyle(typeof(a))
RegularityStyle(::Type{<:Grid}) = Regular()

# Interface

"""
    get_node_coordinates(trian::Grid) -> AbstractArray{<:Point{Dp}}
"""
function get_node_coordinates(trian::Grid)
  @abstractmethod
end

"""
    get_cell_node_ids(trian::Grid)
"""
function get_cell_node_ids(trian::Grid)
  @abstractmethod
end

"""
    get_reffes(trian::Grid) -> Vector{LagrangianRefFE}
"""
function get_reffes(trian::Grid)
  @abstractmethod
end

"""
    get_cell_type(trian::Grid) -> AbstractVector{<:Integer}
"""
function get_cell_type(trian::Grid)
  @abstractmethod
end

"""
    get_facet_normal(trian::Grid)
"""
function get_facet_normal(trian::Grid)
  Dp = num_point_dims(trian)
  Dc = num_cell_dims(trian)
  if Dp == Dc + 1
    @abstractmethod
  else
    @unreachable "get_facet_normal does not make sense for this grid"
  end
end

"""
    test_grid(trian::Grid)
"""
function test_grid(trian::Grid{Dc,Dp}) where {Dc,Dp}
  @test num_cell_dims(trian) == Dc
  @test num_point_dims(trian) == Dp
  @test num_cell_dims(typeof(trian)) == Dc
  @test num_point_dims(typeof(trian)) == Dp
  cell_coords = get_cell_coordinates(trian)
  @test isa(cell_coords,AbstractArray{<:AbstractVector{<:Point}})
  reffes = get_reffes(trian)
  @test isa(reffes,AbstractVector{<:LagrangianRefFE{Dc}})
  cell_types = get_cell_type(trian)
  @test isa(cell_types,AbstractArray{<:Integer})
  ncells = num_cells(trian)
  @test ncells == length(cell_coords)
  @test ncells == length(cell_types)
  nodes_coords = get_node_coordinates(trian)
  @test isa(nodes_coords,AbstractArray{<:Point})
  cell_node_ids = get_cell_node_ids(trian)
  @test isa(cell_node_ids,AbstractArray{<:AbstractArray{<:Integer}})
  @test num_nodes(trian) == length(nodes_coords)
  @test isa(is_oriented(trian),Bool)
  @test isa(is_regular(trian),Bool)
  @test OrientationStyle(trian) in (Oriented(), NonOriented())
  @test RegularityStyle(trian) in (Regular(), Irregular())
  @test is_oriented(trian) == (OrientationStyle(trian) == Oriented())
  @test is_regular(trian) == (RegularityStyle(trian) == Regular())
end

# Some API

"""
    num_cells(trian::Grid) -> Int
"""
num_cells(trian::Grid) = length(get_cell_type(trian))

"""
    num_nodes(trian::Grid) -> Int
"""
num_nodes(trian::Grid) = length(get_node_coordinates(trian))

"""
    num_cell_dims(::Grid) -> Int
    num_cell_dims(::Type{<:Grid}) -> Int
"""
num_cell_dims(::Grid{Dc,Dp}) where {Dc,Dp} = Dc
num_cell_dims(::Type{<:Grid{Dc,Dp}}) where {Dc,Dp} = Dc

"""
    num_point_dims(::Grid) -> Int
    num_point_dims(::Type{<:Grid}) -> Int
"""
num_point_dims(::Grid{Dc,Dp}) where {Dc,Dp} = Dp
num_point_dims(::Type{<:Grid{Dc,Dp}}) where {Dc,Dp} = Dp

"""
    num_dims(::Grid) -> Int
    num_dims(::Type{<:Grid}) -> Int

Equivalent to `num_cell_dims`.
"""
num_dims(g::Grid{Dc}) where Dc = Dc
num_dims(::Type{<:Grid{Dc}}) where Dc = Dc

"""
    is_first_order(trian::Grid) -> Bool
"""
function is_first_order(trian::Grid)
  reffes = get_reffes(trian)
  all(map(is_first_order,reffes))
end

"""
    get_cell_reffe(trian::Grid) -> Vector{<:LagrangianRefFE}

It is not desirable to iterate over the resulting array
for large number of cells if the underlying reference FEs
are of different Julia type.
"""
function get_cell_reffe(trian::Grid)
  type_to_reffe = get_reffes(trian)
  cell_to_type = get_cell_type(trian)
  expand_cell_data(type_to_reffe,cell_to_type)
end

"""
"""
function get_cell_ref_coordinates(trian::Grid)
  type_to_reffe = get_reffes(trian)
  type_to_coords = map(get_node_coordinates,type_to_reffe)
  cell_to_type = get_cell_type(trian)
  expand_cell_data(type_to_coords,cell_to_type)
end

"""
    get_cell_shapefuns(trian::Grid) -> Vector{<:Field}
"""
function get_cell_shapefuns(trian::Grid)
  type_to_reffes = get_reffes(trian)
  cell_to_type = get_cell_type(trian)
  type_to_shapefuns = map(get_shapefuns, type_to_reffes)
  expand_cell_data(type_to_shapefuns,cell_to_type)
end

"""
    get_cell_map(trian::Grid) -> Vector{<:Field}
"""
function get_cell_map(trian::Grid)
  cell_to_coords = get_cell_coordinates(trian)
  cell_to_shapefuns = get_cell_shapefuns(trian)
  lazy_map(linear_combination,cell_to_coords,cell_to_shapefuns)
end

function get_cell_coordinates(trian::Grid)
  node_to_coords = get_node_coordinates(trian)
  cell_to_nodes = get_cell_node_ids(trian)
  lazy_map(Broadcasting(Reindex(node_to_coords)),cell_to_nodes)
end

function Quadrature(trian::Grid,args...;kwargs...)
  cell_ctype = get_cell_type(trian)
  ctype_polytope = map(get_polytope,get_reffes(trian))
  ctype_quad = map(p->Quadrature(p,args...;kwargs...),ctype_polytope)
  cell_quad = expand_cell_data(ctype_quad,cell_ctype)
end

"""
    is_oriented(::Type{<:Grid}) -> Bool
    is_oriented(a::Grid) -> Bool
"""
is_oriented(a::Grid) = is_oriented(typeof(a))
is_oriented(a::Type{T}) where T<:Grid = OrientationStyle(T) == Oriented()

"""
    is_regular(::Type{<:Grid}) -> Bool
    is_regular(a::Grid) -> Bool
"""
is_regular(a::Grid) = is_regular(typeof(a))
is_regular(a::Type{T}) where T<:Grid = RegularityStyle(T) == Regular()

"""
    Grid(reffe::LagrangianRefFE)
"""
function Grid(reffe::LagrangianRefFE)
  UnstructuredGrid(reffe)
end

"""
    compute_linear_grid(reffe::LagrangianRefFE)
"""
function compute_linear_grid(reffe::LagrangianRefFE)
  @notimplementedif ! is_n_cube(get_polytope(reffe)) "linear grid only implemented for n-cubes at the moment"
  if get_order(reffe) == 0
    D = num_cell_dims(reffe)
    partition = tfill(1,Val{D}())
  else
    partition = get_orders(reffe)
  end
  pmin, pmax = get_bounding_box(get_polytope(reffe))
  desc = CartesianDescriptor(pmin,pmax,partition)
  CartesianGrid(desc)
end

"""
    compute_reference_grid(p::LagrangianRefFE, nelems::Integer)
"""
function compute_reference_grid(reffe::LagrangianRefFE, nelems::Integer)
  p = get_polytope(reffe)
  r = LagrangianRefFE(Float64,p,nelems)
  compute_linear_grid(r)
end

"""
    Grid(::Type{ReferenceFE{d}},p::Polytope) where d
"""
function Grid(::Type{ReferenceFE{d}},p::Polytope) where d
  UnstructuredGrid(ReferenceFE{d},p)
end

"""
    simplexify(grid::Grid)
"""
function simplexify(grid::Grid)
  simplexify(UnstructuredGrid(grid))
end


text:
"""
	PipeCell{T,S,U}(x::T, y::S, f::U, layer::Int, tinfo::String, oinfo::String)

Constructs the basic object for complex processing data, the `PipeCell`. The constructor is parametric so it can be
easily extended to new types of processing.

# Main `PipeCell` constructors
  * `PipeCell(x::T where T<:PTuple{Abstractcell}, tinfo::String = "")` creates a `stacked pipe`
  * `PipeCell(x::T where T<:PTuple{AbstractCell}, dispatch::SortedDict, tinfo::String = "")` creates a `parallel pipe`
  * `PipeCell(x::T where T<:PTuple{AbstractCell}, order::Vector{Int}, tinfo::String="")` creates a `serial pipe`

where `const PTuple{T} = Tuple{Vararg{<:T}}`.

The concepts are as follows: when piping data to a stacked pipe, the data is individually piped to each AbstractCell
of the pipe. For parallel pipes, dispatch allows to select which element of the collection (e.g. `Vector`, `Tuple`, data cell etc.)
gets piped where. Serial pipes process data in a serial manner. Pipes can be created also by using the functions
`pipestack`, `pipeparallel` and `pipeserial`.
"""
#############
# PipeCells #
#############
struct PipeCell{T,S,U} <: AbstractCell{T,S,U}
	x::T                                                                                # List of Data/Function/PipeCells
    	y::S                                                                                # Processing information
    	f::U                                                                                # <unused so far>
	layer::Int										# Layer
    	tinfo::String										# Title or name
    	oinfo::String                                                                        	# Other information
end

# General constructors for Stacked pipes (e.g. input routed to all cells)
PipeCell(x::T where T<:PTuple{AbstractCell}, tinfo::String = "") = PipeCell(x, nothing, nothing, countlayers(x)+1, tinfo, oinfoglobal )

# General constructor for Parallel pipes (e.g. input routed according to some dispatch information to cells)
PipeCell(x::T where T<:PTuple{AbstractCell}, dispatch::SortedDict, tinfo::String = "") = begin
	@assert length(x) == length(keys(dispatch)) "Parallel pipe element number does not match dispatch information."
	PipeCell(x, dispatch, nothing, countlayers(x)+1, tinfo, oinfoglobal )
end

# General constructor for Serial pipes (e.g. input passes from element to element of the pipe according to some order)
PipeCell(x::T where T<:PTuple{AbstractCell}, order::Vector{Int}, tinfo::String="") = begin
	@assert length(x) >= length(order) "Serial pipe must contain at least as many elements as order vector."
	PipeCell(x, order, nothing, countlayers(x)+1, tinfo, oinfoglobal )
end

#####################
# PipeCell aliases  #
#####################

# Stacked pipe aliases
const PipeStacked = PipeCell{<:PTuple{AbstractCell}, Void, Void}				# Stacked pipe (generic)
const PipeStackedF = PipeCell{<:PTuple{CellFunF}, Void, Void}  					# Stacked pipe (fixed Cells)
const PipeStackedU = PipeCell{<:PTuple{CellFunU}, Void, Void}  					# Stacked pipe (untrained Cells)
const PipeStackedT =  PipeCell{<:PTuple{CellFunT}, Void, Void} 					# Stacked pipe (trained Cells)

# Parallel pipe aliases
const PipeParallel = PipeCell{<:PTuple{AbstractCell}, <:SortedDict, Void}			#  Parallel pipe (generic)
const PipeParallelF = PipeCell{<:PTuple{CellFunF}, <:SortedDict, Void} 				#  Parallel pipe (fixed Cells)
const PipeParallelU = PipeCell{<:PTuple{CellFunU}, <:SortedDict, Void}				#  Parallel pipe (untrained Cells)
const PipeParallelT = PipeCell{<:PTuple{CellFunT}, <:SortedDict, Void} 				#  Parallel pipe (trained Cells)

# Serial pipe aliases
const PipeSerial = PipeCell{<:PTuple{AbstractCell}, <:Vector{Int}, Void}			# Serial pipe (generic)
const PipeSerialF = PipeCell{<:PTuple{CellFunF}, <:Vector{Int}, Void}				# Serial pipe (fixed Cells)
const PipeSerialU = PipeCell{<:PTuple{CellFunU}, <:Vector{Int}, Void} 				# Serial pipe (untrained Cells)
const PipeSerialT = PipeCell{<:PTuple{CellFunT}, <:Vector{Int}, Void} 				# Serial pipe (trained Cells)

# Pipe alias
const PipeGeneric = PipeCell{<:PTuple{AbstractCell}, <:Any, <:Any}				# Any type of pipe

##############################
# Piping operators for Pipes #
##############################

# Pipe operator (stacked pipes)
|>(x::T where T, p::PipeStacked) = pipestack(ntuple(i->p[i](x), Val{length(p)})) # Pipe operators for generic stacked pipes

# Pipe operator (parallel pipes)
|>(x::T where T<:CellData, p::PipeParallel) = begin	                        		# Pipe operator for generic parallel pipes
	# Read dispatch information
	dispatch = gety!(p)
	K = collect(keys(dispatch))

	# Pass parts from the input data
	pipeparallel(ntuple( (i,K=K,dispatch=dispatch)->p[K[i]](varsubset(x,dispatch[K[i]])), Val{length(p)}), dispatch)
end

# Pipe operator (parallel pipes)
|>(x::T where T, p::PipeParallel) = begin	                        			# Pipe operator for generic parallel pipes
	# Read dispatch information
	dispatch = gety!(p)
	K = collect(keys(dispatch))

	# Pass parts from the input data
	pipeparallel(ntuple( (i,K=K,dispatch=dispatch)->p[K[i]](x[dispatch[K[i]]]), Val{length(p)}), dispatch)
end

# Pipe operator (serial pipes)
|>(x::T where T, p::Union{PipeSerialF, PipeSerialT}) = begin					# Pipe operator for fixed or trained serial pipes (e.g. pass data through the pipe's cells)
    	# Read order information
	order = gety!(p)

	# Start passing data through pipe
	xout = x;                                                                               # Initialize output with the input
    	for i in order		                                                        	# Pass data through each stage of the pipe
        	xout = xout |> p[i]
    	end
    	return xout #data
end

|>(x::T where T, p::PipeSerialU) = begin                                                        # Pipe operator for untrained serial pipes (e.g. train function cells sequentially)
    	# Read order information
	order = gety!(p)

	tp = pipeserial((x |> p[order[1]],))

	# Train first element of the untrained serial pipe and create serial pipe
	for i in order[2:end]		                                                        # Iterate through the rest of the untrained pipe elements and:
        	inp = x |> tp                                                                   # - run data through the trained pipe
        	out = inp |> p[i]                                                               # - train current pipe element
		tp = pipeserial((getx!(tp)...,out))                                             # - add newly trained element in the pipe
    	end
    return tp #trained serial pipe
end

|>(x::T where T, p::PipeSerial) = begin                                                        	# Pipe operator for untrained serial pipes
    	# Read order information
	order = gety!(p)

    	inp = x
	outp = x
	ndo = 0; # (non-data-output) its value indicates if pipe was fully processed or not
	processed_Cells = ()

	# Iterate through the pipe and multiply elements (i order, j just an index in the order vector)
	for (j,i) in enumerate(order)

		# Pass data through each element of the pipe
		outp = inp |> p[i]
		#info("\t\t Pipe element $i done...")

		# Check if output at this stage is data or not
		if ~(outp isa CellFun) && ~(outp isa PipeGeneric)
			processed_Cells = (processed_Cells...,p[i])				# Output is data, pipe element was fixed/trained, insert as is
			inp = outp				  				# former output becomes next cells input
		else 										# Output is not data, store output, the rest of the elements are copied in the output
			processed_Cells = (processed_Cells...,outp)
			ndo = j;
		end

		# Check whether data will not fully pass through the pipe
		if ndo > 0 && ndo < length(p)
			info("[operators] Pipe was partially processed ($(j)/$(length(p)) elements)." )
			# Copy the rest of the pipe unchanged
			for k = j+1:length(order)
				processed_Cells = (processed_Cells..., p[order[k]])
			end
			break;
		end
    	end

	# If at all stages data was returned, return data, else return processed serial pipe
	if ndo == 0 return outp
	else return pipeserial(processed_Cells)
	end
end

##############################################################################################################################
# Functions needed to construct pipes
##############################################################################################################################

# Function that creates stacked pipes
"""
	pipestack(x)

Creates a `stacked pipe`.
"""
pipestack(x::T where T<:DataElement) = x							# leave simple stuff alone
pipestack(x::T where T<:AbstractArray) = x                       				# leave arrays alone
pipestack(x::T where T<:PTuple{DataElement}) = mat(collect(x))					# tuples of elements become 1 column matrices (e.g. 1 obs)
pipestack(x::T where T<:PTuple{AbstractArray}) = dcat(x...) 					# concatenate no matter what type of array
pipestack(x::PTuple{T} where T<:CellData) = vcat(x...) 						# concatenate only datacells of the same type
pipestack(x::T where T<:PTuple{CellData}) = flatten(x, PipeStacked)				# call flatten if the data cell types differ
pipestack(x::T where T<:PTuple{AbstractCell}) = PipeCell(x)              			# tuples of cells call the constructor (e.g. for pipes)
pipestack(x) = pipestack(map(AbstractCell, x)) 							# Make Cells out of arguments, then call pipestack again

# Function that creates parallel pipes
"""
	pipeparallel(x)

Creates a `parallel pipe`.
"""
pipeparallel(x::T where T<:DataElement, dispatch::SortedDict) = x
pipeparallel(x::T where T<:AbstractArray, dispatch::SortedDict) = x
pipeparallel(x::T where T<:PTuple{DataElement}, dispatch::SortedDict) = mat(collect(x))
pipeparallel(x::T where T<:PTuple{AbstractArray}, dispatch::SortedDict) = dcat(x...)
pipeparallel(x::PTuple{T} where T<:CellData, dispatch::SortedDict) = vcat(x...)
pipeparallel(x::T where T<:PTuple{CellData}) = flatten(x, PipeParallel)				# if no dispatch is specified, (e.g. hcat, concatenate data)
pipeparallel(x::T where T<:PTuple{CellData}, dispatch::SortedDict) = flatten(x, PipeStacked)	# if dispatch is specified (e.g. output of pipe, concatenate variables)
pipeparallel(x::T where T<:PTuple{AbstractCell}, dispatch::SortedDict=SortedDict(Dict(i=>i for i in 1:length(x)))) = PipeCell(x, dispatch)
pipeparallel(x, dispatch::SortedDict=SortedDict(Dict(i=>i for i in 1:length(x)))) = pipeparallel(map(AbstractCell, x))

# Function that creates serial pipes
"""
	pipeserial(x)

Creates a `serial pipe`.
"""
#pipeserial{T<:PTuple{CellData}}(x::T, order::Vector{Int}=collect(1:length(x))) = Cell(x, order)# call flatten if the data cell types differ
pipeserial(x::T where T<:PTuple{AbstractCell}, order=collect(1:length(x))) = PipeCell(x, order)
pipeserial(x, order=collect(1:length(x))) = pipeserial(map(AbstractCell, x), order)
pipeserial(x::T where T<:PTuple, order=collect(1:length(x))) = flatten(x, PipeSerial)

#####################################
# Concatenation operators for Pipes #
#####################################
vcat(ac::T where T<:AbstractCell, args...) = pipestack((ac, args...))	                        # Vertical concatenation creates stacked pipes
vcat(ac::T...)  where T<:AbstractCell= pipestack(ac)                   	  			# Vertical concatenation creates stacked pipes
hcat(ac::T where T<:AbstractCell, args...) = pipeparallel((ac, args...))                       	# Horizontal concatenation creates parallel pipes
hcat(ac::T...) where T<:AbstractCell = pipeparallel(ac)                				# Horizontal concatenation creates parallel pipes
+(ac::T where T<:AbstractCell, args...) = pipeserial((ac, args...))                            	# Summation operator creates serial pipes
+(ac::T...) where T<:AbstractCell = pipeserial(ac)             					# Summation operator creates serial pipes

######################
# Indexing for Pipes #
######################
# Generic setindex! for data cells
setindex!(c::T where T<:PipeGeneric, cells, inds...) = setindex!(getx!(c), cells, inds...)

getindex(ac::T where T<:PipeGeneric, i::Int64) = getx!(ac)[i]
getindex(ac::T where T<:PipeStacked, i::Union{Vector{Int},UnitRange{Int}}) = PipeCell(getx!(ac)[i])
getindex(ac::T where T<:PipeParallel, i::Union{Vector{Int}, UnitRange{Int}}) = PipeCell(getx!(ac)[i], SortedDict(Dict(j=>ac.y[j] for j in i)))
getindex(ac::T where T<:PipeSerial, i::Union{Vector{Int}, UnitRange{Int}}) = PipeCell(getx!(ac)[i], collect(1:length(i)) )

#################################
# Iteration interface for pipes #
#################################
start(p::PipeGeneric) = 1
next(p::PipeGeneric, state) = (getx!(p)[state], state+1)
done(p::PipeGeneric, state) = state > length(getx!(p))
endof(p::PipeGeneric) = length(getx!(p))
eltype(p::PipeGeneric) = eltype(getx!(p))
length(p::PipeGeneric) = length(getx!(p))

"""
	flatten(x, t, [tn=t])

Low-level function used for generic pipe creation. `x` should be a `PTuple{<:AbstractCell}` and
`t` one of the aliases `PipeStacked`, `PipeParalel` or `PipeSerial`. Any element of `x` that
is not a `Cell` will be automatically transformed into one. `tn` is used in recursing the pipe
and is used to match the types of parent and nested pipes (e.g. if they match one can flatten them).
"""
flatten(v, op, newop=op) = begin
	#info("!")
	out = ()								# Result of the flattened v
	pipeout = ()								# Result of the flattened Pipe
	for c in v								# Loop through all elements of the Pipe/Tuple
		# Recursion
		if c isa CellFun || c isa CellData
			out=(out..., c)						# If current element is a function or data cell, add to list of cells
		elseif c isa PipeGeneric					# If it is a Pipe, flatten it (if possible)
			# Once entering the nested pipe, change the op
			# so that it matches the pipe type
			arg = ()
			if c isa PipeStacked
				newop = PipeStacked
			elseif c isa PipeParallel
				newop = PipeParallel
				arg = (SortedDict(Dict(i=>i for i in 1:length(c))),)
			elseif c isa PipeSerial
				newop = PipeSerial
				arg = ([i for i in 1:length(c)],)
			else
				error("[flatten]$(typeof(v)) cannot be part of a pipe")
			end
			pipeout= flatten(c,newop)				# Call recursively function, specifying the type

			if op<:newop 						# If Pipe element is the same as its parent Pipe
				out=(out..., pipeout...)			# Concatenate their elements
			else
				out = (out...,newop((pipeout,arg...))) 		# Otherwise, concatenate elements of parent with
			end							# flattened Pipe element
		else								# If element is not a Cell, trasform it into one
			nc = AbstractCell(c)
			out=(out...,nc)
		end
	end

	# Check if we are at the top layer (e.g. v is a Tuple) and if so, construct Pipe;
	# If not at the top layer (e.g. in a pipe, return flattened pipe)
	if v isa Tuple
		try
			if all(c isa CellData for c in v)
				# Check if all label information is identical accross cells
				# - if so, add it to the flattened dataset,
				# - otherwise, concatenate only data content and issue warning
				if !all( isequal(gety!(v[i]),gety!(v[i+1])) for i in 1:length(v)-1 )
				  	labels = nothing
					#warn("[flatten] Label information will be discarded (inconsistent across elements).")
				else
					labels = gety!(v[1])
				end

				# Return concatenated datacell
				if (length(v) > 1)
					if op<:PipeStacked
						return datacell(dcat(getx!.(v)), labels)
					elseif op<:PipeParallel
						return datacell(ocat(getx!.(v)), labels)
					else
						error("[flatten] Could not concatenate data cells, unsupported op=$(op)")
					end
				else
					return datacell(getx!(v[1]), labels)
				end
			end
		catch
			error("[flatten] Flattening failed. Possible reasons:\n
	 			 - inconsistent dimensions or data being concatenated
				 - label size did not match data size. ")
		end

		if op<:PipeStacked arg = ()
		elseif op<:PipeParallel arg = (SortedDict(Dict(i=>i for i in 1:length(out))),)
		elseif op<:PipeSerial arg = ([i for i in 1:length(out)],)
		else error("[flatten] Wrong operation.")
		end
		warn("[flatten] Any custom dispatch/order information for parallel/serial pipes is lost.")
		return op((out,arg...))
	else
		return out
	end
end

ext: jl


text:
function maketempdir()::String
    path::String = mktempdir()
    atexit(() -> rm(path; force = true, recursive = true,))
    return path
end


text:
using Revise
using POMDPs
using Random
using DeepQLearning
using POMDPModels
using POMDPSimulators
using RLInterface
using POMDPPolicies
using Test
using Flux

include("test_env.jl")

mdp = TestMDP((5,5), 4, 6)
mdp = SimpleGridWorld()

env = MDPEnvironment(mdp)
input_dims = obs_dimensions(env)
flattened_input_dims = reduce(*, input_dims)
output_dims = length(actions(env))
model = Chain(x->flattenbatch(x), Dense(flattened_input_dims, 8), Dense(8, output_dims))
# model = Chain(x->flattenbatch(x), LSTM(flattened_input_dims, 32), Dense(32, output_dims))
solver = DeepQLearningSolver(qnetwork = model, prioritized_replay=true, max_steps=3000, learning_rate=0.005,log_freq=500,
                             recurrence=false,trace_length=6, double_q=false, dueling=true)

@time policy = solve(solver, env)

value(policy, [1.,2.])
actionvalues(policy, [1.,2.])

sim = RolloutSimulator(rng=MersenneTwister(1), max_steps=10)
r_drqn =  simulate(sim, mdp, policy, rng)

### Hyperparameters
buffer_size = 3000
batch_size = 32
train_start = 200
learning_rate = 1e-4
trace_length = 6

env = MDPEnvironment(mdp)

### MODEL BUILDING
input_dims = obs_dimensions(env)
flattened_input_dims = reduce(*, input_dims)
output_dims = length(actions(env))
model = Chain(x->flattenbatch(x), LSTM(flattened_input_dims, 32), Dense(32, output_dims))

active_q = model
target_q = deepcopy(active_q)

### Replay Buffer (unchanged from previous implementation)
replay = EpisodeReplayBuffer(env, solver.buffer_size, solver.batch_size, solver.trace_length)
DeepQLearning.populate_replay_buffer!(replay, env, max_pop=solver.train_start)

### Batch Train, in Flux, the batch size is the last dimension
s_batch, a_batch, r_batch, sp_batch, done_batch, trace_mask_batch = DeepQLearning.sample(replay)

s_batch = batch_trajectories(s_batch, solver.trace_length, solver.batch_size)
a_batch = batch_trajectories(a_batch, solver.trace_length, solver.batch_size)
r_batch = batch_trajectories(r_batch, solver.trace_length, solver.batch_size)
sp_batch = batch_trajectories(sp_batch, solver.trace_length, solver.batch_size)
done_batch = batch_trajectories(done_batch, solver.trace_length, solver.batch_size)
trace_mask_batch = batch_trajectories(trace_mask_batch, solver.trace_length, solver.batch_size)

q_values = active_q.(s_batch) # vector of size trace_length nactions x batch_size
q_sa = [zeros(eltype(q_values[1]), solver.batch_size) for i=1:solver.trace_length]
for i=1:solver.trace_length  # there might be a more elegant way of doing this
    for j=1:solver.batch_size
        if a_batch[i][j] != 0
            q_sa[i][j] = q_values[i][a_batch[i][j], j]
        end
    end
end
q_sp_max = vec.(maximum.(target_q.(sp_batch), dims=1))
q_targets = Vector{eltype(q_sa)}(undef, solver.trace_length)
for i=1:solver.trace_length
    q_targets[i] = r_batch[i] .+ (1.0 .- done_batch[i]).*discount(env.problem).*q_sp_max[i]
end

td_tracked = broadcast((x,y) -> x.*y, trace_mask_batch, q_sa - q_targets)
loss_tracked = loss.(td_tracked)
Flux.reset!(active_q)
Flux.truncate!(active_q)
Flux.reset!(target_q)
Flux.truncate!(target_q)

### DQN

active_q = model
target_q = deepcopy(model)

# s_batch is of size (intput_dims..., bs)
o = reset!(env)
model(o[:])

model(flatten(s_batch))

q_values = active_q(flatten_batch(s_batch)) # nactions x batch_size
q_sa = [q_values[a_batch[i], i] for i=1:batch_size]
q_sp_max = @view maximum(target_q(flatten_batch(sp_batch)), dims=1)[:]
q_targets = r_batch .+ (1.0 .- done_batch).*discount(env.problem).*q_sp_max # nactions x batch_size

function loss(q_sa, q_targets)
    mean(huber_loss.(q_sa - q_targets))
end

l, td = loss(q_sa, q_targets)

Flux.data(l)

optimizer = ADAM(Flux.params(active_q), 1e-3)

# use deep copy to update the target network

# use Flux.reset to reset RNN if necessary


text:
@testset "mean" begin
    rng = MersenneTwister(999)
    n = 9

    @testset "Basic" begin
        rrule_test(
            mean,
            randn(rng),
            (randn(rng, n),
            randn(rng, n))
        )
    end

    @testset "with function arg" begin
        rrule_test(
            mean,
            randn(rng),
            (abs2, nothing),
            (randn(rng, n),
            randn(rng, n))
        )
    end

    @testset "with dims kwargs" begin
        X = randn(rng, n, n+1)
        y, mean_pullback = rrule(mean, X; dims=1)
        ȳ = randn(rng, size(y))
        _, dX = mean_pullback(ȳ)
        X̄_ad = extern(dX)
        X̄_fd = j′vp(_fdm, x->mean(x, dims=1), ȳ, X)
        @test X̄_ad ≈ X̄_fd rtol=1e-9 atol=1e-9
    end
end


text:
module BenchmarkLite

import Base: size, length, get, getindex, run, show

export Proc, BenchmarkTable, BenchmarkEntry, cfgname, procname

##### Types

## The abstract type to represent a (typed) procedure
#
# Following methods should be defined on each subtype of Proc
#
# - string(proc):          get a description
# - length(proc, cfg):     the problem size (e.g. the number of elements to process)
# - isvalid(proc, cfg):    test whether a config is valid for the given proc
#
# - s = start(proc, cfg):  set up (not count in the run-time)
# - run(proc, cfg, s):     run the procedure under the given config
# - done(proc, cfg, s):    tear down (not count in the run-time)
#
abstract Proc

##### Results

type BenchmarkTable
    cfgs::Vector        # length = m
    procs::Vector{Proc} # length = n
    plens::Matrix{Int}      # m-by-n vector, proc lengths
    nruns::Matrix{Int}      # m-by-n matrix, number of repeat times (in measuring stage)
    etime::Matrix{Float64}  # m-by-n matrix, elapsed time (in measuring stage)
end

function BenchmarkTable(cfgs::Vector, procs::Vector{Proc})
    m = length(cfgs)
    n = length(procs)
    plens = zeros(m, n)
    for (i, c) in enumerate(cfgs), (j, p) in enumerate(procs)
        plens[i,j] = length(p, c)
    end
    nruns = zeros(Int, m, n)
    etime = zeros(m, n)
    BenchmarkTable(cfgs, procs, plens, nruns, etime)
end

cfgname(bt::BenchmarkTable, i::Integer) = string(bt.cfgs[i])
procname(bt::BenchmarkTable, i::Integer) = string(bt.procs[i])

immutable BenchmarkEntry
    plen::Int       # proc length
    nruns::Int      # number of repeat times
    etime::Float64  # elapsed time (in seconds)
end

size(bt::BenchmarkTable) = (length(bt.cfgs), length(bt.procs))

size(bt::BenchmarkTable, i::Integer) = i == 1 ? length(bt.cfgs) :
                                       i == 2 ? length(bt.procs) : 1

getindex(bt::BenchmarkTable, i::Integer, j::Integer) =
    BenchmarkEntry(bt.plens[i,j], bt.nruns[i,j], bt.etime[i,j])

##### Show Results

abstract BenchmarkUnit

type Sec <: BenchmarkUnit end
type Msec <: BenchmarkUnit end
type Usec <: BenchmarkUnit end
type Nsec <: BenchmarkUnit end
type Ups <: BenchmarkUnit end
type Kps <: BenchmarkUnit end
type Mps <: BenchmarkUnit end
type Gps <: BenchmarkUnit end

get(e::BenchmarkEntry, unit::Sec) = e.etime / e.nruns
get(e::BenchmarkEntry, unit::Msec) = e.etime * 1.0e3 / e.nruns
get(e::BenchmarkEntry, unit::Usec) = e.etime * 1.0e6 / e.nruns
get(e::BenchmarkEntry, unit::Nsec) = e.etime * 1.0e9 / e.nruns
get(e::BenchmarkEntry, unit::Ups) = (1.0 * e.plen * e.nruns) / e.etime
get(e::BenchmarkEntry, unit::Kps) = (1.0e-3 * e.plen * e.nruns) / e.etime
get(e::BenchmarkEntry, unit::Mps) = (1.0e-6 * e.plen * e.nruns) / e.etime
get(e::BenchmarkEntry, unit::Gps) = (1.0e-9 * e.plen * e.nruns) / e.etime

function _show_table(bt::BenchmarkTable, unit::BenchmarkUnit, cfghead::String)
    m = length(bt.cfgs)
    n = length(bt.procs)
    S = Array(String, m+1, n+1)
    S[1,1] = cfghead
    for j=1:n; S[1,j+1] = procname(bt,j); end
    for i=1:m; S[i+1,1] = cfgname(bt,i); end

    for j=1:n, i=1:m
        v = get(bt[i,j], unit)::Float64
        S[i+1,j+1] = @sprintf("%.4f", v)
    end
    return S
end

function show(io::IO, bt::BenchmarkTable; unit::Symbol=:sec, cfghead="config")
    # getting all strings first
    u = unit == :sec ? Sec() :
        unit == :msec ? Msec() :
        unit == :usec ? Usec() :
        unit == :nsec ? Nsec() :
        unit == :kps ? Kps() :
        unit == :mps ? Mps() :
        unit == :gps ? Gps() :
        unit == :ups ? Ups() :
        error("Invalid unit value :$(unit).")
    S = _show_table(bt, u, cfghead)
    nrows, ncols = size(S)

    # calculate the width of each column
    colwids = Array(Int, ncols)
    for j = 1:ncols
        w = 1
        for i = 1:nrows
            wij = length(S[i,j])
            if wij > w; w = wij; end
        end
        colwids[j] = w
    end

    # print the table
    println(io, "$(typeof(bt)) [unit = $unit]")
    for i = 1:nrows
        # first column
        print(io, rpad(S[i,1], colwids[1]))
        print(io, " |  ")
        # remaining columns
        for j = 2:ncols
            print(io, lpad(S[i,j], colwids[j]))
            print(io, "  ")
        end
        println(io)

        if i == 1
            println(io, repeat("-", sum(colwids) + 2 * ncols + 2))
        end
    end
end

show(bt::BenchmarkTable; unit::Symbol=:sec, cfghead="config") =
    show(STDOUT, bt; unit=unit, cfghead=cfghead)

function Base.writecsv(io::IO, bt::BenchmarkTable)
    println(io, "proc,cfg,length,nruns,elapsed")
    for (j, p) in enumerate(bt.procs)
        pname = string(p)
        for (i, c) in enumerate(bt.cfgs)
            e = bt[i,j]
            println(io, join({repr(pname), repr(string(c)), e.plen, e.nruns, e.etime}, ","))
        end
    end
end

##### Run Benchmarks

function run{P<:Proc}(p::P, cfg;
                      nruns::Int = 0, duration::Float64=1.0, allowgc::Bool=true)
    # Run a procedure under a certain config
    #
    # Arguments:
    #
    #   - p:        the procedure to be tested
    #   - cfg:      the configuration to be tested
    #   - nruns:    the number of repeating times (in measuring stage)
    #
    # Keyword arguments:
    #
    #   - duration: the rough duration of the whole multi-run process
    #
    #   When nruns > 0, it runs a fixed number of times in measuring stage,
    #   otherwise, it uses a probing stage to roughly measure the runtime
    #   for each run, and uses it to determine nruns.
    #
    # Returns:  (nruns, etime)
    #
    #   - nruns:    the actual number of running times (in measuring stage)
    #   - etime:    the elapsed seconds (in measuring stage)
    #

    # check validity
    if !isvalid(p, cfg)
        return (0, 0.0)
    end

    # set up
    s = start(p, cfg)

    # warming
    run(p, cfg, s)

    if !allowgc
        gc_disable()
    end

    # probing
    if nruns <= 0
        et = @elapsed run(p, cfg, s)
        # if et is very short, we will perform more accurate probing
        if et < duration / 500
            nr = max(iceil(duration / 500), 2)
            et2 = @elapsed for i=1:nr
                run(p, cfg, s)
            end
            et = et2 / nr
        end
        nruns = iceil(duration / et)
    end

    # measuring
    etime = @elapsed for i = 1:nruns
        run(p, cfg, s)
    end

    # tear down
    done(p, cfg, s)

    if !allowgc
        gc_enable()
    end

    return (nruns, etime)
end

function run(procs::Vector{Proc}, cfgs::Vector;
             duration::Float64=1.0, verbose::Int=2, logger::IO=STDOUT)
    # Run a list of procedures against a list of configurations
    #
    # Arguments:
    #
    #   - procs:    the list of procedures to run
    #   - cfgs:     the list of configurations to run
    #
    # Keyword arguments:
    #
    #   - duration: the duration for each procedure on each config
    #   - verbose:  0 - prints nothing
    #               1 - prints one line for each procedure
    #               2 - prints one line for each procedure on each config
    #   - logger:   the IO to which the information is sent
    #
    # Returns: a BenchmarkTable instance capturing all results
    #

    bt = BenchmarkTable(cfgs, procs)
    m = length(cfgs)
    n = length(procs)
    for (j, p) in enumerate(procs)
        procname = string(p)
        verbose >= 1 && println(logger, "Benchmarking $procname ...")

        for (i, cfg) in enumerate(cfgs)
            cfgname = string(cfg)
            (nr, et) = run(p, cfg; duration=duration)

            bt.nruns[i, j] = nr
            bt.etime[i, j] = et

            verbose >= 2 &&
                println(logger, "  $procname with cfg = $cfgname: nruns = $nr, elapsed = $et secs")
        end
    end
    return bt
end

end # module


text:
module Day21

using Memoize

import ..data_dir # from parent module
input = read(joinpath(data_dir, "day21"), String)

export part1, part2

function parse_player(player_string::AbstractString)
    numbers = [match.match for match in eachmatch(r"[[:digit:]]+", player_string)]
    Player(parse.(Int64, numbers)...)
end

function parse_players(input::AbstractString)
    player_strings = split(input, "\n", keepempty = false)
    parse_player.(player_strings)
end

mutable struct Player
    number::Int64
    position::Int64
    score::Int64
end
Player(number, position) = Player(number, position, 0)
Player(number) = Player(number, 0, 0)

new_position(position::Int64, movement::Int64) = mod(position + movement - 1, 10) + 1

function move!(player::Player, movement::Int64)
    player.position = new_position(player.position, movement)
    player.score += player.position
end

mutable struct DeterministicDice
    sequence::Vector{Int64}
    rolls::Int64
end
DeterministicDice(size::Int64) = DeterministicDice(collect(1:size), 0)

function roll!(dice::DeterministicDice)
    if dice.rolls == length(dice.sequence)
        dice.rolls = 0
    end
    dice.rolls += 1
    return dice.sequence[dice.rolls]
end

function roll!(dice::DeterministicDice, times::Int64)
    sum(roll!(dice) for i = 1:times)
end

function part1(input = input)
    players = parse_players(input)
    sort!(players, by = p -> p.number)
    dice = DeterministicDice(100)
    n_rolls = 0
    while true
        for player in players
            move!(player, roll!(dice, 3))
            n_rolls += 3
            if player.score >= 1000
                losing_player = first(filter(p -> p.number != player.number, players))
                return losing_player.score * n_rolls
            end
        end
    end
end

# enumerating the possible outcomes of rolling 3d3, the frequencies of each total:
const dice_sums = [(3, 1), (4, 3), (5, 6), (6, 7), (7, 6), (8, 3), (9, 1)]

# adapted from code by https://github.com/AxlLind
# https://github.com/AxlLind/AdventOfCode2021/blob/main/src/bin/21.rs
@memoize function dirac_dice_game(
    this_player_position::Int64,
    other_player_position::Int64,
    this_player_score::Int64 = 0,
    other_player_score::Int64 = 0
)
    if other_player_score >= 21
        return 0, 1 # other player wins
    end

    score = (0, 0)
    for (roll, times) in dice_sums
        this_player_new_position = new_position(this_player_position, roll)
        this_player_new_score = this_player_score + this_player_new_position
        # Note that the players alternate below
        other_player_wins, this_player_wins = dirac_dice_game(
            other_player_position,
            this_player_new_position,
            other_player_score,
            this_player_new_score
        )
        score = score .+ times .* (this_player_wins, other_player_wins)
    end

    score
end

function part2(input = input)
    players = parse_players(input)
    sort!(players, by = p -> p.number)
    wins = dirac_dice_game(players[1].position, players[2].position)
    maximum(wins)
end

end #module


text:
using PyPlot, DelimitedFiles, Random, Statistics, Test, LinearAlgebra
using Revise, transD_GP

Random.seed!(23)

zfixed = [-1e5]
ρfixed = [1e12]
nmax = 100

zstart = 0.0
extendfrac, dz = 1.06, 2.
zall, znall, zboundaries = transD_GP.setupz(zstart, extendfrac, dz=dz, n=40)
z, ρ, nfixed = transD_GP.makezρ(zboundaries; zfixed=zfixed, ρfixed=ρfixed)
@info z

## geometry parameters for tempest
zTx = -120
zRx = -80
x_rx = -115.0
y_rx = 0.
rx_roll = 0.
rx_pitch = 0.
rx_yaw = 0.
tx_roll = 0.
tx_pitch = 0.
tx_yaw = 0.
addprimary = true

#current ramp and gate times
ramp =  [  # -0.0400066666667    0.5
			-0.0400000000000    0.0
			-0.0399933333333   -0.5
			-0.0200066666667   -0.5
			-0.0200000000000    0.0
			-0.0199933333333    0.5
			-0.0000066666667    0.5
			 0.0000000000000    0.0
			 0.0000066666667   -0.5]

times = sqrt.([
0.0000066667*0.0000200000
0.0000333333*0.0000466667
0.0000600000*0.0000733333
0.0000866667*0.0001266667
0.0001400000*0.0002066667
0.0002200000*0.0003400000
0.0003533333*0.0005533333
0.0005666667*0.0008733333
0.0008866667*0.0013533333
0.0013666667*0.0021000000
0.0021133333*0.0032733333
0.0032866667*0.0051133333
0.0051266667*0.0079933333
0.0080066667*0.0123933333
0.0124066667*0.0199933333
])

## fill in detail in ohm-m
ρ[(z.>=zstart) .& (z.<50)] .= 20.
ρ[(z.>=50) .&(z.<80)] .= 1
ρ[(z.>=80) .&(z.<100)] .= 20
ρ[(z.>=100) .&(z.<200)] .= 50
ρ[(z.>=200) .&(z.<250)] .= 80
ρ[(z.>=250)]            .= 150

## create tempest operator
tempest = transD_GP.TEMPEST1DInversion.Bfield(
    zTx = zTx, zRx = zRx, x_rx = x_rx, y_rx = y_rx,
    rx_roll = rx_roll, rx_pitch = rx_pitch, rx_yaw = rx_yaw,
    tx_roll = tx_roll, tx_pitch = tx_pitch, tx_yaw = tx_yaw,
	ramp = ramp, times = times,
	z=z,
	ρ=ρ,
	addprimary = addprimary
)

transD_GP.TEMPEST1DInversion.plotmodelfield!(tempest, z, ρ)
## update it
zTx2 = zTx + 20
zRx2 = zRx + 5
x_rx2 = x_rx - 10
y_rx2 = y_rx - 5

rx_roll2 = rx_roll + 3
rx_pitch2 = rx_pitch - 10
rx_yaw2 = rx_yaw + 3

tx_roll2 = tx_roll - 15
tx_pitch2 = tx_pitch - 4
tx_yaw2 = tx_yaw - 3
geovec = [zTx2, zRx2, x_rx2, y_rx2,
		  rx_roll2, rx_pitch2, rx_yaw2,
		  tx_roll2, tx_pitch2, tx_yaw2]

transD_GP.TEMPEST1DInversion.update_geometry(tempest, geovec)
transD_GP.TEMPEST1DInversion.getfieldTD!(tempest, z, ρ)
# make a new operator with the updated values as initial values
tempest2 = transD_GP.TEMPEST1DInversion.Bfield(
    zTx = zTx2, zRx = zRx2, x_rx = x_rx2, y_rx = y_rx2,
    rx_roll = rx_roll2, rx_pitch = rx_pitch2, rx_yaw = rx_yaw2,
    tx_roll = tx_roll2, tx_pitch = tx_pitch2, tx_yaw = tx_yaw2,
	ramp = ramp, times = times,
	z=z,
	ρ=ρ,
	addprimary = addprimary
)

# plot new field
transD_GP.TEMPEST1DInversion.getfieldTD!(tempest2, z, ρ)
μ = transD_GP.AEM_VMD_HMD.μ
ax = gcf().axes
ax[2].loglog(tempest2.F.times, abs.(μ*tempest2.Hz)*1e15, label="Bz_new", "*-")
ax[2].loglog(tempest2.F.times, abs.(μ*tempest2.Hx)*1e15, label="Bx_new", "*-")
# plot updated field
ax[2].loglog(tempest.F.times, abs.(μ*tempest.Hz)*1e15, label="Bz_up", "D--", alpha=0.7)
ax[2].loglog(tempest.F.times, abs.(μ*tempest.Hx)*1e15, label="Bx_up", "D--", alpha=0.7)
legend()

@testset "checking tempest updates" begin
	@test norm(tempest.Hz - tempest2.Hz) < 1e-12
	@test norm(tempest.Hx - tempest2.Hx) < 1e-12
end


text:
# find the maximal end component in an MDP

# find the maximal accepting end components
# algorithm 47 from "Principle of model checking"
# an end component is defined as a set of states

function maximal_end_components(mdp::M; verbose=false) where {M <: Union{MDP, POMDP}}
    verbose ? println("Building graph from mdp ... \n") : nothing
    mdpg = mdp_to_graph(mdp)
    state_space = states(mdp)
    MEC = Vector{Vector{Int64}}(undef, 0)
    MECnew = Vector{Vector{Int64}}(undef, 0)
    push!(MECnew, 1:length(state_space))# initialize
    DEBUG_STEP = 1 #XXX
    DEBUG_MAX_STEPS = 100 # XXX
    verbose ? println("Computing Maximal End Components ... \n") : nothing
    while MEC != MECnew
        MEC = deepcopy(MECnew)
        MECnew = Vector{Vector{Int64}}(undef, 0)
        for sub_state_space in MEC
            # build subgraph that corresponds to sub_state_space, subgraph of mdpgraph, where all actions moving outside of ss are removed
            sub_g = sub_mdp(mdp, sub_state_space, state_space)
            scc = strongly_connected_components(sub_g)
            for component in scc
                # map bag to state indices
                stateindex_component = sub_state_space[component]
                # check that component is not trivial (is at least a cycle)
                if length(component) == 1 && !has_edge(sub_g, component[1], component[1])
                    continue # skip trivial components
                end
                push!(MECnew, stateindex_component)
            end
        end
        verbose ? println("finished $DEBUG_STEP step, old MEC $(length(MEC)) -> new MEC $(length(MECnew))") : nothing
        DEBUG_STEP += 1
    end
    verbose ? println("MECs length: ", length.(MEC)) : nothing
    verbose ? println("MECs computed. \n") : nothing
    return MEC
end

# convert MDP to graph
function mdp_to_graph(mdp::M) where {M <: Union{MDP, POMDP}}
    g = DiGraph(length(states(mdp)))
    for (i, s) in enumerate(ordered_states(mdp))
        si = stateindex(mdp, s)
        if isterminal(mdp, s)
            add_edge!(g, si, si)
        end
        for a in actions(mdp, s)
            d = transition(mdp, s, a)
            for (sp, p) in weighted_iterator(d)
                if !(p ≈ 0.)
                   spi = stateindex(mdp, sp)
                   add_edge!(g, si, spi)
                end
            end
        end
    end
    return g
end

# return a graph of the subMDP defined by ss
function sub_mdp(mdp::M, state_indices::Vector{Int64}, state_space) where {S, M<:Union{MDP,POMDP}}
    g = DiGraph(length(state_indices))
    for (i, si) in enumerate(state_indices)
        s = state_space[si]
        if isterminal(mdp, s)
            add_edge!(g, i, i)
        end
        for a in actions(mdp, s)
            d = transition(mdp, s, a)
            is_action_valid = true
            for (sp, p) in weighted_iterator(d)
                if !(p ≈ 0.)
                    spi = stateindex(mdp, sp)
                    if !in(spi, state_indices)
                        is_action_valid = false
                    end
                end
            end
            if is_action_valid
                for (sp, p) in weighted_iterator(d)
                    if !(p≈0.)
                        spi = stateindex(mdp, sp)
                        j = findfirst(x->x==spi, state_indices)
                        add_edge!(g, i, j)
                    end
                end
            end
        end
    end
    return g
end


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

struct Node
    state::Int
    edges::Vector{Tuple{Edge,Node}}
end

function Node(state::Int)
    return Node(state, Tuple{Edge,Node}[])
end

function Base.show(io::IO, node::Node)
    print(io, summary(node), "(<state=$(node.state),#edges=$(length(node.edges))>)")
end

function findedge(s::Node, b::UInt8)
    for (e, t) in s.edges
        if b ∈ e.labels
            return (e, t)
        end
    end
    error("$(b) ∈ label not found")
end

struct Machine
    start::Node
    states::UnitRange{Int}
    start_state::Int
    final_states::Set{Int}
    eof_actions::Dict{Int,ActionList}
end

function Base.show(io::IO, machine::Machine)
    print(io, summary(machine), "(<states=", machine.states, ",start_state=", machine.start_state, ",final_states=", machine.final_states, ">)")
end

function compile(re::RegExp.RE; optimize::Bool=true)
    dfa = nfa2dfa(remove_dead_nodes(re2nfa(re)))
    if optimize
        dfa = remove_dead_nodes(reduce_nodes(dfa))
    end
    validate(dfa)
    return dfa2machine(dfa)
end

function dfa2machine(dfa::DFA)
    newnodes = Dict{DFANode,Node}()
    new(s) = get!(() -> Node(length(newnodes) + 1), newnodes, s)
    final_states = Set{Int}()
    eof_actions = Dict{Int,ActionList}()
    for s in traverse(dfa.start)
        s′ = new(s)
        if s.final
            push!(final_states, s′.state)
            eof_actions[s′.state] = s.eof_actions
        end
        for (e, t) in s.edges
            push!(s′.edges, (e, new(t)))
        end
    end
    start = new(dfa.start)
    @assert start.state == 1
    return Machine(start, 1:length(newnodes), 1, final_states, eof_actions)
end

function execute(machine::Machine, data::Vector{UInt8})
    s = machine.start
    cs = s.state
    actions = Symbol[]
    for d in data
        try
            e, s = findedge(s, d)
            cs = s.state
            append!(actions, action_names(e.actions))
        catch ex
            if !isa(ex, ErrorException)
                rethrow()
            end
            cs = -cs
        end
    end
    if cs ∈ machine.final_states && haskey(machine.eof_actions, s.state)
        append!(actions, action_names(machine.eof_actions[s.state]))
    end
    if cs > 0
        if cs ∈ machine.final_states
            cs = 0
        else
            cs = -cs
        end
    end
    return cs, actions
end

function execute(machine::Machine, data::String)
    return execute(machine, convert(Vector{UInt8}, codeunits(data)))
end


text:
module FiniteDiff1D

using Symbolics

# Generates a method
# ∇(f_stencil, Δx_stencil, i)
function generate_stencil(n_left, n_right;staggered=false)
    @assert n_left < 0
    @assert n_right > 0
    @assert !(n_left == 0 && n_right == 0)
    staggered && error("staggered not yet supported")

    eqs = []
    k_range = n_left:n_right
    k_eqs = filter(x->x≠0, k_range)
    n_unknowns = length(k_range)
    j_range = 0:(n_unknowns-1)
    k_center = 1-n_left # example: n_left = 2: -2, -1, 0, 1, 2
    @variables α_k[k_eqs], f_k[k_range], fʲ_i[1:n_unknowns-1]
    @info "User input: (n_left, n_right, n_unknowns) = ($n_left, $n_right, $n_unknowns)"
    @info "Unknowns: $fʲ_i"
    @info "Knowns:\n$(join((α_k, f_k), "\n"))"

    # Construct equations:

    # For example: n_left = -2, n_right = 2
    #     k_range:  -2, -1, 0, 1, 2
    #     k_offset:  1,  2, 4, 5
    #     k_eqs:    -2, -1, 1, 2
    #     k_enum:    1,  2, 3, 4
    for (k_enum, k_eq) in enumerate(k_eqs)
        k_offset = k_eq < 0 ? k_enum : k_enum+1
        f_k0 = f_k[k_center]
        ∑_jk = sum(j_range) do j
            fʲ_i_temp = j==0 ? f_k0 : fʲ_i[j]

            α_k[k_enum]^j/factorial(j)*fʲ_i_temp
        end
        push!(eqs, f_k[k_offset] ~ ∑_jk)
    end
    solns = Symbolics.solve_for(eqs, fʲ_i)

    solns = map(solns) do sol
        factor(sol, f_k)
    end

    return solns, f_k, fʲ_i, k_range, n_unknowns
end

# Typically, `f_k` is factoried:
function factor(eq, vars)
    f = []
    for v in vars
        d = Dict(u => (isequal(v,u) ? 1 : 0) for u in vars)
        push!(f, v * substitute(eq, d))
    end
    +(f...)
end

end # module


text:
@testset "Constructors" begin
    @test_nowarn opt = Momentum()
    @test_nowarn opt = Adagrad()
    @test_nowarn opt = Adadelta()
    @test_nowarn opt = RMSprop()
    @test_nowarn opt = Adam()
    @test_nowarn opt = Adamax()
    @test_nowarn opt = Nadam()
    @test_nowarn opt = VanillaGradDescent()
    @test_nowarn opt = Inversedecay()
end


text:
using garw

push!(ARGS, "arg")
garw.julia_main()


text:
@testset "Parallel.Generate Reduce" begin

    function make_vec(g::AbstractGraph{T}) where T<:Integer
        return Vector{T}(undef, nv(g))
    end

    function comp_vec(x::Vector, y::Vector)
        return length(x) < length(y)
    end

    g1 = star_graph(5)

    for parallel in [:distributed, :threads]
        for g in testgraphs(g1)

            s = @inferred(LightGraphs.Parallel.generate_reduce(g, make_vec, comp_vec, 5; parallel=parallel))
            @test length(s) == 5
        end
    end
end


text:
struct RootVolume{S} <: AbstractScalarFunction
    set::S
end
Base.copy(rv::RootVolume) = rv
nth_root(volume::Volume) = RootVolume(volume.set)
Base.show(io::IO, rv::RootVolume) = print(io, "volume^(1/n)(", rv.set, ")")

# Primal:
#   set : x^T Q x ≤ 1
#   volume proportional to 1/det(Q)
#   t ≤ det(Q)^(1/n) <=> 1/t^n ≥ 1/det(Q)
#   volume proportional to 1/t^n
#   For t ≤ det(Q)^(1/n) to be tight we need to maximize `t`
#   hence we need to minimize the volume
# Dual:
#   set : x^T Q^{-1} x ≤ 1
#   volume proportional to det(Q)
#   t ≤ det(Q)^(1/n) <=> t^n ≥ det(Q)
#   volume proportional to t^n
#   For t ≤ det(Q)^(1/n) to be tight we need to maximize `t`
#   hence we need to maximize the volume
function set_space(space::Space, rv::RootVolume, model::JuMP.Model)
    if rv.set isa Ellipsoid
        rv.set.guaranteed_psd = true
    end
    sense = data(model).objective_sense
    if sense == MOI.MIN_SENSE
        return set_space(space, PrimalSpace)
    else
        # The sense cannot be FEASIBILITY_SENSE since the objective function is
        # not nothing
        @assert sense == MOI.MAX_SENSE
        return set_space(space, DualSpace)
    end
end

function ellipsoid_root_volume(model::JuMP.Model, Q::AbstractMatrix)
    n = LinearAlgebra.checksquare(Q)
    t = @variable(model, base_name="t")
    upper_tri = [Q[i, j] for j in 1:n for i in 1:j]
    @constraint(model, [t; upper_tri] in MOI.RootDetConeTriangle(n))
    return t
end

function root_volume(model::JuMP.Model, ell::Union{Sets.PolarOrNot{<:Sets.Ellipsoid},
                                                   Sets.HouseDualOf{<:Sets.AbstractEllipsoid}})
    return ellipsoid_root_volume(model, Sets.convexity_proof(ell))
end

"""
    root_volume(model::JuMP.Model,
                set::Sets.PolarOrNot{<:Sets.ConvexPolySet})

Section IV.A of [MLB05].

[MLB05] A. Magnani, S. Lall and S. Boyd.
*Tractable fitting with convex polynomials via sum-of-squares*.
Proceedings of the 44th IEEE Conference on Decision and Control, and European Control Conference 2005,
**2005**.
"""
function root_volume(model::JuMP.Model,
                     set::Sets.PolarOrNot{<:Sets.ConvexPolySet})
    if Sets.convexity_proof(set) === nothing
        error("Cannot optimize volume of non-convex polynomial sublevel set.",
              " Use PolySet(convex=true, ...)")
    end
    return ellipsoid_root_volume(model, Sets.convexity_proof(set))
end

objective_sense(::JuMP.Model, ::RootVolume) = MOI.MAX_SENSE
function objective_function(model::JuMP.Model, rv::RootVolume)
    return root_volume(model, variablify(rv.set))
end


text:
@testset "mesh" begin
    eM = eMesh{Tri,Tet}()
    push!(eM.point, SVector{3,Float64}(NaN, NaN, NaN))
    @test_throws ErrorException verify_mesh(eM)
    push!(eM.ϵ, NaN)

    append!(eM, output_eMesh_half_plane())
    mesh_remove_unused_points!(eM)
    @test length(eM.point) == 4

    append!(eM, output_eMesh_half_plane())
    mesh_inplace_rekey!(eM)
    @test length(eM.point) == 4

    @test !isempty(eM)
    @test_throws MethodError empty(eM)
    eM = as_tri_eMesh(eM)
    empty!(eM)
    @test isempty(eM)

    eM = as_tri_eMesh(output_eMesh_box())
    @test area(crop_mesh(eM, SMatrix{1,4,Float64,4}(0.0, 0.0, +1.0, 0.0))) ≈ 12.0
    @test area(crop_mesh(eM, SMatrix{1,4,Float64,4}(0.0, 0.0, -1.0, 0.0))) ≈ 12.0
    @test volume(output_eMesh_box()) ≈ 8.0
end

two = 2.0
@testset "half_plane" begin
    eM_hp_f = output_eMesh_half_plane(two, false)
    eM_hp_t = output_eMesh_half_plane(two, true)
    z_val_f = [k[3] for k = get_point(eM_hp_f)]
    z_val_t = [k[3] for k = get_point(eM_hp_t)]
    @test z_val_f == z_val_t
    @test -two == minimum(z_val_f)
    @test 0.0 == maximum(z_val_f)
    @test 4 == n_point(eM_hp_f)
    for k = 1:3
        @test -0.5 ≈ dot(get_point(eM_hp_f)[k], get_point(eM_hp_f)[mod1(k + 1, 3)])
    end
    @test n_tri(eM_hp_f) == 1
    @test n_tri(eM_hp_t) == 4
    @test n_tet(eM_hp_f) == n_tet(eM_hp_t) == 1
end

@testset "sphere" begin
    for k = 1:4
        eM_sphere = output_eMesh_sphere(two, k)
        n_zero = 0
        @test n_tri(eM_sphere) == n_tet(eM_sphere) == 20 * k^2
        for k = 1:n_point(eM_sphere)
            norm_point = norm(get_point(eM_sphere)[k])
            if norm_point == 0.0
                n_zero += 1
            else
                @test two ≈ norm_point
            end
        end
        @test n_zero == 1
        (k == 1) && @test n_point(eM_sphere) == 1 + 12
        (k == 2) && @test n_point(eM_sphere) == 1 + 12 + 30  # each edge bisected
        (k == 3) && @test n_point(eM_sphere) == 1 + 12 + 2 * 30 + 20  # each edge bisected + point in center
        (k == 4) && @test n_point(eM_sphere) == 1 + 12 + 3 * 30 + 3 * 20  # each edge bisected + two more points in center
    end
end


text:
using Pkg
Pkg.build("AprilTags")
using Documenter, AprilTags

makedocs(
    modules = [AprilTags],
    format = :html,
    sitename = "AprilTags.jl",
    pages = Any[
        "Home" => "index.md",
        "Functions" => "func_ref.md"
    ]
    # html_prettyurls = !("local" in ARGS),
    )

deploydocs(
    repo   = "github.com/JuliaRobotics/AprilTags.jl.git",
    target = "build"
)


text:
# mc2e computes energy from mel-cepstrum.
# TODO: theoretical background
function mc2e(mc::AbstractVector, α=0.35, len=256)
    # back to linear frequency domain
    c = SPTK.freqt(mc, len-1, -α)

    # compute impule response from cepsturm
    ir = c2ir(c, len)

    sumabs2(ir)
end

function mc2e(mc::AbstractMatrix, α=0.35, len=256)
    r = Array(eltype(mc), size(mc, 2))
    for i in 1:length(r)
        r[i] = mc2e(view(mc, :, i), α, len)
    end
    r
end

function mc2e(state::SpectralParamState, len=256)
    error("not implemented")
end

function mc2e{T<:Union{MelCepstrum,LinearCepstrum}}(state::SpectralParamState{T},
                                                   len=256)
    assert_not_ready_to_filt(state)

    def = paramdef(state)
    α = allpass_alpha(def)
    data = rawdata(state)
    mc2e(data, α, len)
end


text:
module StochasticNeighborhoodEmbeddings
using Base.Threads
using LinearAlgebra
using SparseArrays
using Statistics
using FillArrays
using ForwardDiff
export tsne, tsne!

# using LazyLinearAlgebra
using CovarianceFunctions
using CovarianceFunctions: Cauchy

# # IDEA: have repulsive_force, attractive_force for SNE
# fast_algorithm_min_size = 2^14 # minimum number of data points for fast algorithm

using NearestNeighbors

const AbstractVecOfVec{T} = AbstractVector{<:AbstractVector{<:T}}
const AbstractMatOrFac = Union{AbstractMatrix, Factorization}

function conditional_neighbor end
function optimize_variances end
function joint_neighbor end

include("util.jl") # definition of statistical functions and kernels
include("parameters.jl")
include("dense.jl") # standard routines based on dense calculations
include("sparse.jl") # accelerated sparse implementations
include("joint.jl")
include("gradient.jl")

# X is the data matrix whose columns are data points
# u is perplexity
# d is embedding dimensionality, usually 2 or 3
# implementation follows "VISUALIZING DATA USING T-SNE" closely
# missing: adaptive learning rate
function tsne(X::AbstractMatrix, u::Real, params::OptimizationParameters = OptimizationParameters())

    n = size(X, 2)
    Y = 1e-2randn(params.d, n) # initialize embedding points
    tsne!(Y, X, u, params)
    return Y
end

# if Y is already allocated
function tsne!(Y::AbstractMatrix, X::AbstractMatrix, u::Real, params::OptimizationParameters = OptimizationParameters())
    P = optimized_symmetrized_neighbor(X, u, params.dofast)
    tsne!(Y, P, params)
    return Y
end

function optimized_symmetrized_neighbor(X::AbstractMatrix, u::Real, dofast::Bool)
    n = size(X, 2)
    Pij = nothing # introduce symbol in scope before if statement
    k = floor(3*u) # k-nearest neighbor cut-off for sparsification of P matrix
    if dofast && k < n
        neighbors = k_nearest_neighbors(X, k)
        σ² = optimize_variances(X, u, neighbors)
        Pij = conditional_neighbor(X, σ², neighbors, gaussian, euclidean)
    else
        σ² = optimize_variances(X, u)
        Pij = conditional_neighbor(X, σ², gaussian, euclidean)
    end
    P = @. (Pij + Pij') / 2n # symmetrized distribution - technically only super-diagonal matters
end

function tsne!(Y::AbstractMatrix, P::AbstractMatrix, params::OptimizationParameters = OptimizationParameters())
    n = size(X, 2)
    G = zeros(params.d, n) # gradient of Y
    G² = zeros(params.d, n) # second moments of gradient
    if params.early_iter > 0 && (params.early_compression > 0 || params.early_exaggeration > 1)
        initialize_embedding!(Y, G, G², P, params)
    end
    @. G = 0 # clearing gradient storage between
    @. G² = 0
    tsne!(Y, G, G², P, params)
    return Y
end

function tsne!(Y::AbstractMatrix, G::AbstractMatrix, G²::AbstractMatrix,
               P::AbstractMatrix, params::OptimizationParameters = OptimizationParameters())
    tsne!(Y, G, G², P, params.max_iter, params.lr, params.beta_1, params.beta_2,
          params.weight_decay, params.eps, params.dofast, params.verbose, params.debias)
    return Y
end

# it's good practice to do a few iterations with weight decay on Y
# see 3.4 Optimization Methods for t-SNE of "VISUALIZING DATA USING T-SNE"
# called "early compression"
# "early exaggeration" multiplies all Pij by a constant,
# e.g. 4 in "VISUALIZING DATA USING T-SNE", 12 in "Accelerating t-SNE using Tree-Based Algorithms"
function initialize_embedding!(Y::AbstractMatrix, G::AbstractMatrix, G²::AbstractMatrix,
                               P::AbstractMatrix, params::OptimizationParameters)
    P .*= params.early_exaggeration # multiplying all Pij by a certain constant
    tsne!(Y, G, G², P, params.early_iter, params.early_lr, params.early_beta_1,
          params.early_beta_2, params.early_compression, params.eps, params.dofast, params.verbose, params.debias)
    P ./= params.early_exaggeration
    return Y
end

function tsne!(Y::AbstractMatrix, G::AbstractMatrix, G²::AbstractMatrix, P::AbstractMatrix,
               max_iter::Int, α::Real, β₁::Real, β₂::Real, γ::Real, ε::Real, dofast::Bool, verbose::Bool, debias::Bool)
    # G² keep track of second moments of gradient
    for t in 1:max_iter
        # NOTE: compared to "Accelerating t-SNE using Tree-Based Algorithms", Q is un-normalized
        Q, Q² = unnormalized_joint_neighbor(Y, cauchy, euclidean, dofast)
        Z = normalization_constant(Q)
        if verbose
            println("iteration ", t)
            println("divergence ", kl(P, Q, Z)) # bottleneck??
        end
        gradient!(G, G², P, Q, Q², Y, Z, β₁, β₂) # gradient + momentum increment
        αₜ = debias ? stepsize(α, β₁, β₂, t) : α
        @. Y = (1 - γ) * Y - αₜ * G / (sqrt(G²) + ε) # increment by objective + weight decay gradients
    end
    return Y
end

end # module


text:
#!/usr/bin/env julia
#=
   Author: Tim Sterne-Weiler, 3/16/2015
   e-mail: tim.sterne.weiler@utoronto.ca

   julia 0.4-dev
=#

## ### ### ### ### ### ### ### ### ### ### ### ### ## #
# ### ### ### ## INITIALIZATION  ## ### ### ### ### ##
## ### ### ### ### ### ### ### ### ### ### ### ### ## #

global head = "[aligater reclass]:"
println(STDERR, "$head Loading Packages..")

using ArgParse
using StatsBase
using Match

function parse_cmd()
  s = ArgParseSettings()

  @add_arg_table s begin
    "--geneid", "-g"
      help = "flag to collapse geneids"
      action = :store_true
    "--biotype", "-b"
      help = "flag to collapse biotypes"
      action = :store_true
    "--uniq", "-u"
      help = "filter for unique chimeras by left/right lengths and 10bp unique match of lig site"
      action = :store_true
    "--load"
      help = "previously stored background file for 'dclass' in HDF5 format"
      arg_type = ASCIIString
    "--save"
      help = "store to file the resulting 'dclass' in jlz format as calculated on the input data"
      arg_type = ASCIIString
    "--barcode"
      help = "random barcode associated with each read in jlz format (optional)"
      arg_type = ASCIIString
    "--size"
      help = "nmer size to use for index-based reclassification"
      arg_type = Int
      default = 32
  end
  return parse_args(s)
end

## ### ### ### ### ### ### ### ### ### ### ### ### ## #
# ### ### ### ### ###  FUNCTIONS ## ### ### ### ### ##
## ### ### ### ### ### ### ### ### ### ### ### ### ## #

include("dictext.jl") #--> savedict, loaddict, dush!, dinc!, dnorm!

function isUniqueJunc!{K <: AbstractString}( used::Dict{K,Bool}, seq, genes, dbar::Dict{K,K}, name )
  barcode = get(dbar, name, "")
  m = match(r"([AGCTUN]{5}_[AGCTUN]{5})", seq)
  cap = m.captures[1]
  m = match(r"([AGCTUN]+)_([AGCTUN]+)", seq)
  lLen,rLen = length(m.captures[1]), length(m.captures[2])
  geneid = sort( genes )
  key = join(geneid, ":") * "_$cap:$lLen:" * barcode
  retbool = false
  if !( haskey( used, key ) )
      used[key] = true
      retbool = true
  end
  retbool
end #--> Bool

# This function can both be used to set, and also get the class associated with two chimeric
# sequences based on window size homology.  First pass should set the database then run a
# second pass to retrieve the final results
function setClass!( dict::Dict{ASCIIString,Char}, class::Char, seqs; size=32, stepsize=floor(size/3), set=false )
  @assert( length(seqs) >= 2 )
  classHeirPair( a::Char, b::Char ) = a < b ? b : a  #--> Char
  classHeirArray( arr::Array{Char,1} ) = max( arr... ) #--> Char
#  retarray = Char[]
  retval = class
  lenA,lenB = map(length, seqs)
  # pad a sequence if it is shorter than our target window size
  seqA = seqs[1] * repeat(".", max( size - lenA, 0 ))
  seqB = seqs[2] * repeat(".", max( size - lenB, 0 ))
  # reset the lengths after padding
  lenA,lenB = length(seqA),length(seqB)
  it=max(1, Int(stepsize)) # calculate the step size, as half the window size
  # iterate through the cartesian product of windows in seqA and seqB by it step size
  for i in 1:it:(lenA-size)+1, j in 1:it:(lenB-size)+1
    winA = seqA[i:(i+size-1)] # access substrings
    winB = seqB[j:(j+size-1)]
    keyA,keyB = winA < winB ? (winA,winB) : (winB,winA)
    jointkey = keyA * "_" * keyB
    # use non-exported ht_keyindex function to avoid over indexing hash
    indx = Base.ht_keyindex( dict, jointkey )
    if indx <= 0
      if set
        dict[jointkey] = class #set val if none exists
        #length(retarray) <= 0 && push!(retarray, class)
      end
    else
      curVal = classHeirPair( dict.vals[indx], class )
      dict.vals[indx] = curVal
      #push!(retarray, curVal)
      retval = classHeirPair( retval, curVal )
    end
  end
#  @assert( length(retarray) > 0 )
#  length(retarray) > 1 ? classHeirArray( retarray ) : retarray[1]
  retval
end #--> Char

function reducegeneid( geneid, biotype, repeatname, repeatclass )
  retval = @match repeatname begin
             r"SSU-rRNA" => "18S_rRNA"
             r"LSU-rRNA" => "28S_rRNA"
             _, if ismatch(r"RNA", repeatclass) end => repeatname
             _ => geneid
           end
   (retval != geneid) && return retval

  retval = @match geneid begin
             r"RNA{0,1}5S|C15orf52" => "5S_rRNA"
             r"RNA{0,1}5-8" => "5.8S_rRNA"
             r"7SL" => "7SLRNA_srpRNA"
             r"7SK" => "7SK_RNA"
             _ => geneid
          end
  (retval != geneid) && return retval

  m = match(r"RNU(\d+)(?=(ATAC|-))", geneid)
  if m != nothing
    retval = @match m.captures[2] begin
               "ATAC" => "RNU"*m.captures[1]*m.captures[2]
               "-" => "U"*m.captures[1]*"_snRNA"
               _ => geneid
             end
  end
  (retval != geneid) && return retval

  m = match(r"RNA{0,1}(\d+)S", geneid)
  if m != nothing
    retval = m.captures[1] * "S_rRNA"
  end
  retval
end #--> ASCIIString

function reducebiotype( geneid, biotype, repeatname, repeatclass )
  biotype = replace(biotype, r"Mt-", "")
  ismatch(r"(srp|sn|r)RNA", biotype) && return(biotype)
  ismatch(r"tRNA", repeatclass) && return("tRNA")
  ismatch(r"SNOR|SCAR|^ACA", geneid) && return("snoRNA")
  ismatch(r"RMRP", geneid) && return("RNase_MRP")
  ismatch(r"RPPH", geneid) && return("RNase_P")
  ismatch(r"VT.*RNA", geneid) && return("vtRNA")
  ismatch(r"7SK_RNA", repeatname) && return(repeatname)
  ismatch(r"U(3|8|9|\d{2,+})|snoRNA", repeatname) && return("snoRNA")
  ismatch(r"snRNA", repeatname) && return("snRNA")
  if repeatclass != "NA"
    return("$biotype\_$repeatclass")
  else
    return(biotype == "NA" ? geneid : biotype)
  end
end #--> ASCIIString

function reducef( func, str )
  const geneInd = 4 # gene-id index
  const biotInd = 7 # biotype index
  const repnInd = 8 # repeat name index
  const repcInd = 9 # repeat class index

  function safesplit( arr, ind; char=':' )
    @assert( length(arr) >= ind )
    split( arr[ind], char )
  end

  gspl = safesplit( str, geneInd )
  bspl = safesplit( str, biotInd )
  nspl = safesplit( str, repnInd )
  cspl = safesplit( str, repcInd )

  lenLim = min( map(length, (gspl, bspl, nspl, cspl))... )
  res = ASCIIString[]
  for i in 1:lenLim
    push!(res, func(gspl[i], bspl[i], nspl[i], cspl[i]))
  end
  res
end #--> ASCIIString[]

function reclassReduceAndPrint( dclass::Dict, dstore::Dict, pargs, seqInd )
  for class in keys(dstore), s in dstore[class]
    seqs = masksplit(s[seqInd], '_')
    if length(seqs) > 1 # try to reclassify
      s[1] = convert(ASCIIString, string( setClass!( dclass, s[1][1], seqs, stepsize=1 ) ))
    end

    # now if collapse flags are true then try to reclassify
    if pargs["biotype"]
      redbiotypes = reducef( reducebiotype, s )
      s = [s; join(redbiotypes, ':')]
    end
    if pargs["geneid"]
      redgenes = reducef( reducegeneid, s )
      if length(redgenes) >= 2
        if length(unique(redgenes)) == 1
          s[1] = replace(s[1], r"[PI]", "S") #re-classify
        end
      end
      s = [s; join(redgenes, ':')]
    end
    # # # print to STDOUT # # #
    println( join( s, '\t' ) )
    # # # # # # # # # # # # # #
  end
end

function masksplit( str, char )
  rep = replace( str, r"[a-z]+", "" )
  split( rep, char )
end

## ### ### ### ### ### ### ### ### ### ### ### ### ## #
# ### ### ### ### ### MAIN  ### ### ### ### ### ### ##
## ### ### ### ### ### ### ### ### ### ### ### ### ## #
function main()
  pargs = parse_cmd()

  # simplify ridiculous split string type
  # in case julia changes by version
  typealias SubType typeof( split("a:b", ':') )

  djunc  = Dict{ASCIIString,Bool}()
  dstore = Dict{Char,Array{SubType,1}}()
  dclass = pargs["load"] == nothing ? Dict{ASCIIString,Char}() : loaddict(pargs["load"])
  dbarcode = pargs["barcode"] == nothing ? Dict{ASCIIString,ASCIIString}() : loaddict(pargs["barcode"])

  const seqInd = 11 # sequence index of .lig
  const geneInd = 4 # gene-id index
  const nameInd = 10 # readName index
  const ligqInd = 13 # LIGQ index

  pargs["load"] == nothing && (sizehint!(dclass, 1000000))
  # first iteration through file, store data, set structures
  for i::ASCIIString in eachline( STDIN )
    s = split(chomp(i), '\t')
    #genes = split(s[geneInd], ':')
    genes = reducef( reducegeneid, s ) # TODO: this gets called twice. fix that.

    # test if this is a unique junction/readset, short circuit by --uniq
    if !pargs["uniq"] || isUniqueJunc!( djunc, s[seqInd], genes, dbarcode, s[nameInd] )
      #set data
      @assert( length(s[1]) == 1 )
      curclass = s[1][1]
      #check ligq for intramolecular matches
      if curclass != 'S' && ismatch(s[ligqInd], r"\[\d+\]")
        s[1] = replace(s[1], r"[PI]", "S")
        curclass = 'S'
      end
      dush!( dstore, curclass, s, arrType=typeof(s) )

      # check if we have already loaded "dclass"
      seqs = masksplit(s[seqInd], '_')
      if length(seqs) > 1 && pargs["load"] == nothing
        setClass!( dclass, curclass, seqs, size=pargs["size"], set=true )
      end
    end
  end

  if pargs["save"] != nothing
#    save(pargs["save"], "dclass", dclass)
    savedict( pargs["save"], dclass )
    quit()
  end

  # finish and print.
  reclassReduceAndPrint( dclass, dstore, pargs, seqInd )
end
#########
main()
########

max_forks_repo_name: timbitz/aligater


text:
function frexp(x::Float64)
    xu = reinterpret(Uint64,x)
    k = int(xu >> 52) & 0x07ff
    if k == 0 # x is subnormal
        x == zero(x) && return x,0
        x *= 0x1p54 # normalise significand
        xu = reinterpret(Uint64,x)
        k = int(xu >> 52) & 0x07ff - 54
    elseif k == 0x07ff
        # NaN or Inf
        return x,0
    end
    k -= 1022
    xu = (xu & 0x800f_ffff_ffff_ffff) | 0x3fe0_0000_0000_0000
    reinterpret(Float64,xu), k
end

function frexp(x::Float32)
    xu = reinterpret(Uint32,x)
    k = int(xu >> 23) & 0x00ff
    if k == 0 # x is subnormal
        x == zero(x) && return x,0
        x *= float32("0x1p25") # no Float32 hex literal
        xu = reinterpret(Uint32,x)
        k = int(xu >> 23) & 0x00ff - 25
    elseif k == 0x00ff
        # NaN or Inf
        return x,0
    end
    k -= 126
    xu = (xu & 0x807f_ffff) | 0x3f00_0000
    reinterpret(Float32,xu), k
end


text:
module NlmeData

using LazyArtifacts
import Arrow

function __init__()
    global Datadir = artifact"Data"
end

const nlmecache = Dict{String, Arrow.Table}()

"""
    dataset(nm)

Return, as an `Arrow.Table`, the data set named `nm`, which can be a `String` or `Symbol`
"""
function dataset(nm::AbstractString)
    get!(nlmecache, nm) do
        path = joinpath(Datadir, nm * ".arrow")
        if !isfile(path)
            throw(ArgumentError(
                "Dataset \"$nm\" is not available.\nUse NlmeData.datasets() for available names."))
        end
        Arrow.Table(path)
    end
end
dataset(nm::Symbol) = dataset(string(nm))

"""
    datasets()

Return a vector of names of the available test data sets
"""
datasets() = first.(Base.Filesystem.splitext.(filter(endswith(".arrow"), readdir(Datadir))))

end


text:
#require(Pkg.dir("MethylUtils", "test","mlml-test-2.jl"))
#include(Pkg.dir("MethylUtils", "src","Features.jl"))

include(Pkg.dir("MethylUtils", "test","mlml-test-2.jl"))
include(Pkg.dir("MethylUtils", "src","mlml-bs-ox-typed.jl"))

require(Pkg.dir("MethylUtils", "src","Features.jl"))
require(Pkg.dir("MethylUtils", "src","mlml-bs-ox-typed.jl"))

println("Get random counts")
meth_counts = rand_bs_ox_counts(10,100,5)

bs_counts = rand_CytosineCounts(1000000,100,5)
ox_counts = rand_CytosineCounts(1000000,100,5)
println("Finished generating random counts")

#println(meth_counts)
#meth_counts
#println(bs_counts)
#println(ox_counts)
#@time a = pmap( EM,  bs_counts,ox_counts, 1:length(bs_counts) )
#println(a)
#println(length(a))

println("run mlml_run" )
@time run_mlml(bs_counts, ox_counts)


text:
function foo(; count::Bool = false)
	foo(:ANY, count=count)	 # Match any DAQ HAT ID in hatlist()
end

function foo(id::Symbol; count::Bool = false)
	number = 5
	if count
		return number
	end
	list = foo(id, number, count=count)
	return list
end

function foo(id::Symbol, number::Integer; count::Bool = false)

	if number < 1
		error("Number  must be >= 1")
	end

	# number of HATS installed
	numbermax = 5

	if count
		return numbermax
	end

	if number <= numbermax
	list = Vector{String}(undef, number)
		for
			list[i] = "$i $id"
		end
		return list
	else
		error("Requesting $number '$id' items $numbermax available")
	end
end

max_stars_repo_name: Spectrum-Tec/MccDaqHats

max_issues_repo_name: Spectrum-Tec/MccDaqHats

max_forks_repo_name: Spectrum-Tec/MccDaqHats


text:
"""
    indices for state and control trajectories over horizon T
"""
struct Indices
    x
    u
end

function indices(n, m, T;
    shift = 0)
    x = []
    u = []

    for t = 1:T
        push!(x, shift + (t - 1) * (n + m) .+ (1:n))

        t == T && continue
        push!(u, shift + (t - 1) * (n + m) + n .+ (1:m))
    end

    return Indices(x,u)
end

max_issues_repo_name: zacmanchester/DirectMotionPlanning

max_forks_repo_name: zacmanchester/DirectMotionPlanning


text:
@testset "univariate" begin
    @testset "Beta" begin
        p = Beta(2, 10)
        q = Beta(3, 5)
        @test KL(p, q) > 0
        @test KL(p, q) ≈ KL(p, q, 100_000) atol = 0.1
    end
    @testset "Exponential" begin
        p = Exponential(2.0)
        q = Exponential(3.0)
        @test KL(p, q) > 0
        @test KL(p, q) ≈ KL(p, q, 100_000) atol = 0.1
    end
    @testset "Gamma" begin
        p = Gamma(2.0, 1.0)
        q = Gamma(3.0, 2.0)
        @test KL(p, q) > 0
        @test KL(p, q) ≈ KL(p, q, 100_000) atol = 0.1
    end
    @testset "InverseGamma" begin
        p = InverseGamma(2.0, 1.0)
        q = InverseGamma(3.0, 2.0)
        @test KL(p, q) > 0
        @test KL(p, q) ≈ KL(p, q, 100_000) atol = 0.1
    end
    @testset "Normal" begin
        p = Normal(0, 1)
        q = Normal(0.5, 0.5)
        @test KL(p, q) > 0
        @test KL(p, q) ≈ KL(p, q, 100_000) atol = 0.1
    end
    @testset "Poisson" begin
        p = Poisson(4.0)
        q = Poisson(3.0)
        @test KL(p, q) > 0
        @test KL(p, q) ≈ KL(p, q, 100_000) atol = 0.1
    end

end


text:
using HypothesisTests, Base.Test
using HypothesisTests: default_tail

# Fisher, 1995 example 4.11
@test abs(pvalue(RayleighTest(0.2370, 60)) - 0.034) <= 0.001

# Fisher, 1995 example 4.12
t = RayleighTest(
    [2, 9, 18, 24, 30, 35, 35, 39, 39, 44,
    44, 49, 56, 70, 76, 76, 81, 86, 91, 112,
    121, 127, 133, 134, 138, 147, 152, 157, 166, 171,
    177, 187, 206, 210, 211, 215, 238, 246, 269, 270,
    285, 292, 305, 315, 325, 328, 329, 343, 354, 359]
    *pi/180)
@test abs(pvalue(t) - 0.20) <= 0.01
@test default_tail(t) == :both
show(IOBuffer(), t)

# Fisher, 1995 example 6.8
wind_direction_6am =
    [356, 97, 211, 232, 343, 292, 157, 302, 335, 302,
    324, 85, 324, 340, 157, 238, 254, 146, 232, 122,
    329]*pi/180
wind_direction_12pm =
    [119, 162, 221, 259, 270, 29, 97, 292, 40, 313,
    94, 45, 47, 108, 221, 270, 119, 248, 270, 45,
    23]*pi/180
t = FisherTLinearAssociation(wind_direction_6am, wind_direction_12pm)
@test abs(t.rho_t- 0.191) < 0.001
@test abs(pvalue(t) - 0.01) < 0.01
@test default_tail(t) == :both
show(IOBuffer(), t)

# Jammaladak, 2001 example 8.1
t = JammalamadakaCircularCorrelation(wind_direction_6am, wind_direction_12pm)
@test abs(t.r - 0.2704648) < 1e-7
@test abs(pvalue(t) - 0.2247383) < 1e-7
@test default_tail(t) == :both
show(IOBuffer(), t)


text:
module VisualizationTests

using Test

@testset "Vtkio" begin include("VtkioTests.jl") end

end # module


text:
using SymPy
using PyCall: PyError
using Test

import PyCall
PyCall.pyimport_conda("sympy.physics.wigner",       "sympy")
PyCall.pyimport_conda("sympy.physics.optics",       "sympy")
PyCall.pyimport_conda("sympy.physics.quantum.spin", "sympy")

import_from(sympy.physics.wigner)
import_from(sympy.physics.optics)
import_from(sympy.physics.quantum.spin)

@testset "Physics" begin
    @test clebsch_gordan(Sym(3)/2, Sym(1)/2, 2, Sym(3)/2, Sym(1)/2, 2) == 1
    @test clebsch_gordan(Sym(3)/2, Sym(1)/2, 1, Sym(3)/2, -Sym(1)/2, 1) == sqrt(Sym(3))/2
    @test clebsch_gordan(Sym(3)/2, Sym(1)/2, 1, -Sym(1)/2, Sym(1)/2, 0) == -sqrt(Sym(2))/2

    θ, ϕ = symbols("theta, phi")
    @test dot_rot_grad_Ynm(Sym(3), 2, 2, 0, θ, ϕ).doit() == 3*sqrt(Sym(55))*Ynm(Sym(5), 2, θ, ϕ)/(11*sqrt(PI))

    @test gaunt(Sym(1),0,1,1,0,-1) == -1/(2*sqrt(PI))
    @test N(gaunt(Sym(1000),1000,1200,9,3,-12)) ≈ 0.00689500421922113448 # takes forever to compute

    @test_throws PyError gaunt(Sym(1),2,0,1.2,0,0,0)
    @test_throws PyError gaunt(Sym(1),0,1,1.1,0,-1.1)

    @test racah(Sym(3),3,3,3,3,3) == -1//14

    @test wigner_3j(Sym(2), 6, 4, 0, 0, 0) == sqrt(Sym(715))/143
    @test wigner_3j(Sym(2), 6, 4, 0, 0, 1) == 0

    @test wigner_6j(Sym(3),3,3,3,3,3) == -1//14
    @test wigner_6j(Sym(5),5,5,5,5,5) == 1//52

    # Optics
    # sympy.physics.optics.RayTransferMatrix not a functino
    @test sympy.physics.optics.RayTransferMatrix(1,2,3,4) == convert(SymMatrix,Sym[1 2; 3 4])
    @test sympy.physics.optics.FlatMirror() == convert(SymMatrix,Sym[1 0; 0 1])

    # Spin
    import_from(sympy.physics.quantum.spin, (:WignerD,), typ=:Any)
    @test WignerD(Sym(3), Sym(2), Sym(1), PI, PI/2, -PI).doit() == -sqrt(Sym(10))/8
    @test WignerD(Sym(1), Sym(1), Sym(0), 0, θ, 0).doit() == -sin(θ)/sqrt(Sym(2))
end


text:
displayinline(xs::List) =
  Collapsible(HTML("List"),
              HTML() do io
                println(io, """
""")
                for x in take(MAX_CELLS, xs)
                  print(io, "
")
                  x = displayinline(x)
                  writemime(io, bestmime(x), x)
                  print(io, "
")
                end
                !isempty(drop(MAX_CELLS, xs)) && println(io, """
⋮
""")
                println(io, """
""")
              end)


text:
# Script used to calculate log-likelihood of multivariate skew t distribution
# for comparison with Julia code. Note that the R script of the same name must be run
# before this, as this script uses its output.

using MultivariateStats, SkewDist

# Load R data
X = readdlm("x.txt", ' ')
Y = readdlm("y.txt", ' ')

n, k = size(Y)
p = size(X,2)

βinit = readdlm("Beta.txt", ' ')
Ωinit = readdlm("Omega.txt", ' ')
αinit = ones(k)
νinit = 4.0

ρ, A, η, logν = SkewDist.dplist2optpar(Ωinit, αinit, νinit)
params = SkewDist.write_params(βinit, ρ, A, η, logν)

println("Parameter vector:", params)
nll, grad = SkewDist.nll_and_grad(params, X, Y)
println("Neg. Log Likelihood: ", nll)
print("Gradient: ", grad, "\n")

βinit = llsq(X,Y; bias=false)
resid = Y - X*βinit


text:
module MedicalModels

using Flux
using Flux: Parallel

include("./unet3D.jl")
include("./utils.jl")

export
    # Export unet3D.jl functions
    unet3D,

    # Export utils.jl functions
    conv,
    tran,
    conv1,
    conv2,
    tran2

end


text:
export Billiard, randominside, randominside_xyφ
#######################################################################################
## Billiard Table
#######################################################################################
struct Billiard{T, D, O<:Tuple, Y}
    obstacles::O
    peridx::Y
end

@inline isperiodic(::Billiard{T,D,<:Tuple,Vector{Int}}) where {T,D} = true
@inline isperiodic(::Billiard{T,D,<:Tuple,Nothing}) where {T,D} = false
@inline isperiodic(i::Int, ::Billiard{T,D,O,Nothing}) where {T,D,O} = false
@inline isperiodic(i::Int, bd::Billiard{T,D,<:Tuple,Vector{Int}}) where {T,D} = i ∈ bd.peridx

#pretty print:
_get_name(o::Obstacle) = :name ∈ fieldnames(typeof(o)) ? o.name : string(typeof(o))
function Base.show(io::IO, bd::Billiard{T,D,BT}) where {T, D, BT}
    s = "Billiard{$T} with $D obstacles:"
    bd.peridx != nothing && (s = "Periodic "*s)
    for i in 1:min(10, length(bd))
        s*="\n  $(_get_name(bd[i]))"
    end
    print(io, s)
    length(bd) > 10 && print(io, "\n  ...")
end

"""
    Billiard(obstacles...)
Construct a `Billiard` from given `obstacles` (tuple, vector, varargs).

For functions like [`boundarymap`](@ref),
it is expected (if possible) that the obstacles of the billiard are sorted,
such that the arc-coordinate `ξ` around the billiard is increasing counter-clockwise.

`ξ` is measured as:
* the distance from start point to end point in `Wall`s
* the arc length measured counterclockwise from the open face in `Semicircle`s
* the arc length measured counterclockwise from the rightmost point
  in `Circular`s
"""
function Billiard(bd::Union{AbstractVector, Tuple})

    T = eltype(bd[1])
    D = length(bd)
    # Assert that all elements of `bd` are of same type:
    for i in 2:D
        eltype(bd[i]) != T && throw(ArgumentError(
        "All obstacles of the billiard must have same type of
        numbers. Found $T and $(eltype(bd[i])) instead."
        ))
    end
    tup = (bd...,)
    peridx = findall(x -> typeof(x) <: PeriodicWall, tup)
    if isodd(length(peridx))
        throw(ArgumentError(
        "A billiard can only have an even number of `PeriodicWall`s, "*
        "since they have to come in pairs."
        ))
    end
    if length(peridx) == 0
        return Billiard{T, D, typeof(tup), Nothing}(tup, nothing)
    else
        return Billiard{T, D, typeof(tup), Vector{Int}}(tup, peridx)
    end
end

function Billiard(bd::Vararg{Obstacle})
    T = eltype(bd[1])
    tup = (bd...,)
    return Billiard(tup)
end

Base.getindex(bd::Billiard, i) = bd.obstacles[i]
Base.lastindex(bd::Billiard) = length(bd.obstacles)
Base.firstindex(bd::Billiard) = 1
# Iteration:
Base.iterate(bd::Billiard) = iterate(bd.obstacles)
Base.iterate(bd::Billiard, state) = iterate(bd.obstacles, state)
Base.length(::Billiard{T, L}) where {T, L} = L
Base.eltype(::Billiard{T}) where {T} = T

#######################################################################################
## Distances
#######################################################################################
for f in (:distance, :distance_init)
    @eval $(f)(p::AbstractParticle, bd::Billiard) = $(f)(p.pos, bd.obstacles)
    @eval $(f)(pos::SV{T}, bd::Billiard) where {T} = $(f)(pos, bd.obstacles)
end

for f in (:distance, :distance_init)
    @eval begin
        function ($f)(p::SV{T}, bd::Tuple)::T where {T}
            dmin::T = T(Inf)
            for obst in bd
                d::T = distance(p, obst)
                d < dmin && (dmin = d)
            end
            return dmin
        end
    end
end

#######################################################################################
## randominside
#######################################################################################
function cellsize(
    bd::Union{Vector{<:Obstacle{T}}, Billiard{T}}) where {T<:AbstractFloat}

    xmin::T = ymin::T = T(Inf)
    xmax::T = ymax::T = T(-Inf)
    for obst ∈ bd
        xs::T, ys::T, xm::T, ym::T = cellsize(obst)
        xmin = xmin > xs ? xs : xmin
        ymin = ymin > ys ? ys : ymin
        xmax = xmax < xm ? xm : xmax
        ymax = ymax < ym ? ym : ymax
    end
    return xmin, ymin, xmax, ymax
end

"""
    randominside(bd::Billiard [, ω]) → p
Return a particle `p` with random allowed initial conditions inside the given
billiard. If supplied with a second argument the
type of the returned particle is `MagneticParticle`, with angular velocity `ω`.

**WARNING** : `randominside` works for any **convex** billiard but it does
not work for non-convex billiards. (this is because it uses `distance`)
"""
randominside(bd::Billiard, ω::Nothing = nothing) = Particle(_randominside(bd)...)
randominside(bd::Billiard{T}, ω::Real) where {T} =
MagneticParticle(_randominside(bd)..., T(ω))

"""
    randominside_xyφ(bd::Billiard) → x, y, φ
Same as [`randominside`](@ref) but only returns position and direction.
"""
function randominside_xyφ(bd::Billiard{T}) where {T<:AbstractFloat}
    xmin::T, ymin::T, xmax::T, ymax::T = cellsize(bd)
    xp = T(rand())*(xmax-xmin) + xmin
    yp = T(rand())*(ymax-ymin) + ymin
    pos = SV{T}(xp, yp)
    dist = distance_init(pos, bd)
    while dist <= sqrt(eps(T))
        xp = T(rand())*(xmax-xmin) + xmin
        yp = T(rand())*(ymax-ymin) + ymin
        pos = SV{T}(xp, yp)
        dist = distance_init(pos, bd)
    end
    return pos[1], pos[2], T(2π*rand())
end

function _randominside(bd::Billiard{T}) where {T<:AbstractFloat}
    x, y, φ = randominside_xyφ(bd)
    pos = SV{T}(x, y)
    vel = SV{T}(sincos(φ)...)
    return pos, vel, zero(SV{T}) # this zero is the `current_cell`
end


text:
module FlightGNC
using LinearAlgebra, Statistics

using Reexport
@reexport using FSimBase
const FSBase = FSimBase
import FSimBase: State, Params, Dynamics!, Dynamics, Command, apply_inputs, Simulator, solve, reinit!, log!, step!, step_until!

using UnPack, ComponentArrays
using Plots                             # required for fig_print
using Convex, Mosek, MosekTools         # required for IASCG
using Roots, OrdinaryDiffEq, Optim      # required for Speed_Predictor in IASCG
using DiffEqFlux                        # required for CTPG_train

# exports -------------------------------------------------------------
## algorithms
export BPNG, BPNG_cmd, Bias_IACG_StationaryTarget, Bias_IACG_StationaryTarget_2D, Bias_Quaternion_IACG, Bias_zero     # defined in BPNG.jl
export IASCG, IASCG_Planner, IASCG_cmd, Path_Generator, Speed_Predictor, Speed_Dynamics         # defined in IASCG.jl
export CTPG_train                                                                               # defined in CTPG.jl

## environments
export PointMassMissile, PursuerEvadorMissile
export VerticalPlaneDynamicMissile
export VerticalPlane_RigidBody_Missile

## utils
export fig_print, equal_AR_3D, view_result

# includes -------------------------------------------------------------
## algorithms
include("algorithms/BPNG.jl")
include("algorithms/IASCG.jl")
include("algorithms/CTPG.jl")

## environments
include("environments/missiles.jl")

## utils
include("utils/fig_utils.jl")
end


text:
using Documenter
using JSExpr

makedocs(
    sitename="JSExpr.jl",
    authors="Travis DePrato and contributors",
    pages=[
        "index.md",
        "basics.md",
        "datastructures.md",
        "juliaisms.md",
        "api.md",
    ],
)

deploydocs(
    repo="github.com/JuliaGizmos/JSExpr.jl.git",
)


text:
@testset "SLIC" begin
  setify(lists) = Set(Set.(lists))

  𝒮 = georef((z=[ones(10,10) 2ones(10,10); 3ones(10,10) 4ones(10,10)],))
  p = partition(𝒮, SLIC(4, 1.0))
  @test length(p) == 4
  @test all(nelements.(p) .== 100)
  p1, p2, p3, p4 = p
  @test mean(coordinates(centroid(p1, ind)) for ind in 1:nelements(p1)) == [5.0,5.0]
  @test mean(coordinates(centroid(p2, ind)) for ind in 1:nelements(p2)) == [15.0,5.0]
  @test mean(coordinates(centroid(p3, ind)) for ind in 1:nelements(p3)) == [5.0,15.0]
  @test mean(coordinates(centroid(p4, ind)) for ind in 1:nelements(p4)) == [15.0,15.0]

  𝒮 = georef((z=[√(i^2+j^2) for i in 1:100, j in 1:100],))
  p = partition(𝒮, SLIC(50, 1.0))
  @test length(p) == 49

  if visualtests
    @test_reference "data/slic.png" plot(p)
  end
end


text:
# test file to test testing
# run this with julia -L extras/test.jl -e "tests(\"test/test_test.jl\")"
require("nearequal")
require("test")

test_context("Testing test tests")
# setup goes here

test_group("string tests")
@test strip("\t  hi   \n") == "hi"
@testfails strip("\t  this should fail   \n") == "hi"

test_group("numeric tests")
@test isapprox(.1+.1+.1, .3)
@testfails isapprox(.1+.1+.1, .4)

test_group("array tests")
a = Array(Float64, 2, 2, 2, 2, 2)
a[1,1,1,1,1] = 10
@test a[1,1,1,1,1] == 10
@testfails a[1,1,1,1,1] == 2

test_group("random tests")
@test rand() != rand() # not very likely to fail
@testfails rand() == rand() # very likely to fail

test_group("exception tests")
@testfails complex(1,2) > 0 # fail
@test throws_exception(complex(1,2) > 0, MethodError)
@testfails throws_exception(complex(1,2) > 0, SystemError)
@testfails throws_exception(2 > 1, DomainError) # would correctly fail

test_group("printing tests")
@test sprint(show, :(1+2)) == ":(+(1,2))"
@test prints(print_joined, ([1,2,3], " : "), "1 : 2 : 3") # prints is a helper
@testfails prints(print_joined, ([1,2,3], " ! "), "1 : 2 : 3")

test_group("performance tests")
fib(n) = n < 2 ? n : fib(n-1) + fib(n-2)
@test fib(20) == 6765
@testfails takes_less_than(fib(20), 1e-6)

# shutdown goes here

max_stars_repo_name: porterjamesj/julia

max_issues_repo_name: porterjamesj/julia

max_forks_repo_name: porterjamesj/julia


text:
module ActionQueue
include("action.jl")

@enum ActionStatus begin
    queued
    approved
    pending
    success
    failed
    canceled
end

@enum ActionType begin
    allocate
    unallocate
    reallocate
    collect
end

abstract type ActionInput end

struct AllocateActionInput <: ActionInput
    status::ActionStatus
    type::ActionType
    deploymentID::AbstractString
    amount::AbstractString
    source::AbstractString
    reason::AbstractString
end

struct UnallocateActionInput <: ActionInput
    status::ActionStatus
    type::ActionType
    allocationID::AbstractString
    source::AbstractString
    reason::AbstractString
end

struct ReallocateActionInput <: ActionInput
    status::ActionStatus
    type::ActionType
    allocationID::AbstractString
    amount::AbstractString
    source::AbstractString
    reason::AbstractString
end

structtodict(x::ActionInput) = Dict(string(k) => getfield(x, k) for k in propertynames(x))

function reallocate_actions(
    proposed_ipfs::Vector{T},
    existing_ipfs::Vector{T},
    proposed_allocations::Dict{T,<:Real},
    existing_allocations::Dict{T,T},
) where {T<:AbstractString}
    ipfses = reallocate_ipfs(existing_ipfs, proposed_ipfs)
    actions = map(
        ipfs -> structtodict(
            ReallocateActionInput(
                queued,
                reallocate,
                existing_allocations[ipfs],
                string(proposed_allocations[ipfs]),
                "AllocationOpt",
                "AllocationOpt",
            ),
        ),
        ipfses,
    )
    return actions, ipfses
end

function allocate_actions(
    proposed_ipfs::Vector{T},
    reallocate_ipfs::Vector{T},
    proposed_allocations::Dict{T,<:Real},
) where {T<:AbstractString}
    ipfses = open_ipfs(proposed_ipfs, reallocate_ipfs)
    actions = map(
        ipfs -> structtodict(
            AllocateActionInput(
                queued,
                allocate,
                ipfs,
                string(proposed_allocations[ipfs]),
                "AllocationOpt",
                "AllocationOpt",
            ),
        ),
        ipfses,
    )
    return actions, ipfses
end

function unallocate_actions(
    existing_allocations::Dict{T,T},
    existing_ipfs::Vector{T},
    reallocate_ipfs::Vector{T},
    frozenlist::Vector{T},
) where {T<:AbstractString}
    ipfses = close_ipfs(existing_ipfs, reallocate_ipfs, frozenlist)
    actions = map(
        ipfs -> structtodict(
            UnallocateActionInput(
                queued,
                unallocate,
                existing_allocations[ipfs],
                "AllocationOpt",
                "AllocationOpt",
            ),
        ),
        ipfses,
    )
    return actions, ipfses
end

end


text:
#!/usr/bin/env julia

using DocOpt
using DataFrames
using RCall
using Glob
import BioBenchmarks

args = docopt("Usage: runbenchmarks.jl <revision>...")
revisions = args["<revision>"]

results = DataFrame()
for revision in revisions
    filepath = BioBenchmarks.run(revision)
    r = readtable(filepath)
    r[:revision] = filepath[9:end-4]
    results = vcat(results, r)
end

# 64-bit integers may overflow in R
for f in [:elapsed, :gc_time, :allocated]
    results[f] = map(Float64, results[f])
end

serial = maximum(vcat(0, [parse(Int, splitext(fn)[1][12:end]) for fn in glob("benchmarks_*.png")]))
plotfile = string("benchmarks_", serial + 1, ".png")
@rput plotfile

R"""
library(ggplot2)

ggplot($(results), aes(x=revision, y=elapsed, color=revision)) +
  geom_boxplot() +
  geom_point(position=position_jitterdodge()) +
  facet_wrap(~name, scales="free") +
  theme(axis.text.x=element_blank()) +
  expand_limits(y=0)

ggsave(plotfile)
"""


text:
# Input constraint is Hyperrectangle (uniform radius)
# Output constriant is HPolytope with one constraint
"""
    FastLin(maxIter::Int64, ϵ0::Float64, accuracy::Float64)

FastLin combines reachability analysis with binary search to find maximum allowable disturbance.

# Problem requirement
1. Network: any depth, ReLU activation
2. Input: hypercube
3. Output: halfspace

# Return
`AdversarialResult`

# Method
Reachability analysis by network approximation and binary search.
- `max_iter` is the maximum iteration in search, default `10`;
- `ϵ0` is the initial search radius, default `100.0`;
- `accuracy` is the stopping criteria, default `0.1`;

# Property
Sound but not complete.

# Reference
[T.-W. Weng, H. Zhang, H. Chen, Z. Song, C.-J. Hsieh, D. Boning, I. S. Dhillon, and L. Daniel,
"Towards Fast Computation of Certified Robustness for ReLU Networks,"
*ArXiv Preprint ArXiv:1804.09699*, 2018.](https://arxiv.org/abs/1804.09699)
"""
@with_kw struct FastLin
    maxIter::Int64    = 10
    ϵ0::Float64       = 100.0
    accuracy::Float64 = 0.1
end

function solve(solver::FastLin, problem::Problem)
    ϵ = fill(solver.ϵ0, solver.maxIter+1)
    ϵ_upper = 2 * max(solver.ϵ0, maximum(problem.input.radius))
    ϵ_lower = 0.0
    n_input = dim(problem.input)
    for i = 1:solver.maxIter
        input_bounds = Hyperrectangle(problem.input.center, fill(ϵ[i], n_input))
        output_bounds = forward_network(solver, problem.network, input_bounds)
        if issubset(output_bounds, problem.output)
            ϵ_lower = ϵ[i]
            ϵ[i+1] = (ϵ[i] + ϵ_upper) / 2
            abs(ϵ[i] - ϵ[i+1]) > solver.accuracy || break
        else
            ϵ_upper = ϵ[i]
            ϵ[i+1] = (ϵ[i] + ϵ_lower) / 2
        end
    end
    if ϵ_lower > minimum(problem.input.radius)
        return AdversarialResult(:SAT, ϵ_lower) # previously :True
    else
        return AdversarialResult(:UNSAT, ϵ_lower)
    end
end

function forward_network(solver::FastLin, nnet::Network, input::Hyperrectangle)
    # use MaxSens to compute bounds:
    L, U = get_bounds(nnet, input.center, input.radius[1])
    output = Hyperrectangle(low = last(L), high = last(U))
    return output
end

# function forward_network(ϵ::Float64, nnet::Network, input::Hyperrectangle)

#     len = length(nnet.layers)

#     act_pattern = Vector{Float64}()
#     A = Vector{Matrix{Float64}}()
#     l = Vector{Vector{Float64}}()
#     u = Vector{Vector{Float64}}()

#     x0 = input.center
#     println("in FN")
#     for m in 1:len
#         println("looping $m")
#         new_l, new_u = get_next_bounds(nnet.layers, x0, ϵ, A, l, act_pattern)

#         update_A!(A, m, new_l, new_u)

#         push!(l, new_l)
#         push!(u, new_u)
#     end

#     return Hyperrectangle(low = last(l), high = last(u))
# end

# #=
# input: input data vector,
# p : ℓp norm,
# ϵ: maximum ℓp-norm perturbation
# l[k], u[k], k∈[m] : layer-wise bounds
# =#
# function get_next_bounds(layers, input_center, ϵ, A, l, act_pattern)
#     W0 = layers[1].weights
#     m = length(A) + 1   # A is of length 0 on the first pass, and is not necessary then.

#     T, H = get_TH(l, A, act_pattern)

#     ## TODO can the order of the for loops be reversed?
#     for j in 1:n_nodes(layers[m])  # for each node `j` in the mth layer:
#         vj = W0[j, :]*input_center + layers[m].bias[j]
#         μ_upper = μ_lower = 0
#         for k in 1:m-1
#             μ_upper = μ_upper - A[k][j, :] * T[k][:, j]     # Eq 6
#             μ_lower = μ_lower - A[k][j, :] * H[k][:, j]
#             vj += A[k][j, :] * layers[m].bias               # Eq 7
#         end

#         γᵁ = vj + ϵ*sum(abs, W0[j, :]) + μ_upper     #||A[0][j, :]||[q] => qnorm of A
#         γᴸ = vj + ϵ*sum(abs, W0[j, :]) + μ_lower
#     end

#     push!(act_pattern, relaxed_ReLU.(γᴸ, γᵁ))

#     return γᴸ, γᵁ
# end

# function get_TH(l, A, act_pattern)
#     T = Vector{SparseMatrixCSC{Float64}}(undef, length(l))
#     H = Vector{SparseMatrixCSC{Float64}}(undef, length(l))

#     for k in 1:length(l)  # for all previous layers, (current one not included)

#         nᵏ = length(act_pattern[k]) # number nodes in the kth layer
#         lᵏ = l[k]

#         T[k] = spzeros(length(lᵏ), nᵏ) # reset/initialize each one
#         H[k] = spzeros(length(lᵏ), nᵏ) # reset/initialize each one

#         for r in 1:length(lᵏ) # NOTE: what does the length of the lower bound represent? The dimensionality of the input space? In the paper they use Iᵏ
#             if 0.0 < act_pattern[k][r] < 1.0
#                 # A⁺ = findall(...)  ## for Julia 0.7+
#                 A⁺ᵣ = find(a->a>0, A[k][:, r])            # positive indices
#                 A⁻ᵣ = setdiff(eachindex(A[k][:, r]), A⁺)  # negative indices
#                 T[k][r, A⁺] = lᵏ[A⁺ᵣ]
#                 H[k][r, A⁻] = lᵏ[A⁻ᵣ]
#             end
#         end
#     end
#     return T, H
# end

# # function update_A!(A, m, l, u)

# #     for k in reverse(1:m-1)
# #         if k == m-1
# #             D = get_D(l[m-1], u[m-1])  # diagonal matrix according to Eq. 5
# #             A[m-1] = layers[m].weights * D
# #         else
# #             A[k] = A[m-1] * A[k]
# #         end
# #     end
# # end
# # Identical to:
# # function update_A!(A, m, l, u)
# #     D = get_D(l, u)  # diagonal matrix according to Eq. 5
# #     A[m-1] = layers[m].weights * D
# #     for k in 1:m-2
# #         A[k] = A[m-1] * A[k]
# #     end
# # end
# function update_A!(A, l, u)
#     D = Diagonal(relaxed_ReLU.(l, u))  # diagonal matrix according to Eq. 5
#     WD = layers[m].weights * D

#     map!(a->WD*a, A, A)
#     push!(A, WD)  # consider pushing I and mapping WD onto it along with everything else
# end