text:
module LReg

"""This module wraps some functionality from [GLMNet.jl](https://github.com/simonster/GLMNet.jl)
for logistic regression."""
LReg

using GLMNet, Distributions

using StatsBase, StatsFuns

import StatsBase: predict
import Distributions: log1pexp

export LR, lreg, predict, linpred

"""The `LR` type contains the coefficent vector of a logistic regression fit, as well as indexes of
included columns in the design matrix.
"""
immutable LR{T<:AbstractFloat}
    β::Vector{T}
    idx::AbstractVector{Int}
    fitwithoffset::Bool
end

"""
Returns the linear component of predicted values from a logistic regression fit with a new X matrix

** Arguments **

* `lr` - an object of type `LR`
* `newx` - matrix of new covariates to predict on. Should be the same dimension as the design matrix used to fit `lr`

** Keyword Arguments **

* `offset` - offsets. The length should be the same as `size(newx, 1)` of `lr` was fit with an offset, or 0 otherwise.
"""
function linpred(lr::LR, newx; offset=Array(eltype(newx), 0))
    if size(newx, 2) == 1
        newx = reshape(newx, size(newx, 1), 1)
    end
    p = newx * lr.β
    if lr.fitwithoffset
        length(offset) == size(newx, 1) || throw(ArgumentError("fit with offset, 'offset' kw arg should have length size(newx, 1)"))
        for i in 1:length(p) #NumericExtensions
          @inbounds p[i] += offset[i]
        end
    else
        length(offset) == 0 || throw(ArgumentError("not fit with offset, 'offset' kw arg should have length 0"))
    end
    p
end

"""
Returns the predicted on the probability scale values from a logistic regression fit with a new X matrix

** Arguments **

* `lr` - an object of type `LR`
* `newx` - matrix of new covariates to predict on. Should be the same dimension as the design matrix used to fit `lr`

** Keyword Arguments

* `offset` - offsets. The length should be the same as `size(newx, 1)` of `lr` was fit with an offset, or 0 otherwise.
"""
predict(lr::LR, newx; offset=Array(eltype(newx), 0)) = logistic(linpred(lr, newx, offset=offset))

"""
Fits a logistic regression model

** Arguments **

* `x` - design matrix
* `y` - response vector. Should have the same length as `size(x, 1)`

** Keyword Arguments **

* `wts` - weight vector. Defaults to all 1s.
* `offset` - offset vector. Defaults to a vector of length 0 for no offset.
* `subset` - column indexes for `x` that should be included in the fit. Defaults to all columns.
* `convTol` - convergence criterion for relative change in deviance. Defaults to 1.0e-8.

** Details **

An intercept is not included by default. If you want one, make sure the first column of your design matrix
is all ones. Because of how `GLMNet` handles intercepts and how `LReg` interfaces with it,
this must be the first column.
If you do that, don't forget to add the column in the same place to `newx` when you call predict.

`subset` is useful if you want to include only some columns of the design matrix but you want to call `predict`
with a `newx` matrix with the same number of columns as `x`. The coefficients corresponding to columns of `x` which
are not used in the fit are set to zero. If you would like to call `predict` with a `newx` matrix that includes
only the columns that you fit on, you should subset `x` yourself before calling `lreg`.
"""
function lreg(x, y; wts=ones(y), offset=similar(y,0), subset=1:size(x,2), convTol=1.0e-8)
    if size(x, 2) == 1
        x = reshape(x, size(x, 1), 1)
    end
    subset=sort(subset)
    fitwithoffset = length(offset) > 0

    if subset == 1:size(x,2) || collect(1:size(x,2)) == subset
        subset = 1:size(x,2)
#         return LR(coef(fit(GeneralizedLinearModel, x, y, Binomial(); wts=wts, offset=offset, convTol=convTol)), subset, fitwithoffset)
        return LR(myfit(x, y, wts=wts, offset=offset, convTol=convTol), subset, fitwithoffset)
    else
        tempx = x[:, subset]
        β = zeros(eltype(x), size(x, 2))
#        β[subset] = coef(fit(GeneralizedLinearModel, tempx, y, Binomial(); wts=wts, offset=offset, convTol=convTol))
        β[subset] = myfit(tempx, y, wts=wts, offset=offset, convTol=convTol)
        return LR(β, subset, fitwithoffset)
    end
end

loss(y, xb) =
    y == one(y)? log1pexp(-xb) :
    y == zero(y)? log1pexp(xb) :
    y * log1pexp(-xb) + (one(y)-y) * log1pexp(xb)

function myfit(x, y; wts=ones(y), offset=similar(y,0), convTol=1.0e-8)
    offsets = length(offset) == 0? nothing: offset
    try
        #annoyingly, glmnet won't do intercept only  models for some reason, even with intercept = true
        # and if there is a constant column in the design matrix, that will get a coefficient of zero
        if size(x, 2) == 1 && all(x .== 1)
            offsets==nothing || error("Fitting an intercept only model with an offset is not supported")
            #intercept only
            [logit(sum(y .* wts)/sum(wts))]
        elseif all(x[:, 1] .== 1)
            #first column is intercept
            est = glmnet(x[:, 2:end], [ones(length(y)) .- y y], Binomial(), weights=wts, offsets=offsets,
                         lambda=[0.0], intercept=true, tol=convTol)
            [est.a0; est.betas[:,1]]
        else
            #no intercept
            glmnet(x, [ones(length(y)) .- y y], Binomial(), weights=wts, offsets=offsets,
                   lambda=[0.0], intercept=false, tol=convTol).betas[:,1]
        end
    catch err
        if isdefined(Main, :LREG_DEBUG) &&
           Main.LREG_DEBUG
            fname = tempname()
            open(fname, "w") do f
                serialize(f, (x, y, wts, offset, convTol))
                info("Error in lreg. x, y, wts, offset and convTol written to $fname")
            end
        end
        rethrow(err)
    end
end

end


text:
#####
##### `frule`/`rrule`
#####

#=
In some weird ideal sense, the fallback for e.g. `frule` should actually be "get
the derivative via forward-mode AD". This is necessary to enable mixed-mode
rules, where e.g. `frule` is used within a `rrule` definition. For example,
broadcasted functions may not themselves be forward-mode *primitives*, but are
often forward-mode *differentiable*.

ChainRulesCore, by design, is decoupled from any specific AD implementation. How,
then, do we know which AD to fall back to when there isn't a primitive defined?

Well, if you're a greedy AD implementation, you can just overload `frule` and/or
`rrule` to use your AD directly. However, this won't play nice with other AD
packages doing the same thing, and thus could cause load-order-dependent
problems for downstream users.

It turns out, Cassette solves this problem nicely by allowing AD authors to
overload the fallbacks w.r.t. their own context. Example using ForwardDiff:

```
using ChainRulesCore, ForwardDiff, Cassette

Cassette.@context MyChainRuleCtx

# ForwardDiff, itself, can call `my_frule` instead of
# `frule` to utilize the ForwardDiff-injected ChainRulesCore
# infrastructure
my_frule(args...) = Cassette.overdub(MyChainRuleCtx(), frule, args...)

function Cassette.execute(::MyChainRuleCtx, ::typeof(frule), f, x::Number)
    r = frule(f, x)
    if isa(r, Nothing)
        fx, df = (f(x), (_, Δx) -> ForwardDiff.derivative(f, x) * Δx)
    else
        fx, df = r
    end
    return fx, df
end
```
=#

"""
    frule(f, x...)

Expressing `x` as the tuple `(x₁, x₂, ...)` and the output tuple of `f(x...)`
as `Ω`, return the tuple:

    (Ω, (ṡelf, ẋ₁, ẋ₂, ...) -> Ω̇₁, Ω̇₂, ...)

The second return value is the propagation rule, or the pushforward.
It takes in differentials corresponding to the inputs (`ẋ₁, ẋ₂, ...`)
and `ṡelf` the internal values of the function (for closures).

If no method matching `frule(f, xs...)` has been defined, then return `nothing`.

Examples:

unary input, unary output scalar function:

```
julia> x = rand();

julia> sinx, sin_pushforward = frule(sin, x);

julia> sinx == sin(x)
true

julia> sin_pushforward(NamedTuple(), 1) == cos(x)
true
```

unary input, binary output scalar function:

```
julia> x = rand();

julia> sincosx, sincos_pushforward = frule(sincos, x);

julia> sincosx == sincos(x)
true

julia> sincos_pushforward(NamedTuple(), 1) == (cos(x), -sin(x))
true
```

See also: [`rrule`](@ref), [`@scalar_rule`](@ref)
"""
frule(::Any, ::Vararg{Any}; kwargs...) = nothing

"""
    rrule(f, x...)

Expressing `x` as the tuple `(x₁, x₂, ...)` and the output tuple of `f(x...)`
as `Ω`, return the tuple:

    (Ω, (Ω̄₁, Ω̄₂, ...) -> (s̄elf, x̄₁, x̄₂, ...))

Where the second return value is the the propagation rule or pullback.
It takes in differentials corresponding to the outputs (`x̄₁, x̄₂, ...`),
and `s̄elf`, the internal values of the function itself (for closures)

If no method matching `rrule(f, xs...)` has been defined, then return `nothing`.

Examples:

unary input, unary output scalar function:

```
julia> x = rand();

julia> sinx, sin_pullback = rrule(sin, x);

julia> sinx == sin(x)
true

julia> sin_pullback(1) == (NO_FIELDS, cos(x))
true
```

binary input, unary output scalar function:

```
julia> x, y = rand(2);

julia> hypotxy, hypot_pullback = rrule(hypot, x, y);

julia> hypotxy == hypot(x, y)
true

julia> hypot_pullback(1) == (NO_FIELDS, (x / hypot(x, y)), (y / hypot(x, y)))
true
```

See also: [`frule`](@ref), [`@scalar_rule`](@ref)
"""
rrule(::Any, ::Vararg{Any}; kwargs...) = nothing

ext: jl


text:
module Onda

using Compat: @compat
using UUIDs, Dates, Random, Mmap
using Compat, TimeSpans, Arrow, Tables, TranscodingStreams, CodecZstd
using ConstructionBase # here for deprecations only
using Legolas
using Legolas: @row

include("utilities.jl")

include("annotations.jl")
export Annotation, write_annotations, validate_annotations, merge_overlapping_annotations

include("signals.jl")
export Signal, SamplesInfo, write_signals, validate_signals,
       channel, channel_count, sample_count, sizeof_samples, sample_type

include("serialization.jl")
export AbstractLPCMFormat, AbstractLPCMStream, LPCMFormat, LPCMZstFormat,
       format, deserialize_lpcm, serialize_lpcm, deserialize_lpcm_callback,
       deserializing_lpcm_stream, serializing_lpcm_stream, finalize_lpcm_stream

include("samples.jl")
export Samples, encode, encode!, decode, decode!, load, store

include("deprecations.jl")

end # module


text:
module LDPCsim

using PyPlot

export
    ldpcMatrix,
    parityCheckMatrix,
    SymErr,
    sumProduct,
    sumProductLog,
    sumProductLogSimple,
    sumProductHard,
    ldpcEncode,
    ldpcSim

#########################################################################
function SymErr(x,y)
    # compare x and y to get # of symbol errors (bit errors for binary data)
    number = length(find(x.!=y));
    ratio= number/length(x);
    return (number,ratio)
end

#########################################################################
function ldpcMatrix(M, N, method, no4cycle, checksPerCol)
    # Create rate 1/2 LDPC adjacency matrix
    #
    # Inputs
    #  M        : num of rows
    #  N        : num of column
    #  method   : technique for distributing non-zero element
    #             0 evencol : For each column, place 1s uniformly at random
    #             1 evenboth: For each column and row, place 1s uniformly at random
    #  no4cyle  : if 1, then eliminate length-4 cycles, otherwise do nothing
    #  checksPerCol: num of ones per column
    #
    # Output
    #  H        : LDPC adjacency matrix

    # num of ones per row (N/M ratio must be 2)
    if N/M != 2
        @printf("Code rate must be 1/2 for now\n");
    end
    ## only used for evenboth
    checksPerRow = int((N/M)*checksPerCol);

    onesInCol = zeros(Int,M,N)
    if method==0   # evencol
        # Distribute 1's uniformly at random within column
        for nn = 1:N
            onesInCol[:, nn] = randperm(M);
        end

        # Create non zero elements (1s) index
        r = reshape(onesInCol[1:checksPerCol, :], N*checksPerCol, 1);
        tmp = repmat([1:N], checksPerCol, 1);
        c = reshape(tmp, N*checksPerCol, 1);

        # Create sparse matrix H
        H = full(sparse(r[:], c[:], ones(Int,N*checksPerCol), M, N,|));

    elseif method==1     # evenboth
        # Distribute 1s uniformly at random within column
        for nn = 1:N
            onesInCol[:, nn] = randperm(M)';
        end

        # Create non zero elements (1s) index
        r = reshape(onesInCol[1:checksPerCol, :], N*checksPerCol, 1);
        tmp = repmat([1:N], checksPerCol, 1);
        c = reshape(tmp, N*checksPerCol, 1);

        # Make the number of 1s between rows as uniform as possible

        # Order row index
        ix = sortperm(r[:]);

        # Order column index based on row index
        cSort = zeros(Int,N*checksPerCol,1)
        for nn = 1:N*checksPerCol
            cSort[nn] = c[ix[nn]];
        end

        # Create new row index with uniform weight
# println("checksPerRow=$(checksPerRow)")
        tmp = repmat([1:M], checksPerRow, 1);
# println("tmp=$(tmp)")
# println("size(tmp)=$(size(tmp))")
# println("N=$(N)")
        r = reshape(tmp, N*checksPerCol, 1);

        # Create  H
        allOnes = ones(Int,N*checksPerCol,1)
        S = sparse(r[:], cSort[:], allOnes[:], M, N,|);
        H = full(S);

    else
        error("unknown method $(method)")
    end

    # Check rows that have no 1 or only have one 1
    for mm = 1:M

        n = randperm(N);
        # Add two 1s if row has no 1
        if length(find(r == mm)) == 0
            H[mm, n[1]] = 1;
            H[mm, n[2]] = 1;
            # Add one 1 if row has only one 1
        elseif length(find(r == mm)) == 1
            H[mm, n[1]] = 1;
        end

    end # for mm

    if no4cycle == 1  # if no4cycle, remove any length-4 cycle

        for mm = 1:M

            # Look for pair of row - column
            for jx = (mm + 1):M
                w = (&)(H[mm, :], H[jx, :]);
                c1 = find(w);
                lc = length(c1);
                if lc > 1

                    # If found, flip one 1 to 0 in the row with less number of 1s
                    if length(find(H[mm, :])) < length(find(H[jx, :]))
                        # Repeat the process until only one column left
                        for cc = 1:lc - 1
                            H[jx, c1[cc]] = 0;
                        end
                    else
                        for cc = 1:lc - 1
                            H[mm, c1[cc]] = 0;
                        end
                    end
                end
            end # for jx
        end # for mm
    end # if no4cycle==1

    return H
end

#########################################################################
function parityCheckMatrix(src, Hin, strategy)
    # Make parity check matrix from adjacency matrix H using LU
    # decomposition.
    #
    # Inputs
    #  src     : Binary source
    #  H       : adjacency matrix
    #  strategy: Strategy for finding the next non-zero diagonal elements
    #            0: First non-zero found by column search
    #            1: Minimum number of non-zeros in later columns
    #            2: Minimum product of:
    #               - num of non-zeros its column minus 1
    #               - num of non-zeros its row minus 1
    # Outputs
    #  c       : Check bits
    #  newH    : Rearrange H

    H = copy(Hin); # copy so as to not overwrite H

    # Get the matric dimension
    (M, N) = size(H);
    # Set a new matrix F for LU decomposition
    F = copy(H);
    # LU matrices
    L = zeros(M, N - M);
    U = zeros(M, N - M);

    # Re-order the M x (N - M) submatrix
    for mx = 1:M

        # strategy {0 = First; 1 = Mincol; 2 = Minprod}
        if strategy==0 # Create diagonal matrix using 'First' strategy

            # Find non-zero elements (1s) for the diagonal
            (r, c) = findn(F[:, mx:end]);

            # Find non-zero diagonal element candidates
            rowIndex = findn(r .== mx)[1];

            # Find the first non-zero column
            chosenCol = c[rowIndex[1]] + (mx - 1);

        elseif strategy==1 # Create diag matrix using 'Mincol' strategy

            # Find non-zero elements (1s) for the diagonal
            (r, c) = findn(F[:, mx:end]);
            colWeight = sum(F[:, mx:end], 1);

            # Find non-zero diagonal element candidates
            rowIndex = findn(r .== mx)[1];

            # Find the minimum column weight
            (x, ix) = findmin(colWeight[c[rowIndex]]);
            # Add offset to the chosen row index to match the dimension of the...
            # original matrix F
            chosenCol = c[rowIndex[ix]] + (mx - 1);

        elseif strategy==2 # Create diag matrix using 'Minprod' strategy

            # Find non-zero elements of F
            (r, c) = findn(F[:, mx:end]);
            colWeight = sum(F[:, mx:end], 1) - 1;
            rowWeight = sum(F[mx, :], 2) - 1;

            # Find non-zero diagonal element candidates
            rowIndex = findn(r .== mx)[1];
# if isempty(rowIndex)
#     println("M=$(M)")
#     println("mx=$(mx)")
#     println("r=$(r)")
# end
            # Find the minimum product
            (x, ix) = findmin(colWeight[c[rowIndex]]*rowWeight);
            # Add offset to the chosen row index to match the dimension of the
            # original matrix F
            chosenCol = c[rowIndex[ix]] + (mx - 1);

        else
            @printf("Please select columns re-ordering strategy!\n");

        end # switch

        # Re-ordering columns of both H and F
        tmp1 = F[:, mx];
        tmp2 = H[:, mx];
        F[:, mx] = F[:, chosenCol];
        H[:, mx] = H[:, chosenCol];
        F[:, chosenCol] = tmp1;
        H[:, chosenCol] = tmp2;

        # Fill the LU matrices column by column
# println("mx=$(mx)")
# println("size(mx)=$(size(mx))")
# println("L=$(L)")
# println("size(L)=$(size(L))")
# println("F=$(F)")
# println("size(F)=$(size(F))")
# println("U=$(U)")
# println("size(U)=$(size(U))")
        L[mx:end, mx] = F[mx:end, mx];
        U[1:mx, mx] = F[1:mx, mx];

        # There will be no rows operation at the last row
        if mx < M

            # Find the later rows with non-zero elements in column mx
            r2 = findn(F[(mx + 1):end, mx])[1];
            # Add current row to the later rows which have a 1 in column mx
            F[(mx + r2), :] = mod(F[(mx + r2), :] + repmat(F[mx, :], length(r2), 1), 2);

        end # if

    end # for mx
# println("H=$(H)")
# println("L=$(L)")
# println("U=$(U)")
# println("L*U=$(L*U)")
    # Find src
    z = mod(H[:, (N - M) + 1:end]*src, 2);

    # Parity check vector found by solving sparse LU
    c = mod(U\(L\z), 2);

    # Return the rearrange H
    newH = H;

    return (c, newH)
end

#########################################################################
function sumProduct(zdata, H, N0, Niter)
    # Probability-domain sum product algorithm LDPC decoder
    #
    # Inputs
    #   zdata     : Received signal vector (column vector)
    #   H         : LDPC matrix
    #   N0        : Noise variance
    #   Niter     : num of iterations
    #
    # Outputs
    #   uhat      : Decoded vector (0/1)
    (M,N) = size(H);

    # priors
    P1 = ones(size(zdata))./(1 + exp(-2*zdata./(N0/2)));
    P0 = 1 - P1;

    # Init
    K0 = zeros(M, N);
    K1 = zeros(M, N);
    R0 = zeros(M, N);
    R1 = zeros(M, N);
    Q0 = H.*repmat(P0', M, 1);
    Q1 = H.*repmat(P1', M, 1);
    uhat = zeros(N,1);

    for n = 1:Niter

        # index over columns
        for mx = 1:M

            # Find non-zeros in the column
            c1 = find(H[mx, :]);

            for k = 1:length(c1)

                # Get column products of dR\c1(l)
                dR = 1.0;
                for l = 1:length(c1)
                    if l!= k
                        dR = dR*(Q0[mx, c1[l]] - Q1[mx, c1[l]]);
                    end
                end # for l

                R0[mx, c1[k]] = (1 + dR)/2;
                R1[mx, c1[k]] = (1 - dR)/2;

            end # for k

        end # for mx

        # index over rows
        for nx = 1:N

            # Find non-zeros in the row
            r1 = find(H[:, nx]);

            for k = 1:length(r1)

                # Get row products of prodOfrij\ri(l)
                prodOfrij0 = 1;
                prodOfrij1 = 1;
                for l = 1:length(r1)
                    if l!= k
                        prodOfrij0 = prodOfrij0*R0[r1[l], nx];
                        prodOfrij1 = prodOfrij1*R1[r1[l], nx];
                    end
                end # for l

                # Update constants
                K0[r1[k], nx] = P0[nx]*prodOfrij0;
                K1[r1[k], nx] = P1[nx]*prodOfrij1;

                # Update Q0 and Q1
                Q0[r1[k], nx] = K0[r1[k], nx]./(K0[r1[k], nx] + K1[r1[k], nx]);
                Q1[r1[k], nx] = K1[r1[k], nx]./(K0[r1[k], nx] + K1[r1[k], nx]);

            end # for k

            # Update constants
            Ki0 = P0[nx]*prod(R0[r1, nx]);
            Ki1 = P1[nx]*prod(R1[r1, nx]);

            # Get Qj
            Qi0 = Ki0/(Ki0 + Ki1);
            Qi1 = Ki1/(Ki0 + Ki1);

            # Decode Qj
            if Qi1 > Qi0
                uhat[nx] = 1;
            else
                uhat[nx] = 0;
            end

        end # for nx

    end # for n

    return uhat
end

#########################################################################
function sumProductLog(zdata, H, N0, Niter)
    # Log-domain sum product algorithm LDPC decoder
    #
    # Inputs
    #  zdata     : Received signal vector (column vector)
    #  H         : LDPC matrix
    #  N0        : Noise variance
    #  Niter     : num of iterations
    #
    # Outputs
    #  uhat      : Decoded vector (0/1)

    (M,N) = size(H);

    # prior log-likelihood. Minus sign is used for 0/1 to -1/1 mapping
    Lci = (-4*zdata./N0)';

    # Init
    LR = zeros(M, N);
    Pibetaij = zeros(M, N);
    uhat = zeros(N,1);

    # Asscociate the L(ci) matrix with non-zero elements of H
    LQ = H.*repmat(Lci, M, 1);

    # Get non-zero elements
    (r, c) = findn(H);

    for n = 1:Niter

        # Get the sign and magnitude of L(Q)
        alphaij = sign(LQ);
        betaij = abs(LQ);

        for l = 1:length(r)
            Pibetaij[r[l], c[l]] = log((exp(betaij[r[l], c[l]]) + 1)/
                                       (exp(betaij[r[l], c[l]]) - 1));
        end

        # index over columns
        for i = 1:M

            # Find non-zeros in the column
            c1 = find(H[i, :]);

            # Get the summation of Pi(betaij))
            for k = 1:length(c1)

                sumOfPibetaij = 0;
                prodOfalphaij = 1;

                # Summation of Pi(betaij)\c1[k]
                Pic1 = Pibetaij[i, c1]
                sumOfPibetaij = sum(Pic1) - Pic1[1,k];
#                sumOfPibetaij = sum(Pibetaij[i, c1]) - Pibetaij[i, c1[k]];

                # Avoid division by zero/very small number, get Pi(sum(Pi[betaij]))
                if sumOfPibetaij < 1e-20
                    sumOfPibetaij = 1e-10;
                end
                PiSumOfPibetaij = log((exp(sumOfPibetaij)+1)/(exp(sumOfPibetaij) - 1));

                # Multiplication of alphaij\c1[k] (use '*' since alphaij are -1/1s)
                prodOfalphaij = prod(alphaij[i, c1])*alphaij[i, c1[k]];

                # Update L[R]
                LR[i, c1[k]] = prodOfalphaij*PiSumOfPibetaij;

            end # for k

        end # for i

        # index over rows
        for j = 1:N

            # Find non-zero in the row
            r1 = find(H[:, j]);

            for k = 1:length(r1)

                # Update L[Q] by summation of L[rij]\r1[k]
                LQ[r1[k], j] = Lci[j] + sum(LR[r1, j]) - LR[r1[k], j];

            end # for k

            # Get L[Qi]
            LQi = Lci[j] + sum(LR[r1, j]);

            # Decode L(Qi)
            if LQi < 0
                uhat[j] = 1;
            else
                uhat[j] = 0;
            end

        end # for j

    end # for n

    return uhat
end

#########################################################################
function sumProductLogSimple(zdata, H, Niter)
    # Simplified log-domain sum product algorithm LDPC decoder
    #
    # Inputs
    #  zdata     : Received signal vector (column vector)
    #  H         : LDPC matrix
    #  Niter     : num of iterations
    #
    # Outputs
    #  uhat      : Decoded vector (0/1)

    (M,N) = size(H);

    # Prior log-likelihood (simplified). Minus sign is used for 0/1 to -1/1 mapping
    Lci = -zdata';

    # Init
    LR = zeros(M, N);
    Pibetaij = zeros(M, N);
    uhat = zeros(N,1);

    # Asscociate the L(ci) matrix with non-zero elements of H
    LQ = H.*repmat(Lci, M, 1);

    for n = 1:Niter

        # Get the sign and magnitude of L(Q)
        alphaij = sign(LQ);
        betaij = abs(LQ);

        # index over columns
        for i = 1:M

            # Find non-zeros in the column
            c1 = find(H[i, :]);

            # Get the minimum of betaij
            for k = 1:length(c1)

                # Minimum of betaij\c1[k]
                minOfbetaij = realmax();
                for l = 1:length(c1)
                    if l != k
                        if betaij[i, c1[l]] < minOfbetaij
                            minOfbetaij = betaij[i, c1[l]];
                        end
                    end
                end # for l

                # Multiplication alphaij\c1(k) (use '*' since alphaij are -1/1s)
                prodOfalphaij = prod(alphaij[i, c1])*alphaij[i, c1[k]];

                # Update L[R]
                LR[i, c1[k]] = prodOfalphaij*minOfbetaij;

            end # for k

        end # for i

        # index over rows
        for j = 1:N

            # Find non-zero in the row
            r1 = find(H[:, j]);

            for k = 1:length(r1)

                # Update L[Q] by summation of L[rij]\r1[k]
                LQ[r1[k], j] = Lci[j] + sum(LR[r1, j]) - LR[r1[k], j];

            end # for k

            # Get L[Qi]
            LQi = Lci[j] + sum(LR[r1, j]);

            # Decode L(Qi)
            if LQi < 0
                uhat[j] = 1;
            else
                uhat[j] = 0;
            end

        end # for j

    end # for n

    return uhat
end

#########################################################################
function sumProductHard(zdata, H, Niter)
    # Hard-decision/bit flipping sum product algorithm LDPC decoder
    #
    # Inputs
    #  zdata     : Received signal vector (column vector)
    #  H         : LDPC matrix
    #  Niter     : num of iterations
    #
    # Outputs
    #  uhat      : Decoded vector (0/1)

    (M,N) = size(H);

    # Prior hard-decision
    ci = 0.5*(sign(zdata') + 1);

    # Init
    R = zeros(M, N);
    uhat = zeros(N,1);

    # Asscociate the ci matrix with non-zero elements of H
    Q = H.*repmat(ci, M, 1);

    for n = 1:Niter

        # index over columns
        for i = 1:M

            # Find non-zeros in the column
            c1 = find(H[i, :]);

            # Get the summation of Q\c1(k)
            for k = 1:length(c1)

                R[i, c1[k]] = mod(sum(Q[i, c1]) + Q[i, c1[k]], 2);

            end # for k

        end # for i

        # index over rows
        for j = 1:N

            # Find non-zero in the row
            r1 = find(H[:, j]);

            # Number of 1s in a row
            numOfOnes = length(find(R[r1, j]));

            for k = 1:length(r1)

                # Update Q, set '1' for majority of 1s else '0', excluding r1(k)
                if numOfOnes + ci[j] >= length(r1) - numOfOnes + R[r1[k], j]
                    Q[r1[k], j] = 1;
                else
                    Q[r1[k], j] = 0;
                end

            end # for k

            # Bit decoding
            if numOfOnes + ci[j] >= length(r1) - numOfOnes
                uhat[j] = 1;
            else
                uhat[j] = 0;
            end

        end # for j

    end # for n
    return uhat
end

#########################################################################
# @doc doc"""ldpcSim(Ns,rho, M,Niter,rate,Mc,mSparse,strategy)
# Simulate transmission of LDPC-coded signals over an AWGN channel
#
# Examples show SER as we deviate away from LTE PUSCH defaults
#  # Performance as a function of SNR
#  ldpcSim(400,[-5.0:5:20],[10])
#  # Performance as a function of M
#  ldpcSim(1000,[3.5],[4,8,16,32,64,128,256])
#  # Performance as a function of SNR for two different M values
#  ldpcSim(400,[-5.0:5:20],[10,100])
# """ ->
function ldpcSim{Ti<:Integer,Tf<:FloatingPoint}(Ns::Ti,rho::Array{Tf,1},
                M::Array{Ti,1};Niter::Ti=5,rate::Tf=.5,Mc::Ti=2,mSparse::Ti=0,
                strategy::Ti=2,Nmat=1)

    println("   Ns   rho     M   Niter rate    Mc mSparse strategy")

    out = cell(0)
    if length(rho)>1 && length(M)==1 && Nmat==1
        for s = 1:length(rho)
            push!(out,ldpcEach(Ns,Mc,M[1],rate,mSparse,strategy,Niter,rho[s]));
        end
        xval = rho;
        xlab = "rho (Eb/N0 in dB)";
        tstr = @sprintf("Mc=%d,Ns=%d,M=%d,rate=%3.2f,mSparse=%d,strategy=%d,Niter=%d",
                        Mc,Ns,M[1],rate,mSparse,strategy,Niter);
    elseif length(rho)>1 && length(M)>1 && Nmat==1
        for mx = 1:length(M)
            for rx = 1:length(rho)
                outi=ldpcEach(Ns,Mc,M[mx],rate,mSparse,strategy,Niter,rho[rx]);
                for a=1:length(outi[2])
                    outi[2][a] = outi[2][a] * ", M=$(M[mx])"
                end
                if mx==1
                    push!(out,outi)
                else
                    t1 = [out[rx][1]; outi[1]]
                    t2 = [out[rx][2]; outi[2]]
                    out[rx] = (t1,t2)
                end
            end
        end
        xval = rho;
        xlab = "rho (Eb/N0 in dB)";
        tstr = @sprintf("Mc=%d,Ns=%d,rate=%3.2f,mSparse=%d,strategy=%d,Niter=%d",
                        Mc,Ns,rate,mSparse,strategy,Niter);
    elseif length(M)>1
        for s = 1:length(M)
            push!(out,ldpcEach(Ns,Mc,M[s],rate,mSparse,strategy,Niter,rho[1]));
        end
        xval = M;
        xlab = "M: message size";
        tstr = @sprintf("Mc=%d,Ns=%d,rate=%3.2f,mSpar=%d,strat=%d,Niter=%d,Eb/N0=%4.1f",
                        Mc,Ns,rate,mSparse,strategy,Niter,rho[1]);
    elseif length(rate)>1
        for s = 1:length(rate)
            push!(out,ldpcEach(Ns,Mc,M[1],rate[s],mSparse,strategy,Niter,rho[1]));
        end
        xval = rate;
        xlab = "rate: M/N";
        tstr = @sprintf("Mc=%d,Ns=%d,M=%d,mSparse=%d,strategy=%d,Niter=%d,Eb/N0=%4.1f",
                        Mc,Ns,M[1],mSparse,strategy,Niter,rho[1]);
    elseif length(Niter)>1
        for s = 1:length(Niter)
            push!(out,ldpcEach(Ns,Mc,M[1],rate,mSparse,strategy,Niter[s],rho[1]));
        end
        xval = Niter;
        xlab = "Niter";
        tstr = @sprintf("Mc=%d,Ns=%d,M=%d,mSparse=%d,strategy=%d,rate=%3.2f,Eb/N0=%4.1f",
                        Mc,Ns,M[1],mSparse,strategy,rate,rho[1]);
    elseif length(Mc)>1
        for s = 1:length(Mc)
            push!(out,ldpcEach(Ns,Mc[s],M[1],rate,mSparse,strategy,Niter,rho[1]));
        end
        xval = log2(Mc);
        xlab = "# constellation bits";
        tstr = "";
    elseif Nmat>1
        for nm = 1:Nmat
            srand(nm)
            for s = 1:length(rho)
                outi = ldpcEach(Ns,Mc,M[1],rate,mSparse,strategy,Niter,rho[s]);
                for a=1:length(outi[2])
                    outi[2][a] = outi[2][a] * " sd$(nm)"
                end
                if nm==1
                    push!(out,outi)
                else
                    for a=1:length(outi[2])
                        t1 = [out[s][1]; outi[1]]
                        t2 = [out[s][2]; outi[2]]
                        out[s] = (t1,t2)
                    end
                end
            end
        end
        xval = rho;
        xlab = "rho (Eb/N0 in dB)";
        tstr = @sprintf("Mc=%d,Ns=%d,M=%d,rate=%3.2f,mSparse=%d,strategy=%d,Niter=%d",
                        Mc,Ns,M[1],rate,mSparse,strategy,Niter);
    else
        @printf("At least one argument must be a vector")
    end

    pColor = {"r>-", "bo--","kx-.","gd-", "c^--","m*-.",
              "rs--","gp-.","bv-", "kh--","c+-.","m.-",};
    pIdx = 1;

    # Find number of samples
    Ns = size(out,1);

    #clf()
    ser = zeros(Ns);
    for a=1:length(out[1][2])
        for k=1:Ns
            ser[k] = out[k][1][a];
        end
        semilogy(xval,ser,pColor[pIdx],label=out[1][2][a]);
        hold(true);
        if pIdx == length(pColor)
            pIdx = 1;
        else
            pIdx = pIdx + 1;
        end
    end

    hold(false);
    xlabel(xlab);
    ylabel("BER");
    legend(loc=0);
    grid(which="both",axis="y")
    grid(which="major",axis="x")
    title(tstr)

    return
end

#######################################################################
function ldpcEach(Ns,Mc,M,rate,mSparse,strategy,Niter,rdb)

    rho     = 10^(rdb/10);

    if rate != 2.0
        @sprintf("Only rate 1/2 tested so far.")
    end

    N = integer(ceil(M/rate));

    # mSparse is method for sparse matrix creation (0 = evencol; 1 = evenboth)
    # Niter is number of iterations
    # strategy is matrix reorder strategy (0 = 1st; 1 = Mincol; 2 = Minprod)
    no4cycle = 1; # if 1, remove length-4 cycles
    checksPerCol = 3; # Number of 1s per column for LDPC matrix

    # # Make the LDPC matrix
    # H = ldpcMatrix(M, N, mSparse, no4cycle, checksPerCol);

    Nalgs = 3;
    ber = zeros(Nalgs);
    times = zeros(Ns);
    algNames = cell(Nalgs)
    u = zeros(N)
    newH = zeros(M,N);
    #Profile.clear()

    if isinteractive(); @printf("     "); end

    # Make random data (0/1)
    dSource = round(rand(M, Ns));
    for jx = 1:Ns
        tic();
        if (Ns<200 || mod(jx,100)==0) && isinteractive()
            @printf("\b\b\b\b\b%5d", jx);
        end

        # Ugly temporary trick to handle parity check errors;
        flag=0
        while flag>-1
            try
                (u,newH) = ldpcEncode(M, N, mSparse, no4cycle,
                             checksPerCol,dSource[:, jx],strategy);
                flag=-1
            catch
                flag+=1
            end
        end
        # # Make the LDPC matrix
        # H = ldpcMatrix(M, N, mSparse, no4cycle, checksPerCol);
        # while rank(H)<M
        #     H = ldpcMatrix(M, N, mSparse, no4cycle, checksPerCol);
        # end
        # # Encoding message
        # (c, newH) = parityCheckMatrix(dSource[:, jx], H, strategy);
        # u = [c; dSource[:, jx]];

        # BPSK modulation
        bpskMod = 2*u - 1;
        # white gaussian noise
        N0 = 1/rho;
        tx = bpskMod + sqrt(N0)*randn(size(bpskMod));

        ax= 0;

        ax = ax+1;
        algNames[ax] = "uncoded";
        uhat = (sign(tx)+1)/2
        (num, rat) = SymErr(uhat[:], u); #  includes parity bits
        ber[ax] += rat;

        ax = ax+1;
        algNames[ax] = "sumProductLog";
        #algNames[ax] = "spl";
        uhat = sumProductLog(tx, newH, N0, Niter);
        (num, rat) = SymErr(uhat[:], u); #  includes parity bits
        ber[ax] += rat;

        ax = ax+1;
        algNames[ax] = "sumProductSimpler";
        #algNames[ax] = "sps";
        uhat = sumProductLogSimple(tx, newH, Niter);
        (num, rat) = SymErr(uhat[:], u);
        ber[ax] += rat;

        # ax = ax+1;
        # algNames[ax] = "sumProduct";
        # uhat = sumProduct(tx, newH, N0, Niter);
        # (num, rat) = SymErr(uhat[:], u);
        # ber[ax] += rat;

        # Better than uncoded only @ high SNR; after the threshold/waterfall
        # ax = ax+1;
        # algNames[ax] = "sumProductHard";
        # uhat = sumProductHard(tx, newH, Niter);
        # (num, rat) = SymErr(uhat[:], u);
        # ber[ax] += rat;

        times[jx] = toq();
        if sum(times)>Ns/40
            if isinteractive(); @printf("\b\b\b\b\b"); end
            Ns = jx
            break
        end
    end # for jx

    ber = ber/Ns;

    # using ProfileView
    # ProfileView.view()

    @printf("%5d %5.1f %5d %5d   %3.2f %5d %5d  %5d  ",
            Ns, rdb,  M, Niter, rate, Mc, mSparse, strategy)

    for ax = 1:Nalgs
        if ax<6
            @printf("    %7.4f",ber[ax])
        end
    end
    @printf("\n")

    return ber, algNames
end

#######################################################################
function ldpcEncode(M, N, mSparse, no4cycle, checksPerCol,src,strategy)
    # Make the LDPC matrix
    H = ldpcMatrix(M, N, mSparse, no4cycle, checksPerCol);
    while rank(H)<M  # Hack to handle low-rank H.
        H = ldpcMatrix(M, N, mSparse, no4cycle, checksPerCol);
    end
    # Encoding message
    (c, newH) = parityCheckMatrix(src, H, strategy);
    u = [c; src];
    return (u,newH)
end

#########################################################################
end # LDPCsim module

max_stars_repo_name: christianpeel/LDPCsim

max_issues_repo_name: christianpeel/LDPCsim

max_forks_repo_name: christianpeel/LDPCsim


text:
# Practical Machine learning
# Reinforcement learning - Q learning example
# Chapter 12

module DeepQLearning

using NNGraph

export DQN, forward, act, learn
include("dqn.jl")

end # module


text:
"The composition of several maps."
struct ComposedMap{T,MAPS} <: CompositeLazyMap{T}
    maps    ::  MAPS
end

ComposedMap(maps...) = ComposedMap{domaintype(maps[1])}(maps...)
ComposedMap{T}(maps...) where {T} = ComposedMap{T,typeof(maps)}(maps)

# TODO: make proper conversion
similarmap(m::ComposedMap, ::Type{T}) where {T} = ComposedMap{T}(m.maps...)

codomaintype(m::ComposedMap) = codomaintype(m.maps[end])

# Maps are applied in the order that they appear in m.maps
applymap(m::ComposedMap, x) = applymap_rec(x, m.maps...)
applymap_rec(x) = x
applymap_rec(x, map1, maps...) = applymap_rec(map1(x), maps...)

# The size of a composite map depends on the first and the last map to be applied
# We check whether they are scalar_to_vector, vector_to_vector, etcetera
mapsize(m::ComposedMap) = _composed_mapsize(m, m.maps[end], m.maps[1], mapsize(m.maps[end]), mapsize(m.maps[1]))
_composed_mapsize(m, m_end, m1, S_end::Tuple{Int,Int}, S1::Tuple{Int,Int}) = (S_end[1],S1[2])
_composed_mapsize(m, m_end, m1, S_end::Tuple{Int,Int}, S1::Tuple{Int}) =
    is_vector_to_scalar(m_end) ? () : (S_end[1],)
_composed_mapsize(m, m_end, m1, S_end::Tuple{Int,Int}, S1::Tuple{}) =
    is_vector_to_scalar(m_end) ? () : (S_end[1],)
_composed_mapsize(m, m_end, m1, S_end::Tuple{Int}, S1::Tuple{Int,Int}) = (S_end[1],S1[2])
_composed_mapsize(m, m_end, m1, S_end::Tuple{Int}, S1::Tuple{Int}) = (S_end[1],)
_composed_mapsize(m, m_end, m1, S_end::Tuple{Int}, S1::Tuple{}) = (S_end[1],)
_composed_mapsize(m, m_end, m1, S_end::Tuple{}, S1::Tuple{Int,Int}) = (1,S1[2])
_composed_mapsize(m, m_end, m1, S_end::Tuple{}, S1::Tuple{Int}) = ()
_composed_mapsize(m, m_end, m1, S_end::Tuple{}, S1::Tuple{}) = ()

function jacobian(m::ComposedMap, x)
    f, fd = backpropagate(x, reverse(components(m))...)
    fd
end
backpropagate(x, m1) = (m1(x), jacobian(m1, x))
function backpropagate(x, m2, ms...)
    f, fd = backpropagate(x, ms...)
    m2(f), jacobian(m2, f) * fd
end

for op in (:inverse, :leftinverse, :rightinverse)
    @eval $op(cmap::ComposedMap) = ComposedMap(reverse(map($op, components(cmap)))...)
end

inverse(m::ComposedMap, x) = inverse_rec(x, reverse(components(m))...)
inverse_rec(x) = x
inverse_rec(x, map1, maps...) = inverse_rec(inverse(map1, x), maps...)

leftinverse(m::ComposedMap, x) = leftinverse_rec(x, reverse(components(m))...)
leftinverse_rec(x) = x
leftinverse_rec(x, map1, maps...) = leftinverse_rec(leftinverse(map1, x), maps...)

rightinverse(m::ComposedMap, x) = rightinverse_rec(x, reverse(components(m))...)
rightinverse_rec(x) = x
rightinverse_rec(x, map1, maps...) = rightinverse_rec(rightinverse(map1, x), maps...)

composedmap() = ()
composedmap(m) = m
composedmap(m1, m2) = composedmap1(m1, m2)
composedmap1(m1, m2) = composedmap2(m1, m2)
composedmap2(m1, m2) = ComposedMap(m1, m2)

composedmap(m1, m2, maps...) = composedmap(composedmap(m1, m2), maps...)

composedmap(m1::ComposedMap, m2::ComposedMap) =
    ComposedMap(components(m1)..., components(m2)...)
composedmap1(m1::ComposedMap, m2) = ComposedMap(components(m1)..., m2)
composedmap2(m1, m2::ComposedMap) = ComposedMap(m1, components(m2)...)

# Arguments to ∘ should be reversed before passing on to mapcompose
(∘)(map1::AbstractMap, map2::AbstractMap) = composedmap(map2, map1)

==(m1::ComposedMap, m2::ComposedMap) =
    ncomponents(m1) == ncomponents(m2) && all(map(isequal, components(m1), components(m2)))

Display.combinationsymbol(m::ComposedMap) = Display.Symbol('∘')
Display.displaystencil(m::ComposedMap) =
    composite_displaystencil(m; reversecomponents=true)
show(io::IO, mime::MIME"text/plain", m::ComposedMap) = composite_show(io, mime, m)
show(io::IO, m::ComposedMap) = composite_show_compact(io, m)

## Lazy multiplication

"The lazy multiplication of one or more maps."
struct MulMap{T,MAPS} <: CompositeLazyMap{T}
    maps    ::  MAPS
end

MulMap(maps::Map{T}...) where {T} = MulMap{T}(maps...)
MulMap{T}(maps::Map{T}...) where {T} = MulMap{T,typeof(maps)}(maps)
MulMap{T}(maps...) where {T} = _mulmap(T, convert.(Map{T}, maps)...)
_mulmap(::Type{T}, maps...) where {T} = MulMap{T,typeof(maps)}(maps)

similarmap(m::MulMap, ::Type{T}) where {T} = MulMap{T}(m.maps...)

applymap(m::MulMap, x) = reduce(*, applymap.(components(m), Ref(x)))

multiply_map() = ()
multiply_map(m) = m
multiply_map(m1, m2) = multiply_map1(m1, m2)
multiply_map1(m1, m2) = multiply_map2(m1, m2)
multiply_map2(m1, m2) = MulMap(m1, m2)

multiply_map(m1, m2, maps...) = multiply_map(multiply_map(m1, m2), maps...)

multiply_map(m1::MulMap, m2::MulMap) =
    MulMap(components(m1)..., components(m2)...)
multiply_map1(m1::MulMap, m2) = MulMap(components(m1)..., m2)
multiply_map2(m1, m2::MulMap) = MulMap(m1, components(m2)...)

function Display.displaystencil(m::MulMap)
    A = Any[]
    list = components(m)
    push!(A, "x -> ")
    push!(A, Display.SymbolObject(list[1]))
    push!(A, "(x)")
    for i in 2:length(list)
        push!(A, " * ")
        push!(A, Display.SymbolObject(list[i]))
        push!(A, "(x)")
    end
    A
end
show(io::IO, mime::MIME"text/plain", m::MulMap) = composite_show(io, mime, m)

## Lazy sum

"The lazy sum of one or more maps."
struct SumMap{T,MAPS} <: CompositeLazyMap{T}
    maps    ::  MAPS
end

SumMap(maps::Map{T}...) where {T} = SumMap{T}(maps...)
SumMap{T}(maps::Map{T}...) where {T} = SumMap{T,typeof(maps)}(maps)
SumMap{T}(maps...) where {T} = _summap(T, convert.(Map{T}, maps)...)
_summap(::Type{T}, maps...) where {T} = SumMap{T,typeof(maps)}(maps)

similarmap(m::SumMap, ::Type{T}) where {T} = SumMap{T}(m.maps...)

applymap(m::SumMap, x) = reduce(+, applymap.(components(m), Ref(x)))

sum_map() = ()
sum_map(m) = m
sum_map(m1, m2) = sum_map1(m1, m2)
sum_map1(m1, m2) = sum_map2(m1, m2)
sum_map2(m1, m2) = SumMap(m1, m2)

sum_map(m1, m2, maps...) = sum_map(sum_map(m1, m2), maps...)

sum_map(m1::SumMap, m2::SumMap) =
    SumMap(components(m1)..., components(m2)...)
sum_map1(m1::SumMap, m2) = SumMap(components(m1)..., m2)
sum_map2(m1, m2::SumMap) = SumMap(m1, components(m2)...)

function Display.displaystencil(m::SumMap)
    A = Any[]
    list = components(m)
    push!(A, "x -> ")
    push!(A, Display.SymbolObject(list[1]))
    push!(A, "(x)")
    for i in 2:length(list)
        push!(A, " + ")
        push!(A, Display.SymbolObject(list[i]))
        push!(A, "(x)")
    end
    A
end
show(io::IO, mime::MIME"text/plain", m::SumMap) = composite_show(io, mime, m)

# Define the jacobian of a composite map
jacobian(m::ComposedMap) = composite_jacobian(reverse(components(m))...)
composite_jacobian(map1) = jacobian(map1)
composite_jacobian(map1, map2) = multiply_map(jacobian(map1) ∘ map2, jacobian(map2))
function composite_jacobian(map1, map2, maps...)
    rest = ComposedMap(reverse(maps)..., map2)
    f1 = jacobian(map1) ∘ rest
    f2 = composite_jacobian(map2, maps...)
    multiply_map(f1, f2)
end

jacobian(m::MulMap) = mul_jacobian(components(m)...)
mul_jacobian() = ()
mul_jacobian(map1) = jacobian(map1)
mul_jacobian(map1, map2) = sum_map(multiply_map(jacobian(map1), map2), multiply_map(map1, jacobian(map2)))
function mul_jacobian(map1, map2, maps...)
    rest = multiply_map(map2, maps...)
    mul_jacobian(map1, rest)
end
function jacobian(m::MulMap, x)
    z = map(t -> applymap(t,x), components(m))
    zd = map(t -> jacobian(t, x), components(m))
    sum(prod(z[1:i-1]) * zd[i] * prod(z[i+1:end]) for i in 1:ncomponents(m))
end

jacobian(m::SumMap) = sum_jacobian(components(m)...)
sum_jacobian() = ()
sum_jacobian(map1) = jacobian(map1)
sum_jacobian(maps...) = sum_map(map(jacobian, maps)...)

jacobian(m::SumMap, x) = sum(jacobian(mc, x) for mc in components(m))


text:
const E_ABSTRACT = "Abstract method. To be overriden by subtype."
const E_NODE_NOT_REACHABLE = "Node could not be reached."


text:
# addition of positive definite matrices

using PDMats

@testset "addition" begin
    for T in (Float64, Float32)
        printstyled("Testing addition with eltype = $T\n"; color=:blue)
        M = convert(Array{T,2}, [4.0 -2.0 -1.0; -2.0 5.0 -1.0; -1.0 -1.0 6.0])
        V = convert(Array{T,1}, [1.5, 2.5, 2.0])
        local X = convert(T, 2.0)

        pm1 = PDMat(M)
        pm2 = PDiagMat(V)
        pm3 = ScalMat(3, X)
        pm4 = X * I
        pm5 = PDSparseMat(sparse(M))

        pmats = Any[pm1, pm2, pm3] #, pm5]

        for p1 in pmats, p2 in pmats
            pr = p1 + p2
            @test size(pr) == size(p1)
            @test Matrix(pr) ≈ Matrix(p1) + Matrix(p2)

            pr = pdadd(p1, p2, convert(T, 1.5))
            @test size(pr) == size(p1)
            @test Matrix(pr) ≈ Matrix(p1) + Matrix(p2) * convert(T, 1.5)
        end

        for p1 in pmats
            pr = p1 + pm4
            @test size(pr) == size(p1)
            @test Matrix(pr) ≈ Matrix(p1) + pm4
        end
    end
end


text:
struct Task
    p::Int
    c::Int
    t::Int
end

mutable struct QuayCrane
    q::Int
    busy::Bool
    work_load::Int
    p::Array{Int, 1}
end

mutable struct LoadingSequence
    order::Array{Task, 1}
    P::Int
    len::Int
    tasks_left::Int
    filled::Array{Int, 1}
    loaded::Array{Int, 1}
end

function horizon_plan(J::Int, C::Int, Q::Int, delta::Int, tt::Int, tasks_by_pos::Dict)
    H=tt*C*(J-1-(delta+1)*(Q-1))
    print("H0 = "); print(H)
    for (key, value) in tasks_by_position
        h=0
        for t in value
            if t.t > h
                h=t.t
            end
        end
        H=H+h
    end
    return(H)
end

function init_ls(P::Int)
    LS = LoadingSequence([], P, 0, P, m, m)
    return(LS)
end

function init_qc(P::Int, Q::Int, qj::Array )
    QC=[]
    for q = 1:Q
        qc_pos=[]
        for j = 1:P
            if qj[j,q]==1
                push!(qc_pos, j)
            end
        end
        push!(QC, QuayCrane(q, false, 0, qc_pos))
    end
    return(QC)
end

max_stars_repo_licenses: Zlib MIT

max_issues_repo_licenses: Zlib MIT

max_forks_repo_licenses: Zlib MIT


text:
#
# ScientificCameras.jl --
#
# Julia abstract layer for scientific cameras.
#
#------------------------------------------------------------------------------
#
# This file is part of the `ScientificCameras.jl` package which is licensed
# under the MIT "Expat" License.
#
# Copyright (C) 2017-2019, Éric Thiébaut.
#

isdefined(Base, :__precompile__) && __precompile__(true)

module ScientificCameras

# Macro for derived modules to re-export the public interface of the
# `ScientificCameras` module so that the end-user does not have to explicitly
# import/use the `ScientificCameras` module.
_publicsymbols = (
    :ScientificCamera,
    :open,
    :close,
    :defaulttimeout,
    :read,
    :start,
    :stop,
    :abort,
    :wait,
    :release,
    :getfullsize,
    :getfullwidth,
    :getfullheight,
    :getroi,
    :setroi!,
    :resetroi!,
    :checkroi,
    :bitsperpixel,
    :equivalentbitstype,
    :supportedpixelformats,
    :getpixelformat,
    :setpixelformat!,
    :getcapturebitstype,
    :getspeed,
    :setspeed!,
    :checkspeed,
    :getgain,
    :setgain!,
    :getbias,
    :setbias!,
    :getgamma,
    :setgamma!,
    :processimages)

macro exportpublicinterface()
    Expr(:export, _publicsymbols...)
end
@exportpublicinterface

import Base: open, read, close, wait
@static if isdefined(Base,:start)
    import Base: start
end
import Statistics: mean, stat

include("types.jl")
include("methods.jl")
include("tools.jl")
include("calibration.jl")
import .Calibration: CalibrationData, fit

end # module


text:
function has_converged(cvchecker::ConvergenceChecker, optimizerstate::OptimizerState)
    return false
end

@with_kw struct FirstOrderOptimality <: ConvergenceChecker
    tol_tangent::Float64 = 1e-11
    tol_normal::Float64 = 1e-11
end

function has_converged(cvchecker::FirstOrderOptimality, pb::CompositeProblem, o::Optimizer, os::OptimizationState)
    return (os.norm_minsubgradient_tangent < cvchecker.tol_tangent) && (os.norm_minsubgradient_normal < cvchecker.tol_normal)
end


text:
using NumPy
using Documenter

DocMeta.setdocmeta!(NumPy, :DocTestSetup, :(using NumPy); recursive=true)

makedocs(;
    modules=[NumPy],
    authors="Mohammad Saad",
    repo="https://github.com/mdsa3d/NumPy.jl/blob/{commit}{path}#{line}",
    sitename="NumPy.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://mdsa3d.github.io/NumPy.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/mdsa3d/NumPy.jl",
    devbranch="master",
)


text:
__precompile__(true)

module FractionalGaussianFields

import LinearAlgebra,
       Random,
       SparseArrays,
       Contour,
       FFTW,
       Interpolations,
       LightGraphs,
       AsyPlots

export torus_gff,
       torus_fgf,
       gff,
       fgf,
       solvelaplace,
       laplacesurface,
       SquareGrid,
       rectangulargrid,
       domainapprox,
       zeroboundary,
       neighbors,
       vertices,
       edges,
       draw,
       boundary,
       flowline,
       interpolate

include("gridgraphs.jl")

"""
    torus_fgf(m,n,s,rng=GLOBAL_RNG)

    torus_fgf(n,s,rng=GLOBAL_RNG)

Return a sample of the FGF on an `m` × `n` square grid
with periodic boundary conditions
"""
function torus_fgf(m::Integer,
                  n::Integer,
                  s::Real,
                  rng=Random.GLOBAL_RNG)
    h = complex(zeros(m,n))
    for j = 1:m
        for k = 1:n
            h[j,k] = 1/sqrt(2.0) *
                (j+k == 2 ? 0 :
                 randn(rng) + im*randn(rng) *
                 1/(sin((j-1)*pi/m)^2+(sin((k-1)*pi/n))^2)^(s/2))
        end
    end
    return real(n/sqrt(2)*FFTW.ifft(h))
end

torus_fgf(n::Integer,s::Real,rng=Random.GLOBAL_RNG) = torus_gff(n,n,s)

torus_gff(m::Integer,n::Integer,rng=Random.GLOBAL_RNG) =
    torus_fgf(m,n,1,rng)

torus_gff(n::Integer,rng=Random.GLOBAL_RNG) = torus_gff(n,n)

using DataStructures

function laplacian(G::SquareGrid,bvertices::Set)
    function newentry!(I,J,V,i,j,v)
        push!(I,i); push!(J,j); push!(V,v);
    end
    m,n = size(G)
    vertexmap = OrderedDict(map(reverse,enumerate([(i,j) for j=1:n
                                        for i=1:m if G.vertices[i,j]])))
    I = Int64[]
    J = Int64[]
    V = Float64[]
    for (v,k) in vertexmap
        if v in bvertices
            newentry!(I,J,V,k,k,1)
        else
            nbs = neighbors(G,v...)
            l = length(nbs)
            for nb in nbs
                j = vertexmap[nb]
                newentry!(I,J,V,k,j,-1)
                newentry!(I,J,V,k,k,1)
            end
        end
    end
    SparseArrays.sparse(I,J,V), vertexmap
end

function fgf(G::SquareGrid,s::Real,bvertices::Set)
    Δ,vertexmap = laplacian(G,bvertices)
    h_vector = real(full(Δ)^(-s/2)*randn(size(Δ,1)))
    h = zeros(size(G))
    for (v,k) in vertexmap
        h[v...] = h_vector[k]
    end
    return h
end

gff(G::SquareGrid,bvertices::Set) =
    fgf(G,1,bvertices)

import Interpolations.interpolate
function interpolate(h::Array{<:Real,2})
    interpolate(h,
                Interpolations.BSpline(Interpolations.Linear()),
                Interpolations.OnGrid())
end

import Base.getindex
getindex(h::Interpolations.BSplineInterpolation,
         z0::Complex) = h[real(z0),imag(z0)]

"""
Find the flow line of angle θ in the field h.
In other words, solve the DE η'(t) = exp(ih(η(t))/χ + iθ)".
"""
function flowline(h::Interpolations.AbstractInterpolation,
                  z0::Complex,
                  χ::Real,
                  θ::Real;
                  δ::Real=0.01,
                  tol::Real=1e-3,
                  S::Set{Complex}=Set{Complex}())
    a,b = size(h.coefs)
    η = [z0]
    while 1.0 ≤ real(η[end]) ≤ a && 1.0 ≤ imag(η[end]) ≤ b
        push!(η, η[end] + δ * exp(im*h[η[end]]/χ + im*θ))
        w = η[end]
        for z in S
            if abs2(z-w) < tol
                return η
            end
        end
    end
    return η
end

flowline(h::Array{<:Real,2},
         z₀::Complex,
         χ::AbstractFloat,
         θ::AbstractFloat;
         kwargs...) = flowline(interpolate(h),z₀,χ,θ;kwargs...)

function zeroboundary(h::Array{<:Real,2})
    m,n = size(h)
    G = rectangulargrid(m,n)
    harmonic = solvelaplace(G,(i,j)->0,(i,j)->h[i,j])
    h - harmonic
end

end # module


text:
using ExtensibleUnions
using Test
using Traceur


text:
using Test
using SparseIR

@testset "basis.jl" begin
    @testset "consistency" begin
        β = 2.0
        wmax = 5.0
        ε = 1e-5
        basis_f, basis_b = SparseIR.finite_temp_bases(β, wmax, ε)
        bs = FiniteTempBasisSet(β, wmax, ε)

        @test length(bs.basis_f) == length(basis_f)
        @test length(bs.basis_b) == length(basis_b)
    end

    @testset "consistency2" begin
        β = 2.0
        wmax = 5.0
        ε = 1e-5

        sve_result = sve_logistic[β * wmax]
        basis_f, basis_b = SparseIR.finite_temp_bases(β, wmax, ε, sve_result)
        smpl_τ_f = TauSampling(basis_f)
        smpl_τ_b = TauSampling(basis_b)
        smpl_wn_f = MatsubaraSampling(basis_f)
        smpl_wn_b = MatsubaraSampling(basis_b)

        bs = FiniteTempBasisSet(β, wmax, ε; sve_result)
        @test smpl_τ_f.sampling_points == smpl_τ_b.sampling_points
        @test bs.smpl_tau_f.matrix == smpl_τ_f.matrix
        @test bs.smpl_tau_b.matrix == smpl_τ_b.matrix

        @test bs.smpl_wn_f.matrix == smpl_wn_f.matrix
        @test bs.smpl_wn_b.matrix == smpl_wn_b.matrix
    end
end


text:
using DSGE
using HDF5
path = dirname(@__FILE__)
writing_output = false

#-----------------------------------------------------------------
# csminwel
#-----------------------------------------------------------------
custom_settings = Dict{Symbol, Setting}(
    :date_forecast_start  => Setting(:date_forecast_start, quartertodate("2015-Q4")))
m = AnSchorfheide(custom_settings = custom_settings, testing = true)

# Load data
file = "$path/../reference/optimize_in.h5"
x0   = h5read(file, "params")
data = h5read(file, "data")'

# For regenerating test file
params_test = deepcopy(x0)
data_test   = Matrix{Float64}(data)

file = "$path/../reference/optimize_out.h5"
minimizer  = h5read(file, "minimizer")
minimum    = h5read(file, "minimum")
H_expected = h5read(file, "H")

# See src/estimate/estimate.jl
DSGE.update!(m, x0)
n_iterations = 3

x0 = Float64[p.value for p in m.parameters]
out, H = optimize!(m, data; iterations = n_iterations)

# Re-generate test file
if writing_output
    minimizer  = out.minimizer
    minimum    = out.minimum
    H_expected = H
    h5open("$path/../reference/optimize_out.h5", "w") do file
        file["minimizer"] = minimizer
        file["minimum"]   = minimum
        file["H"]         = H_expected
    end
end
@testset "Check optimize minimizers are the same [csminwel]" begin
    @test @test_matrix_approx_eq minimizer out.minimizer
    @test minimum ≈ out.minimum atol=1e-6
    @test @test_matrix_approx_eq H_expected H
end

# TODO: fix non-csminwel tests!
#=
#-----------------------------------------------------------------
# Simulated Annealing
#-----------------------------------------------------------------
custom_settings = Dict{Symbol, Setting}(
    :date_forecast_start  => Setting(:date_forecast_start, quartertodate("2015-Q4")))
m = AnSchorfheide(custom_settings = custom_settings, testing = true)

# Load data
file = "$path/../reference/optimize.h5"#_simulated_annealing.h5"
x0   = h5read(file, "params")
data = h5read(file, "data")'

# For regenerating test file
params_test = deepcopy(x0)
data_test   = Matrix{Float64}(data)

minimizer  = h5read(file, "minimizer")
minimum    = h5read(file, "minimum")
H_expected = h5read(file, "H")

# See src/estimate/estimate.jl
DSGE.update!(m, x0)
n_iterations = 1000

x0 = Float64[p.value for p in m.parameters]
out, H = optimize!(m, data; iterations = n_iterations, method = :simulated_annealing)

# Re-generate test file
if writing_output
    h5open("$path/../reference/optimize_simulated_annealing.h5", "w") do file
        file["params"] = params_test
        file["data"] = data_test
        file["minimizer"] = out.minimizer
        file["minimum"] = out.minimum
        file["H"] = H
    end
end
@testset "Check optimize minimizers are the same [simulated_annealing]" begin
    #@test @test_matrix_approx_eq minimizer out.minimizer
    #@test minimum ≈ out.minimum atol=1e-6
    #@test @test_matrix_approx_eq H_expected H
end

#-----------------------------------------------------------------
# Nelder Mead
#-----------------------------------------------------------------
custom_settings = Dict{Symbol, Setting}(
    :date_forecast_start  => Setting(:date_forecast_start, quartertodate("2015-Q4")))
m = AnSchorfheide(custom_settings = custom_settings, testing = true)

# Load data
file = "$path/../reference/optimize_nelder_mead.jld2"
x0   = load(file, "params")
data = load(file, "data")

# For regenerating test file
params_test = deepcopy(x0)
data_test   = data

minimizer  = load(file, "minimizer")
minimum    = load(file, "minimum")
H_expected = load(file, "H")

# See src/estimate/estimate.jl
DSGE.update!(m, x0)
n_iterations = 3

x0 = Float64[p.value for p in m.parameters]
out, H = optimize!(m, data; iterations = n_iterations, method = :nelder_mead)

# Re-generate test file
if writing_output
    JLD2.jldopen("$path/../reference/optimize_nelder_mead.jld2",
                 true, true, true, IOStream) do file
        file["params"] = params_test
        file["data"] = data_test
        file["minimizer"] = out.minimizer
        file["minimum"] = out.minimum
        file["H"] = H
    end
end
@testset "Check optimize minimizers are the same [nelder_mead]" begin
    @test @test_matrix_approx_eq minimizer out.minimizer
    @test minimum ≈ out.minimum atol=1e-6
    @test @test_matrix_approx_eq H_expected H
end

#-----------------------------------------------------------------
# LBFGS
#-----------------------------------------------------------------
custom_settings = Dict{Symbol, Setting}(
    :date_forecast_start  => Setting(:date_forecast_start, quartertodate("2015-Q4")))
m = AnSchorfheide(custom_settings = custom_settings, testing = true)

# Load data
file = "$path/../reference/optimize.h5"#_lbfgs.h5"
x0   = h5read(file, "params")
data = h5read(file, "data")'

# For regenerating test file
params_test = deepcopy(x0)
data_test   = Matrix{Float64}(data)

minimizer  = h5read(file, "minimizer")
minimum    = h5read(file, "minimum")
H_expected = h5read(file, "H")

# See src/estimate/estimate.jl
DSGE.update!(m, x0)
n_iterations = 3

x0 = Float64[p.value for p in m.parameters]
#out, H = optimize!(m, data; iterations = n_iterations, method = :lbfgs)

# Re-generate test file
if writing_output
    h5open("$path/../reference/optimize_lbfgs.h5", "w") do file
        file["params"] = params_test
        file["data"] = data_test
        file["minimizer"] = out.minimizer
        file["minimum"] = out.minimum
        file["H"] = H
    end
end
@testset "Check optimize minimizers are the same [lbfgs]" begin
    #@test @test_matrix_approx_eq minimizer out.minimizer
    #@test minimum ≈ out.minimum atol=1e-6
    #@test @test_matrix_approx_eq H_expected H
end

=#
nothing


text:
using LinearAlgebra
export l1dpαI,multIxAy!,multAxIy!,compute_residual
export thomas_div!,laplace_div!
export l1dpαI_dirichlet,l1dpαI_neuman_left

struct Laplacian1D
    dxm2_::Float64
    n_::Int
    bcr_::BoundaryCondition
    bcl_::BoundaryCondition
    firstdiag_::Float64
    lastdiag_::Float64
    diag_::Float64
    upperdiag_::Float64
    L_::Float64
    function Laplacian1D(n,L,bcleft::BoundaryCondition,bcright::BoundaryCondition)
        dx=L/n
        dxm2=1/(dx^2)
        d=2dxm2
        u=-1dxm2
        (bcright == neumann) ? ld=dxm2 : ld=3dxm2
        (bcleft == neumann) ? fd=dxm2 : fd=3dxm2
        new(dxm2,n,bcleft,bcright,fd,ld,d,u,L)
    end
end

function Base.checkbounds(a::Laplacian1D, i,j)
    ((i>1) && (j>1) && (i<=a.n_) && (j<=a.n_))
end
function Base.getindex(a::Laplacian1D, i,j)
    i==j==1 && return a.firstdiag_
    i==j==a.n_ && return a.lastdiag_
    @boundscheck checkbounds(a,i,j)
    i==j && return a.diag_
    i==(j+1) && return a.upperdiag_
    i==(j-1) && return a.upperdiag_
    0.
end
Base.size(Laplacian1D) = (Laplacian1D.n_,Laplacian1D.n_)
Base.size(Laplacian1D,i) = Laplacian1D.n_

function LinearAlgebra.SymTridiagonal(a::Laplacian1D)
    n=a.n_
    D=ones(n)
    D .*= a.diag_
    D[1]=a.firstdiag_
    D[end]=a.lastdiag_
    U=ones(n-1)
    U .*= a.upperdiag_
    SymTridiagonal(D,U)
end


text:
using ImageSegmentation, Images, Test, SimpleWeightedGraphs, LightGraphs, StaticArrays, RegionTrees

include("core.jl")
include("region_growing.jl")
include("felzenszwalb.jl")
include("fast_scanning.jl")
include("watershed.jl")
include("region_merging.jl")
include("meanshift.jl")


text:
#use basket pricer to price single asset Asians
#supports a term-structure of volatilities and drifts
#past observations may be recorded at time 0 using forward=sum(w*obs) and rescaling the weight.
export priceAsianFixedStrike, priceAsianFloatingStrike

function priceAsianFixedStrike(
    p::DeelstraBasketPricer,
    isCall::Bool,
    strike::AbstractFloat,
    discountFactor::AbstractFloat,
    spot::AbstractFloat,
    forward::AbstractArray{<:Number}, #forward to each Asian observation t_i
    totalVariance::AbstractArray{<:Number}, #vol^2 * t_i
    weight::AbstractArray{<:Number},
)::Number where {TV <: Number}

    nAsset = length(totalVariance)
    spots = fill(spot, nAsset)
    correlation = zeros(eltype(totalVariance), (nAsset, nAsset))
    for (i, vi) in enumerate(totalVariance)
        for j = 1:(i-1)
            vj = totalVariance[j]
            if vi != 0 && vj != 0
                correlation[i, j] = min(vi, vj) / sqrt(vi * vj)
                # else is zero
            end
        end
        correlation[i, i] = one(Float64)
    end
    for i = 1:nAsset
        for j = i+1:nAsset
            correlation[i, j] = correlation[j, i]
        end
    end
    price = priceEuropean(p, isCall, strike, discountFactor, spots, forward, totalVariance, weight, correlation)
    return price
end

function priceAsianFloatingStrike(
    p::DeelstraBasketPricer,
    isCall::Bool,
    strikePercent::Number,
    discountFactor::Number,
    spot::Number,
    forward::Vector{<:Number}, #forward to each Asian observation t_i
    totalVariance::Vector{<:Number}, #vol^2 * t_i
    weight::Vector{<:Number},
)::Number
    #change of numeraire leads to equivalence with fixed strike put
    #floating S,k ,r,mu = k* fixed S/k,S,mu,r
    nAsset = length(totalVariance)
    dfAsset = discountFactor * forward[nAsset]
    forwardAsset = forward ./ forward[nAsset]
    #r,q = q,r => r,q-r = q-r,r
    tvarAsset = zeros(eltype(totalVariance), nAsset)
    correlation = zeros(eltype(totalVariance), (nAsset, nAsset))
    spots = ones(Float64, nAsset)
    stte = totalVariance[nAsset]
    for (i, vi) in enumerate(totalVariance)
        siti = stte - vi
        tvarAsset[i] = siti
        for j = 1:(i-1)
            sjtj = stte - totalVariance[j]
            if siti == zero(Float64) || sjtj == zero(Float64)
                correlation[i, j] = zero(Float64)
            else
                correlation[i, j] = min(siti, sjtj) / sqrt(siti * sjtj) #covar(xi,xj)/voli*volj = E(xi*xj - E(xi)*E(xj))/voli*volj =? min(ti-tj)/sqrt(ti)*sqrt(tj)
            end
        end
        correlation[i, i] = one(Float64)
    end
    for i = 1:nAsset
        for j = i+1:nAsset
            correlation[i, j] = correlation[j, i]
        end
    end

    priceFixed =
        strikePercent *
        priceEuropean(p, !isCall, spot / strikePercent, dfAsset, spots, forwardAsset, tvarAsset, weight, correlation)
    return priceFixed
end


text:
export decimation, samplesPerPeriod, periodsPerFrame, masterTrigger, currentFrame,
     currentPeriod, ramWriterMode, connectADC, startADC, stopADC, readData,
     numSlowDACChan, setSlowDACLUT, enableSlowDAC, currentWP, slowDACInterpolation,
     numSlowADCChan, numLostStepsSlowADC, bufferSize, keepAliveReset, triggerMode,
     slowDACPeriodsPerFrame, enableDACLUT

decimation(rp::RedPitaya) = query(rp,"RP:ADC:DECimation?", Int64)
function decimation(rp::RedPitaya, dec)
  rp.decimation = Int64(dec)
  send(rp, string("RP:ADC:DECimation ", rp.decimation))
end

samplesPerPeriod(rp::RedPitaya) = query(rp,"RP:ADC:PERiod?", Int64)
function samplesPerPeriod(rp::RedPitaya, value)
  rp.samplesPerPeriod = Int64(value)
  send(rp, string("RP:ADC:PERiod ", rp.samplesPerPeriod))
end

numSlowDACChan(rp::RedPitaya) = query(rp,"RP:ADC:SlowDAC?", Int64)
function numSlowDACChan(rp::RedPitaya, value)
  if value <= 0 || value > 4
    error("Num slow DAC channels needs to be between 1 and 4!")
  end
  send(rp, string("RP:ADC:SlowDAC ", Int64(value)))
end

function setSlowDACLUT(rp::RedPitaya, lut::Array)
  lutFloat32 = map(Float32, lut)
  send(rp, string("RP:ADC:SlowDACLUT"))
  @debug "Writing slow DAC LUT"
  write(rp.dataSocket, lutFloat32)
end

function enableDACLUT(rp::RedPitaya, lut::Array)
  lutBool = map(Bool, lut)
  send(rp, string("RP:ADC:EnableDACLUT"))
  @debug "Writing enable DAC LUT"
  write(rp.dataSocket, lutBool)
end

function enableSlowDAC(rp::RedPitaya, enable::Bool, numFrames::Int64=0,
            ffRampUpTime::Float64=0.4, ffRampUpFraction::Float64=0.8)
  enableI = Int32(enable)
  return query(rp, string("RP:ADC:SlowDACEnable ", enableI,
              ",", numFrames, ",", ffRampUpTime, ",", ffRampUpFraction), Int64)
end

function slowDACInterpolation(rp::RedPitaya, enable::Bool)
  enableI = Int32(enable)
  send(rp, string("RP:ADC:SlowDACInterpolation ", enableI))
end

numSlowADCChan(rp::RedPitaya) = query(rp,"RP:ADC:SlowADC?", Int64)
function numSlowADCChan(rp::RedPitaya, value)
  send(rp, string("RP:ADC:SlowADC ", Int64(value)))
end

numLostStepsSlowADC(rp::RedPitaya) = query(rp,"RP:ADC:SlowDACLostSteps?", Int64)

periodsPerFrame(rp::RedPitaya) = query(rp,"RP:ADC:FRAme?", Int64)
function periodsPerFrame(rp::RedPitaya, value)
  rp.periodsPerFrame = Int64(value)
  send(rp, string("RP:ADC:FRAme ", rp.periodsPerFrame))
end

slowDACPeriodsPerFrame(rp::RedPitaya) = query(rp,"RP:ADC:SlowDACPeriodsPerFrame?", Int64)
function slowDACPeriodsPerFrame(rp::RedPitaya, value)
  send(rp, string("RP:ADC:SlowDACPeriodsPerFrame ", value))
end

currentFrame(rp::RedPitaya) = query(rp,"RP:ADC:FRAMES:CURRENT?", Int64)
currentPeriod(rp::RedPitaya) = query(rp,"RP:ADC:PERIODS:CURRENT?", Int64)
currentWP(rp::RedPitaya) = query(rp,"RP:ADC:WP:CURRENT?", Int64)
bufferSize(rp::RedPitaya) = query(rp,"RP:ADC:BUFFER:SIZE?", Int64)

function masterTrigger(rp::RedPitaya, val::Bool)
  valStr = val ? "ON" : "OFF"
  send(rp, string("RP:MasterTrigger ", valStr))
end
masterTrigger(rp::RedPitaya) = occursin("ON", query(rp,"RP:MasterTrigger?"))

function keepAliveReset(rp::RedPitaya, val::Bool)
  valStr = val ? "ON" : "OFF"
  send(rp, string("RP:KeepAliveReset ", valStr))
end
keepAliveReset(rp::RedPitaya) = occursin("ON", query(rp,"RP:KeepAliveReset?"))

# "TRIGGERED" or "CONTINUOUS"
function ramWriterMode(rp::RedPitaya, mode::String)
  send(rp, string("RP:RamWriterMode ", mode))
end

# "INTERNAL" or "EXTERNAL"
function triggerMode(rp::RedPitaya, mode::String)
  send(rp, string("RP:Trigger:Mode ", mode))
end

connectADC(rp::RedPitaya) = send(rp, "RP:ADC:ACQCONNect")
startADC(rp::RedPitaya, wp::Integer = currentWP(rp)) = send(rp, "RP:ADC:ACQSTATUS ON,$wp")
stopADC(rp::RedPitaya) = send(rp, "RP:ADC:ACQSTATUS OFF,0")

# Low level read. One has to take care that the numFrames are available
function readData_(rp::RedPitaya, startFrame, numFrames)
  numSampPerPeriod = rp.samplesPerPeriod
  numPeriods = rp.periodsPerFrame
  numSampPerFrame = numSampPerPeriod * numPeriods

  command = string("RP:ADC:FRAMES:DATA ",Int64(startFrame),",",Int64(numFrames))
  send(rp, command)

  @debug "read data ..."
  u = read!(rp.dataSocket, Array{Int16}(undef, 2 * numFrames * numSampPerFrame))
  @debug "read data!"
  return reshape(u, 2, rp.samplesPerPeriod, numPeriods, numFrames)
end

# Low level read. One has to take care that the numFrames are available
function readDataPeriods_(rp::RedPitaya, startPeriod, numPeriods)
  command = string("RP:ADC:PERiods:DATa ",Int64(startPeriod),",",Int64(numPeriods))
  send(rp, command)

  @debug "read data ..."
  u = read!(rp.dataSocket, Array{Int16}(undef, 2 * numPeriods * rp.samplesPerPeriod))
  @debug "read data!"
  return reshape(u, 2, rp.samplesPerPeriod, numPeriods)
end

# High level read. numFrames can adress a future frame. Data is read in
# chunks
function readData(rp::RedPitaya, startFrame, numFrames, numBlockAverages=1, numPeriodsPerPatch=1)
  dec = rp.decimation
  numSampPerPeriod = rp.samplesPerPeriod
  numSamp = numSampPerPeriod * numFrames
  numPeriods = rp.periodsPerFrame
  numSampPerFrame = numSampPerPeriod * numPeriods

  if rem(numSampPerPeriod,numBlockAverages) != 0
    error("block averages has to be a divider of numSampPerPeriod")
  end

  numTrueSampPerPeriod = div(numSampPerPeriod,numBlockAverages*numPeriodsPerPatch)

  data = zeros(Float32, numTrueSampPerPeriod, 2, numPeriods*numPeriodsPerPatch, numFrames)
  wpRead = startFrame
  l=1

  numFramesInMemoryBuffer = bufferSize(rp) / numSampPerFrame
  @debug "numFramesInMemoryBuffer = $numFramesInMemoryBuffer"

  # This is a wild guess for a good chunk size
  chunkSize = max(1,  round(Int, 1000000 / numSampPerFrame)  )
  @debug "chunkSize = $chunkSize"
  while l<=numFrames
    wpWrite = currentFrame(rp)
    while wpRead >= wpWrite # Wait that startFrame is reached
      wpWrite = currentFrame(rp)
      @debug wpWrite
    end
    chunk = min(wpWrite-wpRead,chunkSize) # Determine how many frames to read
    @debug chunk
    if l+chunk > numFrames
      chunk = numFrames - l + 1
    end

    if wpWrite - numFramesInMemoryBuffer > wpRead
      @error "WARNING: We have lost data !!!!!!!!!!"
    end

    @debug "Read from $wpRead until $(wpRead+chunk-1), WpWrite $(wpWrite), chunk=$(chunk)"

    u = readData_(rp, Int64(wpRead), Int64(chunk))
    utmp1 = reshape(u,2,numTrueSampPerPeriod,numBlockAverages,size(u,3)*numPeriodsPerPatch,size(u,4))
    utmp2 = numBlockAverages > 1 ? mean(utmp1,dims=3) : utmp1

    data[:,1,:,l:(l+chunk-1)] = utmp2[1,:,1,:,:]
    data[:,2,:,l:(l+chunk-1)] = utmp2[2,:,1,:,:]

    l += chunk
    wpRead += chunk
  end

  return data
end

# Low level read. One has to take care that the numFrames are available
function readDataSlow_(rp::RedPitaya, startFrame, numFrames)
  numPeriods = rp.periodsPerFrame
  numChan = numSlowADCChan(rp)

  command = string("RP:ADC:SLOW:FRAMES:DATA ",Int64(startFrame),",",Int64(numFrames))
  #send(rp, command)

  @debug "read data ..."
  #u = read!(rp.dataSocket, Array{Float32}(undef, numChan * numFrames * numPeriods))
  @debug "read data!"
  #return reshape(u, numChan, numPeriods, numFrames)
  return zeros(Float32, numChan, numPeriods, numFrames)
end

max_stars_repo_name: fredrikekre/RedPitayaDAQServer

max_issues_repo_name: fuerpy/RedPitayaDAQServer

max_forks_repo_name: fuerpy/RedPitayaDAQServer


text:
import SeisIO: demean, demean!, taper, taper!,detrend, detrend!
export detrend, detrend!, taper, taper!, demean, demean!, hanningwindow
# Signal processing functions for arrays (rather than SeisData or SeisChannel)

"""
    detrend!(X::AbstractArray{<:AbstractFloat})

Remove linear trend from array `X` using least-squares regression.
"""
function detrend!(X::AbstractArray{<:AbstractFloat})
    T = eltype(X)
    N = size(X,1)

    # create linear trend matrix
    A = similar(X,T,N,2)
    A[:,2] .= T(1)
    A[:,1] .= range(T(0),T(1),length=N)
    # create linear trend matrix
    R = transpose(A) * A

    # do the matrix inverse for 2x2 matrix
    # this is really slow on GPU
    Rinv = inv(Array(R)) |> typeof(R)
    factor = Rinv * transpose(A)

    # remove trend
    X .-= A * (factor * X)
    return nothing
end
detrend(A::AbstractArray{<:AbstractFloat}) = (U = deepcopy(A);
        detrend!(U);return U)
detrend!(R::RawData) = detrend!(R.x)
detrend(R::RawData) = (U = deepcopy(R); detrend!(U.x); return U)
detrend!(C::CorrData) = detrend!(C.corr)
detrend(C::CorrData) = (U = deepcopy(C); detrend!(U.corr); return U)

"""
    demean!(A::AbstractArray{<:AbstractFloat})

Remove mean from array `A`.
"""
function demean!(A::AbstractArray{<:AbstractFloat}; dims=1)
      μ = mean(A,dims=dims)
      A .-= μ
  return nothing
end
demean(A::AbstractArray{<:AbstractFloat}; dims=1) = (U = deepcopy(A);
       demean!(U,dims=dims);return U)
demean!(R::RawData) = demean!(R.x)
demean(R::RawData) = (U = deepcopy(R); demean!(U.x); return U)
demean!(C::CorrData) = demean!(C.corr)
demean(C::CorrData) = (U = deepcopy(C); demean!(U.corr); return U)

"""
   taper!(A,fs; max_percentage=0.05, max_length=20.)

Taper a time series `A` with sampling_rate `fs`.
Defaults to 'hann' window. Uses smallest of `max_percentage` * `fs`
or `max_length`.

# Arguments
- `A::AbstractArray`: Time series.
- `fs::AbstractFloat`: Sampling rate of time series `A`.
- `max_percentage::float`: Decimal percentage of taper at one end (ranging
   from 0. to 0.5).
- `max_length::Real`: Length of taper at one end in seconds.
"""
function taper!(A::AbstractArray{<:AbstractFloat}, fs::Real;
                max_percentage::AbstractFloat=0.05, max_length::Real=20.)
   Nrows = size(A,1)
   T = eltype(A)
   wlen = min(Int(floor(Nrows * max_percentage)), Int(floor(max_length * fs)), Int(
         floor(Nrows/2)))
   taper_sides = hanningwindow(A,2 * wlen)
   A[1:wlen,:] .*= taper_sides[1:wlen]
   A[end-wlen:end,:] .*= taper_sides[wlen:end]
   return nothing
end
taper(A::AbstractArray{<:AbstractFloat}, fs::Real;
       max_percentage::AbstractFloat=0.05,max_length::Real=20.) = (U = deepcopy(A);
       taper!(U,fs,max_percentage=max_percentage,max_length=max_length);return U)
taper!(R::RawData; max_percentage::AbstractFloat=0.05,
       max_length::Real=20.) = taper!(R.x,R.fs,max_percentage=max_percentage,
       max_length=max_length)
taper(R::RawData; max_percentage::AbstractFloat=0.05,
       max_length::Real=20.) = (U = deepcopy(R); taper!(U.x,U.fs,
       max_percentage=max_percentage,max_length=max_length); return U)
taper!(C::CorrData; max_percentage::AbstractFloat=0.05,
      max_length::Real=20.) = taper!(C.corr,C.fs,max_percentage=max_percentage,
      max_length=max_length)
taper(C::CorrData; max_percentage::AbstractFloat=0.05,
      max_length::Real=20.) = (U = deepcopy(C); taper!(U.corr,U.fs,
      max_percentage=max_percentage,max_length=max_length); return U)

"""
  hanningwindow(A,n)

Generate hanning window of length `n`.

Hanning window is sin(n*pi)^2.
"""
function hanningwindow(A::AbstractArray, n::Int)
   T = eltype(A)
   win = similar(A,T,n)
   win .= T(pi) .* range(T(0.),stop=T(n),length=n)
   win .= sin.(win).^2
   return win
end


text:
using Revise
using VectorDataUtils
using Statistics
using BenchmarkTools
using Distances

npts_per_cluster = 4
r = 0.25
x = VectorDataUtils.generate_4_circle_clusters(npts_per_cluster; r)
cntrs = VectorDataUtils.kmeans(x, 4)
cntrs = VectorDataUtils.kmeans(x, 4)
@time cntrs = VectorDataUtils.kmeans(x, 4);

ndim = 20
npts = 1000
x = randn(Float64,ndim,npts)
y = randn(Float64,ndim,npts)

@time VectorDataUtils.max_centre_shift(x,y)
@time VectorDataUtils.max_centre_shift2(x,y)

@benchmark m = VectorDataUtils.mindist2cntrs(x,y)

d1 = Matrix{Float64}(undef,npts,npts)
d2 = Matrix{Float64}(undef,npts,npts)

@benchmark VectorDataUtils.pairwise!(d1,sql2(),x,y)
@benchmark Distances.pairwise!(d2,SqEuclidean(),x,y)

@benchmark d1 = euclidean(x[:,1],y[:,1])
@benchmark d2 = l2dist(x[:,1],y[:,1])

d3 = evaluate(Euclidean(),x[:,1],y[:,1])

tt = pairwise(Euclidean(), x, y, dims=2)


text:
using Dates
using Test
using PowerSystemsUnits

using Unitful: Unitful, @u_str, J, W, hr

@testset "PowerSystemsUnits.jl" begin
    @testset "Units" begin
        @testset "Period Time Types" begin
            time_stamps = collect(now()-Dates.Day(1):Dates.Hour(1):now())
            offsets = time_stamps - time_stamps[1]

            x = PowerSystemsUnits.dt2umin(Dates.Minute.(offsets))
            y = (offsets)u"minute"
            @test x == y
        end
        @testset "Power Related Unit" begin
            @test 1u"Wh" == 1u"W*hr"
            @test 1u"VARh" == 3600u"J"
            @test 1u"USDPerMWh" == 1u"USD/(MW*hr)"
        end
    end
end

max_stars_repo_name: invenia/PowerSystemsUnits

max_issues_repo_name: invenia/PowerSystemsUnits

max_forks_repo_name: invenia/PowerSystemsUnits


text:
"""
    SubsetProperty{N<:Real} <: Property

Type that represents a safety property characterized with a set of safe states.
A safety property is satisfied by a given set of states if the set of states is
fully contained in the set of safe states.

### Fields

- ``safe`` -- convex set representing the safe states
- ``witness`` -- witness point (empty vector if not set)
"""
mutable struct SubsetProperty{N<:Real} <: Property
    safe::LazySet
    witness::Vector{N}

    SubsetProperty{N}(safe::LazySet) where {N<:Real} = new(safe, N[])
end

# type-less convenience constructor
SubsetProperty(safe::LazySet{N}) where {N} = SubsetProperty{N}(safe)

"""
    inout_map_property(prop::SubsetProperty,
                       partition::AbstractVector{<:AbstractVector{Int}},
                       blocks::AbstractVector{Int},
                       n::Int
                      )::SubsetProperty

Map an `SubsetProperty` to the dimensions of analyzed blocks.

### Input

- `prop`      -- property
- `partition` -- block partition; elements are start and end indices of a block
- `blocks`    -- list of all output block indices in the partition
- `n`         -- total number of input dimensions

### Output

A new property of reduced dimension.

### Notes

If the dimension is not reduced, we keep the original set.
Otherwise, the dimension reduction is achieved with a `LinearMap`.
"""
function inout_map_property(prop::SubsetProperty{N},
                            partition::AbstractVector{<:AbstractVector{Int}},
                            blocks::AbstractVector{Int},
                            n::Int
                           )::SubsetProperty{N} where {N<:Real}
    @assert dim(prop.safe) == n "the property has dimension $(dim(prop.safe)) but should have dimension $n"
    proj = projection_map(partition, blocks)
    if length(proj) == n
        # no change in the dimension, copy the old property (keep the set)
        return IntersectionProperty(prop.safe)
    else
        M = sparse(proj, proj, ones(N, length(proj)), n, n)
        return IntersectionProperty(M * prop.safe)
    end
end

"""
    check_property(set::LazySet, prop::SubsetProperty)::Bool

Checks whether a convex set satisfies a safety property.

### Input

- ``set``  -- convex set
- ``prop`` -- safety property with safe states

### Output

`true` iff the safety property is satisfied for the given set of states.
This is the case iff the set of states is a subset of the safe states.
"""
@inline function check_property(set::LazySet, prop::SubsetProperty)::Bool
    is_subset, witness = ⊆(set, prop.safe, true)
    if !is_subset
        # store violation witness
        prop.witness = witness
    end
    return is_subset
end


text:
## External interface

# Main function
"""
    heppcat(Y,k,iters;init=homppca(Y,k))

Estimate probabilistic PCA model for noise that is heteroscedastic across samples.

Inputs are:
+ `Y` : list of matrices (each column is a sample)
+ `k` : number of factors
+ `iters` : number of iterations to run
+ `init`  : initial model (will be modified in-place)
Output is a [`HePPCATModel`](@ref) object.
"""
function heppcat(Y,k,iters;init=homppca(Y,k))
    M = init
    for _ in 1:iters
        updatev!(M,Y,ExpectationMaximization())
        updateF!(M,Y,ExpectationMaximization())
    end
    M
end

# Homoscedastic initialization
"""
    homppca(Y,k)

Estimate probabilistic PCA model for noise that is homoscedastic across samples.

Inputs are:
+ `Y` : list of matrices (each column is a sample)
+ `k` : number of factors
Output is a [`HePPCATModel`](@ref) object.
"""
function homppca(Y,k)
    Yf = reduce(hcat,Y)
    n, L = size(Yf,2), length(Y)
    Uh, s, _ = svd(Yf)
    λh = abs2.(s)./n
    λb = mean(λh[k+1:end])
    HePPCATModel(Uh[:,1:k],λh[1:k] .- λb,I(k),fill(λb,L))
end

# log-likelihood (todo: add constant)
"""
    loglikelihood(M::HePPCATModel,Y)

Log-likelihood of model `M` with respect to data `Y` (dropping constant term).
"""
function loglikelihood(M::HePPCATModel,Y)
    d, k = size(M.U)
    n, L = size.(Y,2), length(Y)
    return 1/2*sum(1:L) do l
        norm(sqrt(Diagonal((M.λ./M.v[l])./(M.λ .+ M.v[l])))*M.U'Y[l])^2 -
        n[l]*sum(log.(M.λ .+ M.v[l])) - n[l]*(d-k)*log(M.v[l]) - norm(Y[l])^2/M.v[l]
    end
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Functions to sample random representations of rotations.
#
# References
# ==============================================================================
#
#   [1] K. Shoemake. Uniform random rotations. In D. Kirk, editor, Graphics Gems
#       III, pages 124-132. Academic, New York, 1992.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

# Euler angle and axis
# ==============================================================================

function Random.rand(rng::AbstractRNG, ::Random.SamplerType{R}) where R <: EulerAngleAxis
    T = eltype(R)
    if T == Any
        T = Float64
    end

    return quat_to_angleaxis(_rand_quat(rng, T))
end

# Euler angles
# ==============================================================================

function Random.rand(rng::AbstractRNG, ::Random.SamplerType{R}) where R <: EulerAngles
    T = eltype(R)
    if T == Any
        T = Float64
    end

    k  = 2 * T(π)
    a₁ = k * rand(rng, T)
    a₂ = k * rand(rng, T)
    a₃ = k * rand(rng, T)

    rot_seq = rand(
        rng,
        (:XYX, :XYZ, :XZX, :XZY, :YXY, :YXZ, :YZX, :YZY, :ZXY, :ZXZ, :ZYX, :ZYZ)
    )

    return EulerAngles{T}(a₁, a₂, a₃, rot_seq)
end

# DCM
# ==============================================================================

function Random.rand(rng::AbstractRNG, ::Random.SamplerType{R}) where R <: DCM
    T = eltype(R)
    if T == Any
        T = Float64
    end

    return quat_to_dcm(_rand_quat(rng, T))
end

# Quaternion
# ==============================================================================

function Random.rand(rng::AbstractRNG, ::Random.SamplerType{R}) where R <: Quaternion
    T = eltype(R)
    if T == Any
        T = Float64
    end

    return _rand_quat(rng, T)
end

#                              Private functions
# ==============================================================================

@inline function _rand_quat(rng::AbstractRNG, T::DataType)
    # A random unit quaternion is sampled using the algorithm in [1].
    k = 2T(π)

    u₁ = rand(rng, T)
    u₂ = k * rand(rng, T)
    u₃ = k * rand(rng, T)

    sin_u₂, cos_u₂ = sincos(u₂)
    sin_u₃, cos_u₃ = sincos(u₃)

    v₁ = √(1 - u₁)
    v₂ = √u₁

    q = Quaternion(
        v₁ * sin_u₂,
        v₁ * cos_u₂,
        v₂ * sin_u₃,
        v₂ * cos_u₃
    )

    return q
end


text:
using CUDAapi

config_path = joinpath(@__DIR__, "ext.jl")
const previous_config_path = config_path * ".bak"

function write_ext(config)
    open(config_path, "w") do io
        println(io, "# autogenerated file, do not edit")
        for (key,val) in config
            println(io, "const $key = $(repr(val))")
        end
    end
end

function main()
    ispath(config_path) && mv(config_path, previous_config_path; force=true)
    config = Dict{Symbol,Any}(:configured => false)
    write_ext(config)

    ## discover stuff
    toolkit_dirs = CUDAapi.find_toolkit()
    # On Windows everything is sad and different
    if Sys.iswindows()
        path = "C:\\Program Files\\NVIDIA Corporation\\NvToolsExt"
        suffix = Sys.ARCH == :x86_64 ? "x64" : "Win32"
        push!(toolkit_dirs, joinpath(path, "bin", suffix))
    end

    config[:libnvtx] = CUDAapi.find_cuda_library("nvToolsExt", toolkit_dirs)
    if config[:libnvtx] == nothing
        error("could not find NVTX")
    end

    ## (re)generate ext.jl

    function globals(mod)
        all_names = names(mod, all=true)
        filter(name-> !any(name .== [nameof(mod), Symbol("#eval"), :eval]), all_names)
    end

    if isfile(previous_config_path)
        @debug("Checking validity of existing ext.jl...")
        @eval module Previous; include($previous_config_path); end
        previous_config = Dict{Symbol,Any}(name => getfield(Previous, name)
                                           for name in globals(Previous))

        if config == previous_config
            info("CuArrays.jl has already been built for this toolchain, no need to rebuild")
            mv(previous_config_path, config_path; remove_destination=true)
            return
        end
    end

    config[:configured] = true
    write_ext(config)

    return
end

main()


text:
function do_ast()
  r = RethinkDB

  #=
  r.db_create("sample_db") |>
    println

  r.db("sample_db") |>
    d -> r.table_create(d, "authors") |>
    println

    cursor = r.db("sample_db") |>
      d -> r.table(d, "authors") |>
      d -> r.insert(d, [
        { "name"=>"William Adama", "tv_show"=>"Battlestar Galactica",
          "posts"=>[
            {"title"=>"Decommissioning speech", "content"=>"The Cylon War is long over..."},
            {"title"=>"We are at war", "content"=>"Moments ago, this ship received..."},
            {"title"=>"The new Earth", "content"=>"The discoveries of the past few days..."}
          ]
        },
        { "name"=>"Laura Roslin", "tv_show"=>"Battlestar Galactica",
          "posts"=>[
            {"title"=>"The oath of office", "content"=>"I, Laura Roslin, ..."},
            {"title"=>"They look like us", "content"=>"The Cylons have the ability..."}
          ]
        },
        { "name"=>"Jean-Luc Picard", "tv_show"=>"Star Trek TNG",
          "posts"=>[
            {"title"=>"Civil rights", "content"=>"There are some words I've known since..."}
          ]
        }
      ]) |> println

  r.db("sample_db") |>
    d -> r.table(d, "authors") |>
    d -> r.get_all(d) |>
    d -> r.filter(d, {"name" => "William Adama"}) |> println

  =#
  #=
  ruby
  [1,
    [39,[ // filter
      [15,["authors"]], // table
      [69,[ // func
        [2,[29]], // Array
        [17,[ // eq
          [170,[ // Bracket
            [10,[29]], // var
          "name"]],
        "William Adama"]]
      ]]
    ]],
  {}]

  js
  [1,
    [39,[ // filter
      [15,["authors"]], // table
      [69,[ // func
        [2,[0]], // Array
        [17,[ // eq
          [170,[ // Bracket
            [13,[]], // IMPLICIT_VAR
          "name"]],
        "William Adama"]]
      ]]
    ]]
  ]

  [39,[15,"authors"],[69,[],[17,[null,"name"],"William Adama"]]]
  =#

  println("")

  r.table("authors") |>
    d -> r.filter(d,
      r.func0(
        r.row("name") |> d -> r.eq(d, "William Adama")
      )
    ) |> to_ast |> JSON.json |> println

  #=
  r.table("marvel").get("IronMan").bracket("firstAppearance")

  [1,
    [170,[ // Bracket
      [16,[ // get
        [15,["marvel"]], // Table
      "IronMan"]],
    "firstAppearance"]]
  ]
  =#
  #r.table("marvel") |>
  #  d -> r.get(d, "IronMan") |>
  #  d -> r.bracket(d, "firstAppearance") |> println

  #=
  r.row("authors")
  [1,
    [170,[ // Bracket
      [13,[]], // IMPLICIT_VAR
    "authors"]]
  ]
  =#

  #r.implicit_var() |>
  #  d -> r.bracket(d, "authors") |> println

  #r.row("authors") |> println

end

function isarray(a::Array)
  true
end

function isarray(a)
  false
end

function isreql(a)
  try
    reql_type(typeof(a))
    true
  catch e
    false
  end
end

#=

[1,
  [39,[ // filter
    [15,["authors"]], // table
    [69,[ // func
      [2,[0]], // Array
      [17,[ // eq
        [170,[ // Bracket
          [13,[]], // IMPLICIT_VAR
        "name"]],
      "William Adama"]]
    ]]
  ]]
]

RethinkDB.ReqlSequence(39, // push()
  Any[
    RethinkDB.ReqlSequence(-1,
      RethinkDB.ReqlTable(15, Any[RethinkDB.ReqlString(-1, "authors")])
    ),

    RethinkDB.ReqlFunction1(-1,
      RethinkDB.ReqlTop(69,
        Any[
          RethinkDB.ReqlArray(2,Any[]),
          RethinkDB.ReqlTop(-1,
            RethinkDB.ReqlBool(17,
              Any[
                RethinkDB.ReqlDatum(170,
                  Any[
                    RethinkDB.ReqlDatum(13,Any[]),
                    RethinkDB.ReqlString(-1, "name")
                  ]),
                RethinkDB.ReqlDatum(-1, "William Adama")
              ]
            )
          )
        ]
      )
    )
  ]
)
=#

function to_ast(ast)
  local ret = []
  push!(ret, 1)
  push!(ret, _to_ast(ast))
  ret
end

function _to_ast(ast)

  if (isarray(ast))
    local tmp = []
    for i in ast
      push!(tmp, _to_ast(i))
    end
    return tmp
  end

  if (!isreql(ast))
    return ast
  end

  if ast.op_code != -1
    local tmp_op_code = []
    push!(tmp_op_code, ast.op_code)
    push!(tmp_op_code, _to_ast(ast.value))
    return tmp_op_code
  else
    return _to_ast(ast.value)
  end
end

max_stars_repo_name: mrcsparker/julia-rethinkdb

max_issues_repo_name: mrcsparker/julia-rethinkdb

max_forks_repo_name: mrcsparker/julia-rethinkdb