text:
"""
    drop1(m::DataFrameRegressionModel, terms::Vector{Symbol}=Symbol[])

Perform a `drop1` analysis of deviance on `terms` in model `m`

Model `m` is modified and refit dropping the terms in `terms` one at a time.
An analysis of deviance is applied to the refit model.
"""
function drop1(m, terms::Vector{Symbol}=Symbol[])
    form = Formula(m.mf.terms)
    basedev = deviance(m)
    deviances, degsfreedom, Χ², pvalues, nms = [basedev], [0], [0.], [NaN], ["<none>"]
    ncoef = length(coef(m))
    dat = m.mf.df
    dinst = m.model.rr.d
    if length(terms) == 0
        terms = filter(x -> isa(x, Symbol) || isa(x, Expr), form.rhs.args[2:end])
    end
    for trm in terms
        newform = dropterm(form, trm)
        newm = fit(GeneralizedLinearModel, newform, dat, dinst, Link(m))
        newdev = deviance(newm)
        newdof = ncoef - length(coef(newm))
        newΧ² = newdev - basedev
        push!(deviances, newdev)
        push!(degsfreedom, newdof)
        push!(Χ², newΧ²)
        push!(pvalues, ccdf(Chisq(newdof), newΧ²))
        push!(nms, string(trm))
    end
    CoefTable([deviances, Χ², degsfreedom, pvalues], ["Deviance", "Χ²", "degfree", "P[> Χ²]"], nms)
end

function StatsBase.counts(v::Union{CategoricalVector,PooledDataVector})
    levs = isa(v, CategoricalArray) ? CategoricalArrays.levels(v) : DataArrays.levels(v)
    counts = zeros(Int, length(levs) + 1)
    for ind in v.refs
        counts[ind + 1] += 1
    end
    counts[1] == 0 ? shift!(counts) : unshift!(levs, "#NULL")
    result = NamedArray(counts)
    NamedArrays.setnames!(result, levs, 1)
    result
end
GLM.Link{T<:GeneralizedLinearModel}(m::DataFrames.DataFrameRegressionModel{T}) = Link(m.model)

"""
    dropterm(f::Formula, trm::Symbol)

Return a copy of `f` without the term `trm`.

# Examples
```jl
julia> dropterm(foo ~ 1 + bar + baz, :bar)
foo ~ 1 + baz
```
"""
function dropterm(f::Formula, trm::Symbol)
    rhs = copy(f.rhs)
    args = rhs.args
    if !(Meta.isexpr(rhs, :call) && args[1] == :+ && (tpos = findlast(args, trm)) > 0)
        throw(ArgumentError("$trm is not a summand in `$rhs`"))
    end
    Formula(f.lhs, Expr(:call, :+, deleteat!(args, [1, tpos])...))
end

"""
    checksingularity(form::Formula, data::DataFrame, tolerance = 1.e-8)

Return a vector of terms in `form` that produce singularity in the model matrix
"""
function checksingularity(form::Formula, data::DataFrame, tolerance = 1.e-8)
    mf = ModelFrame(form, data)
    mm = ModelMatrix(mf)
    qrf = qrfact!(mm.m, Val{true})
    vals = abs.(diag(qrf[:R]))
    firstbad = findfirst(x -> x < min(tolerance, 0.5) * vals[1], vals)
    if firstbad == 0
        return Symbol[]
    end
    mf.terms.terms[view(mm.assign[qrf[:p]], firstbad:length(vals))]
end

"""
    GLMformula(dd::Dict{String,Any})

Create the GLM formula from a configuation dictionary for the model
"""
function GLMformula(dd::Dict)
    haskey(dd, "y_var") || throw(ArgumentError("dd should have a \"y_var\" key"))
    haskey(dd, "finalvars") || throw(ArgumentError("dd should have a \"finalvars\" key"))
    rhs = :(1 + x)
    pop!(rhs.args)
    append!(rhs.args, Symbol.(dd["finalvars"]))
    Formula(Symbol(dd["y_var"]), rhs)
end

"""
    GLMMform(form::Formula, raneff::Symbol)

Replace the last term in the formula by :(1 | ranef)
"""
function GLMMform(form::Formula, raneff::Symbol)
    retrm = :(1 | x)
    retrm.args[end] = raneff
    form.rhs.args[end] = retrm
    form
end

"""
    profileθ(m::GeneralizedLinearMixedModel)

Profile the deviance as a function of θ for models with a single, scalar r.e. term
"""
function profileθ{T}(m::GeneralizedLinearMixedModel{T})
    θ = getθ(m)
    isa(θ, Vector{T}) && length(θ) == 1 ||
        throw(ArgumentError("m must have a single scalar r.e. term"))
    th = zeros(T, 1)
    δ = zeros(T, 1)
    fill!(θ, 0)
    decreasing = false
    basedev = pirls!(setθ!(m, θ), true)
    for t in 0.001:0.001:2.0
        θ[1] = t
        dd = pirls!(setθ!(m, θ), true) - basedev
        sgn = sign(dd - δ[end])
        push!(δ, dd)
        push!(th, t)
        dd ≤ 4.0 || break
        decreasing && sgn > 0 && break
        decreasing = sgn < 0
    end
    th, δ
end

"""
    relevel{T,I}(v::PooledDataVector{T,I}, lev::T)

Create a new `PooledDataVector{T,I}` with `lev` as the first (i.e. reference) level
"""
function relevel{T,I}(v::PooledDataVector{T,I}, lev::T)
    pool = v.pool
    i = convert(I, findfirst(pool, lev))
    i ≠ 0 || throw(ArgumentError("\"$lev\" is not in v.pool"))
    if i == 1
        return copy(v)
    end
    orig = convert(Vector{I}, 1:length(pool))
    perm = append!([i], deleteat!(orig, i))
    PooledDataArray(DataArrays.RefArray(invperm(perm)[v.refs]), v.pool[perm])
end

function cov2cor(m)
    d = Diagonal([inv(sqrt(m[i, i])) for i in 1:size(m, 2)])
    A_mul_B!(d, m * d)
end

"""
    vif(m::DataFrameRegressionModel)

The variance inflation factors for the terms in a model.
"""
function vif(m)
    v = vcov(m)
    assign = m.mm.assign
    nms = coefnames(m.mf)
    if (ipos = findfirst(nms, "(Intercept)")) > 0
        inds = deleteat!(collect(eachindex(nms)), ipos)
        v = view(v, inds, inds)
        assign = view(assign, inds)
    else
        warning("No intercept: vifs may not be sensible.")
    end
    terms = convert(Vector{Symbol}, m.mf.terms.terms)
    n_terms = length(terms)
    if n_terms < 2
        error("model contains fewer than 2 terms")
    end
    R = cov2cor(v)
    detR = det(Symmetric(R))
    result = NamedArray(Array(eltype(R), (length(terms), 3)))
    setnames!(result, string.(terms), 1)
    setnames!(result, ["GVIF", "Df", "GVIF^(1/(2*Df))"], 2)
    inds = collect(eachindex(terms))
    ainds = collect(eachindex(assign))
    for i in inds
        subs = find(x -> x == i, assign)
        notsubs = setdiff(ainds, subs)
        result[i, 1] = det(Symmetric(R[notsubs, notsubs])) / detR
        if (result[i, 2] = length(subs)) > 1
            result[i, 1] *= det(Symmetric(R[subs, subs]))
        end
    end
    if all(result[:, 2] .== 1)
        return result[:, 1:1]
    end
    result
end

max_stars_repo_name: werbac/jlib

max_issues_repo_name: werbac/jlib

max_forks_repo_name: werbac/jlib


text:
export rq_fact, qr_fact

function qr_fact(M::AbstractMatrix, Dcut::Int=typemax(Int), tol::Float64=1E-12, args...)
    F = qr(M, args...)
    q, r = _qr_fix(Array(F.Q), Array(F.R))
    if Dcut > size(q, 2) return q, r end
    U, Σ, V = svd(r, Dcut, tol)
    q * U, Diagonal(Σ) * V'
end

function rq_fact(M::AbstractMatrix, Dcut::Int=typemax(Int), tol::Float64=1E-12, args...)
    q, r = qr_fact(M', Dcut, tol, args...)
    r', q'
end

function _qr_fix(Q::T, R::AbstractMatrix) where {T <: AbstractMatrix}
    d = diag(R)
    for i ∈ eachindex(d)
        @inbounds d[i] = ifelse(isapprox(d[i], 0, atol=1e-14), 1, d[i])
    end
    ph = d ./ abs.(d)
    Q * Diagonal(ph), Diagonal(ph) * R
end

function LinearAlgebra.svd(A::AbstractMatrix, Dcut::Int=typemax(Int), tol::Float64=1E-12, args...)
    #U, Σ, V = psvd(A, rank=Dcut, args...)

    U, Σ, V = svd(A, args...)

    tol = Σ[1] * max(eps(), tol)
    δ = min(Dcut, sum(Σ .> tol))

    U = U[:, 1:δ]
    Σ = Σ[1:δ]
    Σ ./ sum(Σ .^ 2)
    V = V[:, 1:δ]

    d = diag(U)
    for i ∈ eachindex(d)
        @inbounds d[i] = ifelse(isapprox(d[i], 0, atol=1e-14), 1, d[i])
    end
    ph = d ./ abs.(d)
    U * Diagonal(ph), Σ, V * Diagonal(ph)
end

ext: jl


text:
import ArchGDAL
function geolocate(
    slc::AbstractArray,
    lons::AbstractArray,
    lats::AbstractArray,
    step_deg = 0.1,
)
    geolons, geolats = coarse_grid_deg(lons, lats, step_deg)
    geoslc = zeros(eltype(slc), size(geolats, 1), size(geolons, 1))
    # Add counts to take avg when multiple pixels fall in same output bucket
    counts = zeros(size(geoslc))

    println("Using step size $step_deg")
    println("Input/output sizes:")
    @show size(lons), size(lats), size(slc)
    @show size(geolons), size(geolats), size(geoslc)
    for idx in eachindex(lons, lats, slc)
        (lon, lat, val) = lons[idx], lats[idx], slc[idx]
        row = nearest(geolats, lat)
        col = nearest(geolons, lon)
        oob(row, col, geoslc) && continue
        geoslc[row, col] += val
        counts[row, col] += 1
    end

    return (geoslc ./ counts), geolons, geolats
end

function coarse_grid_deg(lons, lats, step)
    lon1, lon2 = extrema(lons)
    lat1, lat2 = extrema(lats)
    return lon1:step:lon2, lat2:(-step):lat1
end

# if we have just an array (with even step size) and point:
nearest(arr::AbstractArray, point) = Int(round((point - first(arr)) / (arr[2] - arr[1])))

# Out of bounds check (if we decide to make output area smaller than input)
oob(row, col, out) = row < 1 || row > size(out, 1) || col < 1 || col > size(out, 2)

function loadfiles()
    rows, cols = 7319, 4720
    slc = Array{ComplexF32,2}(undef, (cols, rows))
    lons = Array{Float32,2}(undef, (cols, rows))
    lats = similar(lons)

    read!("10913_11584.int", slc)
    slc = permutedims(slc)

    read!("lon", lons)
    lons = permutedims(lons)

    read!("lat", lats)
    lats = permutedims(lats)
    return (slc, lats, lons)
end

function write_geotiff(geoslc, geolons, geolats, outname = "geoslc.tif")
    height, width = size(geoslc)

    ArchGDAL.create(
        "geoslc.tif",
        driver = ArchGDAL.getdriver("GTiff"),
        width = width,
        height = height,
        nbands = 1,
        dtype = Float32,
    ) do dataset
        # ArchGDAL expects the dimensions of the buffer to be (cols, rows, bands) or (cols, rows).
        ArchGDAL.write!(dataset, permutedims(Float32.(log10.(abs.(geoslc)))), 1)
    end

    lry, uly = extrema(geolats)
    ulx, lrx = extrema(geolons)
    run(`gdal_edit.py geoslc.tif -a_srs EPSG:4326 -a_ullr $ulx $uly $lrx $lry`)
    # xres = geolons[2] - geolons[1]
    # yres = geolats[2] - geolats[1]
    # run(`gdal_edit.py geoslc.tif -a_srs EPSG:4236 -tr $xres $yres -a_ullr $ulx $uly $lrx $lry
end

function testgeocode(step_deg = 0.0002)
    # step_deg = 0.001  # Leads to 935x892
    slc, lats, lons = loadfiles()
    geoslc, geolons, geolats = geolocate(slc, lons, lats, step_deg)
    write_geotiff(geoslc, geolons, geolats, outname = "geoslc.tif")
end


text:
using DataFrames

df = datatable(; id=1:10, group=repeat([1, 2], 5), age=12)
mat = as_matrix(df)

@testset "cbind" begin
    mat = as_matrix(df)
    @test typeof(cbind(mat, mat, mat)) == Matrix{Int64}
    @test_nowarn cbind(id=1:10, df, df)
    @test_nowarn cbind(id=1:10, df, as_matrix(df))
    @test_nowarn cbind(id=1:10, as_matrix(df), df)
end

@testset "rbind" begin
    @test typeof(rbind(mat, mat, mat)) == Matrix{Int64}
    @test_nowarn rbind(df, mat)
    @test_nowarn rbind(mat, df)
    @test_nowarn rbind(df, mat, mat, df)
end

@testset "melt_list" begin
    probs = [0.90, 0.95, 0.99, 0.999, 0.9999]
    levs = factor(probs)

    d = DataFrame(x=1:40)
    names(d)

    d2 = cbind(d, prob=probs[1])
    @test names(d2) == ["x", "prob"]

    ## for data.frame melt_list
    res = []
    for i = 1:4
        push!(res, d)
    end

    df1 = melt_list(res, id=1:4)
    df2 = melt_list(res)
    # test for empty list
    push!(res, [])
    df3 = melt_list(res)
    @test names(df2) == ["x", "prob", "id", "I"]
end

# @pipe d |> groupby(:y)
# `>(6)` has to in this format
@pipe df |>
      dropmissing |>
      filter(:id => >(6), _) |>
      groupby(:group) |>
      combine(:age => sum)

@pipe df |>
      dropmissing |>
      filter(:id => >(6), _)

@testset "@subset" begin
    x = 1:2
    y = 2
    dt = datatable(; x, y=[2, 3], z=[1, 3])

    @test typeof(dt) == DataFrame
    @test @subset(dt, y == 2) |> nrow == 1
    @test @subset(dt, y == 1) |> nrow == 0
    @test @subset(dt, y == 2 & z == 1) |> nrow == 1

    # works for global variable
    y0 = 2
    @test @pipe(dt |> @subset(y == 2, true)) == @pipe(dt |> @subset(y == y0, true))
end


text:
@doc raw"""
`HODE`: Hamiltonian Ordinary Differential Equation *EXPERIMENTAL*

Defines a Hamiltonian ordinary differential initial value problem, that is
a canonical Hamiltonian system of equations,
```math
\begin{align*}
\dot{q} (t) &= v(t, q(t), p(t)) , & q(t_{0}) &= q_{0} , \\
\dot{p} (t) &= f(t, q(t), p(t)) , & p(t_{0}) &= p_{0} ,
\end{align*}
```
with vector fields ``v`` and ``f``, given by
```math
\begin{align*}
v &=   \frac{\partial H}{\partial p} , &
f &= - \frac{\partial H}{\partial q} ,
\end{align*}
```
initial conditions ``(q_{0}, p_{0})`` and the dynamical variables ``(q,p)``
taking values in ``\mathbb{R}^{d} \times \mathbb{R}^{d}``.

### Fields

* `d`: dimension of dynamical variables ``q`` and ``p`` as well as the vector fields ``v`` and ``f``
* `n`: number of initial conditions
* `v`: function computing the vector field ``v``
* `f`: function computing the vector field ``f``
* `h`: function computing the Hamiltonian ``H``
* `t₀`: initial time
* `q₀`: initial condition for dynamical variable ``q``
* `p₀`: initial condition for dynamical variable ``p``

"""


text:
"""
    causality(x::AbstractUncertainIndexValueDataset,
        y::AbstractUncertainIndexValueDataset,
        test::BinnedDataCausalityTest)

Apply a causality test to `x` and `y`, which are both data series
with uncertainties in both indices and values. To get the data
on an equally-spaced temporal grid, the data are first binned
according to the instructions in `test.binning`.

## Binning methods

- If `test.binning` results in an uncertain value for each bin, then the causality
test is applied `test.n_realizations` times. Examples are [`UncertainData.BinnedResampling`](@ref)
and [`UncertainData.BinnedWeightedResampling`](@ref), which both return a KDE estimate to the
distribution of values in each bin.

- If `test.binning` returns a summary statistic for each bin, then the causality test is
applied to the summarised time series. Examples are [`UncertainData.BinnedMeanResampling`](@ref)
and [`UncertainData.BinnedMeanWeightedResampling`](@ref), which both return the mean of each bin.

## Example

### Some example data

Let's say we have the following datasets `X` and `Y`:

```julia
using CausalityTools, UncertainData

sys = ar1_unidir(uᵢ = [0.1, 0.1], c_xy = 0.41)
vars = (1, 2) # ar1_unidir has only two variables, X and Y
n_steps = 100 # the number of points in the time series
tstep = 10 # the mean of each time value is stepped by `tstep`

X, Y = example_uncertain_indexvalue_datasets(sys, n_steps, vars, tstep = 10,
d_xind = Uniform(7.5, 15.5),
d_yind = Uniform(5.5, 15.5),
d_xval = Uniform(0.1, 0.5));
```

Now we can perform the causality tests on our uncertain data, either over
the bin means, or over an ensemble of independent realisations of
the dataset given the distribution of points in each bin.

### Summary statistic for each bin

#### PredictiveAsymmetryTest

If the binning method returns a summary statistic (e.g. the mean) for each bin,
then the causality test is applied exactly once to the bin summaries.

```julia
# Bin the data by drawing n_draws = 5000 realizations of each uncertain data
# point, then assign the draws to the correct bins and get a distribution
# of values in each bin.
binning = BinnedResampling(0:10:1000, 5000)

# A predictive asymmetry test. For the predictions, we'll use a transfer entropy
# test, which uses the visitation frequency estimator. We'll predict 5 time steps
# forwards and backwards in time (ηs = -5:5). For the state space binning, we'll
# determine how many intervals to split each axis into from the number of points
# in the time series.
k, l, m = 1, 1, 1 # embedding parameters
n_subdivisions = floor(Int, length(0:10:1000)^(1/(k + l + m + 1)))
state_space_binning = RectangularBinning(n_subdivisions)
ηs = -5:5
te_test = VisitationFrequencyTest(k = k, l = l, m = m,
binning = state_space_binning,
ηs = ηs, b = 2) # use base-2 logarithms
pa_test = PredictiveAsymmetryTest(predictive_test = te_test)

# Compute transfer entropy in both directions for the bin means.
# We still have to specify the number of realisations for the test,
# but it is ignored, so it doesn't matter what number we input here.
test = BinnedDataCausalityTest(pa_test, binning, 1)

# Perform the tests
tes_xy = causality(X, Y, test)
tes_yx = causality(Y, X, test)

# Plot the results
plot(xlabel = "Prediction lag (eta)", ylabel = "Predictive asymmetry (bits)")
plot!(ηs[ηs .> 0], tes_xy, label = "x -> y (binned)", c = :black)
plot!(ηs[ηs .> 0], tes_yx, label = "y -> x (binned)", c = :red)
hline!([0], lw = 2, ls = :dot, α = 0.5, label = "", c = :grey)
```

### Uncertain values representing each bin

If the binning method returns an uncertain value for each bin, then the
causality test is applied to `test.n_realizations` independent draws
of the binned dataset.

#### Predictive asymmetry test

```julia
# Bin the data by drawing n_draws = 5000 realizations of each uncertain data
# point, then assign the draws to the correct bins and get a distribution
# of values in each bin.
grid = 0:5:1000
binning = BinnedResampling(grid, 5000)

# A predictive asymmetry test. For the predictions, we'll use a transfer entropy
# test, which uses the visitation frequency estimator. We'll predict 5 time steps
# forwards and backwards in time (ηs = -5:5). For the state space binning, we'll
# determine how many intervals to split each axis into from the number of points
# in the time series.
k, l, m = 1, 1, 1 # embedding parameters
n_subdivisions = floor(Int, length(0:10:1000)^(1/(k + l + m + 1)))
state_space_binning = RectangularBinning(n_subdivisions)
ηs = -5:5
te_test = VisitationFrequencyTest(k = k, l = l, m = m,
binning = state_space_binning,
ηs = ηs, b = 2) # use base-2 logarithms
pa_test = PredictiveAsymmetryTest(predictive_test = te_test)

# Compute transfer entropy in both directions over 50 independent realizations
# of the binned dataset.
n_realizations = 50
test = BinnedDataCausalityTest(pa_test, binning, n_realizations)

# Perform the tests on the bin means. `tes_xy` now contains 50 independent
# computation of the predictive asymmetry at predictive lags 1 to 5,
# and `tes_yx` contains the same but for the other direction
tes_xy = causality(X, Y, test)
tes_yx = causality(Y, X, test)

# Gather results in a matrix and compute means and standard deviations
# for the predictive asymmetries at each prediction lag
M_xy = hcat(tes_xy...,)
M_yx = hcat(tes_yx...,)

means_xy = mean(M_xy, dims = 2)[:, 1]
means_yx = mean(M_yx, dims = 2)[:, 1]
stdevs_xy = std(M_yx, dims = 2)[:, 1]
stdevs_yx = std(M_yx, dims = 2)[:, 1]

# Plot the predictive asymmetry as a function of prediction lag
plot(xlabel = "Prediction lag (eta)", ylabel = "Predictive asymmetry (bits)")
plot!(ηs[ηs .> 0], means_xy, ribbon = stdevs_xy, label = "x -> y (binned)", c = :black)
plot!(ηs[ηs .> 0], means_yx, ribbon = stdevs_yx, label = "y -> x (binned)", c = :red)
hline!([0], lw = 2, ls = :dot, α = 0.5, label = "", c = :grey)
```

#### Convergent cross mapping test

```julia
using Plots

# Bin the data by drawing n_draws = 5000 realizations of each uncertain data
# point, then assign the draws to the correct bins and get a distribution
# of values in each bin.
binning = BinnedResampling(0:10:1000, 5000)

# A convergent cross mapping causality test
ts_lengths = 15:5:100 |> collect
causalitytest = ConvergentCrossMappingTest(
timeseries_lengths = ts_lengths,
n_reps = 100,
libsize = 100,
τ = 1,
replace = true)

# Run ccm in both directions over 50 independent realizations of the binned datasets
n_realizations = 50
test = BinnedDataCausalityTest(causalitytest, binning, n_realizations)

ccms_xy = causality(X, Y, test);
ccms_yx = causality(Y, X, test);

realization_means_xy = [mean.(realization) for realization in ccms_xy]
realization_means_yx = [mean.(realization) for realization in ccms_yx]
M_xy = hcat(realization_means_xy...,)
M_yx = hcat(realization_means_yx...,)

plot(xlabel = "Time series length", ylabel = "Cross map skill")
plot!(ts_lengths, mean(M_xy, dims = 2)[:, 1],
ribbon = std(M_xy, dims = 2)[:, 1],
label = "x -> y (binned)")
plot!(ts_lengths, mean(M_yx, dims = 2)[:, 1],
ribbon = std(M_yx, dims = 2)[:, 1],
label = "y -> x (binned)")
```
"""
function causality(x::AbstractUncertainIndexValueDataset,
    y::AbstractUncertainIndexValueDataset,
    test::BinnedDataCausalityTest)
end

function causality(x::AbstractUncertainIndexValueDataset,
    y::AbstractUncertainIndexValueDataset,
    test::BinnedDataCausalityTest{CT, BR}) where {
        CT <: CausalityTest, BR <: AbstractBinnedUncertainValueResampling}

    binned_x = resample(x, test.binning)
    binned_y = resample(y, test.binning)

    [causality(resample(binned_x.values), resample(binned_y.values),
            test.test) for i = 1:test.n_realizations]
end

function causality(x::AbstractUncertainIndexValueDataset,
    y::AbstractUncertainIndexValueDataset,
    test::BinnedDataCausalityTest{CT, BR}) where {
        CT <: CausalityTest, BR <: AbstractBinnedSummarisedResampling}

    binned_and_summarised_x = resample(x, test.binning)
    binned_and_summarised_y = resample(y, test.binning)

    causality(binned_and_summarised_x, binned_and_summarised_y, test.test)
end


text:
using LuxurySparse, YaoBase, BitBasis
import MLStyle: @match

export @const_gate

"""
    @const_gate <gate name> = <expr>
    @const_gate <gate name>::<type> = <expr>
    @const_gate <gate>::<type>

This macro simplify the definition of a constant gate. It will automatically bind the matrix form
to a constant which will reduce memory allocation in the runtime.

# Example

```julia
@const_gate X = ComplexF64[0 1;1 0]
```

or

```julia
@const_gate X::ComplexF64 = [0 1;1 0]
```

You can bind new element types by simply re-declare with a type annotation.

```julia
@const_gate X::ComplexF32
```
"""
macro const_gate(ex, ex2=:(nlevel=2))
    @match ex2 begin
        :(nlevel = $nlevel) => begin
            if !(nlevel isa Integer)
                error("You should provide an integer as the value of `nlevel`! e.g. `nlevel=3`")
            end
            _const_gate(__module__, __source__, ex, nlevel)
        end
        _ => begin
            error("The second argument should be e.g. `nlevel=3`")
        end
    end
end

function _const_gate(__module__::Module, __source__::LineNumberNode, ex, nlevel)
    @match ex begin
        :($name::$t = $expr) => define_gate(__module__, __source__, name, t, expr, nlevel)
        :($name = $expr) => define_gate(__module__, __source__, name, expr, nlevel)
        :($name::$t) => define_typed_binding(__module__, __source__, name, t)
        _ => throw(Meta.ParseError("use @const_gate <gate name> = <expr> or
            @const_gate <gate name>::<type> = <expr> to define new gate, or
            use @const_gate <gate name>::<new type> to bind new type"))
    end
end

"""
    gatetype_name(name)

Returns the name of a constant gate by adding postfix `Gate`.
"""
gatetype_name(name) = Symbol(name, "Gate")

"""
    constant_name(name)

Return the name of a constant binding of the matrix for a constant gate.
"""
constant_name(name) = gensym(Symbol("Const", "_", name))

"""
    define_gate(__module__::Module, __source__::LineNumberNode, name, type, expr)

Return the expression to define a constant gate with given type (value of `expr`
will be converted to given `type`).
"""
function define_gate(__module__::Module, __source__::LineNumberNode, name, type, expr, nlevel)
    const_binding = constant_name(name)
    gt_name = gatetype_name(name)

    return quote
        $(define_binding(__module__, const_binding, type, expr))
        $(define_struct(__module__, __source__, const_binding, name, nlevel))
        $(define_methods(__module__, const_binding, name))
        $(define_properties(__module__, const_binding, name))
    end
end

"""
    define_gate(__module__::Module, __source__::LineNumberNode, name, expr, nlevel)

Return the expression to define a constant gate, use the inferred type.
"""
function define_gate(__module__::Module, __source__::LineNumberNode, name, expr, nlevel)
    const_binding = constant_name(name)
    gt_name = gatetype_name(name)

    return quote
        $(define_binding(__module__, const_binding, expr))
        $(define_struct(__module__, __source__, const_binding, name, nlevel))
        $(define_methods(__module__, const_binding, name))
        $(define_properties(__module__, const_binding, name))
    end
end

"""
    define_typed_binding(__module__::Module, __source__::LineNumberNode, name, type)

Define a new binding to a new type.
"""
function define_typed_binding(__module__::Module, __source__::LineNumberNode, name, type)
    gt_name = gatetype_name(name)
    const_binding = constant_name(name)
    return quote
        @eval $__module__ begin
            if !isdefined($__module__, $(QuoteNode(name)))
                throw(UndefVarError($(QuoteNode(name))))
            end

            const $const_binding = YaoBlocks.mat($type, $gt_name)
            YaoBlocks.mat(::Type{$type}, ::Type{$gt_name}) = $const_binding
        end
    end
end

function define_struct(__module__::Module, __source__::LineNumberNode, const_binding, name, nlevel)
    gt_name = gatetype_name(name)
    ex = Expr(:block)

    # export bindings if this is not in Main
    if __module__ !== Main
        push!(ex.args, :(export $name, $gt_name))
    end
    # calculate new shape
    N = gensym(:N)
    push!(ex.args, :(@eval $__module__ $N = $logdi(size($(const_binding), 1), $nlevel)))

    # we allow overwrite in order to support the following syntax
    # @const_gate X = BLABLA
    # @const_gate X = Correct_Matrix
    if isdefined(__module__, gt_name)
        msg = "$(string(gt_name)) is already defined, overwritten by new definition at $__source__"
        push!(ex.args, :(@warn $msg))
        push!(
            ex.args,
            :(@eval $__module__ @assert $N == $logdi(size(mat($name), 1), $nlevel) "new constant does not have the same size with previous definitions"),
        )
    else
        push!(
            ex.args,
            :(@eval $__module__ Base.@__doc__ struct $gt_name <:
                                                     YaoBlocks.ConstGate.ConstantGate{$N,$nlevel} end),
        )
    end
    push!(ex.args, :(@eval $__module__ Base.@__doc__ const $(name) = $(gt_name)()))
    return ex
end

function define_binding(__module__::Module, const_binding, type, expr)
    return quote
        @eval $__module__ begin
            const $const_binding = similar($expr, $type)
            if size($(const_binding), 1) != size($(const_binding), 2)
                throw(DimensionMismatch("Quantum Gates must be square matrix."))
            end
            copyto!($const_binding, $(expr))
        end
    end
end

function define_binding(__module__::Module, const_binding, expr)
    return quote
        @eval $__module__ begin
            const $const_binding = $(expr)
            if size($(const_binding), 1) != size($(const_binding), 2)
                throw(DimensionMismatch("Quantum Gates must be square matrix."))
            end
        end
    end
end

# forward instance to const gate type
YaoBlocks.mat(::Type{T}, ::GT) where {T,GT<:ConstantGate} = YaoBlocks.mat(T, GT)

function define_methods(__module__::Module, const_binding, name)
    gt_name = gatetype_name(name)
    return quote
        @eval $__module__ begin
            YaoBlocks.mat(::Type{eltype($const_binding)}, ::Type{$gt_name}) = $const_binding

            function YaoBlocks.mat(::Type{T}, ::Type{$gt_name}) where {T}
                src = YaoBlocks.mat(eltype($const_binding), $gt_name)
                dst = similar(src, T)
                copyto!(dst, src)
                return dst
            end

            function YaoBlocks.print_block(io::IO, ::$gt_name)
                print(io, $(QuoteNode(name)))
            end
        end # eval
    end # quote
end

function define_properties(__module__::Module, const_binding, name)
    gt_name = gatetype_name(name)
    function def(property)
        flag = gensym(:flag)
        return quote
            @eval $__module__ begin
                const $flag = YaoBlocks.$property($const_binding)
                YaoBlocks.$property(::$gt_name) = $flag
                YaoBlocks.$property(::Type{GT}) where {GT<:$gt_name} = $flag
            end
        end
    end

    ex = Expr(:block)
    for each in [:ishermitian, :isreflexive, :isunitary]
        push!(ex.args, def(each))
    end
    return ex
end


text:
module AssociativeWindowAggregation

export FixedWindowAssociativeOp, TimeWindowAssociativeOp, WindowedAssociativeOp
export update_state!, window_full, window_size, window_value

include("base.jl")
include("fixed.jl")
include("time.jl")

end # module


text:
#!/usr/bin/env julia
#=
Common struct
del2z <delta.z@aliyun.com>
=#


text:
if VERSION < v"0.4-"
    using Dates
else
    using Base.Dates
end

typealias URange Union{Union{Range, UnitRange}}

typealias DTArray Union{AbstractMatrix{DateTime}, AbstractMatrix{Date},
                        AbstractVector{DateTime}, AbstractVector{Date}}

typealias GlyphTypes Union{AbstractString, Glyph, Vector{Glyph}}

const EPOCH = DateTime(1970, 1, 1)

unixtime(d::Date) = unixtime(convert(DateTime, d))
unixtime(d::DateTime) = Int(d - EPOCH)

## getglyphs methods
getglyphs(styles::AbstractString, count::Int64) =
    getglyphs(convert(Vector{Glyph}, styles), count)

getglyphs(glyph::Glyph, count::Int64) = getglyphs([glyph], count)

getglyphs(glyphs::Vector{Glyph}, count::Int64) =
    repmat(glyphs, ceil(Int, count / length(glyphs)))

## Plot methods
plot(f::Function, args...;kwargs...) = plot([f], args...; kwargs...)

function plot(fs::Vector{Function}, rng::URange, args...;kwargs...)
    stop = isdefined(rng, :stop) ? rng.stop : rng.len
    plot(fs, rng.start, stop, args...; kwargs...)
end

function plot(fs::Vector{Function}, start::Real=-10, stop::Real=10, args...;
              counteval::Int=COUNTEVAL, kwargs...)
    x = linspace(start, stop, counteval)
    y = Array(Float64, counteval, length(fs))
    for (i, f) in enumerate(fs)
        y[:, i] = map(f, x)
    end
    x = repmat(x, 1, length(fs))
    plot(x, y, args...; kwargs...)
end

plot(x::DTArray, args...; kwargs...) =
    plot(map(unixtime, x), args...; x_axis_type=:datetime, kwargs...)

plot(x::RealArray, y::DTArray, args...; kwargs...) =
    plot(x, map(unixtime, y), args...; y_axis_type=:datetime, kwargs...)

function plot(y::Vector{Vector}, styles::GlyphTypes=DEFAULT_GLYPHS_STR; kwargs...)
    x = [collect(1:length(yv)) for yv in y]
    plot(x, y, styles; kwargs...)
end

function plot(x::RealArray, y::Vector{Vector}, args...; kwargs...)
    x = [x[1:length(yv)] for yv in y]
    plot(x, y, args...; kwargs...)
end

function plot{T<:AbstractVector}(x::Vector{T}, y::Vector{Vector},
                                   styles::GlyphTypes=DEFAULT_GLYPHS_STR;
                                   autoopen::Bool=AUTOOPEN, kwargs...)
    hold_val = HOLD

    # clear CURPLOT if hold was previously false
    hold(true, !hold_val)
    global AUTOOPEN = false

    lx = length(x)

    if lx != length(y)
        error("length of x and y are not equal: x: $(length(x)), y: $(length(y))")
    end

    glyphs = getglyphs(styles, lx)

    for i in 1:lx
        plot(x[i], y[i], glyphs[i]; kwargs...)
    end

    plt = CURPLOT
    hold(hold_val, !hold_val)

    !isinteractive() && autoopen && showplot(plt)

    global AUTOOPEN = autoopen

    return plt
end

function plot(y::RealArray, styles::GlyphTypes=DEFAULT_GLYPHS_STR; kwargs...)
    x = 1:size(y, 1)
    x = ndims(y) > 1 ? repmat(x, 1, size(y, 2)) : x
    plot(x, y, styles; kwargs...)
end

function plot(x::RealVect, y::RealMat, args...; kwargs...)
    x = repmat(x, 1, size(y, 2))
    plot(x, y, args...; kwargs...)
end

function plot(x::RealVect, y::RealVect, args...; kwargs...)
    # is there a more efficient way of forcing a matrix ?
    # this seems ugly but is apparently pretty quick even for big mats
    plot(reshape(x, length(x), 1), reshape(y, length(y), 1), args...; kwargs...)
end

function plot(x::RealMat, y::RealMat, styles::GlyphTypes=DEFAULT_GLYPHS_STR;
              legends::Union{Void, Vector}=nothing, kwargs...)
    if size(x) != size(y)
        error("size of x and y are not equal: x: $(size(x)), y: $(size(y))")
    end

    cols = size(x, 2)
    glyphs = getglyphs(styles, cols)
    legends = legends == nothing ? [nothing for _ in 1:cols] : legends
    dcs = BokehDataSet[BokehDataSet(x[:,i], y[:,i], glyphs[i], legends[i])
                       for i in 1:cols]
    plot(dcs; kwargs...)
end

# there a good if boring reason why we have to use nothing for width, height etc. rather than
# set the default to WIDTH, HEIGHT etc.: its because we wouldn't be able to specify a new width
# or height on an extending plot if the new value happened to be the same as WIDTH or HEIGHT
function plot(columns::Array{BokehDataSet, 1};
              extend::Union{Void, Plot}=nothing,
              title::NullString=nothing, width::NullInt=nothing,
              height::NullInt=nothing, x_axis_type::NullSymbol=nothing,
              y_axis_type::NullSymbol=nothing, legendsgo::NullSymbol=nothing,
              plotfile::NullString=nothing,
              tools::Union{Void, Array{Symbol, 1}}=nothing,
              autoopen::Bool=AUTOOPEN)
    extend == nothing && !HOLD && (global CURPLOT = nothing)

    if extend == nothing && CURPLOT == nothing
        plt = Plot(columns,
                   tools == nothing ? TOOLS : tools,
                   plotfile == nothing ? PLOTFILE : plotfile,
                   title == nothing ? TITLE : title,
                   width == nothing ? WIDTH : width,
                   height == nothing ? HEIGHT : height,
                   x_axis_type == nothing ? X_AXIS_TYPE : x_axis_type,
                   y_axis_type == nothing ? Y_AXIS_TYPE : y_axis_type,
                   legendsgo)
        extend == nothing && (global CURPLOT = plt)
    else
        function add2plot!(p::Plot)
            append!(p.datacolumns, columns)
            tools != nothing && (p.tools = tools)
            plotfile != nothing && (p.filename = plotfile)
            title != nothing && (p.title = title)
            width != nothing && (p.width = width)
            height != nothing && (p.height = height)
            x_axis_type != nothing && (p.x_axis_type = x_axis_type)
            y_axis_type != nothing && (p.y_axis_type = y_axis_type)
            legendsgo != nothing && (p.legendsgo = legendsgo)
        end
        if extend == nothing
            add2plot!(CURPLOT)
            plt = CURPLOT
        else
            add2plot!(extend)
            plt = extend
        end
    end
    !isinteractive() && autoopen && showplot(plt)
    return plt
end


text:
using MPIMeasurements
using Plots
using Unitful
using UnitfulRecipes
using Dates
using UUIDs

#pyplot()
#plotly()
#inspectdr()
default(show = true)

ENV["JULIA_DEBUG"] = "all"

scannerName_ = "MagneticFieldMeasurement"
protocolName_ = "MagneticFieldMeasurement"
operator = "Jonas"
filename_ = joinpath("./field_$(Dates.format(now(), "yyyy-mm-dd_HH-MM-SS")).h5")

protocol = Protocol(protocolName_, scannerName_)
filename(protocol, filename_)
channel = runProtocol(protocol)

measurement_ = MagneticFieldMeasurement(filename_)

result = zeros(shape(measurement_.positions)...)
for pos in measurement_.positions
  idx = posToLinIdx(measurement_.positions, pos)
  #@debug measurement_.fields[idx, :]
  field = norm(measurement_.fields[idx, :])
  result[posToIdx(measurement_.positions, pos)...] = ustrip(u"mT", field)
end

heatmap(result[:,3,:])


text:
module Lexer

using Compat

import Base.UTF8proc
import Base: position
import ..Diagnostics: diag, Incomplete, before
using ..Tokens
using ..Tokens: √

export Token, TokenStream, next_token, set_token!, last_token,
       put_back!, peek_token, take_token, require_token

const SYM_TRUE  = Symbol("true")
const SYM_FALSE = Symbol("false")
const SYM_CTRANSPOSE = Symbol("'")

const EOF_CHAR = convert(Char,typemax(UInt32))

add_dot(s::Symbol) = Symbol(string(".",s))
add_dots!(ops) = append!(ops,add_dot.(ops))

const ALL_OPS = Dict{Symbol,Any}(
:assignment  =>
       [add_dots!([:(=), :(+=), :(-=), :(*=), :(/=), :(//=), :(|=),
            :(^=), :(÷=), :(%=), :(<<=), :(>>=), :(>>>=), :(\=), :(&=)]);
        [:(:=), :(=>), :(~), :($=)]],
:conditional => [:(?)],
:lazy_or     => [:(||)],
:lazy_and    => [:(&&)],
# NOTE: -- is deprecated in 0.4
:arrow       =>
       [Symbol("--"); :(-->);
        add_dots!([:(←), :(→), :(↔), :(↚), :(↛), :(↠), :(↣), :(↦),
        :(↮), :(⇎), :(⇏), :(⇒), :(⇔), :(⇴), :(⇶), :(⇷), :(⇸), :(⇹),
        :(⇺), :(⇻), :(⇼), :(⇽), :(⇾), :(⇿), :(⟵), :(⟶), :(⟷), :(⟷),
        :(⟹), :(⟺), :(⟻), :(⟼), :(⟽), :(⟾), :(⟿), :(⤀), :(⤁), :(⤂),
        :(⤃), :(⤄), :(⤅), :(⤆), :(⤇), :(⤌), :(⤍), :(⤎), :(⤏), :(⤐),
        :(⤑), :(⤔), :(⤕), :(⤖), :(⤗), :(⤘), :(⤝), :(⤞), :(⤟), :(⤠),
        :(⥄), :(⥅), :(⥆), :(⥇), :(⥈), :(⥊), :(⥋), :(⥎), :(⥐), :(⥒),
        :(⥓), :(⥖), :(⥗), :(⥚), :(⥛), :(⥞), :(⥟), :(⥢), :(⥤), :(⥦),
        :(⥧), :(⥨), :(⥩), :(⥪), :(⥫), :(⥬), :(⥭), :(⥰), :(⧴), :(⬱),
        :(⬰), :(⬲), :(⬳), :(⬴), :(⬵), :(⬶), :(⬷), :(⬸), :(⬹), :(⬺),
        :(⬻), :(⬼), :(⬽), :(⬾), :(⬿), :(⭀), :(⭁), :(⭂), :(⭃), :(⭄),
        :(⭇), :(⭈), :(⭉), :(⭊), :(⭋), :(⭌), :(￩), :(￫)])],
:comparison =>
       [Symbol(".!"); :(<:); :(>:);
        add_dots!(
       [:(>),  :(<), :(>=), :(≥), :(<=), :(≤), :(==), :(===), :(≡),
        :(!=), :(≠), :(!==), :(≢), :(∈),
        :(∉), :(∋), :(∌), :(⊆), :(⊈), :(⊂), :(⊄), :(⊊), :(∝), :(∊), :(∍),
        :(∥), :(∦), :(∷), :(∺), :(∻), :(∽), :(∾), :(≁), :(≃), :(≄), :(≅),
        :(≆), :(≇), :(≈), :(≉), :(≊), :(≋), :(≌), :(≍), :(≎), :(≐), :(≑),
        :(≒), :(≓), :(≔), :(≕), :(≖), :(≗), :(≘), :(≙), :(≚), :(≛), :(≜),
        :(≝), :(≞), :(≟), :(≣), :(≦), :(≧), :(≨), :(≩), :(≪), :(≫), :(≬),
        :(≭), :(≮), :(≯), :(≰), :(≱), :(≲), :(≳), :(≴), :(≵), :(≶), :(≷),
        :(≸), :(≹), :(≺), :(≻), :(≼), :(≽), :(≾), :(≿), :(⊀), :(⊁), :(⊃),
        :(⊅), :(⊇), :(⊉), :(⊋), :(⊏), :(⊐), :(⊑), :(⊒), :(⊜), :(⊩), :(⊬),
        :(⊮), :(⊰), :(⊱), :(⊲), :(⊳), :(⊴), :(⊵), :(⊶), :(⊷), :(⋍), :(⋐),
        :(⋑), :(⋕), :(⋖), :(⋗), :(⋘), :(⋙), :(⋚), :(⋛), :(⋜), :(⋝), :(⋞),
        :(⋟), :(⋠), :(⋡), :(⋢), :(⋣), :(⋤), :(⋥), :(⋦), :(⋧), :(⋨), :(⋩),
        :(⋪), :(⋫), :(⋬), :(⋭), :(⋲), :(⋳), :(⋴), :(⋵), :(⋶), :(⋷), :(⋸),
        :(⋹), :(⋺), :(⋻), :(⋼), :(⋽), :(⋾), :(⋿), :(⟈), :(⟉), :(⟒), :(⦷),
        :(⧀), :(⧁), :(⧡), :(⧣), :(⧤), :(⧥), :(⩦), :(⩧), :(⩪), :(⩫), :(⩬),
        :(⩭), :(⩮), :(⩯), :(⩰), :(⩱), :(⩲), :(⩳), :(⩴), :(⩵), :(⩶), :(⩷),
        :(⩸), :(⩹), :(⩺), :(⩻), :(⩼), :(⩽), :(⩾), :(⩿), :(⪀), :(⪁), :(⪂),
        :(⪃), :(⪄), :(⪅), :(⪆), :(⪇), :(⪈), :(⪉), :(⪊), :(⪋), :(⪌), :(⪍),
        :(⪎), :(⪏), :(⪐), :(⪑), :(⪒), :(⪓), :(⪔), :(⪕), :(⪖), :(⪗), :(⪘),
        :(⪙), :(⪚), :(⪛), :(⪜), :(⪝), :(⪞), :(⪟), :(⪠), :(⪡), :(⪢), :(⪣),
        :(⪤), :(⪥), :(⪦), :(⪧), :(⪨), :(⪩), :(⪪), :(⪫), :(⪬), :(⪭), :(⪮),
        :(⪯), :(⪰), :(⪱), :(⪲), :(⪳), :(⪴), :(⪵), :(⪶), :(⪷), :(⪸), :(⪹),
        :(⪺), :(⪻), :(⪼), :(⪽), :(⪾), :(⪿), :(⫀), :(⫁), :(⫂), :(⫃), :(⫄),
        :(⫅), :(⫆), :(⫇), :(⫈), :(⫉), :(⫊), :(⫋), :(⫌), :(⫍), :(⫎), :(⫏),
        :(⫐), :(⫑), :(⫒), :(⫓), :(⫔), :(⫕), :(⫖), :(⫗), :(⫘), :(⫙), :(⫷),
        :(⫸), :(⫹), :(⫺), :(⊢), :(⊣)])],
:pipe       => [:(|>),  :(<|)],
:colon      => [:(:), :(..)],
:plus       => [:($);
        add_dots!([:(+), :(-), :(++), :(⊕), :(⊖), :(⊞), :(⊟), :(|), :(∪), :(∨),
        :(⊔), :(±), :(∓), :(∔), :(∸), :(≂), :(≏), :(⊎), :(⊻), :(⊽),
        :(⋎), :(⋓), :(⧺), :(⧻), :(⨈), :(⨢), :(⨣), :(⨤), :(⨥), :(⨦), :(⨧),
        :(⨨), :(⨩), :(⨪), :(⨫), :(⨬), :(⨭), :(⨮), :(⨹), :(⨺), :(⩁), :(⩂),
        :(⩅), :(⩊), :(⩌), :(⩏), :(⩐), :(⩒), :(⩔), :(⩖), :(⩗), :(⩛), :(⩝),
        :(⩡), :(⩢), :(⩣)])],
:bitshift   => add_dots!([:(<<), :(>>), :(>>>)]),
:times      =>
       add_dots!([:(*), :(/), :(÷), :(%), :(⋅), :(∘), :(×),
        :(\), :(&), :(∩), :(∧), :(⊗), :(⊘), :(⊙), :(⊚), :(⊛), :(⊠), :(⊡),
        :(⊓), :(∗), :(∙), :(∤), :(⅋), :(≀), :(⊼), :(⋄), :(⋆), :(⋇), :(⋉), :(⋊),
        :(⋋), :(⋌), :(⋏), :(⋒), :(⟑), :(⦸), :(⦼), :(⦾), :(⦿), :(⧶), :(⧷), :(⨇),
        :(⨰), :(⨱), :(⨲), :(⨳), :(⨴), :(⨵), :(⨶), :(⨷), :(⨸), :(⨻), :(⨼), :(⨽),
        :(⩀), :(⩃), :(⩄), :(⩋), :(⩍), :(⩎), :(⩑), :(⩓), :(⩕), :(⩘), :(⩚), :(⩜),
        :(⩞), :(⩟), :(⩠), :(⫛), :(⊍)]),
:rational   => add_dots!([:(//)]),
:power      => add_dots!(
       [:(^), :(↑), :(↓), :(⇵), :(⟰), :(⟱), :(⤈), :(⤉), :(⤊), :(⤋),
        :(⤒), :(⤓),  :(⥉), :(⥌), :(⥍), :(⥏), :(⥑), :(⥔), :(⥕), :(⥘), :(⥙),
        :(⥜), :(⥝),  :(⥠), :(⥡), :(⥣), :(⥥), :(⥮), :(⥯), :(￪), :(￬)]),
:decl       => [:(::)],
:dot        => [:(.)])

if VERSION < v"0.4.0-dev+573"
    const OPS_BY_PRECEDENT = Any[
        :assignment, :conditional, :lazy_or, :lazy_and, :arrow, :comparison, :pipe,
        :colon, :plus, :bitshift, :times, :rational, :power, :decl, :dot]
else
    const OPS_BY_PRECEDENT = Any[
        :assignment, :conditional, :arrow, :lazy_or, :lazy_and, :comparison, :pipe,
        :colon, :plus, :bitshift, :times, :rational, :power, :decl, :dot]
end

precedent_ops(n::Integer) = Set{Symbol}(ALL_OPS[OPS_BY_PRECEDENT[n]])
precedent_ops(s::Symbol) = Set{Symbol}(ALL_OPS[s])

const ASSIGNMENT_OPS = precedent_ops(:assignment)

const UNARY_OPS = Set{Symbol}([:(+),  :(-), :(!), :(~), :(<:), :(¬),
                               :(>:), :(√), :(∛), :(∜)])

const UNARY_AND_BINARY_OPS = Set{Symbol}([:(+), :(-), :($), :(&), :(~)])

"Operators are special forms, not function names."
const SYNTACTIC_OPS = Set{Symbol}([:(:=); :(-->); :($=); :(=>); :(||); :(.);
                                    :(...); :(->); add_dots!([
                                    :(=)  , :(+=)  ,  :(-=),
                                    :(/=) , :(//=) ,  :(*=),
                                    :(\=) , :(^=)  ,  :(%=),
                                    :(|=) , :(&=)  ,  :(<<=),
                                    :(>>=), :(>>>=),  :(&&) ])])

const SYNTACTIC_UNARY_OPS = Set{Symbol}([:($), :(&), :(::)])

const OPERATORS = union(Set([:(~), :(!), :(->), :(√), :(∛), :(∜), :(...), :(¬),
                             :(.'), SYM_CTRANSPOSE]),
                             [Set(ALL_OPS[ops]) for ops in OPS_BY_PRECEDENT]...)

const RESERVED_WORDS = Set{Symbol}([:begin,  :while, :if, :for, :try, :return,
                                    :break, :continue, :function, :stagedfunction,
                                    :macro, :quote, :let, :local, :global, :const,
                                    :abstract, :typealias, :type, :bitstype, :immutable,
                                    :ccall, :do, :module, :baremodule, :using, :import,
                                    :export, :importall])
#= Helper functions =#

const PRECEDENCE_MAP = Dict{Symbol, Int}()

let
    for (i, ops) in enumerate(OPS_BY_PRECEDENT)
        for op in ALL_OPS[ops]
            PRECEDENCE_MAP[op] = i
        end
    end
end

operator_precedence(op::Symbol) = PRECEDENCE_MAP[op]

is_syntactic_op(op::Symbol) = in(op, SYNTACTIC_OPS)
is_syntactic_op(op) = false

is_syntactic_unary_op(op::Symbol) = in(op, SYNTACTIC_UNARY_OPS)
is_syntactic_unary_op(op) = false

const SPECIAL_CHARS = Set{Char}("()[]{},;\"`@")

is_special_char(c::Char) = in(c, SPECIAL_CHARS)

isnewline(c::Char) = c === '\n'
isnewline(t::AbstractToken) = isnewline(¬t)
isnewline(c) = false

function isuws(ch::Char)
    c = UInt32(ch)
    return (c==9    || c==10   || c==11   || c==12   || c==13   || c==32 ||
            c==133  || c==160  || c==5760 || c==6158 || c==8192 ||
            c==8193 || c==8194 || c==8195 || c==8196 || c==8197 ||
            c==8198 || c==8199 || c==8200 || c==8201 || c==8202 ||
            c==8232 || c==8233 || c==8239 || c==8287 || c==12288)
end

isbom(c::Char) = UInt32(c) == 0xFEFF

is_zero_width_space(c::Char) = c === '\u200b' || c === '\u2060' || c === '\ufeff'

is_ignorable_char(c::Char) = is_zero_width_space(c) ||
                             ('\u200c' <= c <= '\u200f') ||
                             (c === '\u00ad' || c === '\u2061' || c === '\u115f')

function is_cat_id_start(ch::Char, cat::Integer)
    c = UInt32(ch)
    return (cat == UTF8proc.UTF8PROC_CATEGORY_LU || cat == UTF8proc.UTF8PROC_CATEGORY_LL ||
            cat == UTF8proc.UTF8PROC_CATEGORY_LT || cat == UTF8proc.UTF8PROC_CATEGORY_LM ||
            cat == UTF8proc.UTF8PROC_CATEGORY_LO || cat == UTF8proc.UTF8PROC_CATEGORY_NL ||
            cat == UTF8proc.UTF8PROC_CATEGORY_SC ||  # allow currency symbols
            cat == UTF8proc.UTF8PROC_CATEGORY_SO ||  # other symbols

            # math symbol (category Sm) whitelist
            (c >= 0x2140 && c <= 0x2a1c &&
             ((c >= 0x2140 && c <= 0x2144) || # ⅀, ⅁, ⅂, ⅃, ⅄
              c == 0x223f || c == 0x22be || c == 0x22bf || # ∿, ⊾, ⊿
              c == 0x22a4 || c == 0x22a5 || # ⊤ ⊥
              (c >= 0x22ee && c <= 0x22f1) || # ⋮, ⋯, ⋰, ⋱

              (c >= 0x2202 && c <= 0x2233 &&
               (c == 0x2202 || c == 0x2205 || c == 0x2206 || # ∂, ∅, ∆
                c == 0x2207 || c == 0x220e || c == 0x220f || # ∇, ∎, ∏
                c == 0x2210 || c == 0x2211 || # ∐, ∑
                c == 0x221e || c == 0x221f || # ∞, ∟
                c >= 0x222b)) || # ∫, ∬, ∭, ∮, ∯, ∰, ∱, ∲, ∳

              (c >= 0x22c0 && c <= 0x22c3) ||  # N-ary big ops: ⋀, ⋁, ⋂, ⋃
              (c >= 0x25F8 && c <= 0x25ff) ||  # ◸, ◹, ◺, ◻, ◼, ◽, ◾, ◿

              (c >= 0x266f &&
               (c == 0x266f || c == 0x27d8 || c == 0x27d9 || # ♯, ⟘, ⟙
                (c >= 0x27c0 && c <= 0x27c2) ||  # ⟀, ⟁, ⟂
                (c >= 0x29b0 && c <= 0x29b4) ||  # ⦰, ⦱, ⦲, ⦳, ⦴
                (c >= 0x2a00 && c <= 0x2a06) ||  # ⨀, ⨁, ⨂, ⨃, ⨄, ⨅, ⨆
                (c >= 0x2a09 && c <= 0x2a16) ||  # ⨉, ⨊, ⨋, ⨌, ⨍, ⨎, ⨏, ⨐, ⨑, ⨒,
                                                 # ⨓, ⨔, ⨕, ⨖
                c == 0x2a1b || c == 0x2a1c)))) || # ⨛, ⨜

            (c >= 0x1d6c1 && # variants of \nabla and \partial
             (c == 0x1d6c1 || c == 0x1d6db ||
              c == 0x1d6fb || c == 0x1d715 ||
              c == 0x1d735 || c == 0x1d74f ||
              c == 0x1d76f || c == 0x1d789 ||
              c == 0x1d7a9 || c == 0x1d7c3)) ||

            # super- and subscript +-=()
            (c >= 0x207a && c <= 0x207e) ||
            (c >= 0x208a && c <= 0x208e) ||

            # angle symbols
            (c >= 0x2220 && c <= 0x2222) || # ∠, ∡, ∢
            (c >= 0x299b && c <= 0x29af) || # ⦛, ⦜, ⦝, ⦞, ⦟, ⦠, ⦡, ⦢, ⦣, ⦤, ⦥,
                                            # ⦦, ⦧, ⦨, ⦩, ⦪, ⦫, ⦬, ⦭, ⦮, ⦯
            # Other_ID_Start
            c == 0x2118 || c == 0x212E || # ℘, ℮
            (c >= 0x309B && c <= 0x309C)) # katakana-hiragana sound marks
end

function is_identifier_char(c::Char)
    if ((c >= 'A' && c <= 'Z') ||
        (c >= 'a' && c <= 'z') || c == '_' ||
        (c >= '0' && c <= '9') || c == '!')
        return true
    elseif (UInt32(c) < 0xA1 || UInt32(c) > 0x10ffff)
        return false
    end
    cat = UTF8proc.category_code(c)
    is_cat_id_start(c, cat) && return true
    if cat == UTF8proc.UTF8PROC_CATEGORY_MN || cat == UTF8proc.UTF8PROC_CATEGORY_MC ||
       cat == UTF8proc.UTF8PROC_CATEGORY_ND || cat == UTF8proc.UTF8PROC_CATEGORY_PC ||
       cat == UTF8proc.UTF8PROC_CATEGORY_SK || cat == UTF8proc.UTF8PROC_CATEGORY_ME ||
       cat == UTF8proc.UTF8PROC_CATEGORY_NO ||
       (0x2032 <= UInt32(c) <= 0x2034) || # primes
       UInt32(c) == 0x0387 || UInt32(c) == 0x19da ||
       (0x1369 <= UInt32(c) <= 0x1371)
       return true
    end
    return false
end

function is_identifier_start_char(c::Char)
    if ((c >= 'A' && c <= 'Z') || (c >= 'a' && c <= 'z') || c == '_')
        return true
    elseif (UInt32(c) < 0xA1 || UInt32(c) > 0x10ffff)
        return false
    end
    cat = UTF8proc.category_code(c)
    return is_cat_id_start(c, cat)
end

"Characters that can be in an operator."
const OPERATOR_CHARS = union([Set(string(op)) for op in OPERATORS]...)

is_opchar(c::Char) = in(c, OPERATOR_CHARS)

"Characters that can follow a `.` in an operator."
const DOT_OP_CHARS = unique(map(op->op[2], filter(op->length(op)>=2&&op[1]=='.',map(string,OPERATORS))))

is_dot_opchar(c::Char) = in(c, DOT_OP_CHARS)

is_operator(op::Symbol) = in(op, OPERATORS)
is_operator(op::AbstractToken) = is_operator(¬op)
is_operator(op) = false

"""
`peekchar(io::Bufer)` -> `Char`
`peekchar(s::IOStream)` -> `Char`
`peekchar(ts::TokenStream)` -> `Char`

Version of `peekchar` for `IOBuffer`, `IOStream` and `TokenStream`,
modified from `Base` to give the same semantics as the `IOStream` implementation.
"""
function peekchar end

function peekchar(io::IOBuffer)
    if !io.readable || io.ptr > io.size
        return EOF_CHAR
    end
    ch = convert(UInt8,io.data[io.ptr])
    if ch < 0x80
        return convert(Char,ch)
    end
    # mimic utf8.next function
    trailing = Base.utf8_trailing[ch+1]
    c::UInt32 = 0
    for j = 1:trailing
        c += ch
        c <<= 6
        ch = convert(UInt8,io.data[io.ptr+j])
    end
    c += ch
    c -= Base.utf8_offset[trailing+1]
    return convert(Char,c)
end

# this implementation is copied from Base
const _CHTMP = Array{Char}(1)

peekchar(s::IOStream) = begin
    if ccall(:ios_peekutf8, Int32, (Ptr{Void}, Ptr{Char}), s, _CHTMP) < 0
        return EOF_CHAR
    end
    return _CHTMP[1]
end

eof(io::IO) = Base.eof(io)
eof(t::AbstractToken) = eof(¬t)
eof(c) = c===EOF_CHAR

readchar(io::IO) = eof(io) ? EOF_CHAR : read(io, Char)
takechar(io::IO) = (readchar(io); io)

#= Token Stream =#

type TokenStream{T, I <: IO}
    io::I
    lineno::Int
    lasttoken
    putback
    isspace::Bool
    ateof::Bool
    filename::AbstractString
end

(::Type{TokenStream{T}}){T}(io::IO) = TokenStream{T, typeof(io)}(io, 1, nothing, nothing, false, eof(io), "none")
(::Type{TokenStream{T}}){T}(str::AbstractString) = TokenStream{T}(IOBuffer(str))
TokenStream(x) = TokenStream{Token}(x)

eof(ts::TokenStream) = ts.ateof || eof(ts.io)
skip(ts::TokenStream, n)  = Base.skip(ts.io, n)
position(ts::TokenStream) = Base.position(ts.io)
peekchar(ts::TokenStream) = peekchar(ts.io)

# The last case is the replacement character 0xfffd (3 bytes)
utf8sizeof(c::Char) = UInt32(c) < 0x80 ? 1 : UInt32(c) < 0x800 ? 2 :
    UInt32(c) < 0x10000 ? 3 : UInt32(c) < 0x110000 ? 4 : 3

readchar(ts::TokenStream) = begin
    @assert !eof(ts)
    loc = position(ts.io)
    c = read(ts.io, Char)
    c === '\n' && (ts.lineno += 1)
    return c
end

takechar(ts::TokenStream) = (readchar(ts); ts)

function skipws(ts::TokenStream, newlines::Bool=false)
    nc = peekchar(ts)
    nc === EOF_CHAR && return false
    skipped = false
    while !eof(ts) && (isuws(nc) || isbom(nc)) && (newlines || nc !== '\n')
        takechar(ts)
        nc, skipped = peekchar(ts), true
    end
    return skipped
end

#= Lexer =#
function skip_to_eol(io::IO)
    while !eof(io)
        nc = peekchar(io)
        nc === '\n' && break
        Base.skip(io, utf8sizeof(nc))
    end
    return io
end
skip_to_eol(ts::TokenStream) = skip_to_eol(ts.io)

function read_operator(ts::TokenStream, c::Char, r)
    nc = peekchar(ts)
    if (c === '*') && (nc === '*')
        readchar(ts); loc = Lexer.makerange(ts, r)
        D = diag(loc, "use '^' instead of '**'")
        diag(D, loc, "^", :fixit)
        throw(D)
    end
    # 1 char operator
    if eof(nc) || !is_opchar(nc) || (c === ':' && nc === '-')
        return Symbol(c)
    end
    str, c = [c], nc
    opsym  = Symbol(String(copy(str)))
    while true
        if !eof(c) && is_opchar(c)
            push!(str, c)
            newop = Symbol(String(copy(str)))
            if is_operator(newop)
                skip(ts, utf8sizeof(c))
                c, opsym = peekchar(ts), newop
                continue
            end
        end
        return opsym
    end
end

#=============#
# Read Number
#=============#

function sized_uint_literal(s::AbstractString, b::Integer)
    i = s[1] === '-' ? 3 : 2
    l = (length(s) - i) * b
    l <= 8   && return @compat parse(UInt8,   s)
    l <= 16  && return @compat parse(UInt16,  s)
    l <= 32  && return @compat parse(UInt32,  s)
    l <= 64  && return @compat parse(UInt64,  s)
    l <= 128 && return @compat parse(UInt128, s)
    return @compat parse(BigInt,s)
end

function sized_uint_oct_literal(s::AbstractString)
    contains(s, "o0") && return sized_uint_literal(s, 3)
    len = length(s)
    (len < 5  || (len == 5  && s <= "0o377")) && return @compat parse(UInt8,s)
    (len < 8  || (len == 8  && s <= "0o177777")) && return @compat parse(UInt16, s)
    (len < 13 || (len == 13 && s <= "0o37777777777")) && return @compat parse(UInt32,s)
    (len < 24 || (len == 24 && s <= "0o1777777777777777777777")) && return @compat parse(UInt64,s)
    (len < 45 || (len == 45 && s <= "0o3777777777777777777777777777777777777777777")) && return @compat parse(UInt128,s)
    return @compat parse(BigInt,s)
end

function compare_num_strings(s1::AbstractString, s2::AbstractString)
    s1, s2 = lstrip(s1, '0'), lstrip(s2, '0')
    l1, l2 = length(s1), length(s2)
    return l1 == l2 ? s1 <= s2 : l1 <= l2
end

function is_oct_within_uint128(s::AbstractString)
    len = length(s)
    max = "0o3777777777777777777777777777777777777777777"
    len < 45  && return true
    len > 45  && return false
    len == 45 && s[1] === '-' ? s[2:end] <= max : s <= max
end

function is_within_int128(s::AbstractString)
    len = length(s)
    if s[1] === '-'
        len < 40  && return true
        len > 40  && return false
        len == 40 && return s <= "-170141183460469231731687303715884105728"
    else
        len < 39  && return true
        len > 39  && return false
        len == 39 && s <= "170141183460469231731687303715884105727"
    end
end

# Notes:
# expressions starting with 0x are always hexadecimal literals
# expressions starting with a numeric literal followed by e or E
# are always floating point literals

function string_to_number(str::AbstractString, loc)
    len = length(str)
    len > 0 || throw(ArgumentError("empty string"))
    neg = str[1] === '-'
    # NaN and Infinity
    (str == "NaN" || str == "+NaN" || str == "-NaN") && return NaN
    (str == "Inf" || str == "+Inf" || str == "-Inf") && return Inf
    # floating point literals
    didx, fidx = 0, 0
    isfloat32, isfloat64 = false, false
    for i=1:len
        c = str[i]
        if c === '.'
            @assert isfloat64 == false
            didx, isfloat64 = i, true
        elseif c === 'f'
            @assert i > didx && i != len
            fidx, isfloat32 = i, true
        elseif c === 'e' || c === 'E' || c === 'p' || c === 'P'
            isfloat64 = true
        end
    end
    if isfloat32
        base = @compat parse(Float64,str[1:fidx-1])
        expn = @compat parse(Int,str[fidx+1:end])
        return convert(Float32, base * 10.0 ^ expn)
    elseif isfloat64
        return @compat parse(Float64,str)
    else
        try
            return @compat parse(Int64,str)
        catch ex
            # its better to ask for forgiveness...
            !isa(ex, OverflowError) && rethrow(ex)
            if is_within_int128(str)
                return @compat parse(Int128,str)
            else
                return @compat parse(BigInt,str)
            end
        end
    end
end

# accum digit predicates
is_char_hex(c::Char) = isdigit(c) || ('a' <= c <= 'f')  || ('A' <= c <= 'F')
is_char_oct(c::Char) = '0' <= c <= '7'
is_char_bin(c::Char) = c === '0' || c === '1'

function accum_digits(ts::TokenStream, pred::Function, c::Char, leading_zero::Bool)
    if !leading_zero && c == '_'
        return (Char[], false)
    end
    charr = Char[]
    while true
        if c == '_'
            skip(ts, 1)
            c = peekchar(ts)
            if !eof(c) && pred(c)
                continue
            else
                skip(ts, -1)
                break
            end
        elseif !eof(c) && pred(c)
            skip(ts, 1)
            push!(charr, c)
            c = peekchar(ts)
            continue
        end
        break
    end
    return (charr, true)
end

fix_uint_neg(neg::Bool, n::Number) = neg? -n : n

function disallow_dot!(ts::TokenStream, charr::Vector{Char}, loc)
    if peekchar(ts) === '.'
        skip(ts, 1)
        if is_dot_opchar(peekchar(ts))
            skip(ts, -1)
        else
            throw(diag(loc, "invalid numeric constant \"$(String(charr)).\""))
        end
    end
end

function read_digits!(ts::TokenStream, pred::Function, charr::Vector{Char}, leading_zero::Bool, loc)
    digits, ok = accum_digits(ts, pred, peekchar(ts), leading_zero)
    ok || throw(diag(loc, "invalid use of '_' in numeric constant"))
    isempty(digits) && return false
    append!(charr, digits)
    return true
end

#TODO: try to remove neg as it is not needed for the lexer
function read_number(ts::TokenStream, leading_dot::Bool, neg::Bool)
    r = startrange(ts)
    charr = Char[]
    pred::Function = isdigit

    leading_zero = false
    is_float32_literal  = false
    is_hexfloat_literal = false

    neg && push!(charr, '-')
    if leading_dot
        push!(charr, '.')
    else
        if peekchar(ts) == '0'
            push!(charr, readchar(ts))
            leading_zero = true
            nc = peekchar(ts)
            if nc === 'x'
                skip(ts, 1); push!(charr, nc)
                leading_zero = false
                pred = is_char_hex
            elseif nc === 'o'
                skip(ts, 1); push!(charr, nc)
                leading_zero = false
                pred = is_char_oct
            elseif nc === 'b'
                skip(ts, 1); push!(charr, nc)
                leading_zero = false
                pred = is_char_bin
            end
        else
            nc = peekchar(ts)
            nc === '.' && (skip(ts, 1); push!(charr, nc))
        end
    end
    read_digits!(ts, pred, charr, leading_zero, makerange(ts, r))
    if peekchar(ts) == '.'
        skip(ts, 1)
        if is_dot_opchar(peekchar(ts))
            skip(ts, -1)
        else
            push!(charr, '.')
            read_digits!(ts, pred, charr, false, makerange(ts, r))
            disallow_dot!(ts, charr, Lexer.makerange(ts, r))
        end
    end
    c = peekchar(ts)
    ispP = c === 'p' || c === 'P'
    if (pred === is_char_hex && ispP) || (c === 'e' || c === 'E' || c === 'f')
        skip(ts, 1)
        nc = peekchar(ts)
        if !eof(nc) && (isdigit(nc) || nc === '+' || nc === '-')
            skip(ts, 1)
            is_float32_literal = c === 'f'
            is_hexfloat_literal = ispP
            push!(charr, c)
            push!(charr, nc)
            read_digits!(ts, pred, charr, false, makerange(ts, r))
            disallow_dot!(ts, charr, makerange(ts, r))
        else
            skip(ts, -1)
        end
    # disallow digits after binary or octal literals, e.g. 0b12
    elseif (pred == is_char_bin || pred == is_char_oct) && !eof(c) && isdigit(c)
        push!(charr, c)
        throw(diag(makerange(ts, r), "invalid numeric constant \"$(String(charr))\""))
    end
    str = String(charr)
    base = pred == is_char_hex ? 16 :
           pred == is_char_oct ? 8  :
           pred == is_char_bin ? 2  : 10
    # for an unsigned literal starting with -,
    # remove the - and parse instead as a call to unary -
    (neg && base != 10 && !is_hexfloat_literal) && (str = str[2:end])
    loc = makerange(ts, r)
    try
        if is_hexfloat_literal
            return make_token((@compat parse(Float64,str)), loc)
        elseif pred == is_char_hex
            return make_token(fix_uint_neg(neg, sized_uint_literal(str, 4)), loc)
        elseif pred == is_char_oct
            return make_token(fix_uint_neg(neg, sized_uint_oct_literal(str)), loc)
        elseif pred == is_char_bin
            return make_token(fix_uint_neg(neg, sized_uint_literal(str, 1)), loc)
        elseif is_float32_literal
            return make_token(convert(Float32, string_to_number(str, loc)), loc)
        else
            return make_token(string_to_number(str, loc), loc)
        end
    catch e
        isa(e, ArgumentError) && throw(diag(loc, e.msg))
        rethrow(e)
    end
end

#============================#
# Skip whitespace / comments
#============================#

"""
`skip_multiline_comment(ts::TokenStream, count::Int)` -> `TokenStream`

Skips multiline comments and maintains a count of the number of enclosed `#= =#` pairs
to allow nesting of multi-line comments.

`count`: the number of read `(#=)` tokens, ie:

```
(#= test =#)  (#= test =#)  (#= #= test =# =#)
 ^              ^               ^        ^
cnt 0           cnt 1         cnt 2    cnt 1
```
"""
function skip_multiline_comment(ts::TokenStream, count::Int)
    start, unterminated = -1, true
    startloc = before(here(ts))
    while !eof(ts)
        c = readchar(ts)
        # if "=#" token, decrement the count.
        # The loop is exited when this count goes below zero.
        if c === '='
            start > 0 || (start = position(ts))
            if peekchar(ts) === '#' && position(ts) != start
                takechar(ts)
                count <= 1 || (count -= 1; continue)
                unterminated = false
                break
            end
            continue
        # if "#=" token increase count
        elseif c === '#' && peekchar(ts) === '='
            count += 1
        end
    end
    if unterminated
        D = diag(here(ts),"incomplete: unterminated multi-line comment #= ... =#")
        diag(D, startloc, "starting here", :note)
        throw(Incomplete(:comment, D))
   end
    return ts
end

"""
`skipcomment(ts::TokenStream)` -> `TokenStream`

Skips to the end if this is a mulitiline comment,
otherwise skips to end of line.
"""
function skipcomment(ts::TokenStream)
    @assert readchar(ts) === '#'
    if peekchar(ts) === '='
        skip_multiline_comment(ts, 1)
    else
        skip_to_eol(ts)
    end
    return ts
end

"""
`skipws_and_comments(ts::TokenStream)` -> `TokenStream`

Skips all whitespace before a comment, upon reaching the comment, if it is a
single line comment, skips to the end of the line, otherwise skip to the end
of the multiline comment block.
"""
function skipws_and_comments(ts::TokenStream)
    while !eof(ts)
        skipws(ts, true)
        peekchar(ts) !== '#' && break
        skipcomment(ts)
    end
    return ts
end

function accum_julia_symbol(ts::TokenStream, c::Char)
    nc, charr = c, Char[]
    while is_identifier_char(nc)
        c, nc = readchar(ts), peekchar(ts)
        # make sure that != is always an operator
        if c === '!' && nc === '='
            skip(ts, -1)
            break
        end
        push!(charr, c)
        eof(nc) && break
    end
    str = normalize_string(String(charr), :NFC)
    sym = Symbol(str)
    return sym === SYM_TRUE ? true : sym === SYM_FALSE ? false : sym
end

#= Token stream methods =#

make_token(::Type{Token},val,start,offset) = Token(val)
make_token(::Type{SourceLocToken},val,start,length) =
    SourceLocToken(val,start,length,0)
@noinline make_token(val, r::Void) = Token(val)
make_token(val, r::SourceRange) = SourceLocToken(val,r)

eof_token(ts) = Token(EOF_CHAR)
eof_token(ts::TokenStream{SourceLocToken}) = SourceLocToken(EOF_CHAR,here(ts))

macro tok(val)
    esc(quote
        loc = position(ts.io)
        v = $val
        make_token(T,v,loc,position(ts.io)-loc)
    end)
end

startrange(ts::TokenStream{SourceLocToken}) = position(ts.io)
makerange(ts::TokenStream{SourceLocToken}, r) = SourceRange(r,position(ts.io)-r,0)
nullrange(ts::TokenStream{SourceLocToken}) = SourceRange(position(ts.io),0,0)
function here(ts::TokenStream{SourceLocToken})
    if ts.putback !== nothing
        SourceRange(Tokens.normalize(√(ts.putback)).offset,1,0)
    elseif ts.lasttoken !== nothing
        SourceRange(Tokens.normalize(√(ts.lasttoken)).offset,1,0)
    else
        SourceRange(position(ts.io),1,0)
    end
end

startrange(ts::TokenStream) = nothing
makerange(ts::TokenStream, r) = nothing
nullrange(ts::TokenStream) = nothing
here(ts::TokenStream) = nothing

function next_token{T}(ts::TokenStream{T}, whitespace_newline::Bool)
    ts.ateof && return eof_token(ts)
    tmp = skipws(ts, whitespace_newline)
    eof(ts) && return eof_token(ts)
    ts.isspace = tmp
    while !eof(ts.io)
        c = peekchar(ts)
        if eof(c)
            ts.ateof = true
            return eof_token(T)
        elseif c === ' ' || c === '\t'
            skip(ts, 1)
            continue
        elseif c == '\r'
            skip(ts, 1)
            peekchar(ts) === '\n' || throw(diag(here(ts),"'\\r' not followed by '\\n' is invalid"))
            continue
        elseif c === '#'
            skipcomment(ts)
            if whitespace_newline && peekchar(ts) === '\n'
                takechar(ts)
            end
            continue
        elseif isnewline(c)
            return @tok readchar(ts)
        elseif is_special_char(c)
            return @tok readchar(ts)
        elseif isdigit(c)
            return read_number(ts, false, false)
        elseif c === '.'
            r = startrange(ts)
            skip(ts, 1)
            nc = peekchar(ts)
            if isdigit(nc)
                return read_number(ts, true, false)
            elseif is_opchar(nc)
                return @tok begin
                    op = read_operator(ts, c, r)
                    if op === :(..) && is_opchar(peekchar(ts))
                        throw(diag(makerange(ts, r),
                            string("invalid operator \"", op, peekchar(ts), "\"")))
                    elseif op == Symbol(".!")
                        throw(diag(makerange(ts, r),
                            string("invalid operator \"", op, "\"")))
                    end
                    op
                end
            end
            return @tok :(.)
        elseif is_opchar(c)
            r = startrange(ts)
            return @tok read_operator(ts, readchar(ts), r)
        elseif is_identifier_start_char(c)
            return @tok accum_julia_symbol(ts, c)
        else
            @assert readchar(ts) === c
            if is_ignorable_char(c)
                throw(diag(here(ts),"invisible character \\u$(hex(c))"))
            else
                throw(diag(here(ts),"invalid character \"$c\""))
            end
        end
    end
    ts.ateof = true
    return eof_token(ts)
end

next_token(ts::TokenStream) = next_token(ts, false)
last_token(ts::TokenStream) = ts.lasttoken
set_token!(ts::TokenStream, t::Union{AbstractToken,Void}) = (ts.lasttoken = t; ts)

function put_back!(ts::TokenStream, t::Union{AbstractToken,Void})
    if ts.putback !== nothing
        error("too many pushed back tokens (internal error)")
    end
    ts.putback = t
    return ts
end

function peek_token{T}(ts::TokenStream{T}, whitespace_newline::Bool)
    ts.ateof && return eof_token(ts)
    if ts.putback !== nothing
        return ts.putback::AbstractToken
    end
    lt = last_token(ts)
    if lt !== nothing
        return lt::AbstractToken
    end
    set_token!(ts, next_token(ts, whitespace_newline))
    return last_token(ts)::AbstractToken
end
peek_token(ts::TokenStream) = peek_token(ts, false)::AbstractToken

function take_token{T}(ts::TokenStream{T})
    ts.ateof && return eof_token(ts)
    if ts.putback !== nothing
        t = ts.putback
        ts.putback = nothing
    else
        t = last_token(ts)
        set_token!(ts, nothing)
    end
    return t
end

function require_token(ts::TokenStream, whitespace_newline::Bool)
    if ts.putback !== nothing
        t = ts.putback
    elseif ts.lasttoken !== nothing
        t = ts.lasttoken
    else
        t = next_token(ts, whitespace_newline)
    end
    eof(t) && throw(Incomplete(:other, diag(here(ts), "incomplete: premature end of input")))
    while isnewline(t)
        take_token(ts)
        t = next_token(ts, whitespace_newline)
    end
    eof(t) && throw(Incomplete(:other, diag(here(ts), "incomplete: premature end of input")))
    ts.putback === nothing && set_token!(ts, t)
    return t
end
require_token(ts::TokenStream) = require_token(ts, false)

end


text:
# compute the sparsity pattern of the gradient of an expression
# (i.e., a list of which indices are present)

function compute_gradient_sparsity(nd::Vector{NodeData})

    # this is easy, just build a list of which variable indices appear
    indices = Set{Int}()
    for k in 1:length(nd)
        nod = nd[k]
        if nod.nodetype == VARIABLE
            push!(indices, nod.index)
        end
    end

    return indices
end

function compute_gradient_sparsity!(
    indices::Coloring.IndexedSet,
    nd::Vector{NodeData},
)
    for k in 1:length(nd)
        nod = nd[k]
        if nod.nodetype == VARIABLE
            push!(indices, nod.index)
        elseif nod.nodetype == MOIVARIABLE
            error(
                "Internal error: Invalid to compute sparsity if MOIVARIABLE nodes are present.",
            )
        end
    end
    return nothing
end

export compute_gradient_sparsity, compute_gradient_sparsity!

# compute the sparsity pattern the hessian of an expression

# input_linearity is the linearity with respect to the input
# computed by classify_linearity
# subexpression_edgelist is the edgelist of each subexpression
# subexpression_variables is the list of all variables which appear in
# a subexpression (including recursively)
function compute_hessian_sparsity(
    nd::Vector{NodeData},
    adj,
    input_linearity::Vector{Linearity},
    indexedset::Coloring.IndexedSet,
    subexpression_edgelist::Vector{Set{Tuple{Int,Int}}},
    subexpression_variables::Vector{Vector{Int}},
)

    # idea: consider the linearity/nonlinearity of a node *with respect to the output*
    # The children of any node which is nonlinear with respect to the output
    # should have nonlinear interactions, hence nonzeros in the hessian.
    # This is not true in general but holds for everything we consider.
    # Counterexample: f(x,y,z) = x + y*z, but we don't have any functions like that.
    # By "nonlinear with respect to the output", we mean that the output
    # depends nonlinearly on the value of the node, regardless of
    # how the node itself depends on the input.

    # So start at the root of the tree and classify the linearity wrt the output.
    # For each nonlinear node, do a mini DFS and collect the list of children.
    # Add a nonlinear interaction between all children of a nonlinear node.

    edgelist = Set{Tuple{Int,Int}}()
    nonlinear_wrt_output = fill(false, length(nd))

    children_arr = rowvals(adj)

    stack = Int[]
    stack_ignore = Bool[]
    nonlinear_group = indexedset
    if length(nd) == 1 && nd[1].nodetype == SUBEXPRESSION
        # subexpression comes in linearly, so append edgelist
        for ij in subexpression_edgelist[nd[1].index]
            push!(edgelist, ij)
        end
    end

    for k in 2:length(nd)
        nod = nd[k]
        if nod.nodetype == MOIVARIABLE
            error(
                "Internal error: Invalid to compute sparsity if MOIVARIABLE nodes are present.",
            )
        end
        nonlinear_wrt_output[k] && continue # already seen this node one way or another
        input_linearity[k] == CONSTANT && continue # definitely not nonlinear

        @assert !nonlinear_wrt_output[nod.parent]
        # check if the parent depends nonlinearly on the value of this node
        par = nd[nod.parent]
        if par.nodetype == CALLUNIVAR
            op = par.index
            if op >= USER_UNIVAR_OPERATOR_ID_START ||
               univariate_operators[op] != :+ && univariate_operators[op] != :-
                nonlinear_wrt_output[k] = true
            end
        elseif par.nodetype == CALL
            op = par.index
            if op >= USER_OPERATOR_ID_START
                nonlinear_wrt_output[k] = true
            elseif operators[op] == :+ ||
                   operators[op] == :- ||
                   operators[op] == :ifelse
                # pass
            elseif operators[op] == :*
                # check if all siblings are constant
                sibling_idx = nzrange(adj, nod.parent)
                if !all(
                    i ->
                        input_linearity[children_arr[i]] == CONSTANT ||
                            children_arr[i] == k,
                    sibling_idx,
                )
                    # at least one sibling isn't constant
                    nonlinear_wrt_output[k] = true
                end
            elseif operators[op] == :/
                # check if denominator is nonconstant
                sibling_idx = nzrange(adj, nod.parent)
                if input_linearity[children_arr[last(sibling_idx)]] != CONSTANT
                    nonlinear_wrt_output[k] = true
                end
            else
                nonlinear_wrt_output[k] = true
            end
        end

        if nod.nodetype == SUBEXPRESSION && !nonlinear_wrt_output[k]
            # subexpression comes in linearly, so append edgelist
            for ij in subexpression_edgelist[nod.index]
                push!(edgelist, ij)
            end
        end

        nonlinear_wrt_output[k] || continue

        # do a DFS from here, including all children
        @assert isempty(stack)
        @assert isempty(stack_ignore)
        sibling_idx = nzrange(adj, nod.parent)
        for sidx in sibling_idx
            push!(stack, children_arr[sidx])
            push!(stack_ignore, false)
        end
        empty!(nonlinear_group)
        while length(stack) > 0
            r = pop!(stack)
            should_ignore = pop!(stack_ignore)
            nonlinear_wrt_output[r] = true
            if nd[r].nodetype == LOGIC || nd[r].nodetype == COMPARISON
                # don't count the nonlinear interactions inside
                # logical conditions or comparisons
                should_ignore = true
            end
            children_idx = nzrange(adj, r)
            for cidx in children_idx
                push!(stack, children_arr[cidx])
                push!(stack_ignore, should_ignore)
            end
            should_ignore && continue
            if nd[r].nodetype == VARIABLE
                push!(nonlinear_group, nd[r].index)
            elseif nd[r].nodetype == SUBEXPRESSION
                # append all variables in subexpression
                union!(nonlinear_group, subexpression_variables[nd[r].index])
            end
        end
        for i_ in 1:nonlinear_group.nnz
            i = nonlinear_group.nzidx[i_]
            for j_ in 1:nonlinear_group.nnz
                j = nonlinear_group.nzidx[j_]
                j <= i || continue # only lower triangle
                push!(edgelist, (i, j))
            end
        end
    end

    return edgelist
    #=
    nz = length(edgelist)
    I = Int[]
    J = Int[]
    for (i,j) in edgelist
        push!(I,i)
        push!(J,j)
    end
    return I,J
    =#

end

export compute_hessian_sparsity

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
@named s = sesam3_revMM()
@named pl = plant_const()

# @named s = Sesam.sesam3C()
# sp = compose(ODESystem([
#     s.i_L ~ pl.i_L,
#     # colimitation
#     s.syn_B ~ s.C_synBC,
#   ], t; name=:sp), s, pl)
#sp = structural_simplify(sp)

@named sp = plant_sesam_system(s,pl)
states(sp)

p = pC = Dict(
    s.ϵ_tvr => 0.45,   # carbon use efficiency of microbial tvr (part by predators
        #which respire and corresponding amount of N must be mineralized)
    s.κ_E => 0.8,     ##<< amount of recycling enzyme turnover by biomass (
        # added to assimilable, i.e. uptake instead of R)
    s.a_E => 0.001*365,   ##<< C biomass allocated to enzymes 1/day /microbial biomass
    s.m => 0.005*365,    ##<< maintenance respiration rate   1/day /microbial biomass,
    s.τ => 1/60*365,  ##<< biomass turnover rate (12 days)
    s.k_L => 5.0,       ##<< 1/(x years)   # formerly 1 year
    s.k_R => 1/(20.0),        ##<< 1/(x years) # to demonstrate changes on short time scale
    s.k_mN_L => 0.05 * 60, # enzyme half-saturation constant, in magnitude of enzymes *
        # /yr enzyme turnover 60 times a year
    s.k_mN_R => 0.05 * 60, # enzyme half-saturation constant, in magnitude of enzymes *
    # /yr enzyme turnover 60 times a year
    s.ϵ => 0.5,      ##<< carbon use efficiency for growth respiration
    #i_L => t -> 1 - exp(-t),  # litter input
    pl.i_L0 => 400.0,         # g/m2 input per year (half NPP)
    pl.β_Ni0 => 25,
    pl.i_IN0 => 0,   ##<< input of mineral N,
    #
    # P from plant model parameters not used in CN-Sesam soil model
    pl.β_Pi0 => Inf, #25*20, ## leaf litter N:P ~20(massratio Kang10)
    pl.i_IP0 => Inf, #0.65,   ##<< input of mineral P, weathering: Table3 mixed sedimentary rocks 0.65g/m2/yr Hartmann14 10.1016/j.chemgeo.2013.10.025
    pl.s_EP0 => Inf, # 0.5, # plant 1/20 of typical total microbial enzyme synthesis flux
    pl.u_PlantPmax0 => Inf,
    pl.k_PlantP0 => Inf,
)
pN = Dict(
    s.i_BN => 0.4, ##<< potential immobilization flux rate
    s.β_NEnz => 3.1,     # Sterner02: Protein (Fig. 2.2.), high N investment (low P) need
    s.β_NB => 11.0,
    #s.l_N => 0.96,       #0.00262647*365     ##<< leaching rate of mineralN lN IN
    s.l_N => 0.0,
    s.ν_N =>  0.9,     # microbial N use efficiency accounting for apparent
    ## minceralization of N during uptake in heterogeneous soils
    # for N uptake: take the defaults which take as much N as supplied by litter
    # pl.u_PlantNmax0 => Inf32, # only constrained by k_PlantN in min function
    # pl.k_PlantN0 => 10.57, #0.0289652*365     ##<< plant uptake rate first order of IN
)
p = merge(pC, pN)

u0 = u0C = Dict(
    s.B => 17,
    s.L => 100,
    s.R => 1100,
    #s.cumresp => 0.0,
    s.α_R => 0.1, # enzyme synthesis into L # TODO model by optimality
)
u0C[s.α_L] = 1.0 - u0C[s.α_R]
u0N = Dict(
    s.I_N => 0.04, ##<< inorganic pool gN/m2
    s.L_N => u0[s.L]/p[pl.β_Ni0],
    s.R_N => u0[s.R]/calculate_β_NR_sesam3(p,s) #p[s.β_NB],
    )
u0 = merge(u0C, u0N)
#u0[s.R]/u0[s.R_N] # smaller p[s.β_NB]

tspan = (0.0,100.0)

#prob = ODEProblem(sp,[t for t in u0], tspan, [t for t in p])
#prob = ODEProblem(sp, remove_units(u0), tspan, remove_units(p))
prob = ODEProblem(sp, u0, tspan, p)
#prob = ODEProblem(sp,u0, tspan, p, jac=true)
sol = sol_sesam3 = solve(prob);

i_plot = () -> begin
    #using Plots
    plot(sol)
    plot(sol, vars=[s.R])
end

@testset "non-negative pools" begin
    #st = first(states(sp))
    for st in states(sp)
        @test all(sol[st] .>= 0.0)
    end
end;

@testset "microbial dC balance" begin
    #plot(sol,  vars=[s.u_C, s.r_B+s.syn_B+s.syn_Enz])
    @test all(isapprox.(
        sol[s.u_C],
        sol[s.r_B+s.syn_B+s.syn_Enz],
        rtol = 1e-6))
end;

@testset "system dC balance" begin
    change = sol[s.dB + s.dL + s.dR]
    output = sol[s.r_tot]
    input = sol[s.i_L]
    #plot(sol.t, [input, output, change, change .+ output])
    @test all(isapprox.(input, change .+ output, rtol = 1e-6))
end;

@testset "microbial dN balance" begin
    uptake = p[s.ν_N]*sol[s.u_NOM]
    usage = sol[s.syn_Enz]/p[s.β_NEnz] + sol[s.syn_B]/p[s.β_NB] + sol[s.Φ_NB]
    #plot(sol.t, [uptake, usage])
    @test all(isapprox.(uptake, usage, rtol = 1e-6))
end;

@testset "system dN balance" begin
    change = sol[s.dB]/p[s.β_NB] + sol[s.dL_N + s.dR_N + s.dI_N]
    output = sol[s.u_PlantN + s.leach_N]
    input = sol[s.i_L/s.β_Ni + s.i_IN]
    #sol[s.leach_N]
    #sol[s.i_IN]
    bo = 1:10
    bo = 1:length(input)
    #plot(sol.t[bo], [input[bo], output[bo], change[bo], (change .+ output)[bo]], label = ["input" "output" "change" "change+output"])
    @test all(isapprox.(input, change .+ output, rtol = 1e-6))
end;

@testset "balanced allocation sums to one" begin
    #plot(sol, vars=[s.α_LT,s.α_RT, s.α_LT+s.α_RT])
    @test all(isapprox.(sol[s.α_LT + s.α_RT], 1.0, rtol = 1e-8))
    @test all(isapprox.(sol[s.α_L + s.α_R], 1.0, rtol = 1e-8))
end;


text:
"""
    FunctionTransform(f)

Take a function or object `f` as an argument which is going to act on each vector individually.
Make sure that `f` is supposed to act on a vector.
For example replace `f(x)=sin(x)` by `f(x)=sin.(x)`
```
    f(x) = abs.(x)
    tr = FunctionTransform(f)
```
"""
struct FunctionTransform{F} <: Transform
    f::F
end

@functor FunctionTransform

(t::FunctionTransform)(x) = t.f(x)

_map(t::FunctionTransform, x::AbstractVector{<:Real}) = map(t.f, x)

function _map(t::FunctionTransform, x::ColVecs)
    vals = map(axes(x.X, 2)) do i
        t.f(view(x.X, :, i))
    end
    return ColVecs(reduce(hcat, vals))
end

function _map(t::FunctionTransform, x::RowVecs)
    vals = map(axes(x.X, 1)) do i
        t.f(view(x.X, i, :))
    end
    return RowVecs(reduce(hcat, vals)')
end

duplicate(t::FunctionTransform,f) = FunctionTransform(f)

Base.show(io::IO, t::FunctionTransform) = print(io, "Function Transform: ", t.f)


text:
a = CSV.read("d:/data\\AirOnTimeCSV\\airOT199401.csv")

bb = [Array(b) for b in eachcol(a)]

bbc = Blosc.compress(bb[1])

bbc = Blosc.compress([1,2,3])
io = open("io.test", "w")
write(io, bbc)
write(io, bbc)
close(io)

io = open("io.test", "r")
bbc_new  = Array{UInt8,1}(undef, length(bbc))
bbc_new2  = similar(bbc_new)

read!(io, bbc_new)
read!(io, bbc_new2)

close(io)

Blosc.decompress(Int64, bbc_new)
Blosc.decompress(Int64, bbc_new2)

all(Blosc.decompress(Int64, bbc_new) .== [1,2,3])

using CSV, DataFrames, Bloss, Missing
# a = DataFrame(
# 	a = rand(Int32, 100_000_000), b = rand(Float32, 100_000_000)
# 	)

@time a = CSV.read("d:/data/AirOnTimeCSV/airOT199401.csv")

max_stars_repo_name: xiaodaigh/JDF

max_issues_repo_name: xiaodaigh/JDF

max_forks_repo_name: xiaodaigh/JDF


text:
"""
    AbstractVariableSelection{TakeObjective}

Abstract type for the variable selection. TakeObjective is a boolean saying if one can branch
on the objective or not.
"""
abstract type AbstractVariableSelection{TakeObjective} end

include("mindomain.jl")
include("random.jl")
include("failureBased.jl")


text:
using SafeTestsets

@time begin
    @time @safetestset "FIGTree Tests" begin include("fgt_test.jl") end
    @time @safetestset "Forces Tests" begin include("forces_test.jl") end
    @time @safetestset "Cell List Tests" begin include("celllists_test.jl") end
    @time @safetestset "Check Position Tests" begin include("position_checking_test.jl") end
    @time @safetestset "Full Run Tests" begin include("full_run_test.jl") end
    @time @safetestset "Struct Tests" begin include("struct_test.jl") end
end


text:
import accrue.cryptoerase.runtime.Condition;

class C {
    C(Condition c) { }

    public static void main(String[] args) {}
}

max_stars_repo_name: HarvardPL/cryptoerase

max_issues_repo_name: HarvardPL/cryptoerase

max_forks_repo_name: HarvardPL/cryptoerase


text:
"""
    NoiseSource(eltype, samplerate, nchannels, std=1)

NoiseSource is a multi-channel noise signal generator. The noise on each channel is independent.

Inputs
------
* `samplerate` specifies the sample rate of the signal specified in Hz.
* `nchannels` specifies the number of channels of the signal.
* `std` specifies the desired standard deviation of the signal.

Output
------
* SampleSource object

Example
-------
```julia
source_object = NoiseSource(Float64, 48u"kHz", 2, 0.3)
wn = read(source_object, 480)         # Specify number of samples of signal to generate
wn = read(source_object, 50u"ms")     # Specify length of time of signal to generate
```
"""
mutable struct NoiseSource{T} <: SampleSource
    samplerate::Float64
    nchannels::Int64
    std::Float64

    function NoiseSource(eltype, samplerate::Number, nchannels::Int, std::Number=1)
        new{eltype}(samplerate, nchannels, std)
    end
    function NoiseSource(eltype, samplerate::Unitful.Frequency, nchannels::Int, std::Number=1)
        samplerate = ustrip(uconvert(u"Hz", samplerate))
        NoiseSource(eltype, samplerate, nchannels, std)
    end
end

Base.eltype(::NoiseSource{T}) where T = T
nchannels(source::NoiseSource) = source.nchannels
samplerate(source::NoiseSource) = source.samplerate

function unsafe_read!(source::NoiseSource, buf::Array, frameoffset, framecount)
    buf[1+frameoffset:framecount+frameoffset, 1:source.nchannels] = source.std .* randn(framecount, source.nchannels)
    framecount
end


text:
using CryptoGroups: Group, PGroup, Hash

Base.write(io::IO, tree::Tree) = write(io, encode(tree))
Base.write(f::AbstractString, tree::Tree) = write(f, encode(tree))

Base.string(tree::Tree) = bytes2hex(encode(tree))

function unmarshal_full_public_key(g::Group, tree::Tree)

    G = typeof(g)

    g′, y = convert(Tuple{G, G}, tree)

    @assert g′ == g

    return y
end

marshal_full_public_key(g::G, y::G) where G <: Group = Tree((g, y))

function unmarshal_publickey(tree::Tree; relative::Bool = false)

    g = unmarshal(tree.x[1])
    G = typeof(g)

    g′, y = convert(Tuple{G, G}, tree.x[2])

    if !relative
        @assert g′ == g "Generator does not match specification of the group. Perhaps intentioanl, if so pass `relative=true` as keyword argument."
    end

    return y, g′
end

function marshal_publickey(y::G, g::G) where G <: Group

    group_spec = marshal(g)

    L = bitlength(G) #

    g_tree = Tree(g; L)
    y_tree = Tree(y; L)

    public_key = Tree((g_tree, y_tree))
    tree = Tree((group_spec, public_key))

    return tree
end

function marshal_privatekey(g::Group, s::BigInt)
    group_spec = marshal(g)

    q = order(g)

    @assert s < q "Secret key must be with in the order of the group"

    sleaf = Leaf(s, bytelength(q))

    tree = Tree((group_spec, sleaf))

    return tree
end

function unmarshal_privatekey(tree::Tree)

    g = unmarshal(tree.x[1])
    s = convert(BigInt, tree.x[2])

    return (s, g) ### The group can often be omited when not needed.
end

function map_hash_name(x::AbstractString)
    if x == "SHA-256"
        return "sha256"
    elseif x == "SHA-384"
        return "sha384"
    elseif x == "SHA-512"
        return "sha512"
    else
        error("No corepsonding mapping for $x implemented")
    end
end

function map_hash_name_back(x::AbstractString)
    if x == "sha256"
        return "SHA-256"
    elseif x == "sha384"
        return "SHA-384"
    elseif x == "sha512"
        return "SHA-512"
    else
        error("No corepsonding mapping for $x implemented")
    end
end

map_hash_name_back(x::Hash) = map_hash_name_back(x.spec)

function ro_prefix(protinfo::AbstractDict; auxsid="default")

    version = protinfo["version"]
    sid = protinfo["sid"]

    s_H = protinfo["rohash"]
    s_PRG = protinfo["prg"]
    s_Gq = protinfo["pgroup"]

    nr = parse(Int32, protinfo["statdist"])
    nv = parse(Int32, protinfo["vbitlenro"])
    ne = parse(Int32, protinfo["ebitlenro"])

    data = (version, sid * "." * auxsid, nr, nv, ne, s_PRG, s_Gq, s_H)

    tree = Tree(data)
    binary = encode(tree)

    rohash = Hash(map_hash_name(protinfo["rohash"]))

    ρ = rohash(binary) ### Which hash function shall be used here?

    return ρ
end