text:
"""
BlockDim

An index for a BlockSparseTensor.
"""
const BlockDim = Vector{Int}

# Makes for generic code
dim(d::BlockDim) = sum(d)

"""
BlockDims{N}

Dimensions used for BlockSparse NDTensors.
Each entry lists the block sizes in each dimension.
"""
const BlockDims{N} = NTuple{N,BlockDim}

Base.ndims(ds::Type{<:BlockDims{N}}) where {N} = N

similartype(::Type{<:BlockDims}, ::Type{Val{N}}) where {N} = BlockDims{N}

Base.copy(ds::BlockDims) = ds

"""
dim(::BlockDims,::Integer)

Return the total extent of the specified dimensions.
"""
function dim(ds::BlockDims{N}, i::Integer) where {N}
  return sum(ds[i])
end

"""
dims(::BlockDims)

Return the total extents of the dense space
the block dimensions live in.
"""
function dims(ds::BlockDims{N}) where {N}
  return ntuple(i -> dim(ds, i), Val(N))
end

"""
dim(::BlockDims)

Return the total extent of the dense space
the block dimensions live in.
"""
function dim(ds::BlockDims{N}) where {N}
  return prod(dims(ds))
end

"""
nblocks(::BlockDim)

The number of blocks of the BlockDim.
"""
function nblocks(ind::BlockDim)
  return length(ind)
end

"""
nblocks(::BlockDims,i::Integer)

The number of blocks in the specified dimension.
"""
function nblocks(inds::Tuple, i::Integer)
  return nblocks(inds[i])
end

"""
nblocks(::BlockDims,is)

The number of blocks in the specified dimensions.
"""
function nblocks(inds::Tuple, is::NTuple{N,Int}) where {N}
  return ntuple(i -> nblocks(inds, is[i]), Val(N))
end

"""
nblocks(::BlockDims)

A tuple of the number of blocks in each
dimension.
"""
function nblocks(inds::NTuple{N,<:Any}) where {N}
  return ntuple(i -> nblocks(inds, i), Val(N))
end

"""
blockdim(::BlockDim,::Integer)

The size of the specified block in the specified
dimension.
"""
function blockdim(ind::BlockDim, i::Integer)
  return ind[i]
end

"""
    blockdim(::BlockDims,block,::Integer)

The size of the specified block in the specified
dimension.
"""
function blockdim(inds, block, i::Integer)
  return blockdim(inds[i], block[i])
end

"""
    blockdims(::BlockDims,block)

The size of the specified block.
"""
function blockdims(inds, block)
  return ntuple(i -> blockdim(inds, block, i), ValLength(inds))
end

"""
    blockdim(::BlockDims,block)

The total size of the specified block.
"""
function blockdim(inds, block)
  return prod(blockdims(inds, block))
end

"""
blockdiaglength(inds::BlockDims,block)

The length of the diagonal of the specified block.
"""
function blockdiaglength(inds, block)
  return minimum(blockdims(inds, block))
end

function outer(dim1, dim2, dim3, dims...; kwargs...)
  return outer(outer(dim1, dim2), dim3, dims...; kwargs...)
end

function outer(dim1::BlockDim, dim2::BlockDim)
  dimR = BlockDim(undef, nblocks(dim1) * nblocks(dim2))
  for (i, t) in enumerate(Iterators.product(dim1, dim2))
    dimR[i] = prod(t)
  end
  return dimR
end

function permuteblocks(dim::BlockDim, perm)
  return dim[perm]
end

# Given a CartesianIndex in the range dims(T), get the block it is in
# and the index within that block
function blockindex(T, i::Vararg{Int,N}) where {ElT,N}
  # Start in the (1,1,...,1) block
  current_block_loc = @MVector ones(Int, N)
  current_block_dims = blockdims(T, Tuple(current_block_loc))
  block_index = MVector(i)
  for dim in 1:N
    while block_index[dim] > current_block_dims[dim]
      block_index[dim] -= current_block_dims[dim]
      current_block_loc[dim] += 1
      current_block_dims = blockdims(T, Tuple(current_block_loc))
    end
  end
  return Tuple(block_index), Block{N}(current_block_loc)
end

blockindex(T) = (), Block{0}()

#
# This is to help with ITensor compatibility
#

setblockdim!(dim1::BlockDim, newdim::Int, n::Int) = setindex!(dim1, newdim, n)

sim(dim::BlockDim) = copy(dim)

dir(::BlockDim) = 0

dag(dim::BlockDim) = copy(dim)


text:
"""
Validation PSSE/GENROE:
This case study defines a three bus system with an infinite bus, GENROE and a load.
The fault drop the line connecting the infinite bus and GENROE.
"""

##################################################
############### SOLVE PROBLEM ####################
##################################################

# Define dyr files

names = ["GENROE: Normal Saturation", "GENROE: High Saturation"]

dyr_files = [
    joinpath(TEST_FILES_DIR, "benchmarks/psse/GENROE/ThreeBus_GENROE.dyr"),
    joinpath(TEST_FILES_DIR, "benchmarks/psse/GENROE/ThreeBus_GENROE_HIGH_SAT.dyr"),
]

csv_files = (
    joinpath(TEST_FILES_DIR, "benchmarks/psse/GENROE/TEST_GENROE.csv"),
    joinpath(TEST_FILES_DIR, "benchmarks/psse/GENROE/TEST_GENROE_HIGH_SAT.csv"),
)

init_conditions = [test_psse_genroe_init, test_psse_genroe_high_sat_init]

eigs_values = [test16_eigvals, test16_eigvals_high_sat]

raw_file_dir = joinpath(TEST_FILES_DIR, "benchmarks/psse/GENROE/ThreeBusMulti.raw")
tspan = (0.0, 20.0)

function test_genroe_implicit(dyr_file, csv_file, init_cond, eigs_value)
    path = (joinpath(pwd(), "test-psse-genrou"))
    !isdir(path) && mkdir(path)
    try
        sys = System(raw_file_dir, dyr_file)

        # Define Simulation Problem
        sim = Simulation!(
            ResidualModel,
            sys, #system
            path,
            tspan, #time span
            BranchTrip(1.0, Line, "BUS 1-BUS 2-i_1"), #Type of Fault
        ) #Type of Fault

        # Test Initial Condition
        diff = [0.0]
        res = get_init_values_for_comparison(sim)
        for (k, v) in init_cond
            diff[1] += LinearAlgebra.norm(res[k] - v)
        end

        @test (diff[1] < 1e-3)

        # Obtain small signal results for initial conditions
        small_sig = small_signal_analysis(sim)
        eigs = small_sig.eigenvalues
        @test small_sig.stable

        #Test Eigenvalues
        @test LinearAlgebra.norm(eigs - eigs_value) < 1e-3

        # Solve problem
        @test execute!(sim, IDA(), dtmax = 0.005, saveat = 0.005) ==
              PSID.SIMULATION_FINALIZED
        results = read_results(sim)

        # Obtain data for angles
        series = get_state_series(results, ("generator-102-1", :δ))
        t = series[1]
        δ = series[2]
        series2 = get_voltage_magnitude_series(results, 102)

        t_psse, δ_psse = get_csv_delta(csv_file)

        # Test Transient Simulation Results
        # PSSE results are in Degrees
        @test LinearAlgebra.norm(δ - (δ_psse .* pi / 180), Inf) <= 1e-1
        @test LinearAlgebra.norm(t - round.(t_psse, digits = 3)) == 0.0

        power = PSID.get_activepower_series(results, "generator-102-1")
        rpower = PSID.get_reactivepower_series(results, "generator-102-1")
        @test isa(series2, Tuple{Vector{Float64}, Vector{Float64}})
        @test isa(power, Tuple{Vector{Float64}, Vector{Float64}})
        @test isa(rpower, Tuple{Vector{Float64}, Vector{Float64}})
    finally
        @info("removing test files")
        rm(path, force = true, recursive = true)
    end
end

function test_genroe_mass_matrix(dyr_file, csv_file, init_cond, eigs_value)
    path = (joinpath(pwd(), "test-psse-genrou"))
    !isdir(path) && mkdir(path)
    try
        sys = System(raw_file_dir, dyr_file)

        # Define Simulation Problem
        sim = Simulation!(
            MassMatrixModel,
            sys, #system
            path,
            tspan, #time span
            BranchTrip(1.0, Line, "BUS 1-BUS 2-i_1"), #Type of Fault
        ) #Type of Fault

        # Test Initial Condition
        diff = [0.0]
        res = get_init_values_for_comparison(sim)
        for (k, v) in init_cond
            diff[1] += LinearAlgebra.norm(res[k] - v)
        end

        @test (diff[1] < 1e-3)

        # Obtain small signal results for initial conditions
        small_sig = small_signal_analysis(sim)
        eigs = small_sig.eigenvalues
        @test small_sig.stable

        #Test Eigenvalues
        @test LinearAlgebra.norm(eigs - eigs_value) < 1e-3

        # Solve problem
        @test execute!(sim, Rodas4(), dtmax = 0.005, saveat = 0.005) ==
              PSID.SIMULATION_FINALIZED
        results = read_results(sim)

        #Obtain data for angles
        series = get_state_series(results, ("generator-102-1", :δ))
        t = series[1]
        δ = series[2]
        series2 = get_voltage_magnitude_series(results, 102)

        t_psse, δ_psse = get_csv_delta(csv_file)

        # Test Transient Simulation Results
        # PSSE results are in Degrees
        @test LinearAlgebra.norm(δ - (δ_psse .* pi / 180), Inf) <= 1e-1
        @test LinearAlgebra.norm(t - round.(t_psse, digits = 3)) == 0.0

        power = PSID.get_activepower_series(results, "generator-102-1")
        rpower = PSID.get_reactivepower_series(results, "generator-102-1")
        @test isa(series2, Tuple{Vector{Float64}, Vector{Float64}})
        @test isa(power, Tuple{Vector{Float64}, Vector{Float64}})
        @test isa(rpower, Tuple{Vector{Float64}, Vector{Float64}})
    finally
        @info("removing test files")
        rm(path, force = true, recursive = true)
    end
end

@testset "Test 16 GENROE ResidualModel" begin
    for (ix, name) in enumerate(names)
        @testset "$(name)" begin
            dyr_file = dyr_files[ix]
            csv_file = csv_files[ix]
            init_cond = init_conditions[ix]
            eigs_value = eigs_values[ix]
            test_genroe_implicit(dyr_file, csv_file, init_cond, eigs_value)
        end
    end
end

@testset "Test 16 GENROE MassMatrixModel" begin
    for (ix, name) in enumerate(names)
        @testset "$(name)" begin
            dyr_file = dyr_files[ix]
            csv_file = csv_files[ix]
            init_cond = init_conditions[ix]
            eigs_value = eigs_values[ix]
            test_genroe_mass_matrix(dyr_file, csv_file, init_cond, eigs_value)
        end
    end
end

ext: jl


text:
title = "LyceumMuJoCoViz"
weight = 30


text:
# Concrete implementation ######################################################

"""
    NMRData(data, dims, refdims, name)

The main structure for holding MNR data. A subtype of `AbstractDimensionalArray`, which
maintains and updates its dimensions through transformations and moves dimensions to
`refdims` after reducing operations (like e.g. `mean`).
"""
struct NMRData{T,N,D<:Tuple,R<:Tuple,A<:AbstractArray{T,N},Na<:AbstractString,Me} <: AbstractDimensionalArray{T,N,D,A}
    data::A
    dims::D
    refdims::R
    name::Na
    metadata::Me
    function NMRData(data::A, dims::D, refdims::R, name::Na, metadata::Me
                             ) where {D,R,A<:AbstractArray{T,N},Na,Me} where {T,N}
        map(dims, size(data)) do d, s
            if !(val(d) isa Colon) && length(d) != s
                throw(DimensionMismatch(
                    "dims must have same size as data. This was not true for $dims and size $(size(data)) $(A)."
                ))
            end
        end
        new{T,N,D,R,A,Na,Me}(data, dims, refdims, name, metadata)
    end
end
"""
    NMRData(data, dims::Tuple [, name::String]; refdims=(), metadata=nothing)

Constructor with optional `name`, and keyword `refdims` and `metadata`.

Example:
```julia
using Dates, DimensionalData

ti = (Ti(DateTime(2001):Month(1):DateTime(2001,12)),
x = X(10:10:100))
A = NMRData(rand(12,10), (ti, x), "example")

julia> A[X(Near([12, 35])), Ti(At(DateTime(2001,5)))];

julia> A[Near(DateTime(2001, 5, 4)), Between(20, 50)];
```
"""
NMRData(A::AbstractArray, dims, name::String=""; refdims=(), metadata=Dict()) =
    NMRData(A, DimensionalData.formatdims(A, DimensionalData._to_tuple(dims)), refdims, name, metadata)
#_to_tuple(t::T where T <: Tuple) = t
#_to_tuple(t) = tuple(t)

# Getters
refdims(A::NMRData) = A.refdims
data(A::NMRData) = A.data
name(A::NMRData) = A.name

# get axis values
xval(A::NMRData) = dims(A,X).val
yval(A::NMRData) = dims(A,Y).val
zval(A::NMRData) = dims(A,Z).val
tval(A::NMRData) = dims(A,Ti).val

# metadata accessor functions
metadata(A::NMRData) = A.metadata
metadata(A::NMRData, key::Symbol) = get(metadata(A), key, missing)
metadata(A::NMRData, dim, key::Symbol) = get(metadata(A, dim), key, missing)
getindex(A::NMRData, key::Symbol) = get(metadata(A), key, missing)
getindex(A::NMRData, dim, key::Symbol) = get(metadata(A, dim), key, missing)
setindex!(A::NMRData, v, key::Symbol) = setindex!(metadata(A), v, key)  #(metadata(A)[key] = v)
setindex!(A::NMRData, v, dim, key::Symbol) = setindex!(metadata(A,dim), v, key)  #(metadata(A, dim)[key] = v)

# acqus accessor functions
acqus(A::NMRData) = get(A.metadata, :acqus, missing)
acqus(A::NMRData, key::String) = ismissing(acqus(A)) ? missing : get(acqus(A), uppercase(key), missing)
acqus(A::NMRData, key::Symbol) = acqus(A, string(key))
acqus(A::NMRData, key, index::Int) = acqus(A, key)[index]

# label accessor functions
label(A::NMRData) = get(metadata(A),:label,"")
label(A::NMRData, dim) = get(metadata(A, dim),:label,"")
label!(A::NMRData, labeltext::AbstractString) = (metadata(A)[:label] = labeltext)
label!(A::NMRData, dim, labeltext::AbstractString) = (metadata(A, dim)[:label] = labeltext)

# AbstractDimensionalArray interface
@inline rebuild(A::NMRData, data::AbstractArray, dims::Tuple,
                refdims::Tuple, name::AbstractString, metadata) =
    NMRData(data, dims, refdims, name, metadata)

# Array interface (AbstractDimensionalArray takes care of everything else)
Base.@propagate_inbounds Base.setindex!(A::NMRData, x, I::Vararg{DimensionalData.StandardIndices}) =
    setindex!(data(A), x, I...)

function settval(A::NMRData, tvals)
    haspseudodimension(A) || throw(NMRToolsException("cannot set t values: data does not have a non-frequency dimension."))

    taxis = Ti(tvals, metadata=metadata(A,Ti))
    # two cases to consider: (x, t) or (x, y, t)
    if ndims(A) == 2
        newdims = (dims(A,X), taxis)
    else
        newdims = (dims(A,X), dims(A,Y), taxis)
    end
    NMRData(A.data, newdims, A.name, refdims=A.refdims, metadata=A.metadata)
end

"""
    @traitdef HasPseudoDimension{D}

A trait indicating whether the data object has a non-frequency domain dimension.

# Example

```julia
@traitfn f(x::X) where {X; HasPseudoDimension{X}} = "Fake!"
@traitfn f(x::X) where {X; Not{HasPseudoDimension{X}}} = "The real deal!"
```
"""
@traitdef HasPseudoDimension{D}

# adapted from expansion of @traitimpl
# this is a pure function because it only depends on the type definition of NMRData
Base.@pure function SimpleTraits.trait(t::Type{HasPseudoDimension{D}}) where D<:NMRData{T,N,A} where {T,N,A}
    any(map(dim->(typeintersect(dim, TimeDim) != Union{}), A.parameters)) ? HasPseudoDimension{D} : Not{HasPseudoDimension{D}}
end

"""
    haspseudodimension(spectrum)

Return true if the spectrum contains a non-frequency domain dimension.

# Example

```julia
julia> y2=loadnmr("exampledata/2D_HN/test.ft2");
julia> haspseudodimension(y2)
false
julia> y3=loadnmr("exampledata/pseudo3D_HN_R2/ft/test%03d.ft2");
julia> haspseudodimension(y3)
true
```
"""
function haspseudodimension(spectrum::NMRData{T,N,A}) where {T,N,A}
    any(map(dim->(typeintersect(dim, TimeDim) != Union{}), A.parameters))
end

"""
    scale(d::NMRData)

Return a scaling factor for the data combining the number of scans, receiver gain, and, if
specified, the sample concentration.
"""
function scale(d::NMRData)
    # get ns, rg, and concentration, with safe defaults in case missing
    ns = get(metadata(d), :ns, 1)
    rg = get(metadata(d), :rg, 1)
    conc = get(metadata(d), :conc, 1)
    return ns * rg * conc
end


text:
using ReinforcementLearningBase
using ReinforcementLearningCore
using ReinforcementLearningEnvironments
using ReinforcementLearningTrajectories

using Test
using CUDA
using CircularArrayBuffers
using Flux

@testset "ReinforcementLearningCore.jl" begin
    include("core.jl")
    include("utils/utils.jl")
end


text:
#########################################################################
#
#   rdiff differentiation for functions
#
#########################################################################

### main function

"""
Generates the derivative function for a given function

	rdiff( func::Function, init::Tuple; kwargs...)

Arguments:

- func: is a Julia generic function.
- init: is a tuple containing initial values for each parameter of ``func``. These reference values are needed to to fully evaluate ``ex``, this is a requirement of the derivation algorithm). By default the generated expression will yield the derivative for each variable given unless the variable is listed in the ``ignore`` argument.
- order: (keyword arg, default = 1) is an integer indicating the derivation order (1 for 1st order, etc.). Order 0 is allowed and will produce a function that is a processed version of ``ex`` with some variables names rewritten and possibly some optimizations.
- evalmod: (keyword arg, default=Main) module where the expression is meant to be evaluated. External variables and functions should be evaluable in this module.
- debug: (keyword arg, default=false) if true ``rdiff`` dumps the graph of the generating expression, instead of the expression.
- allorders: (keyword arg, default=true) tells rdiff whether to generate the code for all orders up to ``order`` (true) or only the last order.
- ignore: (keyword arg, default=[]) do not differentiate against the listed variables, useful if you are not interested in having the derivative of one of several variables in ``init``.

```julia
julia> rosenbrock(x) = (1 - x[1])^2 + 100(x[2] - x[1]^2)^2   # function to be derived
julia> rosen2 = rdiff(rosenbrock, (ones(2),), order=2)       # orders up to 2
```

"""
function rdiff(f::Function, sig0::Tuple; args...)
    sig = map( typeof, sig0 )
    fs = methods(f, sig)
    length(fs) == 0 && error("no function '$f' found for signature $sig")
    length(fs) > 1  && error("several functions $f found for signature $sig")  # is that possible ?

	if VERSION >= v"0.5.0-"
		fdef = fs[1].func
		fcode =  Expr(:block, Base.uncompressed_ast(fdef)...)

		# replace slotnames with their initial name
		fcode = deslot(fcode, fdef.slotnames)

		# function parameters
		fargs = map(t -> symbol(t[1]), Base.arg_decl_parts(fs[1])[2][2:end])
	else
		fdef = fs[1].func.code
		flambda = Base.uncompressed_ast(fdef)
		fcode = flambda.args[3]
		fargs = flambda.args[1]  # function parameters
	end

	cargs = [ (fargs[i], sig0[i]) for i in 1:length(sig0) ]
    ex  = transform(fcode) # TODO : add error messages if not parseable

    dex = rdiff(ex; args..., cargs...)

    # Note : new function is created in the same module as original function
    myf = fdef.module.eval( :( $(Expr(:tuple, fargs...)) -> $dex ) )
end

function deslot(ex::Expr, slotnames)
    args = Any[]
    for a in ex.args
        ar = if isa(a,Expr)
                deslot(a, slotnames)
			 elseif isa(a,Slot)
                slotnames[a.id]
			 else
                a
             end
        push!(args, ar)
    end
    Expr(ex.head, args...)
end

### translation functions to recover a workable expression that can be differentiated

# Simplifies expressions for processing
#  - removes Topnodes and linenumbers,
#  - replaces GenSym() with actual symbol
function streamline(ex0::Expr)
    ex = copy(ex0)

    ex.head == :call && isa(ex.args[1], TopNode) && (ex.args[1] = ex.args[1].name)

    args = Any[]
    for a in ex.args
        isa(a, LineNumberNode) && continue
        isa(a, Expr) && a.head==:line && continue

        ar = if isa(a,Expr)
                streamline(a)
             elseif isdefined(:GenSym) && isa(a, GenSym)
                symbol("__gensym$(a.id)")
             else
                a
             end
        push!(args, ar)
    end
    Expr(ex.head, args...)
end

# converts expression to searchable strings
function _e2s(ex::Expr, escape=false)
    ex.head == :macrocall && ex.args[1] == symbol("@rg_str") && return(ex.args[2])

    if ex.head == :call && ex.args[1] == :gotoifnot
        es = "↑gotoifnot"
        ra = 2:length(ex.args)
    else
        es = "↑$(ex.head)"
        ra = 1:length(ex.args)
    end

    for a in ex.args[ra]
        es *= "→" * _e2s(a, escape)
    end
    return es * "↓"
end

function _e2s(thing, escape=false)
    # res = isa(thing, Symbol) ? ":" * string(thing) : repr(thing)
    if isa(thing, Symbol)
        res = ":" * string(thing)
    elseif isdefined(:GlobalRef) && isa(thing, GlobalRef)
        res = _e2s(Expr(:., symbol(thing.mod), QuoteNode(thing.name)))
        # res = string(thing.mod) * "." * string(thing.name)
    else
        res = repr(thing)
    end

    escape || return(res)
    # now escape characters that would otherwise have a meaning in regex
    i = start(res)
    res2 = ""
    while !done(res,i)
        c, j = next(res,i)
        c in "()+*.\$^[]|?" && (res2 *= "\\")
        res2 *= string(c)
        i = j
    end
    res2
end

function e2s(ex::Expr, escape=false)
    if ex.head in [:body, :block]
        return mapreduce(e -> _e2s(e, escape), *, "", ex.args)
    else
        return _e2s(ex, escape)
    end
end

# converts searchable strings back to expressions
function _s2e(s::AbstractString, pos=1)
    cap = match( r"↑([^→↓]*)(.*)", s, pos )
    if cap == nothing # skip junk characters (Labelnodes,..) and return
        cap = match( r".*?↑(.*)", s, pos )
        cap == nothing && return nothing, endof(s)+1
        return nothing, cap.offsets[1]
    end

    he  = symbol(cap.captures[1])
    ar  = Any[]
    pos = cap.offsets[2]
    while s[pos] == '→' && !done(s, pos)
        cap = match( r"→([^→↓]*)(.*)↓$", s, pos )  # s[pos:end]
        cap == nothing && error("[s2e] unexpected string (2)")
        cap1 = cap.captures[1]
        if cap1[1] == '↑'
            ex, pos2 = _s2e(s, cap.offsets[1])
        elseif length(cap1) > 4 && cap1[1:3] == ":(:"    # Quotenodes
            ex = QuoteNode(symbol(cap1[4:end-1]))
            pos2 = cap.offsets[2]
        elseif cap1[1] == ':'        # symbols
            ex = symbol(cap1[2:end])
            pos2 = cap.offsets[2]
        else
            ex = parse(cap1)
            pos2 = cap.offsets[2]
        end
        push!(ar, ex)
        pos = pos2
    end

    c, pos = next(s, pos)
    return Expr(he, ar...), pos
end

function s2e(s::AbstractString)
    res = Expr[]
    pos = 1
    while !done(s, pos)
        ex, pos = _s2e(s, pos)
        ex != nothing && push!(res, ex)
    end
    res
end

if VERSION >= v"0.5.0-"
	# `for` loop search regex string (julia v0.5)
	function formatch(s::AbstractString)
		exreg = quote
		    rg"(?<pre>.*?)"
		    rg"(?<g0>:[#_].+?)" = rg"(?<range>.+?)"
		    rg"(?<iter>.+)" = start(rg"\g{g0}")
			rg":\((?<lab2>\d+): \)"
		    gotoifnot( !(done(rg"\g{g0}", rg"\g{iter}" )) , rg"(?<lab1>\d+)" )
		    rg"(?<g1>.+?)" = next(rg"\g{g0}", rg"\g{iter}")
		    rg"(?<idx>.+?)" = rg":(?:getfield|tupleref)"(rg"\g{g1}", 1)
		    rg"\g{iter}"    = rg":(?:getfield|tupleref)"(rg"\g{g1}", 2)
		    rg"(?<in>.*)"
			rg"(?:\(.*\))?"
		    rg":\((?<lab3>\d+): \)"
		    rg":\(goto \g{lab2}\)"
			rg":\(\g{lab1}: \)"
			rg":\((?<lab4>\d+): \)"
		    rg"(?<post>.*)"
		end
		rexp = Regex(e2s(streamline(exreg), true))

		mm = match(rexp, s)
		if mm != nothing && length(mm.captures) >= 11
			return mm.captures[[1, 3, 8, 9, 12]] # pre, rg, idx, inside, post
		else
			return nothing, nothing, nothing, nothing, nothing
		end
	end
else
	# `for` loop search regex string (julia v0.3.3 + 0.4 latest)
	function formatch(s::AbstractString)
		exreg = quote
		    rg"(?<pre>.*?)"
		    rg"(?<g0>:[#_].+?)" = rg"(?<range>.+?)"
		    rg"(?<iter>.+)" = start(rg"\g{g0}")
		    gotoifnot( !(done(rg"\g{g0}", rg"\g{iter}" )) , rg"(?<lab1>\d+)" )
		    rg":\((?<lab2>\d+): \)"
		    rg"(?<g1>.+?)" = next(rg"\g{g0}", rg"\g{iter}")
		    rg"(?<idx>.+?)" = rg":(?:getfield|tupleref)"(rg"\g{g1}", 1)
		    rg"\g{iter}"    = rg":(?:getfield|tupleref)"(rg"\g{g1}", 2)
		    rg"(?<in>.*)"
		    rg":\((?<lab3>\d+): \)"
		    gotoifnot( !(!(done(rg"\g{g0}", rg"\g{iter}"))) , rg"\g{lab2}" )
		    rg":\(\g{lab1}: \)"
		    rg"(?<post>.*)"
		end
		rexp = Regex(e2s(streamline(exreg), true))

		mm = match(rexp, s)
		if mm != nothing && length(mm.captures) >= 11
	        return mm.captures[[1,3,8,9,11]] # pre, rg, idx, inside, post
		else
			return nothing, nothing, nothing, nothing, nothing
		end
	end
end

function _transform(s::AbstractString)
	pre, rg, idx, inside, post = formatch(s)
    if pre != nothing
        exin = _transform(inside)
        ef = Expr(:for, Expr(:(=), symbol(idx[2:end]), s2e(rg)[1] ), exin)
        return Expr(:block, [ s2e(pre) ; ef ; s2e(post)]...)
    else
        return Expr(:block, s2e(s)...)
    end
end

function transform(ex::Expr)
    s = e2s(streamline(ex))
    tex = _transform(s)

    # remove return statement at the end
    rex = tex.args[end]
    if rex.head == :return
        tex.args[end] = rex.args[1]
    else
        error("[transform] no return statement found at the end")
    end

    tex
end


text:
################
#@require Sundials begin
import Sundials

# note that using a Jacobian is messy in Sundials

# Monkey-patch Sundials v0.2
# TODO: remove this after a new Sundials.jl version gets tagged
eval(Sundials, quote
macro checkflag(ex)
    # Insert a check that the given function call returns 0,
    # throw an error otherwise. Only apply directly to function calls.
    @assert Base.Meta.isexpr(ex, :call)
    fname = ex.args[1]
    quote
        flag = $(esc(ex))
        if flag != 0
            error($(string(fname, " failed with error code = ")), flag)
        end
        flag
    end
end
end)

extraoutputs = 10^6
const sundialssolvers = Dict{Any,Solver}()

# CVode wrapper to high level interface: Sundials.cvode
#####

function wrapper_CVODE_simple{N,T}(tr::TestRun{N,T})
    # Wraps the specific MyPkg.mysolver such that it works within
    # IVPTestSuite setup.
    tc = tr.tc
    so = tr.solver
    ###
    # 0) Wrap tc.fn!, tc.jac!, tc.mass! if necessary
    # not necessary
    ###
    # 1) Make call signature

    tsteps = collect(linspace(tc.tspan[1], tc.tspan[2], extraoutputs)) # TODO: fix by setting mxstep higher
    args = (tc.fn!, tc.ic, tsteps)
    kwargs  = ((:reltol, tr.reltol), (:abstol, tr.abstol))

    ###
    # 2) Call solver, if it does not succeed throw an error (if that
    # is not done anyway)
    out = so.solverfn(args...; kwargs...)
    # (probably no need to modify this section)

    ###
    # 3) Transform output to conform to standard:
    # tend -- end time reached
    # yend -- solution at tend
    # stats -- statistics, if available: (steps_total,steps_accepted, fn_evals, jac_evals, linear_solves)
    #                         otherwise  (-1, -1, -1, -1, -1)
    tend = tc.tspan[2]
    yend = squeeze(out[end,:],1)
    stats = (-1, -1, -1, -1, -1)
    return tend, yend, stats
end
daeindex = 0
adaptive = true
cvode_simple = Solver{:im}(Sundials.cvode, Sundials, wrapper_CVODE_simple, stiff, adaptive, daeindex, explicit_eq)
sundialssolvers[Sundials.cvode] = cvode_simple

# CVode wrapper to low level interface: Sundials.CVode
#####

# Note that this function cannot be nested inside wrapper_CVODE
function cvodefun(t::Float64, y::Sundials.N_Vector, yp::Sundials.N_Vector, fn!::Function)
    y = Sundials.asarray(y)
    yp = Sundials.asarray(yp)
    fn!(t, y, yp)
    return Int32(0)
end

function wrapper_CVODE{N,T}(tr::TestRun{N,T})
    # Adapted from the Sundials.jl wrapper

    # Wraps the specific MyPkg.mysolver such that it works within
    # IVPTestSuite setup.
    tc = tr.tc
    so = tr.solver
    ###
    # 0) Wrap tc.fn!, tc.jac!, tc.mass! if necessary
    # not necessary

    ###
    # 1) Make call signature

    mem = Sundials.CVodeCreate(Sundials.CV_BDF, Sundials.CV_NEWTON)
    if mem == C_NULL
        error("Failed to allocate CVODE solver object")
    end

    flag = Sundials.@checkflag Sundials.CVodeInit(mem,
                              cfunction(cvodefun, Int32, (Sundials.realtype, Sundials.N_Vector, Sundials.N_Vector, Ref{Function})),
                              tc.tspan[1], Sundials.nvector(tc.ic))
    flag = Sundials.@checkflag Sundials.CVodeSetUserData(mem, tc.fn!)
    flag = Sundials.@checkflag Sundials.CVodeSStolerances(mem, tr.reltol, tr.abstol)
    flag = Sundials.@checkflag Sundials.CVodeSetMaxNumSteps(mem, -1)
    flag = Sundials.@checkflag Sundials.CVDense(mem, tc.dof)
    y = copy(tc.ic)
    tout = [0.0]

    ###
    # 2) Call solver, if it does not succeed throw an error (if that
    # is not done anyway)
    for k in 2:length(tc.tspan)
        flag = so.solverfn(mem, tc.tspan[k], y, tout, Sundials.CV_NORMAL) # solverfn==CVode
    end

    ###
    # 3) Transform output to conform to standard:
    # tend -- end time reached
    # yend -- solution at tend
    # stats -- statistics, if available: (steps_total,steps_accepted, fn_evals, jac_evals, linear_solves)
    #                         otherwise  (-1, -1, -1, -1, -1)
    tend = tout[1]
    yend = y
    stats = (-1, -1, -1, -1, -1) # TODO
    return tend, yend, stats
end
daeindex = 0
adaptive = true
cvode = Solver{:im}(Sundials.CVode, Sundials, wrapper_CVODE, stiff, adaptive, daeindex, explicit_eq)
sundialssolvers[Sundials.CVode] = cvode

# IDA wrapper to high-level interface: Sundials.idasol
#####

function wrapper_IDA_simple{N,T}(tr::TestRun{N,T})
    # Wraps the specific MyPkg.mysolver such that it works within
    # IVPTestSuite setup.
    tc = tr.tc
    so = tr.solver
    ###
    # 0) Wrap tc.fn!, tc.jac!, tc.mass! if necessary

    if hasmass(tc)
        function residual!(t, y, dydt, res)
            m = tc.mass!()
            tc.mass!(t,y,m)
            tc.fn!(t,y,res)
            res[:] = res-m*dydt
        end
    else
        function residual!(t, y, dydt, res)
            tc.fn!(t,y,res)
            res[:] = res-dydt
        end
    end

    ###
    # 1) Make call signature
    ic_dydt = tc.fn!()
    tc.fn!(tc.tspan[1], tc.ic, ic_dydt)
    tsteps = collect(linspace(tc.tspan[1], tc.tspan[2], extraoutputs)) # TODO: fix by setting mxstep higher
    args = (residual!, tc.ic, ic_dydt, tsteps)
    kwargs  = ((:reltol, tr.reltol), (:abstol, tr.abstol))

    ###
    # 2) Call solver, if it does not succeed throw an error (if that
    # is not done anyway)
    y,dydt = so.solverfn(args...; kwargs...)
    # (probably no need to modify this section)

    ###
    # 3) Transform output to conform to standard:
    # tend -- end time reached
    # yend -- solution at tend
    # stats -- statistics, if available: (steps_total,steps_accepted, fn_evals, jac_evals, linear_solves)
    #                         otherwise  (-1, -1, -1, -1, -1)
    tend = tc.tspan[2]
    yend = squeeze(y[end,:],1)
    stats = (-1, -1, -1, -1, -1)
    return tend, yend, stats
end
daeindex = 1
adaptive = true
ida_simple = Solver{:im}(Sundials.idasol, Sundials, wrapper_IDA_simple, stiff, adaptive, daeindex, implicit_eq)
sundialssolvers[Sundials.idasol] = ida_simple

# IDA wrapper to low-level interface: Sundials.ISASolve
#####

function idasolfun(t::Float64, y::Sundials.N_Vector, yp::Sundials.N_Vector, r::Sundials.N_Vector, userfun::Function)
    y = Sundials.asarray(y)
    yp = Sundials.asarray(yp)
    r = Sundials.asarray(r)
    userfun(t, y, yp, r)
    return Int32(0)   # indicates normal return
end

function wrapper_IDA{N,T}(tr::TestRun{N,T})
    # Adapted from the Sundials.jl wrapper

    # Wraps the specific MyPkg.mysolver such that it works within
    # IVPTestSuite setup.
    tc = tr.tc
    so = tr.solver
    ###
    # 0) Wrap tc.fn!, tc.jac!, tc.mass! if necessary
    if hasmass(tc)
        function residual!(t, y, dydt, res)
            m = tc.mass!()
            tc.mass!(t,y,m)
            tc.fn!(t,y,res)
            res[:] = res-m*dydt
        end
    else
        function residual!(t, y, dydt, res)
            tc.fn!(t,y,res)
            res[:] = res-dydt
        end
    end

    ###
    # 1) Make call signature
    ic_dydt = tc.fn!()
    tc.fn!(tc.tspan[1], tc.ic, ic_dydt)

    mem = Sundials.IDACreate()
    if mem == C_NULL
        error("Failed to allocate CVODE solver object")
    end

    flag = Sundials.@checkflag Sundials.IDAInit(mem, cfunction(idasolfun, Int32, (Sundials.realtype, Sundials.N_Vector, Sundials.N_Vector, Sundials.N_Vector, Ref{Function})),
                            tc.tspan[1], Sundials.nvector(tc.ic), Sundials.nvector(ic_dydt))
    flag = Sundials.@checkflag Sundials.IDASetUserData(mem, residual!)
    flag = Sundials.@checkflag Sundials.IDASStolerances(mem, tr.reltol, tr.abstol)
    flag = Sundials.@checkflag Sundials.IDADense(mem, tc.dof)
    flag = Sundials.@checkflag Sundials.IDASetMaxNumSteps(mem, -1)

    rtest = zeros(T, tc.dof)
    residual!(tc.tspan[1], tc.ic, ic_dydt, rtest)
    if any(abs(rtest) .>= tr.reltol)
        error("Inconsistent ICs!  This shouldn't happen...")
    end
    yres = zeros(T,length(tc.tspan), length(tc.ic))
    ypres = zeros(T,length(tc.tspan), length(tc.ic))
    yres[1,:] = tc.ic
    ypres[1,:] = ic_dydt
    y = copy(tc.ic)
    yp = copy(ic_dydt)
    tout = [0.0]

    ###
    # 2) Call solver, if it does not succeed throw an error (if that
    # is not done anyway)
    for k in 2:length(tc.tspan)
        flag = so.solverfn(mem, tc.tspan[k], tout, y, yp, Sundials.IDA_NORMAL) # solverfn==IDASolve
    end

    ###
    # 3) Transform output to conform to standard:
    # tend  -- end time reached
    # yend  -- solution at tend
    # stats -- statistics, if available: (steps_total,steps_accepted, fn_evals, jac_evals, linear_solves)
    #                         otherwise  (-1, -1, -1, -1, -1)
    tend = tout[1]
    yend = y
    stats = (-1, -1, -1, -1, -1) # TODO
    return tend, yend, stats
end
daeindex = 1
adaptive = true
ida = Solver{:im}(Sundials.IDASolve, Sundials, wrapper_IDA, stiff, adaptive, daeindex, implicit_eq)

sundialssolvers[Sundials.IDASolve] = ida

merge!(allsolvers, sundialssolvers)


text:
function allexamples()
    include("Convex.jl")
    include("Convex2d.jl")
    include("Grid.jl")
    include("Mesh.jl")
    include("numbering.jl")
    include("Polygon.jl")
    include("Polyhedron.jl")
    include("Text.jl")
end

allexamples()


text:
## splitting up searchcubes

function cuberanges(center, widths)
    t = eltype(center)
    return  NTuple{length(center),Tuple{t,t}}([(center[i]-widths[i], center[i]+widths[i]) for i=1:length(center)])
end

cuberanges(R::RotationVec, T::SVector{3}, σᵣ, σₜ) = cuberanges((R.sx,R.sy,R.sz,T[1],T[2],T[3]), (σᵣ,σᵣ,σᵣ,σₜ,σₜ,σₜ))
cuberanges(R::RotationVec, σᵣ) = cuberanges((R.sx,R.sy,R.sz), (σᵣ,σᵣ,σᵣ))
cuberanges(T::SVector{3}, σₜ) = cuberanges((T[1],T[2],T[3]), (σₜ,σₜ,σₜ))

"""
    sbrngs = subranges(ranges, nsplits)

Takes `ranges`, a nested tuple describing intervals for each dimension in rigid-transformation space
defining a hypercube, and splits the hypercube into `nsplits` even components along each dimension.
If the cube is N-dimensional, the number of returned sub-cubes will be `nsplits^N`.
"""
function subranges(ranges, nsplits::Int)
    t = eltype(eltype(ranges))

    # calculate even splititng points for each dimension
    splitvals = [range(r[1], stop=r[2], length=nsplits+1) |> collect for r in ranges]
    splits = [[(splitvals[i][j], splitvals[i][j+1]) for j=1:nsplits] for i=1:6]
    f(x) = splits[x[1]][x[2]]
    children = fill(ranges, nsplits^6)
    for (i,I) in enumerate(CartesianIndices(NTuple{6,UnitRange{Int}}(fill(1:nsplits, 6))))
        children[i] = NTuple{6,Tuple{t,t}}(map(x->f(x), enumerate(Tuple(I))))
    end
    return children
end

## supertype for search regions
abstract type SearchRegion{T} end
CoordinateTransformations.AffineMap(sr::SearchRegion) = AffineMap(sr.R, sr.T)

## rigid transformation

"""
Describes an transformation uncertainty region centered at rotation R and translation T, with rotation and translation half-widths of σᵣ and σₜ respectively
"""
struct UncertaintyRegion{T<:Real} <: SearchRegion{T}
    R::RotationVec{T}
    T::SVector{3,T}
    σᵣ::T
    σₜ::T
    ranges::NTuple{6,Tuple{T,T}}
end
UncertaintyRegion(R,T,σᵣ,σₜ) = UncertaintyRegion(R, T, σᵣ, σₜ, cuberanges(R, T, σᵣ, σₜ))
UncertaintyRegion(σᵣ, σₜ)    = UncertaintyRegion(one(RotationVec), zero(SVector{3}), σᵣ, σₜ)
UncertaintyRegion(σₜ)        = UncertaintyRegion(one(RotationVec), zero(SVector{3}), π, σₜ)
UncertaintyRegion()         = UncertaintyRegion(one(RotationVec), zero(SVector{3}), π, Inf)

# for speeding up hashing and performance of the priority queue in the branch and bound procedure
const hash_UncertaintyRegion_seed = UInt === UInt64 ? 0x4de49213ae1a23bf : 0xef78ce68
function hash(B::UncertaintyRegion, h::UInt)
    h += hash_UncertaintyRegion_seed
    h = hash(B.ranges, h)
    return h
end

## only rotation
struct RotationRegion{T<:Real} <: SearchRegion{T}
    R::RotationVec{T}
    T::SVector{3,T}
    σᵣ::T
    ranges::NTuple{3,Tuple{T,T}}
end
RotationRegion(R,T,σᵣ)   = RotationRegion(R, T, σᵣ, cuberanges(R, σᵣ))
RotationRegion(σᵣ)       = RotationRegion(one(RotationVec), zero(SVector{3}), σᵣ)
RotationRegion()         = RotationRegion(one(RotationVec), zero(SVector{3}), π)

UncertaintyRegion(rr::RotationRegion) = UncertaintyRegion(rr.R, rr.T, rr.σᵣ, 0)
RotationRegion(ur::UncertaintyRegion) = RotationRegion(ur.R, ur.T, ur.σᵣ)

# for speeding up hashing and performance of the priority queue in the branch and bound procedure
const hash_RotationRegion_seed = UInt === UInt64 ? 0xee63e114344da2b9 : 0xe6cb1eb7
function hash(B::RotationRegion, h::UInt)
    h += hash_RotationRegion_seed
    h = hash(B.ranges, h)
    return h
end

## only translation
struct TranslationRegion{T<:Real} <: SearchRegion{T}
    R::RotationVec{T}
    T::SVector{3,T}
    σₜ::T
    ranges::NTuple{3,Tuple{T,T}}
end
TranslationRegion(R,T,σₜ)   = TranslationRegion(R, T, σₜ, cuberanges(T, σₜ))
TranslationRegion(σₜ)       = TranslationRegion(one(RotationVec), zero(SVector{3}), σₜ)
TranslationRegion()        = TranslationRegion(one(RotationVec), zero(SVector{3}), Inf)

UncertaintyRegion(tr::TranslationRegion) = UncertaintyRegion(tr.R, tr.T, 0, tr.σₜ)
TranslationRegion(ur::UncertaintyRegion) = TranslationRegion(ur.R, ur.T, ur.σₜ)

# for speeding up hashing and performance of the priority queue in the branch and bound procedure
const hash_TranslationRegion_seed = UInt === UInt64 ? 0x24f59aedb6bf903f : 0x76f5f734
function hash(B::TranslationRegion, h::UInt)
    h += hash_TranslationRegion_seed
    h = hash(B.ranges, h)
    return h
end

## initialization of search spaces

# """
#     lim = max_translation(x, y)

# Computes the largest translation needed to ensure that the searchspace contains the best alignment transformation.
# """
# function max_translation(gmmx::AbstractSingleGMM, gmmy::AbstractSingleGMM)
#     trlim = typemin(promote_type(numbertype(gmmx),numbertype(gmmy)))
#     for gaussians in (gmmx.gaussians, gmmy.gaussians)
#         if !isempty(gaussians)
#             trlim = max(trlim, maximum(gaussians) do gauss
#                     maximum(abs, gauss.μ) end)
#         end
#     end
#     return trlim
# end


text:
const PERMS = S_IRUSR | S_IWUSR | S_IRGRP | S_IROTH # 0644

struct NamedSemaphore
    name::String
    handle::Ptr{Nothing}

    function NamedSemaphore(name::String, create::Bool=true, create_exclusive::Bool=false)
        new(name, sem_open(name, create, create_exclusive))
    end
end

lock(sem::NamedSemaphore) = sem_wait(sem.handle)
trylock(sem::NamedSemaphore) = sem_trywait(sem.handle)
unlock(sem::NamedSemaphore) = sem_post(sem.handle)
close(sem::NamedSemaphore) = sem_close(sem.handle)
delete!(sem::NamedSemaphore) = sem_unlink(sem.name)
count(sem::NamedSemaphore) = sem_getvalue(sem.handle)
reset(sem::NamedSemaphore) = error("can not reset a NamedSemaphore")

function serialize(s::AbstractSerializer, sem::NamedSemaphore)
    Serializer.serialize_type(s, typeof(sem))
    serialize(s, sem.name)
end

function deserialize(s::AbstractSerializer, ::Type{NamedSemaphore})
    name = deserialize(s)
    NamedSemaphore(name)
end

# create or open a semaphore
# semaphore is created in unlocked state (i.e. value is 1)
function sem_open(name::String, create::Bool=true, create_exclusive::Bool=false)
    if create
        flags = JL_O_CREAT
        if create_exclusive
            flags |= JL_O_EXCL
        end
        handle = ccall(:sem_open, Ptr{Nothing}, (Cstring, Cint, Base.Cmode_t, Cuint), name, flags, PERMS, Cuint(1))
    else
        handle = ccall(:sem_open, Ptr{Nothing}, (Cstring, Cint), name, JL_O_RDWR)
    end
    systemerror("error creating/opening semaphore", handle == C_NULL)
    handle
end

# unlock a semaphore
function sem_post(handle::Ptr{Nothing})
    ret = ccall(:sem_post, Cint, (Ptr{Nothing},), handle)
    systemerror("error unlocking semaphore", ret != 0)
end

# lock a semaphore
function sem_wait(handle::Ptr{Nothing})
    ret = ccall(:sem_wait, Cint, (Ptr{Nothing},), handle)
    systemerror("error locking semaphore", ret != 0)
end

function sem_trywait(handle::Ptr{Nothing})
    ret = ccall(:sem_trywait, Cint, (Ptr{Nothing},), handle)
    eagain = (Libc.errno() == Libc.EAGAIN)
    locked = (ret == 0)
    systemerror("error locking semaphore", !locked && !eagain)
    locked
end

# close a semaphore (in this process, semaphore will still exist in system)
function sem_close(handle::Ptr{Nothing})
    ret = ccall(:sem_close, Cint, (Ptr{Nothing},), handle)
    systemerror("error closing semaphore", ret != 0)
end

# unlink a semaphore (semaphore will be deleted from system once all opened references to it are closed)
function sem_unlink(name::String)
    ret = ccall(:sem_unlink, Cint, (Cstring,), name)
    systemerror("error unlinking semaphore", ret != 0)
end

# Gets the current value of the semaphore.
# If one or more processes or threads are blocked waiting to lock the semaphore with sem_wait(3),
# POSIX.1 permits two possibilities for the value returned in sval: either 0 is returned; or a
# negative number whose absolute value is the count of the number of processes and threads currently
# blocked in sem_wait(3). Linux adopts the former behavior.
function sem_getvalue(handle::Ptr{Nothing})
    val = Ref{Cint}(0)
    ret = ccall(:sem_getvalue, Cint, (Ptr{Nothing},Ptr{Cint}), handle, val)
    systemerror("error getting semaphore value", ret != 0)
    val[]
end


text:
# This file was generated, do not modify it. # hide
x = 5
println("hello")
x^2

max_stars_repo_name: mcognetta/franklindocs

max_issues_repo_name: mcognetta/franklindocs

max_forks_repo_name: mcognetta/franklindocs


text:
using st16SseqJuliaTools
using ArgParse2
using DataFrames
using Query
using CSV
function julia_main()::Cint
    parser = ArgumentParser( description = "Extract clusters of a proper size and without close chimeras")
    #add_argument!(parser, "-i", "--input-fasta", type = String, help = "Input file")
    add_argument!(parser, "-a", "--cluster-lower-level", type = String, help = "Clustering to group reads for assembly")
    add_argument!(parser, "-b", "--cluster-higher-level", type = String, help = "Clustering to group reads to check 4 chimeras")
    add_argument!(parser, "-t", "--tax-assign", type = String, help = "Taxonomy assignment files - down to genus level. Should match the lower level clustering fasta file ")
    add_argument!(parser, "-m", "--min-cluster-size", type = Float64, help = "Minimum cluster size")
    add_argument!(parser, "-f", "--fasta", type = String, help = "Fasta file matching lower lever of clustering ")
    add_argument!(parser, "-d", "--dir", type = String, help = "Directory for output", default = "grouping_info")
    args = parse_args(parser)
    if !isdir(args.dir)
        mkdir(args.dir)
    end
    aln_all = Alignment(args.fasta)
    centroid_genus = parse_dada_results(args.tax_assign)
    cl1  = read_jc_clusters(args.cluster_lower_level)
    #check which cluster has proper size
    proper_sizes = mark_proper_size_clusters(cl1,args.min_cluster_size)
    #get higher level cluster - it should group such close species as E. Coli & Salmonela and their chimeric reads
    cl2   = read_jc_clusters(args.cluster_higher_level)
    recomb_data = DataFrame()
    groupct = 0
    for cl in keys(cl2)
        ids_for_common_analyses = Array{String,1}()
        for rep in cl2[cl]
                if proper_sizes[rep]
                    push!(ids_for_common_analyses,rep)
                end
        end
        if length(ids_for_common_analyses) > 2
            #a case where chimeric cluster can be formed
            groupct += 1
            #println("------------------")
            #for id in ids_for_common_analyses
            #    println(id," ",length(cl1[id])," ",centroid_genus[id]," ", proper_sizes[id])
            #end
            sub_aln = sub_alignment(aln_all,ids_for_common_analyses)
            recombs = eval_recomb(sub_aln)
            recombs.Group = fill(groupct,nrow(recombs))
            recomb_data = vcat(recomb_data,recombs)
        end
    end
    chimof = args.dir*   "/chimeric.csv"
    println(stderr,"Chimeric reads detection info: ", chimof)
    CSV.write(chimof, recomb_data)
    if nrow(recomb_data) > 0
        chimeric_ids = @from i in recomb_data begin
                                @where i.Unique_snps  < 2
                                @select { i.Read_name}
                                @collect DataFrame
                        end
    else
        chimeric_ids = DataFrame(Read_name = [])
    end
    #filter out detection

    #aln = sanitise.Alignment(args.input_fasta)
    #println(aln_conservation(aln,ignore_amounts=true))
    #construct main output data frame
    finald = DataFrame()
    read_names = collect(keys(cl1))
    finald.Centroids = read_names
    finald.Cluster_name = map( x -> centroid_genus[x], read_names)
    finald.Genus = map( x -> split(centroid_genus[x],"_")[1], read_names)
    finald.Size = map(x -> length(x),values(cl1))
    finald.Proper_size  = collect(values(proper_sizes))
    finald.Chimeric = map(x -> (x in chimeric_ids.Read_name),read_names)
    allof = args.dir* "/all_clusters.csv"
    println(stderr,"Data on all cluster: ", allof)
    CSV.write(allof, finald)
    chosen_clusters = @from i in finald begin
                            @where i.Proper_size  #&& !i.Chimeric
                            @select { i.Cluster_name, i.Genus,i.Centroids , i.Size, i.Chimeric}
                            @collect DataFrame
                        end

    insertcols!(chosen_clusters, 1,:ID => collect(1:nrow(chosen_clusters)))
    chosenof = args.dir*   "/chosen_clusters.csv"
    println(stderr,"Data on chosen cluster: ", chosenof)
    CSV.write(chosenof, chosen_clusters)
    #output contigs
    ct = 0
    for r in eachrow(chosen_clusters)
        ct = r.ID
        f=open(args.dir*"/$ct","w")
        for id in cl1[r.Centroids]
            println(f,id)
        end
        close(f)
        f=open(args.dir*"/$ct."*r.Cluster_name,"w")
        println(f,length(cl1[r.Centroids]))
        close(f)
    end

    return 0
end

julia_main()


text:
"""
    AbstractDataset

Abstract type for dataset objects.
"""
abstract type AbstractDataset end

function download(dataset::AbstractDataset)
    error("download method is not implemented.")
end

function load_dataset(dataset::AbstractDataset)
    error("load_dataset method is not implemented.")
end


text:
#####
##### First derivative operators
#####

for LX in (:ᶜ, :ᶠ), LY in (:ᶜ, :ᶠ), LZ in (:ᶜ, :ᶠ)

    x_derivative = Symbol(:∂x, LX, LY, LZ)
    x_spacing    = Symbol(:Δx, LX, LY, LZ)
    x_difference = Symbol(:δx, LX, :ᵃ, :ᵃ)

    y_derivative = Symbol(:∂y, LX, LY, LZ)
    y_spacing    = Symbol(:Δy, LX, LY, LZ)
    y_difference = Symbol(:δy, :ᵃ, LY, :ᵃ)

    z_derivative = Symbol(:∂z, LX, LY, LZ)
    z_spacing    = Symbol(:Δz, LX, LY, LZ)
    z_difference = Symbol(:δz, :ᵃ, :ᵃ, LZ)

    @eval begin
        @inline $x_derivative(i, j, k, grid, c) = $x_difference(i, j, k, grid, c) / $x_spacing(i, j, k, grid)
        @inline $y_derivative(i, j, k, grid, c) = $y_difference(i, j, k, grid, c) / $y_spacing(i, j, k, grid)
        @inline $z_derivative(i, j, k, grid, c) = $z_difference(i, j, k, grid, c) / $z_spacing(i, j, k, grid)

        @inline $x_derivative(i, j, k, grid, f::Function, args...) = $x_difference(i, j, k, grid, f, args...) / $x_spacing(i, j, k, grid)
        @inline $y_derivative(i, j, k, grid, f::Function, args...) = $y_difference(i, j, k, grid, f, args...) / $y_spacing(i, j, k, grid)
        @inline $z_derivative(i, j, k, grid, f::Function, args...) = $z_difference(i, j, k, grid, f, args...) / $z_spacing(i, j, k, grid)
    end
end

#####
##### Second, Third, and Fourth derivatives
#####

@inline insert_symbol(dir, L, L1, L2) =
                      dir == :x ?
                      (L, L1, L2) :
                      dir == :y ?
                      (L1, L, L2) :
                      (L1, L2, L)

for dir in (:x, :y, :z), L1 in (:ᶜ, :ᶠ), L2 in (:ᶜ, :ᶠ)

     first_order_face = Symbol(:∂,  dir, insert_symbol(dir, :ᶠ, L1, L2)...)
    second_order_face = Symbol(:∂², dir, insert_symbol(dir, :ᶠ, L1, L2)...)
     third_order_face = Symbol(:∂³, dir, insert_symbol(dir, :ᶠ, L1, L2)...)
    fourth_order_face = Symbol(:∂⁴, dir, insert_symbol(dir, :ᶠ, L1, L2)...)

     first_order_center = Symbol(:∂,  dir, insert_symbol(dir, :ᶜ, L1, L2)...)
    second_order_center = Symbol(:∂², dir, insert_symbol(dir, :ᶜ, L1, L2)...)
     third_order_center = Symbol(:∂³, dir, insert_symbol(dir, :ᶜ, L1, L2)...)
    fourth_order_center = Symbol(:∂⁴, dir, insert_symbol(dir, :ᶜ, L1, L2)...)

    @eval begin
        $second_order_face(i, j, k, grid, c) =  $first_order_face(i, j, k, grid, $first_order_center,  c)
         $third_order_face(i, j, k, grid, c) =  $first_order_face(i, j, k, grid, $second_order_center, c)
        $fourth_order_face(i, j, k, grid, c) = $second_order_face(i, j, k, grid, $second_order_face,   c)

        $second_order_center(i, j, k, grid, c) =  $first_order_center(i, j, k, grid, $first_order_face,    c)
         $third_order_center(i, j, k, grid, c) =  $first_order_center(i, j, k, grid, $second_order_face,   c)
        $fourth_order_center(i, j, k, grid, c) = $second_order_center(i, j, k, grid, $second_order_center, c)

        $second_order_face(i, j, k, grid, f::Function, args...) =  $first_order_face(i, j, k, grid, $first_order_center,  f::Function, args...)
         $third_order_face(i, j, k, grid, f::Function, args...) =  $first_order_face(i, j, k, grid, $second_order_center, f::Function, args...)
        $fourth_order_face(i, j, k, grid, f::Function, args...) = $second_order_face(i, j, k, grid, $second_order_face,   f::Function, args...)

        $second_order_center(i, j, k, grid, f::Function, args...) =  $first_order_center(i, j, k, grid, $first_order_face,    f::Function, args...)
         $third_order_center(i, j, k, grid, f::Function, args...) =  $first_order_center(i, j, k, grid, $second_order_face,   f::Function, args...)
        $fourth_order_center(i, j, k, grid, f::Function, args...) = $second_order_center(i, j, k, grid, $second_order_center, f::Function, args...)
    end
end

#####
##### Operators of the form A*∂(q) where A is an area and q is some quantity.
#####

for dir in (:x, :y, :z), LX in (:ᶜ, :ᶠ), LY in (:ᶜ, :ᶠ), LZ in (:ᶜ, :ᶠ)

    operator   = Symbol(:A, dir, :_∂, dir, LX, LY, LZ)
    area       = Symbol(:A, dir, LX, LY, LZ)
    derivative = Symbol(:∂, dir, LX, LY, LZ)

    @eval begin
        $operator(i, j, k, grid, c) = $area(i, j, k, grid) * $derivative(i, j, k, grid, c)
    end
end


text:
module Links
#################### TransformDistribution ####################

import Distributions: Continuous, ContinuousDistribution,
ContinuousUnivariateDistribution, ContinuousMultivariateDistribution,
## Univariate Distributions
Truncated, Cauchy, Gumbel, Laplace, Logistic, NoncentralT, Normal,
NormalCanon, TDist, Beta, BetaPrime, Chi, Chisq, Erlang, Exponential,
FDist, Frechet, Gamma, InverseGamma, InverseGaussian, Kolmogorov,
LogNormal, NoncentralChisq, NoncentralF, Rayleigh, Weibull, KSOneSided,
NoncentralBeta

import StatsFuns: logit

using Reexport
@reexport using Distributions

using BayesianTools
@reexport using BayesianTools

#################### Bounded Distribution ####################
const TransformDistribution{T<:ContinuousUnivariateDistribution} = Union{T,Truncated{T},Improper}

function link(d::TransformDistribution, x::Real)
  a, b = minimum(d), maximum(d)
  lowerbounded, upperbounded = isfinite(a), isfinite(b)
  if lowerbounded && upperbounded
    logit((x - a) / (b - a))
  elseif lowerbounded
    log(x - a)
  elseif upperbounded
    log(b - x)
  else
    x
  end
end

function invlink(d::TransformDistribution, x::Real)
  a, b = minimum(d), maximum(d)
  lowerbounded, upperbounded = isfinite(a), isfinite(b)
  if lowerbounded && upperbounded
    (b - a) * invlogit(x) + a
  elseif lowerbounded
    exp(x) + a
  elseif upperbounded
    b - exp(x)
  else
    x
  end
end

function logjacobian(d::TransformDistribution, x::Real)
    # x - a >= 0
    # b - x >= 0
    a, b = minimum(d), maximum(d)
    lowerbounded, upperbounded = isfinite(a), isfinite(b)
    if lowerbounded && upperbounded
      log((b-a)/((b-x)*(x-a)))
    elseif lowerbounded
      -log(x - a)
    elseif upperbounded
      -log(x - b)
    else
      0.0
    end
end
#################### RealDistribution ####################

const RealDistribution =Union{Cauchy, Gumbel, Laplace, Logistic, NoncentralT, Normal,
                              NormalCanon, TDist}

link(d::RealDistribution, x::Real) = x
invlink(d::RealDistribution, x::Real) = x
logjacobian(d::RealDistribution, x::Real) = zero(eltype(x))
#################### PositiveDistribution ####################

const PositiveDistribution = Union{BetaPrime, Chi, Chisq, Erlang, Exponential, FDist, Frechet,
                                   Gamma, InverseGamma, InverseGaussian, Kolmogorov, LogNormal,
                                   NoncentralChisq, NoncentralF, Rayleigh, Weibull}

link(d::PositiveDistribution, x::Real) = log(x)
invlink(d::PositiveDistribution, x::Real) = exp(x)
logjacobian(d::PositiveDistribution, x::Real) = log(x)
#################### UnitDistribution ####################

const UnitDistribution = Union{Beta, KSOneSided, NoncentralBeta}

link(d::UnitDistribution, x::Real) = logit(x)
invlink(d::UnitDistribution, x::Real) = invlogit(x)
## logjacobian is d link / d x = log (1/(1-x))
logjacobian(d::UnitDistribution, x::Real) = -log(x)-log(1.0-x)

invlogit(x::Real) = 1.0 / (exp(-x) + 1.0)
invlogit(x::AbstractArray) = map(invlogit, x)

export link, invlink, logjacobian

end

max_stars_repo_name: gragusa/ProductDistributions

max_issues_repo_name: gragusa/ProductDistributions

max_forks_repo_name: gragusa/ProductDistributions


text:
"""
    resize(files,dest,size)
Resize all the images in the image files named in the string vector `files`, using 2-vector `size` for [height,width]. Rewrite them as PNG files in directory `dest`.

Also computes the same data as `summary` and returns a named tuple of the result.
"""
function resize(folder::AbstractVector{String},dest,sz)
	N = length(folder)
	summary = (
		size = fill([0,0],N),
		avgG = zeros(N),
		purkrow = fill(NaN,N),
		purkcol = fill(NaN,N),
		freqR = fill([NaN],N),
		freqG = fill([NaN],N),
		freqB = fill([NaN],N)
	)
	@showprogress for (idx,fname) in enumerate(folder)
		X = load(fname)
		X = imresize(X,sz)
		pngname = splitext(basename(fname))[1]*".png"
		save(joinpath(dest,pngname),X)
		update_summary!(summary,idx,X)
	end
	return summary
end

"""
    summarize(files)
Get summary data for all the images in the image files named in the string vector `files`. Returns named tuple with fields:

* `size` [height,width] in pixels
* `purkrow` row index of the purkinje (float), `NaN` if none found
* `purkcol` column index of the purkinje (float), `NaN` if none found
* `avgG` average intensity in the green channel
* `freqR` histogram bin values for intensity in the red channel
* `freqG` histogram bin values for intensity in the green channel
* `freqB` histogram bin values for intensity in the blue channel
"""
function summarize(folder::AbstractVector{String})
	N = length(folder)
	summary = (
		size = fill([0,0],N),
		purkrow = fill(NaN,N),
		purkcol = fill(NaN,N),
		avgG = zeros(N),
		freqR = fill([NaN],N),
		freqG = fill([NaN],N),
		freqB = fill([NaN],N)
	)
	@showprogress for (idx,fname) in enumerate(folder)
		update_summary!(summary,idx,load(fname))
	end
	return summary
end

# utility for summary
function update_summary!(summary,idx,X)
	sz = size(X)
	summary.size[idx] = [sz...]

	purk = findpurkinje(X,0.6)
	if length(purk) > 50
		summary.purkrow[idx] = median(i[1] for i in purk)
		summary.purkcol[idx] = median(i[2] for i in purk)
	end
	edges = collect(LinRange(0,1,33))
	edges[1] = -eps()

	summary.avgG[idx] = mean(green.(X))

	w = fit(Histogram,vec(red.(X)),edges,closed=:right).weights
	summary.freqR[idx] = w/prod(sz)
	w = fit(Histogram,vec(green.(X)),edges,closed=:right).weights
	summary.freqG[idx] = w/prod(sz)
	w = fit(Histogram,vec(blue.(X)),edges,closed=:right).weights
	summary.freqB[idx] = w/prod(sz)

end

"""
	intensity(data)
Return a vector of an intensity measure in the green channel of each image of a movie. The input should be a DataFrame of the data returned by `summary`.
"""
function intensity(data::DataFrame)
	# Construct an average from the intensity histogram, but leaving out the last bin.
	# The reason is that most of the pixels in that bin are in the sclera, and the
	# average value becomes sensitive to small motions of the eyelids.
	# Might be better to use a median...
	bins = (0.5:31.5)/33
	intensity = Float64[]
	for (i,distr) in enumerate(data.freqG)
		# data might be in string form, if it was loaded from CSV–so parse it
		if !(eltype(distr) <: Number)
			distr = parse.(Float64,split(replace(distr[2:end-1],","=>" ")))
		end
		push!(intensity,sum(bins[j]*distr[j] for j in 1:31))
	end
	return intensity
end

"""
	darkframes(data)
Return a boolean vector indicating which images in a folder are likely too dark to be used. The input should be a DataFrame of the data returned by `summary`.
"""
function darkframes(data::DataFrame)
	q = intensity(data)
	ip,jp = data.purkrow,data.purkcol
	nopurk = isnan.(ip)
	isdark = @. ( q < 0.15 ) | ( nopurk & (q < 0.25) )

	# Some "bright" frames might still need to get marked dark, if they are a good
	# bit below the "typical bright" frame. Iterate until this stabilizes.
	ct = 3
	while ct < count(isdark)
		ct = count(isdark)
		μ = median(q[.!isdark])
		gap = 0.3*(μ - median(q[isdark]))
		isdark = @. isdark | (q < μ - gap)
	end
	return isdark
end

"""
	goodframes(data)
Return a boolean vector indicating which images in a folder appear to be in the first illuminated interblink period. The input should be a DataFrame of the data returned by `summary`.
"""
function goodframes(data::DataFrame)
	isdark = darkframes(data)
	N = length(isdark)
	# discard where at least 2 frames within each window of three are marked as dark
	w = [ sum(isdark[i:i+2]) for i in 1:N-2 ]
	fd = findall(w.>1)

	# Look for a jump in the indices of dark frames. The period within that jump is our
	# target.
	idx = findfirst(diff(fd).>2)
	!isnothing(idx) && (isdark[1:fd[idx]] .= true)
	!isnothing(idx) && (isdark[fd[idx+1]:end] .= true)

	return .!isdark
end

max_stars_repo_name: tobydriscoll/cornea-detection

max_issues_repo_name: tobydriscoll/cornea-detection

max_forks_repo_name: tobydriscoll/cornea-detection


text:
module DistancePlan

using StaticArrays

function norm(v)
  return sqrt(sum(v .^ 2))
end

function distance(center, obstacle::AbstractVector{T}) where {T <: AbstractFloat}
  return norm(center .- obstacle)
end

abstract type Volume end

struct Box <: Volume
  lows
  highs
end

export Box

function distance(center, obstacle::Box)
  diffs = [(obstacle.lows .- center) (center .- obstacle.highs)]
  return maximum(diffs)
end

struct Sphere <: Volume
  center
  radius
end

export Sphere

function distance(center, obstacle::Sphere)
  return norm(center .- obstacle.center) - obstacle.radius
end

function distance(center, obstacle::AbstractVector)
  return minimum([distance(center, obs) for obs=obstacle])
end

struct SphereTree
  centers::Vector
  squared_radii::Vector
  children::Vector{Vector{Int}}
  parent::Vector{Int}
  SphereTree() = new([], [], [], [])
end

export distance

function nearest_sphere(tree::SphereTree, point)::Int
  return argmin([sum((center .- point) .^ 2) - srad
                 for (center, srad) in zip(tree.centers, tree.squared_radii)])
end

function project(point, sphere::Sphere)
  v = point .- sphere.center
  v /= norm(v)
  v = v .* sphere.radius .+ sphere.center
  draw!(v, color="orange")
  return v
end

function sample(bounds::Box)::SVector{N} where {N}
  r = rand(N)
  x = SVector{N}((bounds.highs .- bounds.lows) .* r .+ bounds.lows)
  draw!(x, color="blue")
  return x
end

function push_sphere!(tree::SphereTree, parent::Int, sphere::Sphere)
  child = 1 + length(tree.parent)
  push!(tree.centers, sphere.center)
  push!(tree.squared_radii, sphere.radius ^ 2)
  push!(tree.parent, parent)
  push!(tree.children, [])
  push!(tree.children[parent], child)
end

function get_sphere(tree::SphereTree, sphere::Int)
  return Sphere(tree.centers[sphere], sqrt(tree.squared_radii[sphere]))
end

function evt_step!(tree::SphereTree, bounds::Box, obstacles::Vector; min_radius=10.) where {N}
  x = sample{N}(bounds)
  x_near = nearest_sphere(tree, x)
  x_new = project(x, get_sphere(tree, x_near))
  x_new_rad = distance(x_new, obstacles)
  # Sphere is of minimum radius and center is in bounds.
  # Center can be out of bounds due to projection from inside sphere.
  if x_new_rad > min_radius && distance(x_new, bounds) < 0
    push_sphere!(tree, x_near, Sphere(x_new, x_new_rad))
  end
end

function evt(origin::SVector, bounds::Box, obstacles::Vector, steps::Int) where {N}
  tree = SphereTree()
  radius = distance(origin, obstacles)
  origin_sphere = Sphere(origin, radius)
  push_sphere!(tree, 1, origin_sphere)
  for iter in 1:steps
    evt_step!{N}(tree, bounds, obstacles)
  end
  return tree
end

export evt

function get_path(tree, last_idx::Int)::Vector
  path = [tree.centers[last_idx]]
  while last_idx != 1
    last_idx = tree.parent[last_idx]
    pushfirst!(path, tree.centers[last_idx])
  end
  return path
end

function evt_goal(origin::SVector, bounds::Box, obstacles::Vector, goal::SVector; min_radius=10.)
  tree = SphereTree()
  radius = distance(origin, obstacles)
  origin_sphere = Sphere(origin, radius)
  push_sphere!(tree, 1, origin_sphere)
  while true
    gidx = nearest_sphere(tree, goal)
    goal_sphere = Sphere(tree.centers[gidx], sqrt(tree.squared_radii[gidx]))
    if distance(goal, goal_sphere) <= 0
      return tree, get_path(tree, gidx)
    end
    evt_step!(tree, bounds, obstacles, min_radius=min_radius)
  end
end

export evt_goal

import Makie

function draw!(box::Box; color="red")
  lx, ly = box.lows
  hx, hy = box.highs
  points = [
      lx ly
      hx ly
      hx hy
      lx hy
      lx ly
      NaN NaN
  ]
  Makie.plot!(points[:, 1], points[:, 2], color=color)
end

function draw!(sphere::Sphere; points=64, color="red")
  cx, cy = sphere.center
  r = sphere.radius
  t = range(0, 2pi, length=points)
  x = cx .+ r .* cos.(t)
  y = cy .+ r .* sin.(t)
  Makie.plot!(x, y, color=color)
end

function draw!(obstacle::SVector; delta=0.1, color="red")
  x, y = obstacle
  points = [
      x - delta y - delta
      x + delta y + delta
      NaN NaN
      x - delta y + delta
      x + delta y - delta
      NaN NaN
  ]
  Makie.plot!(points[:, 1], points[:, 2], color=color)
end

function draw!(tree::SphereTree; color="green")
  for (center, srad) in zip(tree.centers, tree.squared_radii)
    draw!(Sphere(center, sqrt(srad)), color=color)
  end
end

export draw!

end # module

max_stars_repo_name: krzentner/DistancePlan

max_issues_repo_name: krzentner/DistancePlan

max_forks_repo_name: krzentner/DistancePlan


text:
using GLVisualize, GeometryTypes, Reactive, GLAbstraction

if !isdefined(:runtests)
    window = glscreen()
    timesignal = loop(linspace(0f0, 1f0, 360))
end

description = """
Efficiently animated 3D unicode arrow field
"""

# let the visualization rotate later on
rotation = map(rotationmatrix_z, const_lift(*, timesignal, 2f0 * pi))
# create some random 3D vectors
vectors3d = rand(Vec3f0, 5, 5, 5)
# this is not the best way to use 2D sprites, but this will space them on
# a 3D grid and use the rotation from `vectors3d` and the length of them
# to look up the a color from the optional keyword argument `color_map`.
arrows = visualize(('➤', vectors3d), scale = Vec2f0(0.1), model = rotation)

_view(arrows, camera = :perspective)

if !isdefined(:runtests)
    renderloop(window)
end


text:
using Distributions
using FreqTables
using Gadfly
using Cairo

#This Julia script can be used to generate random (binomially)
#distributed codon usage bias data, which can then be compared to
#the actual sequences.
#It contains a function (getSample) to produce the random data.
#The cub can be adjusted in the function body
#So can the length and the noise around the cub.
#The script also contains a few lines to evaluate the slope

#############################################
#Calculates artificial (random) Entropy data for a
#probability "bas"
#Returns the data that can then be used to calulcated the
#temperature
#If basenorm is true, then entropy is calculated
#relative to equal opportunity model
#############################################
function getSample(bas=0.5,basenorm=true)
#tmp is the success probability
#rt is the length we consider
local tmp5,tmp3,tmp4,tbl,outp,v,d,based
 tmp5,tmp3,tmp4,tbl,outp=-1,-1,-1,-1,-1

v=0;
for i in 1:1

#first decide on the size N-> rt and the probablity p->tmp
rt = rand([5,6,7,8,9,10])
rt=5
#tmp = 0.5 + 0.1*randn()
tmp =bas + 0.*randn()

if (tmp > 1) tmp =1 end
if(tmp <0) tmp = 0 end

#Now define your distributions d is the biased one, based is the one with 0.5
d=Binomial(rt,tmp)

#if basenorm is true, then use as base probability the 0.5 prob otherwise the same as the input prob
if(basenorm)
based= Binomial(rt,0.5)
else
based= Binomial(rt,tmp)
end

#DELETEME DELTE NEXT LINE
#based= Binomial(rt,0.2)

#inner loop
#generate samples with the chosen probability and size

for j in 1:1

#generate the actual sample
tmp10=rand(d,10000)

#calculate the probability
a=round.(1000000000*pdf.(based,tmp10))/1000000000

if ( ((rt % 2) == 0) & (tmp10 ==  Int64(round(rt/2))  ))
#println("Twice")
v=cat(1,v,a)
else

if(basenorm)
v=cat(1,v,2*a)
else # don't mulyoply by 2
v=cat(1,v,a)
end#if

end#if

end#innerfor

#now produce the freqtable
#This freqtable has the frequency as col2
# and the probability as col1
deleteat!(v,1)
a= v
tbl=freqtable(a)
tmp4=convert(Array, tbl)
en=sum(tmp4)
rel=names(tbl)[1]
tmp4=tmp4/sum(tmp4)# sampled frequencies (empirical distribution)
#tmp6= log.(tmp4)
tmp6= tmp4
#tmp3=-log.(names(tbl)[1])
tmp3=names(tbl)[1]#theoretical distributions
#println(tbl)
println(tmp4)
println( exp(-en*sum(tmp4.*log.(tmp4./rel)) ))
println( sum(tmp4.*log.(tmp4./rel)) )
tmp5=[tmp3 tmp6]
if(outp != -1)
outp = [outp; tmp5]
else
outp = tmp5
end

v=0
#Now start the whole thing again, by choosing a new number and length
end#outerfor

return(outp)
end#function
#############################################
#############################################

#if you pass false to this function, then the probability is evaluated
#according to the same bias as the sample sequence.
#Otheriwise to the binomial with p=1/2
a= getSample(0.2,false)

#writedlm("/tmp/dta.txt",a)

writecsv("/tmp/dta.txt",a)
tmp=linreg(a[:,1],a[:,2])
print(tmp)
#myplot=plot(x=a[:,1],y=a[:,2],Geom.point,Geom.smooth(method=:lm))
#draw(PS("myplot.eps", 8inch, 6inch), myplot)

max_stars_repo_name: jkobject/CodonUsageBias

max_issues_repo_name: jkobject/CodonUsageBias

max_forks_repo_name: jkobject/CodonUsageBias


text:
"""
    symbolize(x)

Generalized surrogate function of `Symbol`
"""
symbolize(x::AbstractString) = Symbol(x)
symbolize(x::AbstractArray{T}) where {T<:AbstractString} = [symbolize(y) for y in x]

"""
    static_array(x)

Transform to static array
"""
function static_array(x::AbstractVector)
    y = MVector{length(x)}(collect(x))

    if x isa OffsetArray
        idx0 = eachindex(x) |> first
        idx1 = eachindex(x) |> last

        y = OffsetArray(y, idx0:idx1)
    end

    return y
end

function static_array(x::AbstractMatrix)
    y = MMatrix{size(x, 1),size(x, 2)}(collect(x))

    if x isa OffsetArray
        idx0 = axes(x, 1) |> first
        idx1 = axes(x, 1) |> last
        idy0 = axes(x, 2) |> first
        idy1 = axes(x, 2) |> last

        y = OffsetArray(y, idx0:idx1, idy0:idy1)
    end

    return y
end

function static_array(x::AbstractArray{T,3}) where T
    y = MArray{Tuple{size(x, 1),size(x, 2),size(x, 3)}}(collect(x))

    if x isa OffsetArray
        idx0 = axes(x, 1) |> first
        idx1 = axes(x, 1) |> last
        idy0 = axes(x, 2) |> first
        idy1 = axes(x, 2) |> last
        idz0 = axes(x, 3) |> first
        idz1 = axes(x, 3) |> last

        y = OffsetArray(y, idx0:idx1, idy0:idy1, idz0:idz1)
    end

    return y
end

function static_array(x::AbstractArray{T,4}) where T
    y = MArray{Tuple{size(x, 1),size(x, 2),size(x, 3),size(x, 4)}}(collect(x))

    if x isa OffsetArray
        ida0 = axes(x, 1) |> first
        ida1 = axes(x, 1) |> last
        idb0 = axes(x, 2) |> first
        idb1 = axes(x, 2) |> last
        idc0 = axes(x, 3) |> first
        idc1 = axes(x, 3) |> last
        idd0 = axes(x, 4) |> first
        idd1 = axes(x, 4) |> last

        y = OffsetArray(y, ida0:ida1, idb0:idb1, idc0:idc1, idd0:idd1)
    end

    return y
end

"""
    dynamic_array(x)

Transform to dynamic array
"""
dynamic_array(x) = x
dynamic_array(x::StaticArray) = Array(x)


text:
# test deterministic vs. stochastic
using Distributed;
addprocs(20);
@everywhere include("loadMod.jl");
@everywhere const GUROBI_ENV = Gurobi.Env();

caseList = [13,33,123];
NN = 1000;
Δt = 0.25;
N = 5;
iterMax = 20;
pathListDRaw = load("pathHist_1000.jld");
TList = [24,36,48,72,96];

for ci in 1:length(caseList)
    pathDictA = pathListDRaw["pathDict"][ci];
    detOut = Dict();
    stochOut = Dict();
    detCostDict = Dict();
    for T in TList
        τ = Int64(1/6*T);
        for j in procs()
            remotecall_fetch(readInData,j,ci,caseList,T,τ);
            #remotecall_fetch(readInData_old,j,T,ωSet0,10000,0);
        end

        # select a preset pathDict
        pathDict = pathDictA[T];
        solDet,costDet = exeDet(T, Δt, fData, bData, dData, pDistr, NN, pathDict);
        listDet = [costDet[i] for i in 1:NN];
        meanDet = mean(listDet);
        sigmaDet = std(listDet);
        println(round(meanDet,digits = 2)," ",
                round(meanDet - 1.96*sigmaDet,digits = 2)," ",
                round(meanDet + 1.96*sigmaDet,digits = 2));
        detOut[T] = [costDet,listDet,meanDet,sigmaDet];

        # train the stochastic programming strategy
        cutDictPG = preGen(T, Δt, N, iterMax);
        # cutDict,LBHist,UBHist,UBuHist,UBlHist,timeHist = solveMain(τ, T, Δt, N, true, false, 10, 10, cutDictPG);
        cutDict,LBHist,UBHist,UBuHist,UBlHist,timeHist = solveMain(T, Δt, N, false, false,
            max(Int64(round(500/N)),20), max(Int64(round(500/N)),20),cutDictPG);
        for j in procs()
            remotecall_fetch(cutIni,j,cutDict);
        end

        solSDDP, LBSDDP, costSDDP = exeForward(T, Δt, NN, false, pathDict);
        listSDDP = [costSDDP[i] for i in 1:NN];
        meanSDDP = mean(listSDDP);
        sigmaSDDP = std(listSDDP);
        println(round(meanSDDP,digits = 2)," ",
                round(meanSDDP - 1.96*sigmaSDDP,digits = 2)," ",
                round(meanSDDP + 1.96*sigmaSDDP,digits = 2));
        stochOut[T] = [LBSDDP, costSDDP, listSDDP, meanSDDP, sigmaSDDP, LBHist];

        detSol,detObj = detBuild(Δt, T, fData, bData, dData);
        detCostDict[T] = detObj;

        #println("Output: ",round(meanDet,1)," & ",round(meanSDDP,1)," & ", round(meanDet - meanSDDP,1)," & ", round((meanDet - meanSDDP)/meanDet*100,1));
        #save("detResults_old.jld", "detOut", detOut, "stochOut", stochOut, "nomOut", detCostDict);
        save("detResults_$(ci).jld", "detOut", detOut, "stochOut", stochOut, "nomOut", detCostDict);
    end
end

# for T in TList
#     println(T," & ",round(detOut[T][3],2)," & ",round(stochOut[T][4],2)," & ",round((detOut[T][3] - stochOut[T][4])/detOut[T][3]*100,2)," & ");
# end

# pathListData = pmap(i -> simuPath(τ,T,pDistr), 1:NN);
# pathDict = Dict();
# for i in 1:length(pathListData)
#     pathDict[i] = pathListData[i];
# end
# solDet,costDet = exeDet(τ, T, Δt, fData, bData, dData, pDistr, NN, pathDict);
# listDet = [costDet[i] for i in 1:NN];
# meanDet = mean(listDet);
# sigmaDet = std(listDet);
# println(round(meanDet,2)," ",round(meanDet - 1.96*sigmaDet,2)," ",round(meanDet + 1.96*sigmaDet,2));
#
# cutDict,LBHist,UBHist,UBuHist,UBlHist = solveMain(τ, T, Δt, fData, pDistr, bData, dData, N, 20, 20);
# solSDDP, LBSDDP, costSDDP = exeForward(τ, T, Δt, fData, bData, dData, pDistr, NN, cutDict, pathDict);
# listSDDP = [costSDDP[i] for i in 1:NN];
# meanSDDP = mean(listSDDP);
# sigmaSDDP = std(listSDDP);
# println(round(meanSDDP,2)," ",round(meanSDDP - 1.96*sigmaSDDP,2)," ",round(meanSDDP + 1.96*sigmaSDDP,2));
#
# outputData = [solDet,costDet,meanDet,sigmaDet,
#     cutDict,LBHist,UBHist,UBuHist,UBlHist,solSDDP,LBSDDP,costSDDP,meanSDDP,sigmaSDDP];
# save("detOut.jld","data",outputData);


text:
# This file was generated, do not modify it. # hide
n = 8
k = 5

max_stars_repo_name: cscherrer/InformativePrior

max_issues_repo_name: cscherrer/InformativePrior

max_forks_repo_name: cscherrer/InformativePrior


text:
using InteractiveUtils: code_llvm
using Test

include("../benchmark/map.jl")

"""
    llvm_ir(f, args) :: String

Get LLVM IR of `f(args...)` as a string.
"""
llvm_ir(f, args) = sprint(code_llvm, f, Base.typesof(args...))

nmatches(r, s) = count(_ -> true, eachmatch(r, s))

twice(x) = 2x

@testset "setindex!! (@inbounds)" begin
    @test mapbb(twice, [1:3;], Val(false)) == [2, 4, 6]

    ir_noinbounds = llvm_ir(mapbb, (twice, Float64[], Val(false)))
    if Base.JLOptions().check_bounds == 2  # --check-bounds=no
        @test nmatches(r"fmul <4 x double>", ir_noinbounds) >= 4
    else
        @test nmatches(r"fmul <4 x double>", ir_noinbounds) == 0
    end
    ir_inbounds = llvm_ir(mapbb, (twice, Float64[], Val(true)))
    @test nmatches(r"fmul <4 x double>", ir_inbounds) >= 4
end


text:
import DASSL  # note this requires current master (Feb 2015)

function dassl_wrapper{N,T}(tr::TestRun{N,T})
    # Wraps DASSL.dasslSolve
    tc = tr.tc
    so = tr.solver

    ###
    # 0) Wrap tc.fn!, tc.jac!, tc.mass! if necessary
    if hasmass(tc)
        fn(t,y,dydt) = (m = tc.mass!(); tc.mass!(t,y,m); out = tc.fn!(); tc.fn!(t,y,out); out-m*dydt)
    else
        fn(t,y,dydt) = (out = tc.fn!(); tc.fn!(t,y,out); out-dydt)
    end
    # TODO once DASSL updates interface:
    if hasjacobian(tc)
        Fy(t,y,dydt) = (out = tc.jac!(); tc.jac!(t,y,out); out)
        if hasmass(tc)
            mass = nonmod_mass(tc)
            Fdy(t,y,dydt) = -mass(t,y)
        else
            if issparse(tc)
                Fdy(t,y,dydt) = -speye(T, tc.dof)
            else
                Fdy(t,y,dydt) = -eye(T, tc.dof)
            end
        end
        jac(t,y,dydt,a) = Fy(t,y,dydt) + a*Fdy(t,y,dydt)
    end
    ###
    # 1) Make call signature
    args = (fn, tc.ic, tc.tspan)
    if hasjacobian(tc)
        kwargs  = ((:reltol, tr.reltol), (:abstol, tr.abstol), (:jacobian, jac))
    else
        kwargs  = ((:reltol, tr.reltol), (:abstol, tr.abstol))
    end

    ###
    # 2) Call solver, if it does not succeed throw an error (if that
    # is not done anyway)
    (t,y) = so.solverfn(args...; kwargs...)
    # (probably no need to modify this section)

    ###
    # 3) Transform output to conform to standard:
    # tend -- end time reached
    # yend -- solution at tend
    # stats -- statistics, if available: (steps_total,steps_accepted, fn_evals, jac_evals, linear_solves)
    #                         otherwise  (-1, -1, -1, -1, -1)
    tend = t[end]
    yend = y[end]
    stats = (-1, -1, -1, -1, -1)
    return tend, yend, stats
end

daeindex = 1
adaptive = true
dassl = Solver{:im}(DASSL.dasslSolve, DASSL, dassl_wrapper, stiff, adaptive, daeindex, explicit_mass_eq)

DASSLsolvers = Dict{Any,Solver}()
DASSLsolvers[DASSL.dasslSolve] = dassl

merge!(allsolvers, DASSLsolvers)


text:
"""
## Extract epoch data from SSR

#### Arguments
* `a`: A SSR object
* `valid_triggers`: Trigger numbers that are considered valid ([1,2])
* `remove_first`: Remove the first n triggers (0)
* `remove_last`: Remove the last n triggers (0)

#### Example

```julia
epochs = extract_epochs(SSR, valid_triggers=[1,2])
```
"""
function extract_epochs(
    a::SSR;
    valid_triggers::Union{AbstractArray,Int} = [1, 2],
    remove_first::Int = 0,
    remove_last::Int = 0,
    kwargs...,
)

    merge!(
        a.processing,
        Dict(
            "epochs" => extract_epochs(
                a.data,
                a.triggers,
                valid_triggers,
                remove_first,
                remove_last,
            ),
        ),
    )

    return a
end

function epoch_rejection(a::SSR; retain_percentage::Number = 0.95, kwargs...)

    a.processing["epochs"] =
        epoch_rejection(a.processing["epochs"], retain_percentage; kwargs...)

    return a
end

function create_sweeps(a::SSR; epochsPerSweep::Int = 64, kwargs...)

    if epochsPerSweep > size(a.processing["epochs"], 2)
        error("Sweep length is longer than number of epochs will allow")
    end

    merge!(
        a.processing,
        Dict("sweeps" => create_sweeps(a.processing["epochs"], epochsPerSweep)),
    )

    return a
end

#######################################
#
# Add triggers for more epochs
#
#######################################

function add_triggers(a::SSR; kwargs...)

    @debug("Adding triggers to reduce SSR. Using SSR modulation frequency")

    add_triggers(a, modulationrate(a); kwargs...)
end

function add_triggers(a::SSR, mod_freq::Number; kwargs...)

    @debug("Adding triggers to reduce SSR. Using $(mod_freq)Hz")

    epochIndex =
        DataFrame(Code = a.triggers["Code"], Index = [Int(i) for i in a.triggers["Index"]])
    epochIndex[:Code] = epochIndex[:Code] .- 252

    add_triggers(a, mod_freq, epochIndex; kwargs...)
end

function add_triggers(
    a::SSR,
    mod_freq::Number,
    epochIndex;
    cycle_per_epoch::Int = 1,
    kwargs...,
)

    @info(
        "Adding triggers to reduce SSR. Reducing $(mod_freq)Hz to $cycle_per_epoch cycle(s)."
    )

    # Existing epochs
    existing_epoch_length = median(diff(epochIndex[:Index]))     # samples
    existing_epoch_length_s = existing_epoch_length / samplingrate(a)
    @debug("Existing epoch length: $(existing_epoch_length_s)s")

    # New epochs
    new_epoch_length_s = cycle_per_epoch / mod_freq
    new_epochs_num = round.(existing_epoch_length_s / new_epoch_length_s) - 2
    new_epoch_times = collect(1:new_epochs_num) * new_epoch_length_s
    new_epoch_indx = [0; round.(new_epoch_times * samplingrate(a))]
    @debug("New epoch length = $new_epoch_length_s")
    @debug("New # epochs     = $new_epochs_num")

    # Place new epoch indices
    @debug("Was $(length(epochIndex[:Index])) indices")
    new_indx = epochIndex[:Index][1:end-1] .+ new_epoch_indx'
    new_indx = reshape(new_indx', length(new_indx), 1)[1:end-1]
    @debug("Now $(length(new_indx)) indices")

    # Place in dict
    new_code = round.(Int, ones(1, length(new_indx))) .+ 252
    a.triggers = Dict(
        "Index" => vec((new_indx)'),
        "Code" => vec(new_code),
        "Duration" => ones(length(new_code), 1)',
    )
    #TODO Possible the trigger duration of one is not long enough

    return a
end

#######################################
#
# Rejection channels
#
#######################################

function channel_rejection(
    a::SSR;
    threshold_abs::Number = 1000,
    threshold_var::Number = 2,
    kwargs...,
)

    # Run on epochs if available, else run on the raw data
    if haskey(a.processing, "epochs")
        data = reshape(
            a.processing["epochs"],
            size(a.processing["epochs"], 1) * size(a.processing["epochs"], 2),
            size(a.processing["epochs"], 3),
        )
    else
        data = a.data
    end

    valid = vec(channel_rejection(data, threshold_abs, threshold_var))

    @info(
        "Rejected $(sum(.!valid)) channels $(join(channelnames(a)[findall(.!valid)], " "))"
    )

    remove_channel!(a, channelnames(a)[.!valid])

    return a
end


text:
using Base.Test
using PETSc2

@testset "  ---Checking Petsc data types---" begin

  @test ( PetscScalar )== Float64
  @test ( PetscReal )== Float64
  @test ( PetscInt )== Int32

end


text:
using EVQueues, Plots

lambda=120.0;
mu=1.0;
gamma=1.0;
C=60.0;

Tfinal=1000.0;

sim = ev_edf(lambda,mu,gamma,Tfinal,C,snapshots=[Tfinal])
compute_statistics!(sim)

#TIme plot of occupation and power
p=plot(sim.timetrace)
display(p)

#CDF of departure attained workloads
p=servicecdf(sim.EVs)
display(p)

#State space of the last snapshot
snap = sim.snapshots[end];

p=stateplot(snap.charging)
display(p)

#fairness measure over time

t=collect(0:10:Tfinal);
h=20.0;
J=compute_fairness(sim,t,h)

p = plot(   xlabel = "Time",
            ylabel = "J",
            title = "Fairness measure",
            ylim = (0,1)
            )

plot!(p,t,J,lw=2,legend=:none)

display(p)


text:
using Gadfly

plot(layer(y=20*rand(20), Geom.bar), layer(x=20*rand(20), Geom.line))