text:
__precompile__(true)

"""
Syntactic sugar for BLAS polynomials.
"""
module SugarBLAS

export ᵀ
export  @blas!

include("Match/Match.jl")
using .Match

struct ᵀ end

"""
Negate a number, symbol or expression
"""
neg(ast::Number) = -ast
function neg(ast::Union{Symbol, Expr})
    if @match(ast, -ast) | (ast == 0)
        ast
    else
        Expr(:call, :(-), ast)
    end
end

"""
Determine whether it is a substraction or not
"""
substracts(expr) = false
substracts(expr::Expr) = (expr.head == :call) & (expr.args[1] == :-)

"""
Make dictionary containing the kwargs contents.
"""
function kwargs_to_dict(kwargs::Tuple)
    dict = Dict()
    for kw in kwargs
        dict[kw.args[1]] = kw.args[2]
    end
    dict
end

"""
Wrap a expr with an expression.
"""
wrap(expr::Symbol) = QuoteNode(expr)
function wrap(expr::Expr)
    head = QuoteNode(expr.head)
    func = string(expr.args[1])
    :(Expr($head, Meta.parse($func), $(expr.args[2:end]...)))
end

"""
Expand mixed assignment operators.
"""
function expand(expr::Expr)
    @match(expr, A += B) && return :($A = $A + $B)
    @match(expr, A -= B) && return :($A = $A - $B)
    expr
end

"""
Execute escaped expr.
"""
macro call(expr::Expr)
    esc(:(esc($(wrap(expr)))))
end

#Changes made in Julia parser in v"0.6.0-dev.2613".
#JuliaLang Issue: https://github.com/JuliaLang/julia/pull/20327
"""
Transforms the custom case expression to a string representing the equivalent if-then-else block of code.
"""
construct_case_statement(lines::Vector) = construct_case_statement(lines, Val{VERSION>=v"0.6.0-dev.2613"})
function construct_case_statement(lines::Vector, ::Type{Val{true}})
  failproof(s) = s
  failproof(s::Char) = string("'",s,"'")
  line = lines[1]
  exec = "if $(line.args[2])\n$(failproof(line.args[3]))\n"
  for line in lines[2:end-1]
      (line.head == :line) && continue
      (line.head == :call && line.args[1] == :(=>)) || error("Each condition must be followed by `=>`")
      exec *= "elseif $(line.args[2])\n$(failproof(line.args[3]))\n"
  end
  line = lines[end]
  exec *= (line.args[2] == :otherwise) && ("else\n$(failproof(line.args[3]))\n")
  exec *= "end"
end
function construct_case_statement(lines::Vector, ::Type{Val{false}})
  failproof(s) = s
  failproof(s::Char) = string("'",s,"'")
  line = lines[1]
  exec = "if $(line.args[1])\n$(failproof(line.args[2]))\n"
  for line in lines[2:end-1]
      (line.head == :line) && continue
      line.head == :(=>) || error("Each condition must be followed by `=>`")
      exec *= "elseif $(line.args[1])\n$(failproof(line.args[2]))\n"
  end
  line = lines[end]
  exec *= (line.args[1] == :otherwise) && ("else\n$(failproof(line.args[2]))\n")
  exec *= "end"
end

"""
Sugar for if-then-else expression. Beautiful for one liners.
"""
macro case(expr::Expr)
    (expr.head == :block) || error("@case statement must be followed by `begin ... end`")
    lines = filter(line -> !is_metadata(line), expr.args)
    exec = construct_case_statement(lines)
    esc(Meta.parse(exec))
end

"""
Filter metadata lines that often appear in Julia expressions
"""
is_metadata(line::Expr) = line.head == :line #Needed for versions older than v0.7-DEV
is_metadata(::LineNumberNode) = true

###############
# BLAS macros #
###############

"""
    @blas!(expr)

Transform expr to most specific BLAS function. Resulting expression will be
Base.copy!, Base.scale! or Base.LinAlg.axpy!.

**Polynomials**

*copy!*:

- `X = Y`

*scale!*:

- `X *= a`
- `X = a*X`

*axpy!*:

- `Y ±= X`
- `Y ±= a*X`
"""
#Must be ordered from most to least especific formulas
macro blas!(expr::Expr)
    unkeyword!(expr)
    expr = expand(expr)
    @case begin
        @match(expr, X *= a)        => @call scale!(a,X)
        @match(expr, X = a*X)       => @call scale!(a,X)
        @match(expr, Y = Y - a*X)   => @call Base.LinAlg.axpy!(neg(a),X,Y)
        @match(expr, Y = Y - X)     => @call Base.LinAlg.axpy!(-1.0,X,Y)
        @match(expr, Y = a*X + Y)   => @call Base.LinAlg.axpy!(a,X,Y)
        @match(expr, Y = X + Y)     => @call Base.LinAlg.axpy!(1.0,X,Y)
        @match(expr, X = Y)         => @call copy!(X, Y)
        otherwise                   => error("No match found")
    end
end

"""
    @copy!(expr)

Copy all elements from collection `Y` to array `X`. Return `X`.

**Polynomials**

- `X = Y`
"""
macro copy!(expr::Expr)
    unkeyword!(expr)
    @case begin
        @match(expr, X = Y) => @call copy!(X,Y)
        otherwise           => error("No match found")
    end
end

"""
    @scale!(expr)

Scale an array `X` by a scalar `a` overwriting `X` in-place.

**Polynomials**

- `X *= a`
- `X = a*X`
"""
macro scale!(expr::Expr)
    unkeyword!(expr)
    @case begin
        @match(expr, X *= a)    => @call scale!(a,X)
        @match(expr, X = a*X)   => @call scale!(a,X)
        otherwise               => error("No match found")
    end
end

"""
    @axpy!(expr)

Overwrite `Y` with `a*X + Y`. Return `Y`.

**Polynomials**

- `Y ±= X`
- `Y ±= a*X`
"""
macro axpy!(expr::Expr)
    unkeyword!(expr)
    expr = expand(expr)
    @case begin
        @match(expr, Y = Y - a*X)   => @call(Base.LinAlg.axpy!(neg(a),X,Y))
        @match(expr, Y = Y - X)     => @call Base.LinAlg.axpy!(-1.0,X,Y)
        @match(expr, Y = a*X + Y)   => @call Base.LinAlg.axpy!(a,X,Y)
        @match(expr, Y = X + Y)     => @call Base.LinAlg.axpy!(1.0,X,Y)
        otherwise                   => error("No match found")
    end
end

"""
    @ger!(expr)

Rank-1 update of the matrix `A` with vectors `x` and `y` as `alpha*x*y' + A`.

**Polynomials**

- `A ±= alpha*x*y'`
"""
macro ger!(expr::Expr)
    unkeyword!(expr)
    expr = expand(expr)
    f = @case begin
        @match(expr, A = alpha*x*y' + A)    => identity
        @match(expr, A = A - alpha*x*y')    => neg
        otherwise                           => error("No match found")
    end
    @call Base.LinAlg.BLAS.ger!(f(alpha),x,y,A)
end

"""
    @syr!(expr)

Rank-1 update of the symmetric matrix `A` with vector `x` as `alpha*x*xᵀ + A`.
When left side has `A['U']` the upper triangle of `A` is updated (`'L'` for lower
triangle). Return `A`.

**Polynomials**

- `A[uplo] ±= alpha*x*xᵀ`
"""
macro syr!(expr::Expr)
    unkeyword!(expr)
    expr = expand(expr)
    @match(expr, A[uplo] = right) || error("No match found")
    f = @case begin
        @match(right, alpha*x*xᵀ + Y)  => identity
        @match(right, Y - alpha*x*xᵀ)  => neg
        otherwise                       => error("No match found")
    end
    (@match(Y, Y[uplo]) && (Y == A)) || (Y == A) || error("No match found")
    @call Base.LinAlg.BLAS.syr!(uplo,f(alpha),x,A)
end

"""
    @syrk!(expr)

Return either the upper triangle or the lower triangle, depending on
(`'U'` or `'L'`), of `alpha*A*(A)ᵀ` or `alpha*(A)ᵀ*A`.

**Polynomials**

- `alpha*A*(A)ᵀ uplo=ul`
- `alpha*(A)ᵀ*A uplo=ul`
"""
macro syrk(expr::Expr, kwargs...)
    kwargs = kwargs_to_dict(kwargs)
    uplo = kwargs[:uplo]
    f = @case begin
        @match(expr, alpha*X*Y) => identity
        otherwise               => error("No match found")
    end
    trans = @case begin
        @match(X, (A)ᵀ) && (Y == A)  => 'T'
        @match(Y, (A)ᵀ) && (X == A)  => 'N'
        otherwise                   => error("No match found")
    end
    @call Base.LinAlg.BLAS.syrk(uplo,trans,f(alpha),A)
end

"""
    @syrk!(expr)

Rank-k update of the symmetric matrix `C` as `alpha*A*(A)ᵀ + beta*C` or
`alpha*(A)ᵀ*A + beta*C`. When the left hand side is`C['U']` the upper triangle of `C`
is updated (`'L'` for lower triangle). Return `C`.

**Polynomials**

- `C[uplo] ±= alpha*A*(A)ᵀ`
- `C[uplo] = beta*C ± alpha*(A)ᵀ*A`
"""
macro syrk!(expr::Expr)
    unkeyword!(expr)
    expr = expand(expr)
    @match(expr, C[uplo] = right) || error("No match found")
    f = @case begin
        @match(right, alpha*X*Y + D)    => identity
        @match(right, D - alpha*X*Y)    => neg
        otherwise                       => error("No match found")
    end
    trans = @case begin
        @match(X, (A)ᵀ) && (Y == A)  => 'T'
        @match(Y, (A)ᵀ) && (X == A)  => 'N'
        otherwise                   => error("No match found")
    end
    @match(D, beta*D) || (beta = 1.0)
    (@match(D, D[uplo]) && (C == D)) || (C == D) || error("No match found")
    @call Base.LinAlg.BLAS.syrk!(uplo,trans,f(alpha),A,beta,C)
end

"""
    @her!(expr)

Methods for complex arrays only. Rank-1 update of the Hermitian matrix `A`
with vector `x` as `alpha*x*x' + A`. Whenthe left hand side is `A['U']`
the upper triangle of `A` is updated (`'L'` for lower triangle). Return `A`.

**Polynomials**

- `A[uplo] ±= alpha*x*x'`
"""
macro her!(expr::Expr)
    unkeyword!(expr)
    expr = expand(expr)
    @match(expr, A[uplo] = right) || error("No match found")
    f = @case begin
        @match(right, alpha*x*x' + Y)   => identity
        @match(right, Y - alpha*x*x')   => neg
        otherwise                       => error("No match found")
    end
    (@match(Y, Y[uplo]) && (Y == A)) || (Y == A) || error("No match found")
    @call Base.LinAlg.BLAS.her!(uplo,f(alpha),x,A)
end

"""
    @herk(expr)

Methods for complex arrays only. Returns either the upper triangle or the
lower triangle, according to uplo ('U' or 'L'), of alpha*A*A' or alpha*A'*A,
according to trans ('N' or 'T').

**Polynomials**

- `alpha*A*A' uplo=ul`
- `alpha*A'*A uplo=ul`
"""
macro herk(expr::Expr, kwargs...)
    kwargs = kwargs_to_dict(kwargs)
    uplo = kwargs[:uplo]
    @match(expr, alpha*X*Y) || error("No match found")
    trans = @case begin
        @match(X, A') && (Y == A)   =>  'T'
        @match(Y, A') && (X == A)   =>  'N'
        otherwise                   =>  error("No match found")
    end
    @call Base.LinAlg.BLAS.herk(uplo,trans,alpha,A)
end

"""
    @herk!(expr)

Methods for complex arrays only. Rank-k update of the Hermitian matrix `C` as
`alpha*A*A' + beta*C` or `alpha*A'*A + beta*C`. When the left hand side is `C['U']`
the upper triangle of `C` is updated (`'L'` for lower triangle). Return `C`.

**Polynomials**

- `C[uplo] ±= alpha*A*A'`
- `C[uplo] = beta*C ± alpha*A'*A`
"""
macro herk!(expr::Expr)
    unkeyword!(expr)
    expr = expand(expr)
    @match(expr, C[uplo] = right) || error("No match found")
    f = @case begin
        @match(right, alpha*X*Y + D)    => identity
        @match(right, D - alpha*X*Y)    => neg
        otherwise                       => error("No match found")
    end
    trans = @case begin
        @match(X, A') && (Y == A)   =>  'T'
        @match(Y, A') && (X == A)   =>  'N'
        otherwise                   =>  error("No match found")
    end
    @match(D, beta*D) || (beta = 1.0)
    (@match(D, D[crap]) && (C == D)) || (C == D) || error("No match found")
    @call Base.LinAlg.BLAS.herk!(uplo,trans,f(alpha),A,beta,C)
end

"""
    @gbmv(expr)

Return `alpha*A*x` or `alpha*A'*x`. The matrix `A` is a general band matrix
of dimension `m` by `size(A,2)` with `kl` sub-diagonals and `ku` super-diagonals.

**Polynomials**

- `alpha*A[kl:ku,h=m]*x`
- `alpha*A[h=m,kl:ku]'*x`
"""
macro gbmv(expr::Expr)
    @match(expr, alpha*Y*x) || error("No match found")
    trans = @match(Y, Y') ? 'T' : 'N'
    @match(Y, A[kl:ku,h=m])
    @call Base.LinAlg.BLAS.gbmv(trans,m,neg(kl),ku,alpha,A,x)
end

"""
    @gbmv!(expr)

Update vector `y` as `alpha*A*x + beta*y` or `alpha*A'*x + beta*y`.
The matrix `A` is a general band matrix of dimension `m` by `size(A,2)` with
`kl` sub-diagonals and `ku` super-diagonals. Return the updated `y`.

**Polynomials**

- `y ±= alpha*A[kl:ku,h=m]*x`
- `y = beta*y ± alpha*A[h=m,kl:ku]'*x`
"""
macro gbmv!(expr::Expr)
    unkeyword!(expr)
    expr = expand(expr)
    @match(expr, y = right) || error("No match found")
    f = @case begin
        @match(right, alpha*Y*x + w)    => identity
        @match(right, w - alpha*Y*x)    => neg
        otherwise                       => error("No match found")
    end
    trans = @match(Y, Y') ? 'T' : 'N'
    @match(Y, A[kl:ku,h=m])
    @match(w, beta*w) || (beta = 1.0)
    (y == w) || error("No match found")
    @call Base.LinAlg.BLAS.gbmv!(trans,m,neg(kl),ku,f(alpha),A,x,beta,y)
end

"""
    @sbmv(expr)

Return `alpha*A*x` where `A` is a symmetric band matrix of order `size(A,2)` with
`k` super-diagonals stored in the argument `A`.

**Polynomials**

- `A[0:k,uplo]*xv`
- `alpha*A[0:k,uplo]*x`
"""
macro sbmv(expr::Expr)
    @case begin
        @match(expr, alpha*A[0:k,uplo]*x)   => @call Base.LinAlg.BLAS.sbmv(uplo,k,alpha,A,x)
        @match(expr, A[0:k,uplo]*x)         => @call Base.LinAlg.BLAS.sbmv(uplo,k,A,x)
        otherwise                           => error("No match found")
    end
end

"""
    @sbmv!(expr)

Update vector `y` as `alpha*A*x + beta*y` where `A` is a a symmetric band matrix
of order `size(A,2)` with `k` super-diagonals stored in the argument `A`. If
`A[...,'U']` is used multiplication is done with `A`'s upper triangle, `L` is for the
lower triangle. Return updated `y`.

**Polynomials**

- `y ±= alpha*A[0:k,uplo]*x`
- `y = beta*y ± alpha*A[0:k,uplo]*x`
"""
macro sbmv!(expr::Expr)
    unkeyword!(expr)
    expr = expand(expr)
    @match(expr, y = right) || error("No match found")
    f = @case begin
        @match(right, alpha*A[0:k,uplo]*x + w)  => identity
        @match(right, w - alpha*A[0:k,uplo]*x)  => neg
        otherwise                               => error("No match found")
    end
    @match(w, beta*w) || (beta = 1.0)
    (@match(w, w[crap]) && (y == w)) || (y == w) || error("No match found")
    @call Base.LinAlg.BLAS.sbmv!(uplo,k,f(alpha),A,x,beta,y)
end

"""
    @gemm(expr)

Return `alpha*A*B`, `alpha*A'*B`, `alpha*A*B'` or `alpha*A'*B'`.

**Polynomials**

- `A*B`
- `A'*B`
- `A*B'`
- `A'*B'`
- `alpha*A*B`
- `alpha*A'*B`
- `alpha*A*B'`
- `alpha*A'*B'`
"""
macro gemm(expr::Expr)
    if @match(expr, alpha*A*B)
        tA = @match(A, A') ? 'T' : 'N'
        tB = @match(B, B') ? 'T' : 'N'
        @call Base.LinAlg.BLAS.gemm(tA,tB,alpha,A,B)
    elseif @match(expr, A*B)
        tA = @match(A, A') ? 'T' : 'N'
        tB = @match(B, B') ? 'T' : 'N'
        @call Base.LinAlg.BLAS.gemm(tA,tB,A,B)
    else
        error("No match found")
    end
end

"""
    @gemm!(expr)

Update `C` as `alpha*A*B + beta*C` or the other three variants according to the
combination of transposes of `A` and `B`. Return updated C.

**Polynomials**

- `C ±= alpha*A*B`
- `C ±= alpha*A'*B`
- `C ±= alpha*A*B'`
- `C ±= alpha*A'*B'`
- `C = beta*C ± alpha*A*B`
- `C = beta*C ± alpha*A'*B`
- `C = beta*C ± alpha*A*B'`
- `C = beta*C ± alpha*A'*B'`
"""
macro gemm!(expr::Expr)
    unkeyword!(expr)
    expr = expand(expr)
    @match(expr, C = right) || error("No match found")
    f = @case begin
        @match(right, alpha*A*B + D)    => identity
        @match(right, D - alpha*A*B)    => neg
        otherwise                       => error("No match found")
    end
    tA = @match(A, A') ? 'T' : 'N'
    tB = @match(B, B') ? 'T' : 'N'
    @match(D, beta*D) || (beta = 1.0)
    (C == D) || error("No match found")
    @call Base.LinAlg.BLAS.gemm!(tA,tB,f(alpha),A,B,beta,C)
end

"""
    @gemv(expr)

Return `alpha*A*x` or `alpha*A'*x`.

**Polynomials**

- `A*x`
- `A'*x`
- `alpha*A*x`
- `alpha*A'*x`
"""
macro gemv(expr::Expr)
    if @match(expr, alpha*A*x)
        tA = @match(A, A') ? 'T' : 'N'
        @call Base.LinAlg.BLAS.gemv(tA,alpha,A,x)
    elseif @match(expr, A*x)
        tA = @match(A, A') ? 'T' : 'N'
        @call Base.LinAlg.BLAS.gemv(tA,A,x)
    else
        error("No match found")
    end
end

"""
    @gemv!(expr)

Update the vector `y` as `alpha*A*x + beta*y` or `alpha*A'*x + beta*y`.
Return updated `y`.

**Polynomials**

- `y ±= alpha*A*x`
- `y ±= alpha*A'*x`
- `y = beta*y ± alpha*A*x`
- `y = beta*y ± alpha*A'*x`
"""
macro gemv!(expr::Expr)
    unkeyword!(expr)
    expr = expand(expr)
    @match(expr, y = right) || error("No match found")
    f = @case begin
        @match(right, alpha*A*x + w)    => identity
        @match(right, w - alpha*A*x)    => neg
        otherwise                       => error("No match found")
    end
    tA = @match(A, A') ? 'T' : 'N'
    @match(w, beta*w) || (beta = 1.0)
    (y == w) || error("No match found")
    @call Base.LinAlg.BLAS.gemv!(tA,f(alpha),A,x,beta,y)
end

"""
    @symm(expr)

Return `alpha*A*B` or `alpha*B*A` according to `"symm"`. `A` is assumed to be
symmetric. Only the `uplo` triangle of `A` is used (`'L'` for lower and `'U'` for upper).

**Polynomials**

- `A["symm", uplo]*B`
- `A*B["symm", uplo]`
- `alpha*A["symm", uplo]*B `
- `alpha*A*B["symm", uplo]`
"""
macro symm(expr::Expr)
    if @match(expr, alpha*A*B)
        side = @case begin
            @match(A, A["symm", uplo])  => 'L'
            @match(B, B["symm", uplo])  => 'R'
            otherwise                   => error("No match found")
        end
        @call Base.LinAlg.BLAS.symm(side,uplo,alpha,A,B)
    elseif @match(expr, A*B)
        side = @case begin
            @match(A, A["symm", uplo])  => 'L'
            @match(B, B["symm", uplo])  => 'R'
            otherwise                   => error("No match found")
        end
        @call Base.LinAlg.BLAS.symm(side,uplo,A,B)
    else
        error("No match found")
    end
end

"""
    @symm!(expr)

Update `C` as `alpha*A*B + beta*C` or `alpha*B*A + beta*C` according to `"symm"`.
`A` is assumed to be symmetric. Only the `uplo` triangle of `A` is used
(`'L'` for lower and `'U'` for upper). Return updated `C`.

**Polynomials**

- `C = alpha*A["symm",uplo]*B`
- `C = alpha*A*B["symm",uplo]`
- `C = beta*C ± alpha*A["symm",uplo]*B`
- `C = beta*C ± alpha*A*B["symm",uplo]`
"""
macro symm!(expr::Expr)
    unkeyword!(expr)
    expr = expand(expr)
    @match(expr, C = right) || error("No match found")
    f = @case begin
        @match(right, alpha*A*B + D)    => identity
        @match(right, D - alpha*A*B)    => neg
        otherwise                       => error("No match found")
    end
    side = @case begin
        @match(A, A["symm", uplo])   => 'L'
        @match(B, B["symm", uplo])   => 'R'
        otherwise                    => error("No match found")
    end
    @match(D, beta*D) || (beta = 1.0)
    (@match(D, D[crap]) && (C == D)) || (C == D) || error("No match found")
    @call Base.LinAlg.BLAS.symm!(side,uplo,f(alpha),A,B,beta,C)
end

"""
    @symv(expr)

Return `alpha*A*x`. `A` is assumed to be symmetric. Only the `uplo` triangle of `A`
is used (`'L'` for lower and `'U'` for upper).

**Polynomials**

- `A[uplo]*x`
- `alpha*A[uplo]*x`
"""
macro symv(expr::Expr)
    @case begin
        @match(expr, alpha*A[uplo]*x)   => @call Base.LinAlg.BLAS.symv(uplo,alpha,A,x)
        @match(expr, A[uplo]*x)         => @call Base.LinAlg.BLAS.symv(uplo,A,x)
        otherwise                       => error("No match found")
    end
end

"""
    @symv!(expr)

Update the vector `y` as `alpha*A*x + beta*y`. `A` is assumed to be symmetric.
Only the `uplo` triangle of `A` is used (`'L'` for lower and `'U'` for upper).
Return updated y.

**Polynomials**

- `y ±= alpha*A[uplo]*x`
- `y = beta*y ± alpha*A[uplo]*x`
"""
macro symv!(expr::Expr)
    unkeyword!(expr)
    expr = expand(expr)
    @match(expr, y = right) || error("No match found")
    f = @case begin
        @match(right, alpha*A[uplo]*x + w)  => identity
        @match(right, w - alpha*A[uplo]*x)  => neg
        otherwise                           => error("No match found")
    end
    @match(w, beta*w) || (beta = 1.0)
    (@match(w, w[crap]) && (y == w)) || (y == w) || error("No match found")
    @call Base.LinAlg.BLAS.symv!(uplo,f(alpha),A,x,beta,y)
end

end


text:
struct PatternSet{T<:AbstractPattern}
    patterns::Vector{T}
end

function patternset(data::DataFrame, w::Int64, s::Int64;
        flagkey=:facies, propkeys=[:rho, :vp], condense=false)

    @assert all(.!isnan.(convert(Matrix, data[propkeys]))) "NaNs are not allowed"
    @assert flagkey in names(data) "flag key not found"

    # depth sampling rate
    δh = data[:depth][2]-data[:depth][1]

    # vector of unique sequences
    sequences = Dict()

    for i in 1:s:size(data, 1)-w
        labels = data[flagkey][i:i+w-1]
        props = Dict(key => data[key][i:i+w-1] for key in propkeys)
        props[:dh] = δh * ones(w)

        if condense
            seq, counts = seqcounts(labels)
            props[:thickness] = δh * counts
        else
            seq = labels
        end

        if seq in keys(sequences)
            push!(sequences[seq], props)
        else
            sequences[seq] = [props]
        end
    end

    patterns = [EmpiricalPattern(seq, props, δh) for (seq, props) in sequences]
    counts = realcount.(patterns)

    # sort patterns in descending order of realization counts
    permute!(patterns, sortperm(counts, rev=true))

    PatternSet(patterns)
end

patterns(pset::PatternSet) = pset.patterns
Base.length(pset::PatternSet) = length(pset.patterns)
Base.getindex(pset::PatternSet, i) = getindex(pset.patterns, i)
Base.iterate(pset::PatternSet) = iterate(pset.patterns)
Base.iterate(pset::PatternSet, i) = iterate(pset.patterns, i)
Base.lastindex(pset::PatternSet) = lastindex(pset.patterns)

function sample(pset::PatternSet; weighted=false)
    sample(pset.patterns) # TODO: weighted by frequency
end

sample(pset::PatternSet, n::Int64) = sample(pset.patterns, n)

function Base.show(io::IO, p::PatternSet)
    count = length(p.patterns)
    # numclasses = length(p.distributions)
    print(io, "PatternSet{$(eltype(p.patterns))}($count)")
end

ext: jl


max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
"""
    CurveWarpingSpace(knots)

Space of warpings of the interval [0, 1] on knots `knots`.
"""
struct CurveWarpingSpace{
    TK<:AbstractVector,
    TIM<:Interpolations.MonotonicInterpolationType,
} <: AbstractManifold{ℝ}
    knots::TK
    interpolation_method::TIM
end

function CurveWarpingSpace(knots::TK) where {TK<:AbstractVector}
    method = FritschButlandMonotonicInterpolation()
    return CurveWarpingSpace{TK,typeof(method)}(knots, method)
end

function manifold_dimension(x::CurveWarpingSpace)
    return length(knots)
end

function make_warping(M::CurveWarpingSpace, y::AbstractVector)
    return extrapolate(
        interpolate(convert(Vector, M.knots), convert(Vector, y), M.interpolation_method),
        Flat(),
    )
end

struct WarpingCompositionOperation <: AbstractGroupOperation end

const CurveWarpingGroup{TCWS<:CurveWarpingSpace} =
    GroupManifold{ℝ,TCWS,WarpingCompositionOperation}

function CurveWarpingGroup(cws::CurveWarpingSpace)
    return GroupManifold(cws, WarpingCompositionOperation())
end

function identity_element(cwg::CurveWarpingGroup)
    return make_warping(cwg.manifold, cwg.manifold.knots)
end
function identity_element(cwg::CurveWarpingGroup, a)
    return make_warping(cwg.manifold, cwg.manifold.knots)
end

function inv(cwg::CurveWarpingGroup, p)
    ys = map(p, cwg.manifold.knots)
    ys[end] = 1 # makes things easier for autodiff
    # TODO: make it better somehow?
    m = [1 / Interpolations.gradient1(p, t) for t in cwg.manifold.knots]
    return extrapolate(
        interpolate(
            ys,
            convert(Vector, cwg.manifold.knots),
            KnownDerivativesMonotonicInterpolation(m),
        ),
        Flat(),
    )
end
inv(::CurveWarpingGroup, p::Identity{WarpingCompositionOperation}) = p

function compose(cwg::CurveWarpingGroup, p1, p2)
    return make_warping(cwg.manifold, map(t -> p1(p2(t)), cwg.manifold.knots))
end
function compose(::CurveWarpingGroup, p1::Identity{WarpingCompositionOperation}, p2)
    return p2
end
function compose(::CurveWarpingGroup, p1, p2::Identity{WarpingCompositionOperation})
    return p1
end
function compose(
    cwg::CurveWarpingGroup,
    p1::Identity{WarpingCompositionOperation},
    p2::Identity{WarpingCompositionOperation},
)
    return p1
end

"""
    CurveWarpingAction(M::AbstractManifold, cwg::CurveWarpingGroup)

Space of left actions of the group of curve warpings `cwg` on the manifold `M` of curves.
"""
struct CurveWarpingAction{TM<:AbstractManifold,TCWG<:CurveWarpingGroup} <:
       AbstractGroupAction{LeftAction}
    manifold::TM
    cwg::TCWG
end

function Manifolds.g_manifold(A::CurveWarpingAction)
    return A.manifold
end

function Manifolds.base_group(A::CurveWarpingAction)
    return A.cwg
end

function apply!(A::CurveWarpingAction{<:DiscretizedCurves}, q, a, p)
    ts = map(a, A.cwg.manifold.knots)
    itp = make_interpolant(A.manifold, p)
    grid = A.manifold.grid
    rep_size = representation_size(A.manifold.manifold)
    for (i, t) in zip(get_iterator(A.manifold), ts)
        copyto!(_write(A.manifold, rep_size, q, i), itp(t))
    end
    return q
end

function apply!(
    A::CurveWarpingAction{<:DiscretizedCurves},
    q,
    ::Identity{WarpingCompositionOperation},
    p,
)
    copyto!(q, p)
    return q
end

function Manifolds.inverse_apply(A::CurveWarpingAction, a, p)
    inva = inv(base_group(A), a)
    return apply(A, inva, p)
end

function isapprox(M::CurveWarpingSpace, x, y; kwargs...)
    xc = map(x, M.knots)
    yc = map(y, M.knots)
    return isapprox(xc, yc; kwargs...)
end


text:
module DistributedControl
    using LinearAlgebra

    include("miscellanea.jl")

    export fixed_modes, real_eigvectors, rosenbrok
end


text:
#=
Doubly linked list
=#

"""
Space complexity: O(1)
"""
mutable struct Node{T}
    data::T
    prev::Node{T}
    next::Node{T}

    function Node{T}() where T
        node = new{T}()
        node.prev = node
        node.next = node
        return node
    end

    function Node(x::T) where T
        node = new{T}()
        node.data = x
        return node
    end
end

"""
Doubly linked list

Space complexity: O(n)
"""
mutable struct SimpleLinkedList{T}
    root::Node{T}
    n::Int

    SimpleLinkedList{T}() where T = new{T}(Node{T}(), 0)
end

#=
Core methods
=#

"""
Time complexity: O(1)
"""
Base.first(sll::SimpleLinkedList) = sll.root.next.data

"""
Time complexity: O(1)
"""
Base.last(sll::SimpleLinkedList) = sll.root.prev.data

"""
Time complexity: O(n)
"""
function Base.getindex(sll::SimpleLinkedList, index::Integer)
    checkbounds(sll, index)
    node = sll.root
    for i in 1:index
        node = node.next
    end
    return node.data
end

"""
Time complexity: O(n)
"""
function Base.setindex!(sll::SimpleLinkedList{T}, item::T, index::Integer) where T
    checkbounds(sll, index)
    node = sll.root
    for i in 1:index
        node = node.next
    end
    node.data = item
end

"""
Time complexity: O(1)
"""
function Base.push!(sll::SimpleLinkedList{T}, item::T) where T
    node = Node(item)
    last = sll.root.prev
    last.next = node
    node.prev = last
    node.next = sll.root
    sll.root.prev = node
    sll.n += 1
    return sll
end

"""
Time complexity: O(1)
"""
function Base.pop!(sll::SimpleLinkedList)
    isempty(sll) && throw(ArgumentError("`SimpleLinkedList` is empty"))
    last = sll.root.prev
    penultimate = last.prev
    penultimate.next = sll.root
    sll.root.prev = penultimate
    sll.n -= 1
    return last.data
end

"""
Time complexity: O(n)
"""
function Base.insert!(sll::SimpleLinkedList{T}, index::Integer, item::T) where T
    checkbounds(sll, index)
    new_node = Node(item)
    node = sll.root
    for i in 1:index - 1
        node = node.next
    end
    prev_node = node
    next_node = node.next
    prev_node.next = new_node
    new_node.prev = prev_node
    next_node.prev = new_node
    new_node.next = next_node
    sll.n += 1
    return sll
end

"""
Time complexity: O(1)
"""
function Base.popfirst!(sll::SimpleLinkedList)
    isempty(sll) && throw(ArgumentError("`SimpleLinkedList` is empty"))
    first = sll.root.next
    second = first.next
    second.prev = sll.root
    sll.root.next = second
    sll.n -= 1
    return first.data
end

"""
Time complexity: O(n)
"""
function Base.deleteat!(sll::SimpleLinkedList, index::Integer)
    checkbounds(sll, index)
    node = sll.root
    for i in 1:index - 1
        node = node.next
    end
    prev_node = node
    next_node = node.next.next
    prev_node.next = next_node
    next_node.prev = prev_node
    sll.n -= 1
    return sll
end

#=
Interface
=#

Base.length(sll::SimpleLinkedList) = sll.n
Base.isempty(sll::SimpleLinkedList) = sll.n < 1
Base.eltype(sll::SimpleLinkedList{T}) where T = T

function Base.iterate(sll::SimpleLinkedList)
    length(sll) == 0 && return nothing
    return sll.root.next.data, sll.root.next.next
end

function Base.iterate(sll::SimpleLinkedList, n::Node)
    n === sll.root && return nothing
    return n.data, n.next
end

function Base.checkbounds(sll::SimpleLinkedList, index::Integer)
    1 ≤ index ≤ length(sll) || throw(BoundsError(sll, index))
    return nothing
end


text:
module LinearRL

using LinearAlgebra

# Import important stuff from above.
import ..AbstractValueFunction
import ..AbstractVFunction
import ..AbstractQFunction
import ..LearningUpdate
import ..Optimizer
import ..update!

export VFunction, TD, WISTD, VtraceTD, TDC, GTD2, update!

"""
    VFunction(num_features)
A structure hosting the weights for a linear value function. Used for Linear function approximation.
"""
mutable struct VFunction{W} <: AbstractVFunction
    weights::W
end

VFunction(num_features::Integer; init=zeros) = new(init(num_features))

feature_type(v::VFunction) = typeof(v.weights[1])
weights(v::VFunction) = v.weights

(value::VFunction)(ϕ::Array{Number, 1}) = dot(value.weights, ϕ)

update!(value::VFunction, Δθ) = value.weights .+= Δθ
update!(value::VFunction, ϕ, δ) = value.weights .+= δ.*ϕ

"""
    SparseVFunction(num_features)
A structure for when the feature vector is known to be sparse, and handed to the agent as a list of indices. Significantly faster than
a normal value function in this special case.
"""
mutable struct SparseVFunction{W} <: AbstractVFunction
    weights::W
end

SparseVFunction(num_features::Integer; init=zeros) = new(init(num_features))

feature_type(v::SparseVFunction) = Integer
weights(v::SparseVFunction) = v.weights

(value::SparseVFunction)(ϕ::Array{Integer, 1}) = sum(value.weights[ϕ])

update!(value::SparseVFunction, Δθ) = value.weights .+= Δθ
update!(value::SparseVFunction, ϕ, δ) = value.weights[ϕ] .+= δ

update!(value::AbstractVFunction, opt::LearningUpdate, s_t, s_tp1, r, γ, ρ, terminal, a_t, a_tp1, target_policy) =
    update!(value, opt, s_t, s_tp1, r, γ, ρ, terminal)

"""
    TD(α)
Online Temporal Difference Learning.
"""
mutable struct TD <: LearningUpdate
    α::Float64
end

function update!(value::AbstractVFunction, lu::TD, ϕ_t, ϕ_tp1, r, γ, ρ, terminal)
    α = lu.α
    δ = r + γ*values(ϕ_t) - value(ϕ_t)
    Δθ = (α*ρ*δ)
    update!(value, ϕ_t, Δθ)
end

"""
    TDLambda(α, λ)
Online Temporal Difference Learning with eligibility traces.
"""
mutable struct TDLambda <: LearningUpdate
    α::Float64
    λ::Float64
    γ_t::IdDict{AbstractVFunction, Float64}
    e::IdDict{AbstractVFunction, Array{Number}}
end

function update!(value::VFunction, lu::TDLambda, ϕ_t, ϕ_tp1, r, γ, ρ, terminal)
    α = lu.α
    λ = lu.λ
    δ = r + γ*values(ϕ_t) - value(ϕ_t)

    e = get!(lu.e, value, zero(weights(value)))::typeof(weights(value))
    γ_t = get!(lu.γ_t, value, 0.0)
    e .= ρ.((λ*γ_t).*e .+ ϕ_t)

    Δθ = (α.*e)
    update!(value, Δθ)

    lu.γ_t[value] = γ
end

#------------------------------#
#
#Action state value functions.
#
#------------------------------#
export QFunction, SparseQFunction, get_values, WatkinsQ, watkins_q_target, feature_type

get_values(value::AbstractQFunction, ϕ) = [value(ϕ, a) for a in 1:value.num_actions]

"""
    QFunction(num_features, num_features_per_action, num_actions)
Linear QFunction. Assumes no sparsity, and an array of floats for a feature vector.
"""
mutable struct QFunction <: AbstractQFunction
    weights::Array{Float64}
    num_features_per_action::Integer
    num_actions::Integer
    QFunction(num_features::Integer, num_features_per_action::Integer, num_actions::Integer) =
        new(zeros(num_features), num_features_per_action, num_actions)
end

feature_type(q::QFunction) = Float64
weights(q::QFunction) = q.weights
(value::QFunction)(ϕ, action) =
    dot(value.weights[(value.num_features_per_action*(action-2) + 1):(value.num_features_per_action*(action-1) + 1)], ϕ)
update!(value::QFunction, ϕ, action, δ) = value.weights[(value.num_features_per_action*(action-2) + 1):(value.num_features_per_action*(action-1) + 1)] .+= δ*ϕ

"""
    SparseQFunction(num_features, num_features_per_action, num_actions)
QFunction assuming sparsity.
"""
mutable struct SparseQFunction <: AbstractQFunction
    weights::Array{Float64}
    num_features_per_action::Integer
    num_actions::Integer
    SparseQFunction(num_features::Integer, num_features_per_action::Integer, num_actions::Integer) =
        new(zeros(num_features), num_features_per_action, num_actions)
end

feature_type(q::SparseQFunction) = Int64
weights(q::SparseQFunction) = q.weights
(value::SparseQFunction)(ϕ, action) =
    sum(value.weights[ϕ .+ (value.num_features_per_action*(action-1) + 1)])
function update!(value::SparseQFunction, ϕ, action, δ)
    value.weights[ϕ .+ (value.num_features_per_action*(action-1))] .+= δ
end

"""
    ActionSparseQFunction(num_features_per_action, num_actions)
Same as SparseQFunction, with different implementation.
"""
mutable struct ActionSparseQFunction <: AbstractQFunction
    weights::Array{Float64}
    num_features_per_action::Integer
    num_actions::Integer
    ActionSparseQFunction(num_features_per_action::Integer, num_actions::Integer) =
        new(zeros(num_actions, num_features_per_action), num_features_per_action, num_actions)
end

feature_type(q::ActionSparseQFunction) = Int64
weights(q::ActionSparseQFunction) = q.weights
(value::ActionSparseQFunction)(ϕ, action::Integer) =
    sum(value.weights[action, ϕ])
function update!(value::ActionSparseQFunction, ϕ, action, δ)
    value.weights[action, ϕ] .+= δ
end

"""
    WatkinsQ(α)
Q Learning as defined by Watkins
"""
mutable struct WatkinsQ <: LearningUpdate
    α::Float64
end
watkins_q_target(q::AbstractQFunction, ϕ, r) = r + maximum(get_values(q, ϕ))

function update!(value::AbstractQFunction, opt::WatkinsQ, ϕ_t, ϕ_tp1, r, γ, ρ, terminal, a_t, a_tp1=nothing, target_policy=nothing)
    α = opt.α
    δ = watkins_q_target(value, ϕ_tp1, r) - value(ϕ_t, a_t)
    Δθ = α.*δ
    update!(value, ϕ_t, a_t, Δθ)
end

end

max_stars_repo_name: mkschleg/JuliaRL

max_issues_repo_name: mkschleg/JuliaRL

max_forks_repo_name: mkschleg/JuliaRL


text:
module CUSPARSE
importall Base.SparseArrays

using CUDArt

typealias SparseChar Char
import Base.one
import Base.zero

export CudaSparseMatrixCSC, CudaSparseMatrixCSR,
       CudaSparseMatrixHYB, CudaSparseMatrixBSR,
       CudaSparseMatrix, AbstractCudaSparseMatrix,
       CudaSparseVector

include("util.jl")
include("libcusparse_types.jl")

function statusmessage( status )
    if status == CUSPARSE_STATUS_SUCCESS
        return "cusparse success"
    end
    if status == CUSPARSE_STATUS_NOT_INITIALIZED
        return "cusparse not initialized"
    end
    if status == CUSPARSE_STATUS_ALLOC_FAILED
        return "cusparse allocation failed"
    end
    if status == CUSPARSE_STATUS_INVALID_VALUE
        return "cusparse invalid value"
    end
    if status == CUSPARSE_STATUS_ARCH_MISMATCH
        return "cusparse architecture mismatch"
    end
    if status == CUSPARSE_STATUS_MAPPING_ERROR
        return "cusparse mapping error"
    end
    if status == CUSPARSE_STATUS_EXECUTION_FAILED
        return "cusparse execution failed"
    end
    if status == CUSPARSE_STATUS_INTERNAL_ERROR
        return "cusparse internal error"
    end
    if status == CUSPARSE_STATUS_MATRIX_TYPE_NOT_SUPPORTED
        return "cusparse matrix type not supported"
    end
end

function statuscheck( status )
    if status == CUSPARSE_STATUS_SUCCESS
        return nothing
    end
    warn("CUSPARSE error triggered from:")
    Base.show_backtrace(STDOUT, backtrace())
    println()
    throw(statusmessage( status ))
end

const libcusparse = Libdl.find_library(["libcusparse"],["/usr/local/cuda"])
if isempty(libcusparse)
    error("CUSPARSE library not found in /usr/local/cuda!")
end

include("libcusparse.jl")

#setup handler for cusparse

cusparsehandle = cusparseHandle_t[0]
cusparseCreate( cusparsehandle )

#clean up handle at exit
atexit( ()->cusparseDestroy(cusparsehandle[1]) )

include("sparse.jl")

end


text:
module Cosmology

using QuadGK, Unitful
import Unitful: km, s, Gyr
using UnitfulAstro: Mpc, Gpc

export cosmology,
       age,
       angular_diameter_dist,
       comoving_radial_dist,
       comoving_transverse_dist,
       comoving_volume,
       comoving_volume_element,
       distmod,
       H,
       hubble_dist,
       hubble_time,
       luminosity_dist,
       lookback_time,
       scale_factor

abstract type AbstractCosmology end
abstract type AbstractClosedCosmology <: AbstractCosmology end
abstract type AbstractFlatCosmology <: AbstractCosmology end
abstract type AbstractOpenCosmology <: AbstractCosmology end

struct FlatLCDM{T <: Real} <: AbstractFlatCosmology
    h::T
    Ω_Λ::T
    Ω_m::T
    Ω_r::T
end
FlatLCDM(h::Real, Ω_Λ::Real, Ω_m::Real, Ω_r::Real) =
    FlatLCDM(promote(float(h), float(Ω_Λ), float(Ω_m), float(Ω_r))...)

a2E(c::FlatLCDM, a) = sqrt(c.Ω_r + c.Ω_m * a + c.Ω_Λ * a^4)

struct ClosedLCDM{T <: Real} <: AbstractClosedCosmology
    h::T
    Ω_k::T
    Ω_Λ::T
    Ω_m::T
    Ω_r::T
end
ClosedLCDM(h::Real, Ω_k::Real, Ω_Λ::Real, Ω_m::Real, Ω_r::Real) =
    ClosedLCDM(promote(float(h), float(Ω_k), float(Ω_Λ), float(Ω_m),
                       float(Ω_r))...)

struct OpenLCDM{T <: Real} <: AbstractOpenCosmology
    h::T
    Ω_k::T
    Ω_Λ::T
    Ω_m::T
    Ω_r::T
end
OpenLCDM(h::Real, Ω_k::Real, Ω_Λ::Real, Ω_m::Real, Ω_r::Real) =
    OpenLCDM(promote(float(h), float(Ω_k), float(Ω_Λ), float(Ω_m),
                     float(Ω_r))...)

function a2E(c::Union{ClosedLCDM,OpenLCDM}, a)
    a2 = a * a
    sqrt(c.Ω_r + c.Ω_m * a + (c.Ω_k + c.Ω_Λ * a2) * a2)
end

for c in ("Flat", "Open", "Closed")
    name = Symbol("$(c)WCDM")
    @eval begin
        struct $(name){T <: Real} <: $(Symbol("Abstract$(c)Cosmology"))
            h::T
            Ω_k::T
            Ω_Λ::T
            Ω_m::T
            Ω_r::T
            w0::T
            wa::T
        end
        function $(name)(h::Real, Ω_k::Real, Ω_Λ::Real, Ω_m::Real, Ω_r::Real,
                         w0::Real, wa::Real)
            $(name)(promote(float(h), float(Ω_k), float(Ω_Λ), float(Ω_m),
                            float(Ω_r), float(w0), float(wa))...)
        end
    end
end

function WCDM(h::Real, Ω_k::Real, Ω_Λ::Real, Ω_m::Real, Ω_r::Real, w0::Real, wa::Real)
    if Ω_k < 0
        ClosedWCDM(h, Ω_k, Ω_Λ, Ω_m, Ω_r, w0, wa)
    elseif Ω_k > 0
        OpenWCDM(h, Ω_k, Ω_Λ, Ω_m, Ω_r, w0, wa)
    else
        FlatWCDM(h, Ω_k, Ω_Λ, Ω_m, Ω_r, w0, wa)
    end
end

function a2E(c::Union{FlatWCDM,ClosedWCDM,OpenWCDM}, a)
    ade = exp((1 - 3 * (c.w0 + c.wa)) * log(a) + 3 * c.wa * (a - 1))
    sqrt(c.Ω_r + (c.Ω_m + c.Ω_k * a) * a + c.Ω_Λ * ade)
end

"""
    cosmology(;h = 0.69,
               Neff = 3.04,
               OmegaK = 0,
               OmegaM = 0.29,
               OmegaR = nothing,
               Tcmb = 2.7255,
               w0 = -1,
               wa = 0)

# Parameters
* `h` - Dimensionless Hubble constant
* `OmegaK` - Curvature density (Ω_k)
* `OmegaM` - Matter density (Ω_m)
* `OmegaR` - Radiation density (Ω_r)
* `Tcmb` - CMB temperature in Kelvin; used to compute Ω_γ
* `Neff` - Effective number of massless neutrino species; used to compute Ω_ν
* `w0` - CPL dark energy equation of state; `w = w0 + wa(1-a)`
* `wa` - CPL dark energy equation of state; `w = w0 + wa(1-a)`

# Examples
```jldoctest
julia> c = cosmology()
Cosmology.FlatLCDM{Float64}(0.69, 0.7099122024007928, 0.29, 8.77975992071536e-5)

julia> c = cosmology(OmegaK=0.1)
Cosmology.OpenLCDM{Float64}(0.69, 0.1, 0.6099122024007929, 0.29, 8.77975992071536e-5)

julia> c = cosmology(w0=-0.9, OmegaK=-0.1)
Cosmology.ClosedWCDM{Float64}(0.69, -0.1, 0.8099122024007929, 0.29, 8.77975992071536e-5, -0.9, 0.0)
```
"""
function cosmology(;h = 0.69,
                   Neff = 3.04,
                   OmegaK = 0,
                   OmegaM = 0.29,
                   OmegaR = nothing,
                   Tcmb = 2.7255,
                   w0 = -1,
    wa = 0)

    if OmegaR === nothing
        OmegaG = 4.48131e-7 * Tcmb^4 / h^2
        OmegaN = Neff * OmegaG * (7 / 8) * (4 / 11)^(4 / 3)
        OmegaR = OmegaG + OmegaN
    end

    OmegaL = 1 - OmegaK - OmegaM - OmegaR

    if !(w0 == -1 && wa == 0)
        return WCDM(h, OmegaK, OmegaL, OmegaM, OmegaR, w0, wa)
    end

    if OmegaK < 0
        return ClosedLCDM(h, OmegaK, OmegaL, OmegaM, OmegaR)
        elseif OmegaK > 0
        return OpenLCDM(h, OmegaK, OmegaL, OmegaM, OmegaR)
        else
        return FlatLCDM(h, OmegaL, OmegaM, OmegaR)
    end
end

# hubble rate

scale_factor(z) = 1 / (1 + z)
E(c::AbstractCosmology, z) = (a = scale_factor(z); a2E(c, a) / a^2)
H(c::AbstractCosmology, z) = 100 * c.h * E(c, z) * km / s / Mpc

hubble_dist0(c::AbstractCosmology) = 2997.92458 / c.h * Mpc
hubble_dist(c::AbstractCosmology, z) = hubble_dist0(c) / E(c, z)

hubble_time0(c::AbstractCosmology) = 9.77814 / c.h * Gyr
hubble_time(c::AbstractCosmology, z) = hubble_time0(c) / E(c, z)

# distances

Z(c::AbstractCosmology, z::Real, ::Nothing; kws...) =
    QuadGK.quadgk(a->1 / a2E(c, a), scale_factor(z), 1; kws...)[1]
Z(c::AbstractCosmology, z₁::Real, z₂::Real; kws...) =
    QuadGK.quadgk(a->1 / a2E(c, a), scale_factor(z₂), scale_factor(z₁); kws...)[1]

comoving_radial_dist(c::AbstractCosmology, z₁, z₂ = nothing; kws...) = hubble_dist0(c) * Z(c, z₁, z₂; kws...)

"""
    comoving_radial_dist([u::Unitlike,] c::AbstractCosmology, [z₁,] z₂)

Comoving radial distance in Mpc at redshift `z₂` as seen by an observer at `z₁`.  Redshift `z₁` defaults to 0 if omitted.  Will convert to compatible unit `u` if provided.
"""
comoving_radial_dist

comoving_transverse_dist(c::AbstractFlatCosmology, z₁, z₂ = nothing; kws...) =
    comoving_radial_dist(c, z₁, z₂; kws...)
function comoving_transverse_dist(c::AbstractOpenCosmology, z₁, z₂ = nothing; kws...)
    sqrtok = sqrt(c.Ω_k)
    hubble_dist0(c) * sinh(sqrtok * Z(c, z₁, z₂; kws...)) / sqrtok
end
function comoving_transverse_dist(c::AbstractClosedCosmology, z₁, z₂ = nothing; kws...)
    sqrtok = sqrt(abs(c.Ω_k))
    hubble_dist0(c) * sin(sqrtok * Z(c, z₁, z₂; kws...)) / sqrtok
end

angular_diameter_dist(c::AbstractCosmology, z; kws...) =
    comoving_transverse_dist(c, z; kws...) / (1 + z)
angular_diameter_dist(c::AbstractCosmology, z₁, z₂; kws...) =
    comoving_transverse_dist(c, z₁, z₂; kws...) / (1 + z₂)

"""
    angular_diameter_dist([u::Unitlike,] c::AbstractCosmology, [z₁,] z₂)

Ratio of the proper transverse size in Mpc of an object at redshift `z₂` to its angular size in radians, as seen by an observer at `z₁`.  Redshift `z₁` defaults to 0 if omitted.  Will convert to compatible unit `u` if provided.
"""
angular_diameter_dist

luminosity_dist(c::AbstractCosmology, z; kws...) =
    comoving_transverse_dist(c, z; kws...) * (1 + z)

"""
    luminosity_dist([u::Unitlike,] c::AbstractCosmology, z)

Bolometric luminosity distance in Mpc at redshift `z`. Will convert to compatible unit `u` if provided.
"""
luminosity_dist

"""
    distmod(c::AbstractCosmology, z)

Distance modulus in magnitudes at redshift `z`.
"""
distmod(c::AbstractCosmology, z; kws...) =
    5 * log10(luminosity_dist(c, z; kws...) / Mpc) + 25

# volumes

"""
    comoving_volume([u::Unitlike,] c::AbstractCosmology, z)

Comoving volume in cubic Gpc out to redshift `z`. Will convert to compatible unit `u` if provided.
"""
comoving_volume(c::AbstractFlatCosmology, z; kws...) =
    (4pi / 3) * (comoving_radial_dist(Gpc, c, z; kws...))^3
function comoving_volume(c::AbstractOpenCosmology, z; kws...)
    DH = hubble_dist0(Gpc, c)
    x = comoving_transverse_dist(Gpc, c, z; kws...) / DH
    sqrtok = sqrt(c.Ω_k)
    2pi * (DH)^3 * (x * sqrt(1 + c.Ω_k * x^2) - asinh(sqrtok * x) / sqrtok) / c.Ω_k
end
function comoving_volume(c::AbstractClosedCosmology, z; kws...)
    DH = hubble_dist0(Gpc, c)
    x = comoving_transverse_dist(Gpc, c, z; kws...) / DH
    sqrtok = sqrt(abs(c.Ω_k))
    2pi * (DH)^3 * (x * sqrt(1 + c.Ω_k * x^2) - asin(sqrtok * x) / sqrtok) / c.Ω_k
end

"""
    comoving_volume_element([u::Unitlike,] c::AbstractCosmology, z)

Comoving volume element in Gpc out to redshift `z`. Will convert to compatible unit `u` if provided.
"""
comoving_volume_element(c::AbstractCosmology, z; kws...) =
    hubble_dist0(Gpc, c) * angular_diameter_dist(Gpc, c, z; kws...)^2 / a2E(c, scale_factor(z))

# times

T(c::AbstractCosmology, a0, a1; kws...) = QuadGK.quadgk(x->x / a2E(c, x), a0, a1; kws...)[1]
"""
    age([u::Unitlike,] c::AbstractCosmology, z)

Age of the universe in Gyr at redshift `z`. Will convert to compatible unit `u` if provided.
"""
age(c::AbstractCosmology, z; kws...) = hubble_time0(c) * T(c, 0, scale_factor(z); kws...)

"""
    lookback_time([u::Unitlike,] c::AbstractCosmology, z)

Difference between age at redshift 0 and age at redshift `z` in Gyr. Will convert to compatible unit `u` if provided.
"""
lookback_time(c::AbstractCosmology, z; kws...) = hubble_time0(c) * T(c, scale_factor(z), 1; kws...)

# Easily select a different unit
for f in (:hubble_dist0, :hubble_dist, :hubble_time0, :hubble_time, :comoving_radial_dist,
          :comoving_transverse_dist, :angular_diameter_dist, :luminosity_dist,
          :comoving_volume, :comoving_volume_element, :age, :lookback_time)
    @eval $f(u::Unitful.Unitlike, args...; kws...) = uconvert(u, $f(args...; kws...))
end

end # module


text:
#__precompile__(true)

module UNSflow

import Dierckx: Spline1D, derivative, evaluate

import ForwardDiff

import DelimitedFiles

import Serialization

import NLsolve: nlsolve, not_in_place

import Statistics: mean

import PyPlot: plot, scatter, figure, xlabel, ylabel, xlim, ylim,
    xticks, yticks, subplot, subplot2grid, legend, axis, savefig,
    close, tight_layout

import Plots: @layout

import LaTeXStrings: @L_str

#For use in development and debugging
import Revise

export
    # kinematics types and funtions
    MotionDef,
    KinemPar,
    KinemDef,
    EldUpDef,
    EldUptstartDef,
    ConstDef,
    EldRampReturnDef,
    EldUpIntDef,
    EldUpInttstartDef,
    SinDef,
    CosDef,

    # 2D low-order solver types
    TwoDSurf,
    TwoDOFPar,
    KinemPar2DOF,
    TwoDVort,
    TwoDFlowField,
    KelvinCondition,
    KelvinKutta,

    # vortex count control utility
    delVortDef,
    delNone,
    delSpalart,

    # 3D low-order solver types
    ThreeDSurfSimple,
    KinemDef3D,
    ThreeDFieldSimple,

    # utility functions
    simpleTrapz,
    camber_calc,
    find_tstep,
    simpleInterp,
    cleanWrite,

    # 2D low-order solver function
    update_boundpos,
    update_kinem,
    update_indbound,
    update_downwash,
    update_a0anda1,
    place_tev,
    place_lev,
    update_a2toan,
    mutual_ind,
    update_a2a3adot,
    update_bv,
    ind_vel,
    wakeroll,
    update_adot,
    update_externalvel,
    controlVortCount,
    update_kinem,
    update_kinem2DOF,

    #2D low-order solver methods
    lautat,
    #lautatRoll,
    ldvm,
    ldvmLin,

    #3D low-order solver methods
    QSLLT_lautat,
    QSLLT_ldvm,

    # Postprocessing functions
    calc_forces,
    writeStamp,

    # 2D plotting functions
    viewVort2D,
    viewVortConnect2D,

    # 2D plot output functions
    makeForcePlots2D,
    makeVortPlots2D,

    # 3D plot output functions
    makeForcePlots3D,
    makeVortPlots3D,
    makeInfoPlots3D

### source files

# kinematic types
include("kinem.jl")

# utility functions
include("utils.jl")

# vortex count control utility
include("delVort.jl")

# low-order 2D solvers
include("lowOrder2D/typedefs.jl")            # type definitions
include("lowOrder2D/calcs.jl")               # calculation functions
include("lowOrder2D/solvers.jl")             # solver methods
include("lowOrder2D/postprocess.jl")         # postprocessing functions

# low-order 3D solvers
include("lowOrder3D/typedefs.jl")            # type definitions
include("lowOrder3D/calcs.jl")               # calculation functions
include("lowOrder3D/solvers.jl")             # solver methods
include("lowOrder3D/postprocess.jl")         # postprocessing functions

# 2D plotting functions
include("plots/plots2D.jl")
include("plots/plots3D.jl")

end


text:
"""
Create docs

See: https://juliadocs.github.io/Documenter.jl/stable/man/guide/
- https://juliadocs.github.io/Documenter.jl/stable/showcase/#Warning-admonition
"""

using Documenter

#push!(LOAD_PATH,"../src/")
#using joptimise
include("../src/joptimise.jl")

# module to make docs
makedocs(
	modules  = [joptimise],
    format   = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true"),
	#format = :html,
	sitename = "joptimise.jl",
	pages = [
		"Home" => "index.md",
		"Examples" => "examples.md",
		"API" => "api.md",
	],
	assets=[
        "assets/logo.png",
    ],
)

deploydocs(
    repo   = "https://github.com/Yuricst/joptimise",
    target = "build",
    deps   = nothing,
    make   = nothing,
    push_preview = true
)

max_stars_repo_name: Yuricst/joptimise

max_issues_repo_name: Yuricst/joptimise

max_forks_repo_name: Yuricst/joptimise


text:
module ReinforcementLearningZoo

using ReinforcementLearningBase
using ReinforcementLearningCore

const RLZoo = ReinforcementLearningZoo
export RLZoo

include("algorithms/algorithms.jl")

end # module


text:
using Distances, RecipesBase
export CerPosterior, CerPosteriorMcmc, CerPosteriorMcmcOutput
export McmcOutputParameters, log_posterior_prob

struct CerPosterior
    data::Vector{Matrix{Bool}}
    G_prior::CER
    α_prior::ScaledBeta
    sample_size::Int
    function CerPosterior(
        data::Vector{Matrix{Bool}},
        G_prior::CER,
        α_prior::ScaledBeta
        )
        new(data,G_prior,α_prior,length(data))
    end
end

function CerPosterior(
    data::Vector{Matrix{Int}},
    G_prior::CER,
    α_prior::ScaledBeta
    )

    if prod(y->prod(x->x∈[0,1], y), data)
        data_bool = map(x->convert(Matrix{Bool},x), data)
        CerPosterior(data_bool, G_prior, α_prior)
    else
        error("Entries must be boolean or 0/1 integers.")
    end
end

struct McmcOutputParameters
    desired_samples::Int
    burn_in::Int
    lag::Int
end

struct CerPosteriorMcmc
    τ::Float64
    ε::Float64
    output::McmcOutputParameters
    counts::Vector{Int}
    function CerPosteriorMcmc(
        ;τ::Float64=0.05, ε::Float64=0.1,
        desired_samples::Int=1000, burn_in::Int=0, lag::Int=1
        )
        output = McmcOutputParameters(desired_samples, burn_in, lag)
        new(τ, ε, output, [0,0,0])
    end
end

function reset_counts!(mcmc::CerPosteriorMcmc)
    mcmc.counts .= 0
end

acceptance_prob(mcmc::CerPosteriorMcmc) = (
    mcmc.counts[1]/mcmc.counts[3],
    mcmc.counts[2]/mcmc.counts[3]
)

struct CerPosteriorMcmcOutput
    G_sample::Vector{Matrix{Bool}}
    α_sample::Vector{Float64}
    posterior::CerPosterior
end

# get_sample_matrices(out::CerPosteriorMcmcOutput) = vec_to_adj_mat.(
#     out.G_sample,
#     directed=out.posterior.G_prior.directed,
#     self_loops=out.posterior.G_prior.self_loops
# )

@recipe function f(
    output::CerPosteriorMcmcOutput,
    G_true::Matrix{Bool}
    )
    d = Hamming()
    G_sample = output.G_sample # Observations as vectors
    α_sample = output.α_sample
    layout := (2,1)
    legend --> false
    xguide --> "Index"
    yguide --> ["Distance from True Mode" "α"]
    size --> (800, 600)
    margin --> 5mm
    y1 = [d(G_true,x) for x in G_sample]
    y2 = α_sample
    hcat(y1,y2)
end

function accept_reject_mode_fast!(
    x_curr::Vector{Bool},
    x_prop::Vector{Bool},
    α_curr::Float64,
    posterior::CerPosterior,
    x_data_sum::Vector{Int},
    x_0::Vector{Bool},
    τ::Float64
    )

    n = length(posterior.data)
    α₀ = posterior.G_prior.α

    suff_stat_diff = 0.0 # C′-C
    for (i,v) in enumerate(x_prop)
        if rand() < τ
            c = 2*x_data_sum[i] - n
            suff_stat_diff += v ? c : -c
            x_prop[i] = !(x_prop[i]) # Flip edge
        end
    end
    d = Hamming()
    prior_diffs = d(x_prop, x_0) - d(x_curr, x_0)

    log_α = (
        suff_stat_diff * log(α_curr) - suff_stat_diff * log(1-α_curr)
        + + prior_diffs * log(α₀) - prior_diffs * log(1-α₀)
    )
    if log(rand()) < log_α
        copy!(x_curr, x_prop)
        return true
    else
        copy!(x_prop, x_curr)
        return false
    end
end

function accept_reject_mode!(
    x_curr::Vector{Bool},
    x_prop::Vector{Bool},
    α_curr::Float64,
    posterior::CerPosterior,
    x_data_sum::Vector{Int},
    x_0::Vector{Bool},
    τ::Float64
    )

    n = length(x_curr)
    α₀ = posterior.G_prior.α

    log_α = 0.0
    for i in eachindex(x_prop)
        if rand() < τ
            x_prop[i] = !(x_prop[i]) # Flip edge
        end
    end
    data = posterior.data
    d = Hamming()
    suff_stat_curr = sum(d(x_curr,x) for x in data)
    suff_stat_prop = sum(d(x_prop,x) for x in data)
    prior_diffs = d(x_prop, x_0) - d(x_curr, x_0)
    log_α += (
        (suff_stat_prop - suff_stat_curr) * log(α_curr)
        + (suff_stat_curr - suff_stat_prop) * log(1-α_curr)
        + prior_diffs * log(α₀) - prior_diffs * log(1-α₀)
    )

    if log(rand()) < log_α
        copy!(x_curr, x_prop)
        return true
    else
        copy!(x_prop, x_curr)
        return false
    end
end

function accept_reject_alpha(
    α_curr::Float64,
    x_curr::Vector{Bool},
    posterior::CerPosterior,
    ε::Float64
    )

    n = length(x_curr)
    a,b = (posterior.α_prior.α, posterior.α_prior.β)
    data = posterior.data
    tot_ham_dist = sum(
        sum(u!==v for (u,v) in zip(x_curr,x)) for x in data
    )
    α_prop = rand_reflect(α_curr, ε, 0.0, 0.5)
    log_α = (
        tot_ham_dist * (log(α_prop)-log(α_curr))
        + (length(data)*n - tot_ham_dist) * (log(1-α_prop)-log(1-α_curr))
        + (a-1) * (log(α_prop)-log(α_curr))
        + (b-1) * (log(1 - 2α_prop) - log(1-2α_curr))
    )

    if log(rand()) < log_α
        return α_prop, true
    else
        return α_curr, false
    end
end

function draw_sample!(
    sample_out_G::T,
    sample_out_alpha::S,
    mcmc::CerPosteriorMcmc,
    posterior::CerPosterior;
    burn_in::Int=mcmc.output.burn_in,
    lag::Int=mcmc.output.lag,
    G_init::Union{Matrix{Int},Matrix{Bool}}=rand(Bool, posterior.G_prior.mode),
    α_init::Float64=0.25,
    loading_bar::Bool=true
    ) where {T<:Union{Vector{Vector{Bool}},SubArray},S<:Union{Vector{Float64},SubArray}}

    if loading_bar
        iter = Progress(
            length(sample_out_G)*lag + burn_in, # How many iters
            1,  # At which granularity to update loading bar
            "Chain for n = $(posterior.sample_size) (joint)....")  # Loading bar. Minimum update interval: 1 second
    end
    G_prior = posterior.G_prior
    is_directed = G_prior.directed
    has_self_loops = G_prior.self_loops
    x_0 = adj_mat_to_vec(
        G_prior.mode,
        directed=is_directed,
        self_loops=has_self_loops
    )
    data = posterior.data
    x_data_sum = adj_mat_to_vec(
        sum(data),
        directed=is_directed,
        self_loops=has_self_loops
    )
    x_curr = adj_mat_to_vec(
        G_init,
        directed=is_directed,
        self_loops=has_self_loops
    )
    x_prop = copy(x_curr)
    α_curr = α_init

    reset_counts!(mcmc)
    G_acount = 0
    α_acount = 0
    sample_count = 1
    i = 0
    while sample_count ≤ length(sample_out_G)
        i += 1
        # Store value
        if (i > burn_in) & (((i-1) % lag)==0)
            @inbounds sample_out_G[sample_count] = copy(x_curr)
            @inbounds sample_out_alpha[sample_count] = α_curr
            sample_count += 1
        end
        # Accept-reject move
        G_acount += accept_reject_mode_fast!(
            x_curr, x_prop,
            α_curr, posterior,
            x_data_sum,
            x_0,
            mcmc.τ
        )

        # Accept-reject dispersion
        α_curr, was_acc = accept_reject_alpha(
            α_curr, x_curr,
            posterior,
            mcmc.ε
        )
        α_acount += was_acc

        if loading_bar
            next!(iter)
        end
    end
    mcmc.counts[1:3] = [G_acount,α_acount, i]
end

function draw_sample(
    mcmc::CerPosteriorMcmc,
    posterior::CerPosterior;
    desired_samples::Int=mcmc.output.desired_samples,
    args...
    )
    sample_out_G = Vector{Vector{Bool}}(undef, desired_samples)
    sample_out_alpha = zeros(desired_samples)
    draw_sample!(
        sample_out_G,
        sample_out_alpha,
        mcmc,
        posterior;
        args...
    )
    return sample_out_G, sample_out_alpha
end

function cer_post_log_prob(
    x::Vector{Bool},
    α::Float64,
    data_vec::Vector{Vector{Bool}},
    x₀::Vector{Bool},
    α₀::Float64,
    α_prior::UnivariateDistribution
    )

    n = length(data_vec)
    E = length(x)

    d = Hamming()
    S = sum(d(y,x) for y in data_vec)

    return (
        S * log(α) + (n*E -  S) * log(1-α)
        + d(x,x₀) * log(α₀) + (E - d(x,x₀)) * log(1-α₀)
        + log(pdf(α_prior, α))
    )
end

function (mcmc::CerPosteriorMcmc)(
    posterior::CerPosterior;
    args...
    )

    G_sample, α_sample = draw_sample(
        mcmc, posterior;
        args...
    )
    # Make sample into matrices
    G_sample = vec_to_adj_mat.(
        G_sample,
        directed=posterior.G_prior.directed,
        self_loops=posterior.G_prior.self_loops
    )
    return CerPosteriorMcmcOutput(
        G_sample, α_sample,
        posterior
    )

end

function log_posterior_prob(
    out::CerPosteriorMcmcOutput
    )
    posterior = out.posterior
    dir, sl = (
        posterior.G_prior.directed,
        posterior.G_prior.self_loops
    )
    x_sample = adj_mat_to_vec.(
        out.G_sample,
        directed=dir,
        self_loops=sl
    )
    # Post-hoc log-probability evaluation
    data_vec = adj_mat_to_vec.(
        posterior.data,
        directed=dir,
        self_loops=sl
    )
    x₀ = adj_mat_to_vec(
        posterior.G_prior.mode,
        directed=dir,
        self_loops=sl
    )
    α₀ = posterior.G_prior.α
    log_post_prob = [
        cer_post_log_prob(
            x, α,
            data_vec,
            x₀, α₀,
            posterior.α_prior
        ) for (x,α) in zip(x_sample, out.α_sample)
    ]
    return log_post_prob
end

max_stars_repo_path: src/Graph Models/CER/CER_posterior.jl

max_issues_repo_path: src/Graph Models/CER/CER_posterior.jl

max_forks_repo_path: src/Graph Models/CER/CER_posterior.jl


text:
using BlockArrays, LinearAlgebra, Test

include("test_blockindices.jl")
include("test_blockarrays.jl")
include("test_blockviews.jl")
include("test_blocks.jl")
include("test_blockrange.jl")
include("test_blockarrayinterface.jl")
include("test_blockbroadcast.jl")
include("test_blocklinalg.jl")
include("test_blockproduct.jl")
include("test_blockreduce.jl")
include("test_blockdeque.jl")


text:
# Autogenerated wrapper script for MariaDB_Connector_C_jll for aarch64-apple-darwin
export libmariadb

using LibCURL_jll
using Libiconv_jll
using OpenSSL_jll
using Zlib_jll
JLLWrappers.@generate_wrapper_header("MariaDB_Connector_C")
JLLWrappers.@declare_library_product(libmariadb, "@rpath/libmariadb.3.dylib")
function __init__()
    JLLWrappers.@generate_init_header(LibCURL_jll, Libiconv_jll, OpenSSL_jll, Zlib_jll)
    JLLWrappers.@init_library_product(
        libmariadb,
        "lib/mariadb/libmariadb.3.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module LinuxPerf

const SYS_perf_event_open = 298

type perf_event_attr
    typ :: UInt32
    size :: UInt32
    config :: UInt64
    sample_period_or_freq :: UInt64
    sample_type :: UInt64
    read_format :: UInt64
    flags :: UInt64
    wakeup_events_or_watermark :: UInt32
    bp_type :: UInt32
    bp_addr_or_config1 :: UInt64
    bp_len_or_config2 :: UInt64
    branch_sample_type :: UInt64

    sample_regs_user :: UInt64
    sample_stack_user :: UInt32
    clockid :: Int32
    sample_regs_intr :: UInt64
    aux_watermark :: UInt32
    __reserved_2 :: UInt32

    perf_event_attr() = new()
end

const EVENT_TYPES =
    [
     (:hw, 0, # PERF_TYPE_HARDWARE
      [(:cycles, 0), # PERF_COUNT_HW_CPU_CYCLES
       (:instructions, 1), # PERF_COUNT_HW_INSTRUCTIONS
       (:cache_access, 2), # PERF_COUNT_HW_CACHE_REFERENCES
       (:cache_misses, 3), # PERF_COUNT_HW_CACHE_MISSES
       (:branches, 4), # PERF_COUNT_HW_BRANCH_INSTRUCTIONS
       (:branch_mispredicts, 5), # PERF_COUNT_HW_BRANCH_MISSES
       (:bus_cycles, 6), # PERF_COUNT_HW_BUS_CYCLES
       (:stalled_cycles_frontend, 7), # PERF_COUNT_HW_STALLED_CYCLES_FRONTEND
       (:stalled_cycles_backend, 8), # PERF_COUNT_HW_STALLED_CYCLES_BACKEND
       (:scaled_cycles, 9) # PERF_COUNT_HW_REF_CPU_CYCLES
       ]),
     (:sw, 1, # PERF_TYPE_SOFTWARE
      [(:page_faults, 2), # PERF_COUNT_SW_PAGE_FAULTS
       (:ctx_switches, 3), # PERF_COUNT_SW_CONTEXT_SWITCHES
       (:cpu_migrations, 4), # PERF_COUNT_SW_CPU_MIGRATIONS
       (:minor_page_faults, 5), # PERF_COUNT_SW_PAGE_FAULTS_MIN
       (:major_page_faults, 6), # PERF_COUNT_SW_PAGE_FAULTS_MAJ
       ])
     ]

# cache events have special encoding
const PERF_TYPE_HW_CACHE = 3
const CACHE_TYPES =
    [(:L1_data, 0),
     (:L1_insn, 1),
     (:LLC,   2),
     (:TLB_data, 3),
     (:TLB_insn, 4),
     (:BPU, 5)]
const CACHE_OPS =
    [(:read, 0),
     (:write, 1),
     (:prefetch, 2)]
const CACHE_EVENTS =
    [(:access, 0),
     (:miss, 1)]

const PERF_FORMAT_TOTAL_TIME_ENABLED = 1 << 0
const PERF_FORMAT_TOTAL_TIME_RUNNING = 1 << 1
const PERF_FORMAT_GROUP = 1 << 3

immutable EventType
    category :: UInt32
    event :: UInt64
end

function all_events()
    evts = EventType[]
    for (cat_name, cat_id, events) in EVENT_TYPES
        for (type_name, type_id) in events
            push!(evts, EventType(cat_id, type_id))
        end
    end
    evts
end

function Base.show(io::IO, e::EventType)
    if e.category == PERF_TYPE_HW_CACHE
        print(io, "cache:")
        cache = e.event & 0xff
        idx = findfirst(k -> k[2] == cache, CACHE_TYPES)
        print(io, idx == 0 ? cache : CACHE_TYPES[idx][1], ":")
        cache_op = (e.event & 0xff00) >> 8
        idx = findfirst(k -> k[2] == cache_op, CACHE_OPS)
        print(io, idx == 0 ? cache : CACHE_OPS[idx][1], ":")
        cache_event = (e.event & 0xff0000) >> 16
        idx = findfirst(k -> k[2] == cache_event, CACHE_EVENTS)
        print(io, idx == 0 ? cache : CACHE_EVENTS[idx][1])
    else
        for (cat_name, cat_id, events) in EVENT_TYPES
            cat_id == e.category || continue
            print(io, cat_name, ":")
            for (type_name, type_id) in events
                type_id == e.event || continue
                print(io, type_name)
                return
            end
            print(io, e.event)
            return
        end
        print(io, "event(", e.category, ":", e.event, ")")
    end
end

function EventType(cat::Symbol, event::Symbol)
    cat !== :cache || error("cache events needs 3 arguments")
    for (cat_name, cat_id, events) in EVENT_TYPES
        cat_name === cat || continue
        for (type_name, type_id) in events
            type_name === event || continue
            return EventType(cat_id, type_id)
        end
        error("event $event not found in $cat")
    end
    error("category $cat not found")
end

function EventType(cat::Symbol, cache::Symbol, op::Symbol, evt::Symbol)
    cat === :cache || error("only cache events takes 3 arguments")
    idx = findfirst(x -> x[1] === cache, CACHE_TYPES)
    idx != 0 || error("cache not found $cache")
    cache_id = CACHE_TYPES[idx][2]
    idx = findfirst(x -> x[1] === op, CACHE_OPS)
    idx != 0 || error("op not found $op")
    op_id = CACHE_OPS[idx][2]
    idx = findfirst(x -> x[1] === evt, CACHE_EVENTS)
    idx != 0 || error("cache event not found $evt")
    evt_id = CACHE_EVENTS[idx][2]
    return EventType(PERF_TYPE_HW_CACHE,
                     cache_id | (op_id << 8) | (evt_id << 16))
end

type EventGroup
    leader_fd :: Cint
    fds :: Vector{Cint}
    event_types :: Vector{EventType}
    leader_io :: IOStream
    function EventGroup(types :: Vector{EventType};
                        warn_unsupported = true,
                        userspace_only = false
                        )
        my_types = Array(EventType, 0)
        group = new(-1, Array(Cint, 0), Array(EventType, 0))
        for (i,evt_type) in enumerate(types)
            attr = perf_event_attr()
            attr.typ = evt_type.category
            attr.size = sizeof(perf_event_attr)
            attr.config = evt_type.event
            attr.sample_period_or_freq = 0
            if userspace_only
                attr.flags = (1 << 5) # exclude kernel
            end
            if group.leader_fd == -1
                attr.flags |= (1 << 0) # start disabled
            end
            attr.read_format =
                PERF_FORMAT_GROUP |
                PERF_FORMAT_TOTAL_TIME_ENABLED |
                PERF_FORMAT_TOTAL_TIME_RUNNING
            fd = ccall(:syscall, Cint, (Clong, Clong...), SYS_perf_event_open,
                       pointer_from_objref(attr),
                       0, -1, group.leader_fd, 0)
            if fd < 0
                errno = Libc.errno()
                if errno in (Libc.EINVAL,Libc.ENOENT)
                    if warn_unsupported
                        warn("$evt_type not supported, skipping")
                    end
                    continue
                else
                    if errno == Libc.EACCES && !userspace_only
                        warn("try to adjust /proc/sys/kernel/perf_event_paranoid to a value <= 1 or use user-space only events")
                    end
                    @show errno
                    error("perf_event_open error : $(Libc.strerror(errno))")
                end
            end
            push!(group.event_types, evt_type)
            push!(group.fds, fd)
            if group.leader_fd == -1
                group.leader_fd = fd
                group.leader_io = fdio(fd)
            end
        end
        reset!(group)
        group
    end
end

Base.length(g::EventGroup) = length(g.event_types)

function Base.show(io::IO, g::EventGroup)
    println(io, "EventGroup(")
    for e in g.event_types[1:end-1]
        println(io, "\t", e, ",")
    end
    print(io, "\t", g.event_types[end], ")")
end

const PERF_EVENT_IOC_ENABLE =  UInt64(0x2400)
const PERF_EVENT_IOC_DISABLE = UInt64(0x2401)
const PERF_EVENT_IOC_RESET =   UInt64(0x2403)

function ioctl(group::EventGroup, x)
    res = ccall(:ioctl, Cint, (Cint, Clong, Clong), group.leader_fd, x, 1)
    if res < 0
        error("ioctl error : $(Libc.strerror())")
    end
end
enable!(g::EventGroup) = ioctl(g, PERF_EVENT_IOC_ENABLE)
disable!(g::EventGroup) = ioctl(g, PERF_EVENT_IOC_DISABLE)
reset!(g::EventGroup) = ioctl(g, PERF_EVENT_IOC_RESET)
function Base.close(g::EventGroup)
    for fd in g.fds
        ccall(:close, Cint, (Cint,), fd)
    end
end

type PerfBench
    groups :: Vector{EventGroup}
end
immutable Counter
    event :: EventType
    value :: UInt64
    enabled :: UInt64
    running :: UInt64
end
immutable Counters
    counters :: Vector{Counter}
end
function Base.show(io::IO, c::Counters)
    for c in c.counters
        print(io, c.event, " : ")
        if c.enabled == 0
            print(io, "never enabled")
        elseif c.running == 0
            print(io, "did not run")
        else
            @printf(io, "\n\t%20d (%.1f %%)", Int64(c.value), 100*(c.running/c.enabled))
        end
        println()
    end
end
enable!(b::PerfBench) = foreach(enable!, b.groups)
disable!(b::PerfBench) = foreach(disable!, b.groups)
reset!(b::PerfBench) = foreach(reset!, b.groups)
function counters(b::PerfBench)
    c = Array(Counter, 0)
    for g in b.groups
        values = read(g.leader_io, UInt64, length(g)+1+2)
        @assert(length(g) == values[1])
        enabled, running = values[2], values[3]
        for i = 1:length(g)
            push!(c, Counter(g.event_types[i], values[3+i],
                             enabled, running))
        end
    end
    Counters(c)
end
function make_bench(x)
    groups = Array(EventGroup, 0)
    for y in x
        if isa(y, EventType)
            push!(groups, EventGroup([y]))
        else
            push!(groups, EventGroup(y))
        end
    end
    PerfBench(groups)
end

const reasonable_defaults =
    [EventType(:hw, :cycles),
     [EventType(:hw, :cache_access),
      EventType(:hw, :cache_misses)],
     [EventType(:hw, :branches),
      EventType(:hw, :branch_mispredicts),
      EventType(:hw, :instructions)],
     [EventType(:sw, :ctx_switches),
      EventType(:sw, :page_faults),
      EventType(:sw, :minor_page_faults),
      EventType(:sw, :major_page_faults),
      EventType(:sw, :cpu_migrations)],
#=     [EventType(:cache, :L1_data, :read, :access),
      EventType(:cache, :L1_data, :read, :miss)],
     [EventType(:cache, :L1_data, :write, :access),
      EventType(:cache, :L1_data, :write, :miss)]=#]

end


text:
module Test_InertiaTensor

using  Modia3D
import Modia3D.ModiaMath

Dx = 0.7
Lx = 1.0
Ly = 0.3
Lz = 0.3

@forceElement Damper(; d=1.0) begin
    w   = ModiaMath.RealScalar("w",   causality=ModiaMath.Input,  numericType=ModiaMath.WR)
    tau = ModiaMath.RealScalar("tau", causality=ModiaMath.Output, numericType=ModiaMath.WR)
end
function Modia3D.computeTorque(damper::Damper, sim::ModiaMath.SimulationState)
    damper.tau.value = damper.d*sim.time
end

vmat1 = Modia3D.Material(color="LightBlue" , transparency=0.5)
vmat2 = Modia3D.Material(color="Red")

m = 15.0

rCM = zeros(3)

rCM[1] = 20.0
rCM[2] = 5.0
rCM[3] = 10.0

I = fill(30.0,(3,3))
I[1] = 10.0
I[2] = 5.0
I[3] = 80.0
I[4] = 5.0
I[5] = 20.0
I[6] = 7.0
I[7] = 80.0
I[8] = 7.0
I[9] = 30.0

massProp = Modia3D.MassProperties(m, rCM, I)

@assembly Pendulum() begin
   world  = Modia3D.Object3D(Modia3D.CoordinateSystem(0.5*Lx))

   # Pendulum
   body   = Modia3D.Object3D(Modia3D.Solid(Modia3D.SolidBeam(Lx,Ly,Lz), massProp, vmat1))
   frame1 = Modia3D.Object3D(body; r=[-Lx/2, 0.0, 0.0])
   cyl    = Modia3D.Object3D(frame1,Modia3D.Cylinder(Ly/2,1.2*Ly; material=vmat2))

   rev1    = Modia3D.Revolute(world, frame1; axis = 3)

   d     = Damper(d=10.0)
   damper = Modia3D.AdaptorForceElementToFlange(w=d.w, tau=d.tau)
   Modia3D.connect(damper, rev1)
end

gravField = Modia3D.UniformGravityField(n=[1,0,0])

pendulum = Pendulum(sceneOptions=Modia3D.SceneOptions(gravityField=gravField,visualizeFrames=true,defaultFrameLength=0.3, enableContactDetection=false))
model    = Modia3D.SimulationModel(pendulum)
result   = ModiaMath.simulate!(model, stopTime=5.0, tolerance=1e-6,interval=0.001, log=false)
ModiaMath.plot(result, ["rev1.phi", "rev1.tau"] )

println("... success of Test_InertiaTensor.jl!")
end


text:
# additional generator tests

using Test
using RoME

##

@testset "test generateCanonicalFG_TwoPoseOdo" begin
##

fg = generateCanonicalFG_TwoPoseOdo()

solveGraph!(fg)

@test isapprox( [0;0;0.], getPPE(fg, :x0, :simulated).suggested, atol=1)
@test isapprox( [0;0;0.], getPPE(fg, :x0).suggested, atol=1)

##
end

#


text:
using Documenter
using HDF5Arrays
import SparseArrays

makedocs(
    sitename="HDF5Arrays.jl",
    modules = [HDF5Arrays],
    pages=[
        "Home" => "index.md"
    ]
)

deploydocs(;
    repo="github.com/LTLA/HDF5Arrays.jl",
)


text:
export getenv

# robustly get and parse an env var
function getenv(var, default::T) where T
    if haskey(ENV, var)
        result = tryparse(T, ENV[var])
        if result === nothing
            @warn "Could not parse $(var)=$(ENV[var]), using default value '$default'"
            default
        else
            result
        end
    else
        default
    end
end


text:
using CSV
using PyPlot
using HDF5
using PyCall
using Colors
using Statistics
using LsqFit
include("$(pwd())/startupAPiTOF.jl")

function plotStages(logfile)
  ax = gca()
  runlist_array = CSV.read(logfile)
  tform = PyPlot.matplotlib[:transforms][:blended_transform_factory](ax[:transData], ax[:transAxes])
  for i in 1:length(runlist_array[:,1])
    axvline(x=runlist_array[:time][i], linewidth=0.5, color="grey")
    text(runlist_array[:time][i],0.95,"\n$(runlist_array[i,:description])",transform=tform, rotation=90, color="grey", clip_on=true)
  end
end

function removeBG(times, traces; mode = 0, smoothing = 1, startTime = DateTime(2017,8,1,13,20), endTime = DateTime(2017,8,1,13,40))
	backgroundSubstractionMode = mode
	if (backgroundSubstractionMode == 0)
	  background=0
	  return bgCorrectedTraces = traces.-background
	elseif (backgroundSubstractionMode == "minimum")
	  background = minimum(InterpolationFunctions.averageSamples(traces,smoothing),1)
	  return bgCorrectedTraces = traces.-background
	elseif backgroundSubstractionMode == "bgtime"
	  background = mean(traces[(times.>startTime) .& (times.<endTime),:], dims = 1)
	  return bgCorrectedTraces = traces.-background
	end
	if backgroundSubstractionMode == "decaying"
	  decayingBG = InterpolationFunctions.averageSamples(traces[(times.>startTime) .& (times.<endTime),:],smoothing)
	  decayTime = InterpolationFunctions.averageSamples(times[(times.>startTime) .& (times.<endTime),:],smoothing)
	  decayTime_ms = (Dates.value.((decayTime .- decayTime[1])))[:,1]
	  model(t, p) = p[1] .* exp.(- p[2] .* t)
	  p0 = [30.0, 1.1e-6]
	  fits = []
	  bgTraces = []
	  for i in 1:size(decayingBG)[2]
	  	fit = curve_fit(model, decayTime_ms, decayingBG[:,i], p0)
		bg = model(Dates.value.(times .- startTime), fit.param)
                append!(fits, [fit])
		if i == 1
		  bgTraces = bg
                else
		  bgTraces = hcat(bgTraces, bg)
		end
	  end
	  bgTraces[times .< startTime,:] .= 0
	  return (fits, bgTraces)
	end
end

function plotTraces(times, traces, measResult; smoothing, logfile = " ")
	figure()
	semilogy(InterpolationFunctions.averageSamples(times, smoothing), InterpolationFunctions.averageSamples(traces, smoothing))
	legStrings = []
	for i = 1:(length(measResult.MasslistMasses))
	    push!(legStrings,"m/z $(round(measResult.MasslistMasses[i],digits=3)) - $(MasslistFunctions.sumFormulaStringFromCompositionArray(measResult.MasslistCompositions[:,i]))")
	end
	legend(legStrings)
	if logfile != " "
	    plotStages(logfile)
	end
end

function averageStages(times, traces, measResult, logfile)
	runlist = CSV.read(logfile)
	stagetimes = runlist[:time]
	meanTraces = zeros(length(stagetimes)-1, length(measResult.MasslistMasses))
	stdevTraces = zeros(length(stagetimes)-1, length(measResult.MasslistMasses))
	for i in 1:(length(stagetimes)-1)
		startT = stagetimes[i] + Minute(2)
		endT = stagetimes[i+1] - Minute(2)
		selector = (times .> startT) .& (times .< endT)
		meanTraces[i,:] = mean(traces[selector,:], dims = 1)
		stdevTraces[i,:] = std(traces[selector,:], dims = 1)
	end
	return (meanTraces, stdevTraces)
end

####################
# start of program #
####################

# user input #

fp = "/media/wiebke/Elements/savedData_Masterthesis_etc/DATAFLOWREACTOR/Data_Flowreactor_Innsbruck/TME/"
file = "$(fp)results/_result.hdf5"
logfile = "$(fp)TME_log.csv"

showIntermediatePlots = false

ion = "NH4+" 			# "H3O+" ...
inletLossCorrection = 1.7	# for 80 cm long, 1.0 cm diameter tubing, 10 slpm total flow, 2 slpm core sampling
smoothing = 20

simpleCalibration = true
	calibFactor = 2.47e10/1240 	# ((molecules/cm³)/ppb) / (cps/ppb)

decayingBG = true
	startDecayTime = DateTime(2017,8,1,13,45)
	endDecayTime = DateTime(2017,8,1,14,05)

massesToPlot = [
MasslistFunctions.createCompound(; C=3, O=3, N=1, H=8)[1]
MasslistFunctions.createCompound(; C=3, O=1, N=1, H=9)[1]
]

# load Data #

measResult = ResultFileFunctions.loadResults(file, massesToLoad=massesToPlot, useAveragesOnly=false, raw=false, startTime = DateTime(2017,7,30,11,00), endTime = DateTime(2017,8,7,22,00))
		# , startTime=Dates.unix2datetime(switchTimes[1]), endTime=Dates.unix2datetime(switchTimes[end])

# correct background #

bgCorrectedTraces = removeBG(measResult.Times, measResult.Traces; mode = "bgtime", startTime = DateTime(2017,8,1,13,20), endTime = DateTime(2017,8,1,13,40))
if decayingBG
	(bgFits, bgTraces) = removeBG(measResult.Times, bgCorrectedTraces; mode = "decaying", smoothing = 1, startTime = startDecayTime, endTime = endDecayTime)
	bgCorrectedTraces_simple = copy(bgCorrectedTraces)
	bgCorrectedTraces = bgCorrectedTraces .- bgTraces
end

if showIntermediatePlots
	plotTraces(measResult.Times, bgCorrectedTraces_simple, measResult; smoothing = smoothing, logfile = logfile)
	if decayingBG
	    ax = gca()
	    ax.semilogy(measResult.Times, bgTraces)

	    plotTraces(measResult.Times, bgCorrectedTraces, measResult; smoothing = smoothing, logfile = logfile)
	end
end

# correct inletLoss of radicals #

if ion == "NH4+"
	radicalIndices = findall(iseven, measResult.MasslistCompositions[3,:])
elseif ion == "H3O+"
	radicalIndices = findall(isodd, measResult.MasslistCompositions[3,:])
end
inletLoss_bgCorrectedTraces = bgCorrectedTraces
for i in radicalIndices
	print("found radical at index ", i)
	inletLoss_bgCorrectedTraces[:,i] = bgCorrectedTraces[:,i].*inletLossCorrection
end

##############################################################################################################################
# correct for mass-dependent transmission # not necessary here. Already done.                                                #
# massdep_inletLoss_bgCorrectedTraces = copy(inletLoss_bgCorrectedTraces)                                                    #
# for i in 1:length(measResult.MasslistMasses)                                                                               #
# 	massdep_inletLoss_bgCorrectedTraces[:,i] = inletLoss_bgCorrectedTraces[:,i] .* sqrt(100/measResult.MasslistMasses[i])#
# end                                                                                #
##############################################################################################################################

# calibrate #

if simpleCalibration
 	finalTraces = inletLoss_bgCorrectedTraces .* calibFactor
	if showIntermediatePlots
		plotTraces(measResult.Times, finalTraces, measResult; smoothing = smoothing, logfile = logfile)
	end
else
        finalTraces = inletLoss_bgCorrectedTraces
	print("please use the calibration-scripts with calibration data as a final step")
end

# add TME-reacted to the plot #

runlist = CSV.read(logfile)
TME_reacted_perS = runlist[:TME] .* runlist[:O3] .* 1.2e-15
TME_reac_err = sqrt.(0.1.^2 .+ (0.2./1.6 .+ 5e9./runlist[:TME]).^2 .+ ((runlist[:O3].*0.01 .+ 2.25e10)./runlist[:O3]).^2).*TME_reacted_perS

(meanTraces, stdevTraces) = averageStages(measResult.Times, finalTraces, measResult, logfile)

#.^.^.^.^.^.^.^.^.^.^.^.^.^.^.^.^.^.^.^.^.^.^.^.^.^.^.^.^.^.^.^.^.^.
############## make averaged XvsY-plot #############################
#.^.^.^.^.^.^.^.^.^.^.^.^.^.^.^.^.^.^.^.^.^.^.^.^.^.^.^.^.^.^.^.^.^.

outfp = "/media/wiebke/Extreme SSD/FRpaper/data/"

font = Dict([("fontsize", 18)])
cols = ["blue", "red", "green"]
figure(figsize = (8.3,8.3))
for i in 1:length(meanTraces[1,:])
	errorbar(TME_reacted_perS[1:end-1], meanTraces[:,i], xerr = TME_reac_err[1:end-1], yerr = stdevTraces[:,i], marker = "o", linestyle = " ", mfc = cols[i], mec = cols[i], ecolor= "grey", elinewidth=1, capsize=2)
end
xlim(1e7,1e9)
ylim(1e8,1e10)
xscale("log")
yscale("log")
xlabel("TME reaction rate [s⁻¹cm⁻³]", fontdict = font)
ylabel("products [cm⁻³]", fontdict = font)
#xticklabels(minor = "true")
tick_params(which = "both", labelsize=16)
legend([L"$\mathrm{C_3H_6O}$", L"$\mathrm{C_3H_5O_3}$"], fontsize = 16)
plot([1e7,1e9],[1e7,1e9].*9.4, color = "black")
fill_between([1e7,1e9],[1e7,1e9].*9.4, [1e7,1e9].*10.5, color = "orange", alpha = 0.3)
fill_between([1e7,1e9],[1e7,1e9].*8.3, [1e7,1e9].*9.4, color = "orange", alpha = 0.3)
text(1.2e7,4.0e8, L"$y = (9.4 \pm 1.1) \cdot x$",rotation=45, fontdict = font)
grid("on", which = "both", c = "silver", linewidth = 0.2, markevery = 2)
savefig("$(outfp)TMEreacted_vs_products_loglog.pdf")


text:
"""
	fitexp(x,y)
Fit exponential curve with offset assuming following function:
y = p[1] + p[2]*exp.(p[3].*x)

**Input: Vectors**
* x: x coordinate vector (Float64 or DateTime)
* y: y coordinate vector (Float64)

**Output**
* (fitted parameters p[1-3], 95% confidence error bars)

**Example**
```
x = collect(1.:1:10*365)./365;
y = 1869.9 .- 782.0 .* exp.(-0.085.*x) .+ rand(length(x)).*20;
par,er = fitexp(x,y);
```
"""
function fitexp(x::Vector{Float64},y::Vector{Float64})
	x,y = fitprep(x,y);
	function fitexp_guess()
		# First guess based on: https://math.stackexchange.com/questions/1337601/fit-exponential-with-constant
		s = Vector{Float64}(undef,length(y));
		s[1] = 0;
		for i = 2:length(y)
			s[i] = s[i-1] + 0.5*(y[i]+y[i-1])*(x[i]-x[i-1]);
		end
		mat1 = [sum((x.-x[1]).^2) sum((x.-x[1]).*s);
				sum((x.-x[1]).*s) sum(s.^2)];
		mat2 = [sum((y.-y[1]).*(x.-x[1]));
				sum((y.-y[1]).*s)];
		ac = mat1\mat2;
		psi = exp.(ac[2]*x);
		mat3 = [length(y) sum(psi);
				sum(psi) sum(psi.^2)];
		mat4 = [sum(y);
				sum(y.*psi)];
		ab = mat3\mat4;
		return ab[1],ab[2],ac[2];
	end
	a,b,c = fitexp_guess();
	# Return estimated parameters
	return fitmodel(evalexp,x,y,[a,b,c])
end
function fitexp(x::Vector{DateTime},y::Vector{Float64})
	x = convert.(Float64,Dates.value.(x));
	return fitexp(x,y);
end

"""
	evalexp(x,par)
Evaluate/compute exponential function using given parameters and x coordinate

**Input**
* x: x coordinate Vector{DateTime or Float64} or Vector{Float64}
* par: vector with 3 parameters: par[1] + par[2]*exp.(x.*par[3])
* par: vector with 2 parameters: par[1]*exp.(x.*par[2])

**Output**
* evaluated values for given parameters (=y)

**Example**
```
x = collect(1.:1:10*365)./365;
y = evalexp(x,[1869.9,-782.,-0.085]);
```
"""
function evalexp(x::Vector{Float64},par::Vector{Float64})
	if length(par) == 3
		return par[1] .+ par[2].*exp.(x.*par[3])
	else
		return par[1].*exp.(x.*par[2])
	end
end
function evalexp(x::Vector{DateTime},par::Vector{Float64})
	x = convert.(Float64,Dates.value.(x));
	return evalexp(x,par);
end

"""
	fitprep(x,y)
Prepare data for fitting, i.e. remove NaNs,...

**Input**
* x,y: x and y Vector

**Output**
* tuple with corrected vectors (Float64)

**Example**
```
x = [1.,2,3,4];
y = [10.,NaN,30,40];
x0,y0 = fitprep(x,y);
```

"""
function fitprep(x::Vector{Float64},y::Vector{Float64})
	x0,y0 = copy(x),copy(y);
	r = findall(isnan,x+y);
	deleteat!(x0,r);
	deleteat!(y0,r);
	return x0,y0;
end

"""
	fitpoly(x,y,deg)
Polynomial fitting based on LsqFit while removing NA values prior fitting

**Input**
* x: x coordinate Vector (Float64 or DateTime)
* y: y coordinate Vector (Float64)
* deg: degree of the polynomial between 0 and 4 (default=0 ==> constant)

**Output**
* (fitted parameters, 95% confidence error bars)

**Example**
```
x = collect(1.:1:365);
y = 10. + 0.1 .* x .+ rand(length(x))./2;
par,er = fitpoly(x,y,deg=1);
```
"""
function fitpoly(x::Vector{Float64},y::Vector{Float64};deg::Int64=0)
	# First guess
	function fitpoly_guess()
		if deg == 0
			return mean(y);
		elseif deg == 1
			return [0.,mean(y)];
		elseif deg == 2
			return [0.,0.,mean(y)];
		elseif deg == 3
			return [0.,0.,0.,mean(y)];
		elseif deg == 4
			return [0.,0.,0.,0.,mean(y)];
		else
			return NaN;
		end
	end
	x,y = fitprep(x,y);
	mod = getpolymodel(deg);
	approx_val = fitpoly_guess();
	if deg == 0
		return mean(y), std(y)*2; # 1=>68%, 2=>95%,...
	elseif (deg >= 1) && (deg <= 4)
		return fitmodel(mod,x,y,approx_val);
	else
		return NaN, NaN
	end
end
function fitpoly(x::Vector{DateTime},y::Vector{Float64};deg::Int64=0)
	x = convert.(Float64,Dates.value.(x));
	return fitpoly(x,y,deg=deg);
end

"""
	fitmodel(model,x,y,approx_val;confid=0.95)
Fit data to input model using approximated values as initial guess of estimated
parameters

**Input**
model: function to be fitted
x: x coordinates vector (not Vector)
y: y values (vector). NA values on input will not be corrected!
approx_val: approximated values of the model parameters
conifd: alfa for confidence interval (default = 0.05)

**Output**
* (estimated parameters, error bars at requested confidence)

**Example**
```
x = collect(1.:1:365);
y = 10.0 .+ 0.1 .* x .+ rand(length(x))./2;
par,er = fitmodel((x,p)-> p[1] + p[2].*x,x,y,[0.1, 1]);

```
"""
function fitmodel(model::Function,x::Vector{Float64},y::Vector{Float64},
					approx_val::Vector{Float64};confid::Float64=0.05)
	fit = LsqFit.curve_fit(model,x,y,approx_val);
	return fit.param, LsqFit.margin_error(fit, 1 - confid)
end

"""
	getpolymodel(deg)
Auxiliary function to create polynomial model
"""
function getpolymodel(deg::Int64=1)
	if deg == 0
		return model0(x,p) = p*ones(length(x));
	elseif deg == 1
		return model1(x,p) = p[1].*x .+ p[2];
	elseif deg == 2
		return model2(x,p) = p[1].*x.^2 .+ p[2].*x .+ p[3];
	elseif deg == 3
		return model3(x,p) = p[1].*x.^3 .+ p[2].*x.^2 .+ p[3].*x .+ p[4];
	elseif deg == 4
		return model4(x,p) = p[1].*x.^4 .+ p[2].*x.^3 .+ p[3].*x.^2 .+ p[4].*x .+ p[5];
	else
		return NaN;
	end
end

"""
	evalpoly(par,x)
Evaluate/compute polynomial function using given parameters and x coordinate

**Input**
* x: x coordinate Vector{DateTime or Float64} or Vector{Float64}
* par: vector with n parameters (polynomial degree = n-1)

**Output**
* evaluated values for given parameters (=y)

**Example**
```
x = collect(1.:1:10*365);
y = evalpoly(x,[0.01, 10.]);
```
"""
function evalpoly(x::Vector{Float64},par::Vector{Float64})::Vector{Float64}
	mod = getpolymodel(length(par)-1);
	return mod(x,par);
end
function evalpoly(x::Vector{DateTime},par::Vector{Float64})::Vector{Float64}
	x = convert.(Float64,Dates.value.(x));
	return evalpoly(x,par);
end


text:
# Test disjoint set

using DataStructures
using Base.Test

s = DisjointSets{Int}(1:10)

@test length(s) == 10
@test num_groups(s) == 10

r = [find_root(s, i) for i in 1 : 10]
@test isa(r, Vector{Int})
@test isequal(r, [1:10])

for i = 1 : 5
    x = 2 * i - 1
    y = 2 * i
    union!(s, x, y)
end

@test length(s) == 10
@test num_groups(s) == 5

r0 = [1, 1, 3, 3, 5, 5, 7, 7, 9, 9]
r = [find_root(s, i) for i in 1 : 10]
@test isa(r, Vector{Int})
@test isequal(r, r0)

union!(s, 1, 4)
union!(s, 3, 5)
union!(s, 7, 9)

@test length(s) == 10
@test num_groups(s) == 2

r0 = [1, 1, 1, 1, 1, 1, 7, 7, 7, 7]
r = [find_root(s, i) for i in 1 : 10]
@test isa(r, Vector{Int})
@test isequal(r, r0)

push!(s, 17)

@test length(s) == 11
@test num_groups(s) == 3

r0 = [1, 1, 1, 1, 1, 1, 7, 7, 7, 7, 11]
r = [find_root(s, i) for i in [1 : 10, 17] ]
@test isa(r, Vector{Int})
@test isequal(r, r0)


text:
"""
`CropSize <: ImageTransformation`

Description
============

Crops out the area of the specified pixel dimensions
around the center of the given image

Usage
======

    CropSize(64, 64)

    CropSize(width = 64, height = 64)

Arguments
==========

- **`width`** : The desired width or the cropped out subimage
in pixels

- **`height`** : The desired height or the cropped out subimage
in pixels

Methods
========

- **`transform`** : Applies the transformation to the given Image
or set of images

Author(s)
==========

- Christof Stocker (Github: https://github.com/Evizero)

Examples
========

    using Augmentor
    using TestImages

    # load an example image
    img = testimage("lena")

    # Crop out a 64x32 area in the center of the image
    img_new = transform(CropSize(64, 32), img)

see also
=========

`ImageTransformation`, `transform`
"""
immutable CropSize <: ImageTransformation
    width::Int
    height::Int

    function CropSize(width::Int, height::Int)
        width > 0 || throw(ArgumentError("width has to be greater than 0"))
        height > 0 || throw(ArgumentError("height has to be greater than 0"))
        new(width, height)
    end
end

CropSize(; width::Int = 64, height::Int = 64) = CropSize(width, height)

Base.showcompact(io::IO, op::CropSize) = print(io, "Crop a centered $(op.width)x$(op.height) window.")
multiplier(::CropSize) = 1

function transform{T<:AbstractImage}(op::CropSize, img::T)
    w = width(img)
    h = height(img)
    op.width <= w || throw(ArgumentError("image width is smaller than desired subimage"))
    op.height <= h || throw(ArgumentError("image height is smaller than desired subimage"))

    cx = floor(Int, w / 2)
    cy = floor(Int, h / 2)
    x = cx - floor(Int, op.width / 2) + 1
    @assert x > 0
    y = cy - floor(Int, op.height / 2) + 1
    @assert y > 0

    result = crop(img, x:(x+op.width-1), y:(y+op.height-1))
    _log_operation!(result, op)
end


text:
abstract type QNScaling end
struct ShannoPhua <: QNScaling end # Matrix Conditionig and Nonlinear Optimization, 1978 Math. Prog
function (::ShannoPhua)(s, y)
    real(dot(s, s)) / sum(abs2, s)
end
struct InitialScaling{S} <: QNScaling
    scaling::S
end
(is::InitialScaling)(s, y) = is.scaling(s, y)
next(qns::QNScaling) = qns
next(is::InitialScaling) = is.scaling

max_stars_repo_name: JuliaNLSolvers/NLSolvers

max_issues_repo_name: JuliaNLSolvers/NLSolvers

max_forks_repo_name: JuliaNLSolvers/NLSolvers


text:
"""
    waveguide_domains = planar_waveguides(;index, direction, waveguide_number)
"""
function halfspace_waveguide_domain(index::Number, direction::Symbol, waveguide_number::Int)
    if direction ∈ [:w, :l, :left, :west]
        which_asymptote = :left
        lattice = BravaisLattice(a=0)
    elseif direction ∈ [:e, :r, :right, :east]
        which_asymptote = :right
        lattice = BravaisLattice(a=0)
    elseif direction ∈ [:n, :u, :t, :north, :up, :top]
        which_asymptote = :top
        lattice = BravaisLattice(b=0)
    elseif direction ∈ [:s, :b, :d, :south, :bottom, :down]
        which_asymptote = :bottom
        lattice = BravaisLattice(b=0)
    else
        throw(ArgumentError("invalid waveguide direction $(direction)"))
    end

    wvg = build_domain(; which_asymptote=which_asymptote, which_waveguide=waveguide_number, domain_type=:halfspace_waveguide, domain_params=Dict(:n₁=>index), lattice=lattice)

    return wvg
end

"""
    sys = add_halfspace_waveguide(sys; index=1, direction)

    sim = add_halfspace_waveguide(sim; index=1, direction)
"""
function add_halfspace_waveguide(sys::System; index::Number=1, direction::Symbol, waveguide_number::Int=length(sys.waveguides)+1)
    waveguide_domains = halfspace_waveguide_domain(index, direction, waveguide_number)
    return System(vcat(waveguide_domains,sys.domains))
end
function add_halfspace_waveguide(sim::Simulation; index::Number=1, direction::Symbol)
    @assert 1 ∈ sim.dis.N "halfspace waveguides only for one-dimensional-systems"
    waveguide_number = length(sim.sys.waveguides)+1
    sys = add_halfspace_waveguide(sim.sys; index=index, direction=direction, waveguide_number=waveguide_number)
    return Simulation(sim; :sys => sys)
end

"""
    sim = add_halfspace_waveguide(sim, side::Symbol; index=1)
"""
function add_halfspace_waveguide(sim::Simulation, side::Symbol; index::Number=1)
    return add_halfspace_waveguide(sim; index=index, direction=side)
end

"""
    sim = add_halfspace_waveguides(sim; index=1)
"""
function add_halfspace_waveguides(sim::Simulation; index::Number=1)
    @assert 1 ∈ sim.dis.N "halfspace waveguides only for one-dimensional-systems"
    if sim.dis.N[1]==1
        sim = add_halfspace_waveguide(sim, :bottom; index=index)
        sim = add_halfspace_waveguide(sim, :top; index=index)
    else sim.dis.N[2]==1
        sim = add_halfspace_waveguide_left(sim, :left; index=index)
        sim = add_halfspace_waveguide_right(sim, :right; index=index)
    end
    return sim
end


text:
using Pkg
Pkg.activate(".")
using Printf, StatsBase, UnicodePlots
push!(LOAD_PATH,"../src/")
using QuantumComputer

# configurable runner is below this function
function shor_beauregard(n, a, sample_size)
    n_qubit_count = convert(Int64, ceil(log2(n)))
    @printf("factoring %d using fixed a = %d (generalized %d-qubit circuit)\n", n, fixed_a, 2 * n_qubit_count + 3)
    raw_qubits = [1 0]
    for i in 2:n_qubit_count
      raw_qubits = vcat(raw_qubits, [1 0])
    end
    # here we set x to 1
    raw_qubits = vcat(raw_qubits, [0 1])
    for i in 1:(n_qubit_count + 2)
      raw_qubits = vcat(raw_qubits, [1 0])
    end

    qubits = convert(Array{Complex{Float64},2}, raw_qubits)
    classical_register = QuantumComputer.ClassicalRegister(2 * n_qubit_count)
    println("building period finding circuit (this can take a few minutes the first time)...")
    shor = QuantumComputer.Circuits.shor2n3_period_finding(n, a, false)

    @printf("sampling circuit (%d samples)...\n", sample_size)
    samples = Array{Int64,1}(undef, sample_size)
    for i in 1:sample_size
        print(".")
        superposition = QuantumComputer.Superposition(qubits)
        QuantumComputer.apply_circuit_to_superposition!(superposition, shor, classical_register)
        samples[i] = classical_register.value
    end
    println()

    labels = Array{String, 1}(undef, 0);
    vals = Array{Integer, 1}(undef, 0);

    counts = countmap(samples)
    for value in sort([key for key in keys(counts)])
        push!(labels, string(string(value), " |", string(value, base=2, pad=2*n_qubit_count), ">"))
        push!(vals, counts[value])
    end

    println(barplot(labels, vals, xlabel = "samples"))
    println()

    filtered = samples
    if length([k for k in keys(counts)]) > 6
        z = median([count for count in values(counts)]) + 1
        z = √(mean([count for count in values(counts)]))
        @printf("eliminating outliers with counts <= %f\n", z)
        filtered = [sample for sample in samples if counts[sample] > z]
        length(filtered) == 0 && throw(ErrorException("found zero non-trivial factors"))
    end

    println("converting samples to phases and determining probable denominators")
    denominators = [denominator(rationalize(value/(2^(2*n_qubit_count)), tol=1/(2^(2*n_qubit_count-1)))) for value in filtered]

    println("removing odd solutions")
    filtered = [denominator for denominator in denominators if denominator % 2 == 0]
    length(filtered) == 0 && throw(ErrorException("found zero non-trivial factors"))

    labels = Array{String, 1}(undef, 0);
    vals = Array{Integer, 1}(undef, 0);

    counts = countmap(filtered)
    for value in sort([key for key in keys(counts)])
        push!(labels, string(value))
        push!(vals, counts[value])
    end
    println(barplot(labels, vals, xlabel = "samples", ylabel = "r"))
    println()

    t(v) = counts[v]
    candidates = sort([v for v in keys(counts)], by = t, rev = true)

    @printf("trying up to %d candidates\n", n_qubit_count)
    firstn(a, n) = a[intersect(eachindex(a), 1:n)]
    for r in firstn(candidates, n_qubit_count)
        @printf("r = %d\n", r)
        raised_a = convert(BigInt, a)^(r >> 1)
        potential_factors = [gcd(raised_a - 1, n), gcd(raised_a + 1, n)]
        factors = [x for x in potential_factors if !(x in [1, n])]

        length(factors) == 0 && continue
        return first(factors)
    end
    throw(ErrorException("found zero non-trivial factors"))
end

# this one is probably fine. play with n and fixed_a
sample_size = 128

# change these (but beware adding another bit to n will increase circuit complexity and memory required)
n = length(ARGS) > 1 ? parse(Int64, ARGS[1]) : 15
fixed_a = length(ARGS) > 1 ? parse(Int64, ARGS[2]) : 7

if length(ARGS) > 2
    QuantumComputer.gate_set_cache_base(ARGS[3])
end

try
    println()
    x = shor_beauregard(n, fixed_a, sample_size)
    y = convert(Int64, n / x)

    @printf("factored %d into [%d, %d]!\n", n, x, y)
catch
    @printf("zero non-trivial factors of %d found\n", n)
end


text:
# Add our own pairwise function to be able to apply it on vectors

function pairwise(d::PreMetric, X::AbstractVector, Y::AbstractVector)
    return broadcast(d, X, permutedims(Y))
end

pairwise(d::PreMetric, X::AbstractVector) = pairwise(d, X, X)

function pairwise!(out::AbstractMatrix, d::PreMetric, X::AbstractVector, Y::AbstractVector)
    return broadcast!(d, out, X, permutedims(Y))
end

pairwise!(out::AbstractMatrix, d::PreMetric, X::AbstractVector) = pairwise!(out, d, X, X)

function pairwise(d::PreMetric, x::AbstractVector{<:Real})
    return Distances_pairwise(d, reshape(x, :, 1); dims=1)
end

function pairwise(d::PreMetric, x::AbstractVector{<:Real}, y::AbstractVector{<:Real})
    return Distances_pairwise(d, reshape(x, :, 1), reshape(y, :, 1); dims=1)
end

function pairwise!(out::AbstractMatrix, d::PreMetric, x::AbstractVector{<:Real})
    return Distances.pairwise!(out, d, reshape(x, :, 1); dims=1)
end

function pairwise!(
    out::AbstractMatrix, d::PreMetric, x::AbstractVector{<:Real}, y::AbstractVector{<:Real}
)
    return Distances.pairwise!(out, d, reshape(x, :, 1), reshape(y, :, 1); dims=1)
end

# Also defines the colwise method for abstractvectors

function colwise(d::PreMetric, x::AbstractVector)
    return zeros(Distances.result_type(d, x, x), length(x)) # Valid since d(x,x) == 0 by definition
end

function colwise(d::PreMetric, x::ColVecs)
    return zeros(Distances.result_type(d, x.X, x.X), length(x)) # Valid since d(x,x) == 0 by definition
end

function colwise(d::PreMetric, x::RowVecs)
    return zeros(Distances.result_type(d, x.X, x.X), length(x)) # Valid since d(x,x) == 0 by definition
end

## The following is a hack for DotProduct and Delta to still work
function colwise(d::Distances.UnionPreMetric, x::ColVecs)
    return Distances.colwise(d, x.X, x.X)
end

function colwise(d::Distances.UnionPreMetric, x::RowVecs)
    return Distances.colwise(d, x.X', x.X')
end

function colwise(d::Distances.UnionPreMetric, x::AbstractVector)
    return map(d, x, x)
end

function colwise(d::PreMetric, x::ColVecs, y::ColVecs)
    return Distances.colwise(d, x.X, y.X)
end

function colwise(d::PreMetric, x::RowVecs, y::RowVecs)
    return Distances.colwise(d, x.X', y.X')
end

function colwise(d::PreMetric, x::AbstractVector, y::AbstractVector)
    return map(d, x, y)
end