text:
module TestImages

export indicator_matrix, indicator_matrix_2
export AbstractImageParams, ImageParams, CircleParams
export circle_position, square_position
export calibration_position, background_position
export gray_scale_indices
export circle_image, gray_scale_image
export square_image, circle_polar_image
export combined_images_size, combine_images
export create_image, plateau_length
export AbstractTestImage, AbstractGrayScale
export GrayScaleLine, GrayScalePyramid, WhiteRect
export CircleImage, SquareImage

using IntervalSets, SimpleTraits
using Statistics: mean
import ..Monads: mjoin, mreturn
using ..Monads
import ..CTImages: rescale, rescale!
using ..CTImages
using ..Interpolation: NoInterpolation
import ..Geometry: ParallelBeamGeometry, FanBeamGeometry, AbstractParallelBeamGeometry
import ..Geometry
using ..AbstractAlgorithms: AbstractProjectionAlgorithm
import ..Utils: _atype, yaml_repr, struct2dict, linspace
import ..Info: CTInfo
import ..AbstractAlgorithms: project_image
import ..CalibrationBase:
    calibrate_image,
    calibrate_image!,
    calibrate_tomogram,
    calibrate_tomogram!,
    calibration_data
import Base: getproperty, size, show, getindex, setindex!, IndexStyle, eltype
using Base: @propagate_inbounds

abstract type AbstractImageParams{T<:Real} end

eltype(::Type{P}) where {P<:AbstractImageParams{T}} where T = T
eltype(imp::AbstractImageParams) = eltype(typeof(imp))
size(p::AbstractImageParams) = p.rows, p.cols
function size(p::AbstractImageParams, dim::Integer)
    dim == 1 && return p.rows
    dim == 2 && return p.cols
    return 1
end

yaml_repr(imp::AbstractImageParams) = struct2dict(imp)
CTInfo(imp::AbstractImageParams) = CTInfo(pairs(struct2dict(imp))...)

function indicator_matrix(::Type{T}, sq::Int, pad::Int = 0) where {T<:Real}
    if pad == 0
        pad = sq + 1
    end
    rows = 2 * (sq + pad)
    cols = 2 * (sq + pad)
    matrix = zeros(T, rows, cols)
    matrix[pad+sq+1:pad+2sq, pad+1:pad+sq] .= 1.0
    matrix[pad+1:pad+sq, pad+sq+1:pad+2sq] .= 1.0
    return matrix
end

indicator_matrix(sq::Int) = indicator_matrix(Float32, sq)

function indicator_matrix_2(
    ::Type{T},
    sq1::Int,
    sq2::Int,
    pad::Int = 0,
) where {T<:Real}
    sqmin, sqmax = minmax(sq1, sq2)
    if pad == 0
        pad = sqmin + 1
    end
    rows = 2 * (sqmax + pad)
    cols = 2 * (sqmax + pad)
    matrix = zeros(T, rows, cols)
    matrix[pad+1:pad+sqmin, pad+1:pad+sqmin] .= 1.0
    matrix[pad+sqmax+1:pad+2sqmax, pad+sqmax+1:pad+2sqmax] .= 1.0
    return matrix
end

indicator_matrix_2(sq1::Int, sq2::Int, pad::Int = 0) =
    indicator_matrix_2(Float32, sq1, sq2, pad)

const _default_gray_scale_params = Dict(
    :swidth => 200, # Gray scale width
    :sheight => 40, # Gray scale height
    :pad => 30,     # Padding around gray scale and calibration circle
    :dist => 10,    # Distance between gray scale and calibration circle
    :radius => 15,  # Calibration circle radius
)

"""
    struct ImageParams{T<:Real}

Hold information to construct the input test image.
"""
struct ImageParams{T} <: AbstractImageParams{T}
    rows::Int # Test image rows
    cols::Int # Test image cols
    gray_scale::ClosedInterval{T} # Gray scale interval.
    calibration_value::T # Value for algorithm calibration.
    background::T
    kw::Dict{Symbol,Any}

    function ImageParams{T}(
        rows::Optional{Integer} = nothing,
        cols::Optional{Integer} = nothing;
        gray_scale::Optional{ClosedInterval} = nothing,
        calibration_value::Optional{Real} = nothing,
        background::Optional{Real} = nothing,
        hounsfield::Bool = true,
        kwargs...
    ) where {T<:Real}
        gray_scale = maybe(ifelse(hounsfield, -1000..1000, 0..1), gray_scale)
        kw = isempty(kwargs) ? _default_gray_scale_params : Dict(kwargs)
        srows, scols = combined_images_size(; kw...)
        rows, cols = maybe(srows, rows), maybe(scols, cols)
        @assert(
            rows >= srows && cols >= scols,
            "Requested image size is too small: " *
            "rows($rows) ≥ $srows && columns($cols) ≥ $scols"
        )
        background = maybe(minimum(gray_scale), background)
        calibration_value = isnothing(calibration_value) ?
            mean(gray_scale) : calibration_value

        if calibration_value ∉ background..maximum(gray_scale)
            midpoint = mean(gray_scale)
            @warn(
                "Calibration value should be in the interval of the gray scale ($calibration_value ∉ $gray_scale), taking midpoint: $midpoint"
            )
            calibration_value = midpoint
        end

        new(
            rows,
            cols,
            gray_scale,
            calibration_value,
            background,
            kw,
        )
    end
end

function getproperty(p::ImageParams, s::Symbol)
    s ∈ fieldnames(ImageParams) && return getfield(p, s)
    s === :width && return p.cols
    s === :height && return p.rows
    p.kw[s]
end

"""
    ImageParams([T=Float32]; <keyword arguments>) where {T<:Real}

Construct ImageParams object.

# Arguments:
- `width=200`: width of the rectangle gray scale image.
- `height=40`: height of the rectangle gray scale image.
- `pad=30`: padding around image.
- `dist=10`: distance between images.
- `radius=20`: radius of calibration circle.
- `gray_scale=1000..1000`: interval of gray scale values.
- `calibration_value=nothing`: value of the calibration circle. If not specified, is `gray_scale[2]`.
- `background=nothing`: value to be used as background.
- `hounsfield=false`: whether the gray scale should be in Hounsfield units.
"""
function ImageParams(
    ::Type{T} = Float32;
    rows::Optional{Integer} = nothing,
    cols::Optional{Integer} = nothing,
    width::Optional{Integer} = nothing,
    height::Optional{Integer} = nothing,
    kwargs...
) where {T<:Real}
    rows = maybe(rows, height)
    cols = maybe(cols, width)
    ImageParams{T}(rows, cols; kwargs...)
end

show(io::IO, p::AbstractImageParams) = print(
    io,
    """Image Params:
- size (W×H): $(p.width) × $(p.height)
- gray scale: $(p.gray_scale)
- calibration value: $(p.calibration_value)
- background value: $(p.background)""",
)

struct CircleParams{T} <: AbstractImageParams{T}
    radius::Int # Calibration circle radius
    rows::Int # Test image rows
    cols::Int # Test image cols
    nϕ::Int # Number of angles to compute the circle
    gray_scale::ClosedInterval{T} # Gray scale interval.
    calibration_value::T # Value for algorithm calibration.
    background::T

    function CircleParams{T}(
        radius::Integer,
        rows::Integer,
        cols::Integer,
        nϕ::Optional{Integer} = nothing;
        gray_scale::Optional{ClosedInterval} = nothing,
        calibration_value::Optional{Real} = nothing,
        background::Optional{Real} = nothing,
        hounsfield::Bool = false,
    ) where {T<:Real}
        gray_scale = maybe(ifelse(hounsfield, -1000..0, 0..1), gray_scale)
        background = maybe(minimum(gray_scale), background)
        calibration_value = isnothing(calibration_value) ?
            maximum(gray_scale) : calibration_value

        if calibration_value ∉ background..maximum(gray_scale)
            value = maximum(gray_scale)
            @warn "Calibration value should be in the interval of the gray scale ($calibration_value ∉ $gray_scale), taking maximum: $value"
            calibration_value = value
        end

        new(
            radius,
            rows,
            cols,
            maybe(max(rows,cols), nϕ),
            gray_scale,
            calibration_value,
            background,
        )
    end
end

function getproperty(p::CircleParams, s::Symbol)
    s ≡ :width && return p.cols
    s ≡ :height && return p.rows
    s ≡ :nphi && return getfield(p, nϕ)
    getfield(p, s)
end

function CircleParams(
    ::Type{T} = Float32;
    radius::Optional{Integer} = nothing,
    width::Optional{Integer} = nothing,
    height::Optional{Integer} = nothing,
    rows::Optional{Integer} = nothing,
    cols::Optional{Integer} = nothing,
    nϕ::Optional{Integer} = nothing,
    nphi::Optional{Integer} = nothing,
    kwargs...
) where {T <: Real}
    nϕ = maybe(nϕ, nphi)
    #factor = 31 / 50 # magic number!
    factor = 1
    #circ_zoom = 1.3 / 17
    circ_zoom = 1//5
    #=
        If width or height are present, use them, otherwise use rows and cols.
        Finally, if rows or cols are not present, compute defaults.
    =#
    rows, cols = maybe(rows, height), maybe(cols, width)
    width = maybe(512, cols)
    height = isnothing(rows) ? round(Int, width * factor) : rows
    radius = isnothing(radius) ? round(Int, width * circ_zoom) : radius
    rows, cols = height, width
    CircleParams{T}(radius, rows, cols, nϕ; kwargs...)
end

struct SquareParams{T} <: AbstractImageParams{T}
    l::Int
    rows::Int
    cols::Int
    gray_scale::ClosedInterval{T}
    calibration_value::T
    background::T

    function SquareParams{T}(
        l::Int,
        rows::Integer,
        cols::Integer;
        gray_scale::Optional{ClosedInterval} = nothing,
        calibration_value::Optional{Real} = nothing,
        background::Optional{Real} = nothing,
        hounsfield::Bool = false,
    ) where {T <: Real}
        gray_scale = maybe(ifelse(hounsfield, -1000..0, 0..1), gray_scale)
        background = maybe(minimum(gray_scale), background)
        calibration_value = isnothing(calibration_value) ?
            maximum(gray_scale) : calibration_value

        if calibration_value ∉ background..maximum(gray_scale)
            value = maximum(gray_scale)
            @warn "Calibration value should be in the interval of the gray scale ($calibration_value ∉ $gray_scale), taking maximum: $value"
            calibration_value = value
        end
        new(l, rows, cols, gray_scale, calibration_value, background)
    end
end

function getproperty(p::SquareParams, s::Symbol)
    s ≡ :width && return p.cols
    s ≡ :height && return p.rows
    getfield(p, s)
end

function SquareParams(
    ::Type{T} = Float32;
    l::Optional{Integer} = nothing,
    rows::Optional{Integer} = nothing,
    cols::Optional{Integer} = nothing,
    width::Optional{Integer} = nothing,
    height::Optional{Integer} = nothing,
    kwargs...
) where {T <: Real}
    rows, cols = maybe(rows, height), maybe(cols, width)
    width = maybe(maybe(512, rows), cols)
    height = maybe(width, rows)
    rows, cols = height, width
    l = maybe(min(rows, cols) ÷ 2, l)
    SquareParams{T}(l, rows, cols; kwargs...)
end

"""
    circle_position(imp::AbstractImageParams)

Return circle position inside image given image parameters.
"""
function circle_position end

function circle_position(imp::ImageParams)
    srows, _ = combined_images_size(imp)
    rows, cols = size(imp)
    r = (rows - srows) ÷ 2 + imp.pad + imp.radius + 1
    c = cols ÷ 2
    r, c
end
circle_position(imp::CircleParams) = imp.rows ÷ 2, imp.cols ÷ 2

"""
    square_position(imp::SquareParams)

Return square position inside image given image parameters.
"""
function square_position(imp::SquareParams)
    imp.rows ÷ 2, imp.cols ÷ 2
end

"""
    background_position(imp::AbstractImageParams)

Return suitable position to calibrate background.
"""
function background_position end

function background_position(imp::ImageParams)
    cr, cc = circle_position(imp)
    cr, cc ÷ 2
end

function background_position(imp::CircleParams)
    cr, cc = circle_position(imp)
    max(1, cr - imp.radius), max(1, cc - imp.radius)
end

function background_position(imp::SquareParams)
    imp.l == imp.rows == imp.cols && @warn "Cannot calibrate square image if square side equals image size"
    cr, cc = square_position(imp)
    max(1, cr), max(1, (1 + cc - imp.l) ÷ 2)
end

"""
    gray_scale_indices(imp::ImageParams)

Return a tuple `(row, column_range)` where `row` is the row index of the scale and `column_range` is the range of columns.
"""
function gray_scale_indices(imp::ImageParams)
    rows, cols = size(imp)
    srows, scols = combined_images_size(imp)
    min_row = (rows - srows) ÷ 2 + imp.pad + 2 * imp.radius + imp.dist + 1
    max_row = min_row + imp.sheight - 1
    min_col = (cols - scols) ÷ 2 + imp.pad + 1
    max_col = min_col + imp.swidth - 1
    min_row:max_row, min_col:max_col
end

"""
    circle_image([T=Float32]; radius=20, value=1) where {T<:Real}

Create a square image with a circle of given value.

# Examples
```julia-repl
julia> circle_image(30, 0.8)
30×30 Array{Float32,2}:
[...]
```

See also: [`gray_scale_image`](@ref), [`combine_images`](@ref)
"""
function circle_image(
    ::Type{T} = Float32;
    radius::Integer = 20,
    calibration_value::Real = one(T),
    background::Real = 0,
    rows::Optional{Integer} = nothing,
    cols::Optional{Integer} = nothing,
    swidth::Optional{Integer} = nothing,
    sheight::Optional{Integer} = nothing,
    nϕ::Optional{Integer} = nothing,
    nphi::Optional{Integer} = nothing,
    ν::Real = 1,
    scale::Optional{Real} = nothing,
    kwargs...,
) where {T<:Real}
    nϕ = maybe(nϕ, nphi)
    ν = maybe(ν, scale)
    rows = isnothing(sheight) ? maybe(2radius, rows) : sheight
    cols = isnothing(swidth) ? maybe(rows, cols) : swidth
    nr = min(rows, cols) ÷ 2
    polar2cart(
        circle_polar_image(T; nr, nϕ, radius, calibration_value, background);
        rows, cols, background, ν, kwargs...)
end

"""
    circle_image(imp::AbstractImageParams)

Create a square image with a circle of given value from parameters `imp`.
"""
function circle_image end

function circle_image(imp::ImageParams; kwargs...)
    circle_image(
        eltype(imp);
        imp.radius,
        imp.calibration_value,
        imp.background,
        rows = 2imp.radius,
        kwargs...,
    )
end

function circle_image(imp::CircleParams; kwargs...)
    circle_image(
        eltype(imp);
        imp.radius,
        imp.calibration_value,
        imp.background,
        imp.rows,
        imp.cols,
        imp.nϕ,
        kwargs...,
    )
end

"""
    gray_scale_image(
        [T=Float32];
        rows=40,
        cols=200,
        swidth=nothing,
        sheight=nothing,
        gray_scale=-1000..1000,
    ) where {T <: Real}

Create an image with a gray scale rectangle with given scale gray_scale.

# Examples
```julia-repl
julia> gray_scale_image(swidth=80, sheight=40, gray_scale=-1000..1000)
40×80 Array{Float32,2}:
[...]
```

See also: [`circle_image`](@ref), [`combine_images`](@ref)
"""
function gray_scale_image(
    ::Type{T} = Float32;
    rows::Integer = 200,
    cols::Integer = 40,
    swidth::Optional{Integer} = nothing,
    sheight::Optional{Integer} = nothing,
    gray_scale::ClosedInterval = -1000..1000,
    background::Optional{Real} = nothing,
) where {T<:Real}
    rows = maybe(rows, sheight)
    cols = maybe(cols, swidth)
    minv, maxv = endpoints(gray_scale)
    z::T = maybe(infimum(gray_scale), background)
    val_range = minv == maxv ? minv : linspace(T, minv..maxv, cols)
    image = fill(z, rows, cols)
    image[1:rows,1:cols] .= val_range'
    image
end

"""
pyramid_gray_scale_image([T=Float32]; swidth=200, sheight=40, gray_scale=-1000..1000, plateau = 0) where {T <: Real}

Create an image with a pyramid gray scale rectangle with given scale gray_scale.

# Examples
```julia-repl
julia> pyramid_gray_scale_image(swidth=80, sheight=40, gray_scale=-1000..1000)
40×80 Array{Float32,2}:
[...]
```

See also: [`gray_scale_image`](@ref), [`circle_image`](@ref), [`combine_images`](@ref)
"""
function pyramid_gray_scale_image(
    ::Type{T} = Float32;
    swidth::Integer = 200,
    sheight::Integer = 40,
    gray_scale::ClosedInterval = -1000..1000,
    background::Optional{Real} = nothing,
    plateau::Real = 0,
) where {T<:Real}
    @assert 0 ≤ plateau ≤ 1 "Plateau is a fraction of the gray scale width!"
    minv, maxv = endpoints(gray_scale)
    val_range = fill(maxv, swidth)
    plateau_len = round(Int, swidth * plateau)
    len = (swidth - plateau_len) ÷ 2
    if len ≠ 0
        step::T = (maxv - minv) / (swidth - plateau_len - len - 1)
        half_scale = range(minv; step, length = len)
        val_range[1:len] .= half_scale
        val_range[end:-1:(swidth-len+1)] .= half_scale # Odd case included.
    end
    z::T = maybe(minv, background)
    image = fill(z, sheight, swidth)
    image .= val_range'
    image
end

"""
    gray_scale_image(imp::ImageParams)

Create an image with a gray scale rectangle with given scale gray_scale
from parameters `imp`.
"""
function gray_scale_image(imp::ImageParams{T}) where {T<:Real}
    gray_scale_image(
        T; imp.swidth, imp.sheight, imp.gray_scale, imp.background)
end

"""
    pyramid_gray_scale_image(imp::ImageParams)

Create an image with a pyramid gray scale rectangle with given scale gray_scale
from parameters `imp`.

Create an image with a pyramid gray scale rectangle with given scale gray_scale
from parameters `imp`.
"""
function pyramid_gray_scale_image(
    imp::ImageParams{T};
    plateau::Real = 0,
) where {T<:Real}
    pyramid_gray_scale_image(
        T;
        imp.swidth,
        imp.sheight,
        imp.gray_scale,
        imp.background,
        plateau,
    )
end

function combined_images_size(; swidth, sheight, radius, pad, dist)
    sheight + 2radius + 2pad + dist, swidth + 2pad
end

function combined_images_size(p::ImageParams)
    combined_images_size(;
        p.swidth,
        p.sheight,
        p.radius,
        p.pad,
        p.dist,
    )
end

"""
    combine_images(imp::ImageParams{T}, rect::AbstractMatrix{T}, circle::AbstractMatrix{T}) where {T <: Real}

Helper function to create a gray scale image.

Combine the `rect` image with the `circle` image for calibration.
The circle image should be smaller than the gray scale image!

See also: [`gray_scale_image`](@ref), [`circle_image`](@ref)
"""
function combine_images(
    imp::ImageParams{T},
    rect::AbstractMatrix{T},
    circle::AbstractMatrix{T},
) where {T<:Real}
    rect_rows, rect_cols = size(rect)
    circ_rows, circ_cols = size(circle)

    @assert((rect_rows, rect_cols) == (imp.sheight, imp.swidth),
        "Expected gray scale size $((imp.sheight, imp.swidth)), got $((srows, scols)) instead")

    @assert(2imp.radius == circ_rows,
        "Expected radius $(imp.radius), got $(circ_rows ÷ 2)")

    if circ_rows >= rect_rows || circ_cols >= rect_cols
        @warn(
            "Circle image should be smaller than gray scale rectangle",
            (circ_rows, circ_cols),
            (rect_rows, rect_cols)
        )
    end

    srows, scols = combined_images_size(imp)
    rows, cols = size(imp)
    @assert(
        rows >= srows && cols >= scols,
        "Requested image size too small: " *
        "rows($rows) ≥ $srows && columns($cols) ≥ $scols"
    )

    simage = fill(imp.background, srows, scols)
    image = fill(imp.background, rows, cols)

    crow_beg = imp.pad + 1
    crow_end = crow_beg + circ_rows - 1
    ccol_beg = scols ÷ 2 - imp.radius + 1
    ccol_end = ccol_beg + circ_cols - 1

    simage[crow_beg:crow_end, ccol_beg:ccol_end] .= circle

    rrow_beg = crow_end + imp.dist + 1
    rrow_end = rrow_beg + rect_rows - 1
    rcol_beg = imp.pad + 1
    rcol_end = rcol_beg + rect_cols - 1

    simage[rrow_beg:rrow_end, rcol_beg:rcol_end] .= rect

    ri = (rows - srows) ÷ 2 + 1
    rf = ri + srows - 1
    ci = (cols - scols) ÷ 2 + 1
    cf = ci + scols - 1
    image[ri:rf,ci:cf] .= simage
    image
end

"""
    create_image(par::ImageParams)

Create gray scale image.
"""
function create_image(par::ImageParams)
    circ_img = circle_image(par)
    rect_img = gray_scale_image(par)
    combine_images(par, rect_img, circ_img) |> CTImage
end

"""
    create_pyramid_image(par::ImageParams)

Create pyramid gray scale image.
"""
function create_pyramid_image(par::ImageParams; plateau::Real = 0)
    circ_img = circle_image(par)
    rect_img = pyramid_gray_scale_image(par; plateau)
    combine_images(par, rect_img, circ_img) |> CTImage
end

abstract type AbstractTestImage{T<:Real} <: AbstractMatrix{T} end
abstract type AbstractGrayScale{T} <: AbstractTestImage{T} end

eltype(::Type{A}) where {A <: AbstractTestImage{T}} where T = T
@inline _get_image(gs::AbstractTestImage) = getfield(gs, :image)
_atype(gs::AbstractTestImage) = _atype(_get_image(gs))
size(gs::AbstractTestImage) = size(_get_image(gs))
size(gs::AbstractTestImage, d::Integer) = size(_get_image(gs), d)
@propagate_inbounds function getindex(gs::AbstractTestImage, i::Int)
    @boundscheck checkbounds(_get_image(gs), i)
    @inbounds getindex(_get_image(gs), i)
end
@propagate_inbounds function getindex(
    gs::AbstractTestImage, I::Vararg{Int,N}
) where N
    @boundscheck checkbounds(_get_image(gs), I...)
    @inbounds getindex(_get_image(gs), I...)
end
@propagate_inbounds function setindex!(gs::AbstractTestImage, v, i::Int)
    @boundscheck checkbounds(_get_image(gs), i)
    @inbounds setindex!(_get_image(gs), v, i)
end
@propagate_inbounds function setindex!(
    gs::AbstractTestImage, v, I::Vararg{Int,N}
) where N
    @boundscheck checkbounds(_get_image(gs), I...)
    @inbounds setindex!(_get_image(gs), v, I...)
end
IndexStyle(::Type{T}) where {T<:AbstractTestImage} = IndexLinear()

struct CircleImage{T<:Real} <: AbstractTestImage{T}
    params::CircleParams{T}
    image::CTImageMat{T}
end

function getproperty(grimg::CircleImage, s::Symbol)
    s ∈ fieldnames(CircleImage) && return getfield(grimg, s)
    getproperty(grimg.params, s)
end

function CircleImage(imp::CircleParams; kwargs...)
    image = circle_image(imp; kwargs...)
    CircleImage(imp, CTImage(image))
end

function CircleImage(
    ::Type{T} = Float32; interpolation = nothing, kwargs...) where{T <: Real}
    interpolation = maybe(NoInterpolation(), interpolation)
    CircleImage(CircleParams(T; kwargs...); interpolation)
end

function CircleImage(imp::ImageParams; kwargs...)
    CircleImage(
        eltype(imp);
        imp.radius,
        imp.rows,
        imp.cols,
        imp.gray_scale,
        imp.calibration_value,
        imp.background,
        kwargs...
    )
end

function CircleImage(pbg::AbstractParallelBeamGeometry; kwargs...)
    CircleImage(eltype(pbg); pbg.width, pbg.height, kwargs...)
end

struct GrayScaleLine{T} <: AbstractGrayScale{T}
    params::ImageParams{T}
    image::CTImageMat{T}
end

function getproperty(grimg::GrayScaleLine, s::Symbol)
    s ∈ fieldnames(GrayScaleLine) && return getfield(grimg, s)
    getproperty(grimg.params, s)
end

function GrayScaleLine(imp::ImageParams; kwargs...)
    GrayScaleLine(imp, create_image(imp); kwargs...)
end

function GrayScaleLine(
    ::Type{T} = Float32;
    radius::Optional{Integer} = nothing,
    width::Optional{Integer} = nothing,
    height::Optional{Integer} = nothing,
    rows::Optional{Integer} = nothing,
    cols::Optional{Integer} = nothing,
    gray_scale::ClosedInterval = -1000..1000,
    calibration_value::Real = zero(T),
    background::Real = -1000,
    kwargs...,
) where {T<:Real}
    factor = 31 / 50 # magic number!
    #=
        If width or height are present, use them, otherwise use rows and cols.
        Finally, if rows or cols are not present, compute defaults.
    =#
    rows, cols = maybe(rows, height), maybe(cols, width)
    width = maybe(500, cols)
    height = isnothing(rows) ? round(Int, width * factor) : rows
    circ_zoom = 1
    swidth = round(Int, 17 / 25 * width)
    sheight = round(Int, 6 / 31 * height)
    pad = round(Int, 3 / 25 * width / circ_zoom)
    dist = 2pad ÷ 3
    radius = isnothing(radius) ? round(Int, width * circ_zoom / 17) : radius
    rows, cols = height, width
    imp = ImageParams(
        T;
        rows,
        cols,
        swidth,
        sheight,
        pad,
        dist,
        radius,
        gray_scale,
        calibration_value,
        background,
    )
    GrayScaleLine(imp; kwargs...)
end

function GrayScaleLine(geometry::AbstractParallelBeamGeometry; kwargs...)
    GrayScaleLine(
        eltype(geometry);
        geometry.width,
        geometry.height,
        kwargs...
    )
end

struct GrayScalePyramid{T} <: AbstractGrayScale{T}
    params::ImageParams{T}
    plateau::T
    image::CTImageMat{T}
end

function getproperty(grimg::GrayScalePyramid, s::Symbol)
    s ∈ fieldnames(GrayScalePyramid) && return getfield(grimg, s)
    getproperty(grimg.params, s)
end

function GrayScalePyramid(
    imp::ImageParams{T}; plateau::Real = zero(T), kwargs...) where {T<:Real}
    image = create_pyramid_image(imp; plateau)
    GrayScalePyramid(imp, plateau, image)
end

function GrayScalePyramid(
    ::Type{T} = Float32;
    radius::Optional{Integer} = nothing,
    width::Optional{Integer} = nothing,
    height::Optional{Integer} = nothing,
    rows::Optional{Integer} = nothing,
    cols::Optional{Integer} = nothing,
    gray_scale::ClosedInterval = -1000..1000,
    calibration_value::Real = zero(T),
    background::Real = -1000,
    plateau::Real = zero(T),
    kwargs...,
) where {T<:Real}
    factor = 31 / 50 # magic number!
    #=
        If width or height are present, use them, otherwise use rows and cols.
        Finally, if rows or cols are not present, compute defaults.
    =#
    rows, cols = maybe(rows, height), maybe(cols, width)
    width = maybe(500, cols)
    height = isnothing(rows) ? round(Int, width * factor) : rows
    circ_zoom = 1
    swidth = round(Int, 17 / 25 * width)
    sheight = round(Int, 6 / 31 * height)
    pad = round(Int, 3 / 25 * width / circ_zoom)
    dist = 2pad ÷ 3
    radius = isnothing(radius) ? round(Int, width * circ_zoom / 17) : radius
    rows, cols = height, width
    imp = ImageParams(
        T;
        rows,
        cols,
        swidth,
        sheight,
        pad,
        dist,
        radius,
        gray_scale,
        calibration_value,
        background,
    )
    GrayScalePyramid(imp; plateau, kwargs...)
end

function GrayScalePyramid(geometry::AbstractParallelBeamGeometry; kwargs...)
    GrayScalePyramid(
        eltype(geometry);
        geometry.width,
        geometry.height,
        kwargs...
    )
end

"""
    circle_polar_image([T=Float32] nr, nϕ, radius; value=1) where {T<:Real}

Create a circle of radius `radius` inside a `nr×nϕ` image
in polar coordinates .

# Examples
```julia-repl
julia> circle_polar_image(30, 30, 10)
30×30 Array{Float32,2}:
[...]
```
"""
function circle_polar_image(
    ::Type{T},
    nr::Integer,
    nϕ::Integer,
    radius::Integer;
    calibration_value::Real = one(T),
    background::Real = 0,
) where {T<:Real}
    @assert radius <= nr "Radius should be less than number of radial points"
    mat = fill(T(background), nϕ, nr)
    mat[:, 1:radius] .= calibration_value
    mat
end
circle_polar_image(::Type{T}, nr::Integer, radius::Integer; kwargs...) where {T} =
    circle_polar_image(T, nr, nr, radius; kwargs...)
circle_polar_image(nr::Integer, nϕ::Integer, radius::Integer; kwargs...) =
    circle_polar_image(Float32, nr, nϕ, radius; kwargs...)
circle_polar_image(nr::Integer, radius::Integer; kwargs...) =
    circle_polar_image(nr, nr, radius; kwargs...)
circle_polar_image(::Type{T}, radius::Integer; kwargs...) where {T} =
    circle_polar_image(T, radius, radius; kwargs...)
circle_polar_image(radius::Integer; kwargs...) =
    circle_polar_image(radius, radius; kwargs...)
function circle_polar_image(
    ::Type{T} = Float32;
    radius::Integer = 60,
    nr::Optional{Integer} = nothing,
    nϕ::Optional{Integer} = nothing,
    nphi::Optional{Integer} = nothing,
    kwargs...) where {T<:Real}
    nr = maybe(2radius, nr)
    nϕ = maybe(nr, maybe(nϕ, nphi))
    circle_polar_image(T, nr, nϕ, radius; kwargs...)
end

"""
    square_image([T=Float32] r, c; l=nothing) where {T<:Real}

Create a `l×l` square inside a `r×c` image.

# Examples
```julia-repl
julia> square_image(30, 30; l=10)
30×30 Array{Float32,2}:
[...]
```
"""
function square_image(
    ::Type{T},
    r::Integer,
    c::Integer,
    l::Optional{Integer} = nothing;
    calibration_value::Real = one(T),
    background::Real = 0,
) where {T<:Real}
    matrix = fill(T(background), r, c)
    if isnothing(l)
        l = min(r, c)
    end
    @assert l <= min(r, c)
    rin = r ÷ 2 - l ÷ 2 + 1
    rfi = rin + l - 1
    cin = c ÷ 2 - l ÷ 2 + 1
    cfi = cin + l - 1
    matrix[rin:rfi, cin:cfi] .= calibration_value
    matrix
end

square_image(r, c, l; kwargs...) = square_image(Float32, r, c, l; kwargs...)

function square_image(
    ::Type{T} = Float32;
    rows::Integer,
    cols::Integer,
    l::Optional{Integer} = nothing,
    kwargs...
) where {T <: Real}
    square_image(T, rows, cols, l; kwargs...)
end

function square_image(imp::SquareParams)
    square_image(
        eltype(imp);
        imp.l,
        imp.rows,
        imp.cols,
        imp.calibration_value,
        imp.background,
    )
end

struct SquareImage{T} <: AbstractTestImage{T}
    params::SquareParams{T}
    image::CTImageMat{T}
end

function getproperty(grimg::SquareImage, s::Symbol)
    s ∈ fieldnames(SquareImage) && return getfield(grimg, s)
    getproperty(grimg.params, s)
end

function SquareImage(imp::SquareParams)
    img = square_image(imp)
    SquareImage(imp, CTImage(img))
end

function SquareImage(
    ::Type{T} = Float32;
    l::Optional{Integer} = nothing,
    kwargs...
) where {T <: Real}
    imp = SquareParams(T; l, kwargs...)
    SquareImage(imp)
end

struct WhiteRect{T} <: AbstractGrayScale{T}
    params::ImageParams{T}
    image::CTImageMat{T}
end

function getproperty(grimg::WhiteRect, s::Symbol)
    s ∈ fieldnames(WhiteRect) && return getfield(grimg, s)
    getproperty(grimg.params, s)
end

function WhiteRect(imp::ImageParams; kwargs...)
    grimg = GrayScaleLine(imp; kwargs...)
    WhiteRect(grimg.params, grimg.image)
end

function WhiteRect(
    ::Type{T} = Float32;
    radius::Optional{Integer} = nothing,
    width::Optional{Integer} = nothing,
    height::Optional{Integer} = nothing,
    rows::Optional{Integer} = nothing,
    cols::Optional{Integer} = nothing,
    gray_scale::Union{ClosedInterval,U} = 1000,
    calibration_value::Real = zero(T),
    background::Real = -1000,
    kwargs...,
) where {T<:Real,U<:Real}
    gray_scale = gray_scale isa Real ? gray_scale..gray_scale : gray_scale
    factor = 31 / 50 # magic number!
    #=
        If width or height are present, use them, otherwise use rows and cols.
        Finally, if rows or cols are not present, compute defaults.
    =#
    rows, cols = maybe(rows, height), maybe(cols, width)
    width = maybe(500, cols)
    height = isnothing(rows) ? round(Int, width * factor) : rows
    circ_zoom = 1
    swidth = round(Int, 17 / 25 * width)
    sheight = round(Int, 6 / 31 * height)
    pad = round(Int, 3 / 25 * width / circ_zoom)
    dist = 2pad ÷ 3
    radius = isnothing(radius) ? round(Int, width * circ_zoom / 17) : radius
    rows, cols = height, width
    imp = ImageParams(
        T;
        rows,
        cols,
        swidth,
        sheight,
        pad,
        dist,
        radius,
        gray_scale,
        calibration_value,
        background,
    )
    WhiteRect(imp; kwargs...)
end

function WhiteRect(geometry::AbstractParallelBeamGeometry; kwargs...)
    WhiteRect(
        eltype(geometry);
        geometry.width,
        geometry.height,
        kwargs...
    )
end

for nm ∈ Geometry._geometry_names
    @eval begin
        $nm(imp::ImageParams; kwargs...) =
            $nm(eltype(imp); imp.width, imp.height, kwargs...)
        $nm(gs::AbstractTestImage; kwargs...) =
            $nm(eltype(gs); gs.width, gs.height, kwargs...)
    end
end

plateau_length(grsc::GrayScalePyramid) = round(Int, grsc.swidth * grsc.plateau)

for nm ∈ (:calibrate_image, :calibrate_image!, :calibrate_tomogram,:calibrate_tomogram!)
    @eval begin
        $nm(imp::AbstractImageParams; kwargs...) = x -> $nm(x, imp; kwargs...)
        $nm(x, gs::AbstractTestImage; kwargs...) = $nm(x, gs.params; kwargs...)
    end
end

"""
    calibration_position(imp::AbstractImageParams)

Get the position for calibration (max value).
"""
function calibration_position end

calibration_position(imp::Union{ImageParams,CircleParams}) = circle_position(imp)
calibration_position(imp::SquareParams) = square_position(imp)

function calibration_data(imp::AbstractImageParams)
    min_pos = background_position(imp)
    max_pos = calibration_position(imp)
    min_pos, max_pos
end

const _gs_images =
    (:CircleImage, :WhiteRect, :GrayScaleLine, :GrayScalePyramid, :SquareImage)

end # module


text:
using Base.Test
using QuantumOptics

@testset "stochastic_master" begin

b_spin = SpinBasis(1//2)
sz = sigmaz(b_spin)
sm = sigmam(b_spin)
sp = sigmap(b_spin)
zero_op = 0*sz
γ = 0.1
noise_op = 0.5γ*sz

H = γ*(sp + sm)
Hs = [noise_op]

ψ0 = spindown(b_spin)

ρ0 = dm(ψ0)
rates = [0.1]
J = [sm]
Js = [sm]
Jdagger = dagger.(J)
Js .*= rates
Jsdagger = dagger.(Js)

rates_mat = [0.1 0.05; 0.05 0.1]

dt = 1/30.0
T = [0:0.1:1;]
T_short = [0:dt:dt;]

# Function definitions for master_dynamic
function fdeterm_master(t, rho)
    H, J, Jdagger
end
function fstoch1_master(t, rho)
    [zero_op], [zero_op], rates
end
function fstoch2_master(t, rho)
    Js, Jsdagger
end
function fstoch3_master(t, rho)
    Hs
end
function fstoch4_master(t, rho)
    J, Jdagger, rates
end

# Test master
tout, ρt_det = timeevolution.master(T, ψ0, H, J; rates=rates)
tout, ρt1 = stochastic.master(T, ψ0, H, J, 0 .*J; rates=rates, dt=dt)
tout, ρt2 = stochastic.master(T, ρ0, LazyProduct(H, one(H)), sqrt.(rates).*J, 0 .* J; Hs=0 .* Hs, dt=dt)
tout, ρt3 = stochastic.master_dynamic(T, ρ0, fdeterm_master, fstoch1_master; rates=rates, dt=dt)
for i=1:length(tout)
    @test tracedistance(ρt1[i], ρt_det[i]) < dt
    @test tracedistance(ρt2[i], ρt_det[i]) < dt
    @test tracedistance(ρt3[i], ρt_det[i]) < dt
end

# Test remaining function calls for short times to test whether they work in principle
tout, ρt = stochastic.master(T_short, ρ0, H, J, 0.*J; dt=dt)
tout, ρt = stochastic.master(T_short, ψ0, H, [sm, sm], [sm, sm]; rates=rates_mat, dt=dt)
tout, ρt = stochastic.master(T_short, ρ0, H, J, J; rates_s=[0.1], Hs=Hs, dt=dt)

# Test master dynamic
tout, ρt = stochastic.master_dynamic(T_short, ψ0, fdeterm_master, fstoch1_master; dt=dt, noise_processes=1)
tout, ρt = stochastic.master_dynamic(T_short, ρ0, fdeterm_master, fstoch2_master; dt=dt)
tout, ρt = stochastic.master_dynamic(T_short, ρ0, fdeterm_master, fstoch2_master; fstoch_H=fstoch3_master, dt=dt)
tout, ρt = stochastic.master_dynamic(T_short, ρ0, fdeterm_master, fstoch2_master; fstoch_J=fstoch4_master, dt=dt)
tout, ρt = stochastic.master_dynamic(T_short, ρ0, fdeterm_master, fstoch2_master; fstoch_J=fstoch2_master, dt=dt)
tout, ρt = stochastic.master_dynamic(T_short, ρ0, fdeterm_master, fstoch2_master;
            fstoch_H=fstoch3_master, fstoch_J=fstoch4_master, dt=dt)
tout, ρt = stochastic.master_dynamic(T_short, ρ0, fdeterm_master, fstoch2_master;
            fstoch_H=fstoch3_master, fstoch_J=fstoch2_master, dt=dt)

# Test error messages
@test_throws ArgumentError stochastic.master(T, ψ0, H, [sm, sm], [sm, sm]; rates_s=rates_mat, dt=dt)
@test_throws ArgumentError stochastic.master_dynamic(T, ψ0, fdeterm_master, fstoch1_master; rates_s=rates_mat, dt=dt)

end # testset

ext: jl


text:
using Test, InnovationPrivacy, StatsBase

@testset "Check Transition Probabilities sum up to one" begin
    for N in 5:1:15, n in 0:1:10
       for i in 1:N
           @test sum(Transition_Probability(n, ParsGrowth(N = N), normalize = :no)[i, :]) ≈ 1.0
       end
    end
end

@testset "Check Bayesian Identities" begin
    for σ²_θ in 0.1:0.2:1.5, σ²_ϵ in 0.1:0.2:1.5, n in 0:1:10
        bpar = BPars(n, ParsGrowth(σ²_θ = σ²_θ, σ²_ϵ = σ²_ϵ))
        @test bpar.w_an + bpar.w_pn ≈ 1.0
    end
end

@testset "Check Expected Profit Identity" begin
    for σ²_θ in 0.1:0.5:1.6, σ²_ϵ in 0.1:0.5:1.6, n in 0:2:10, a_bar in -2:0.2:2
        par  = ParsGrowth(σ²_θ = σ²_θ, σ²_ϵ = σ²_ϵ)
        bpar = BPars(n, par)

        e_profits = Eπ(a_bar; par, bpar)
        e_revenue = Er(a_bar; par, bpar)
        wage_bill = labor_costs(a_bar; par, bpar)

        @test e_profits - (e_revenue - wage_bill - par.w*par.f) ≈ 0.
    end
end

@testset "Check Transition Matrix" begin
    for σ²_θ in 0.1:0.5:1.1, σ²_ϵ in 0.1:0.5:1.1, n in 0:5:10

        par      = ParsGrowth(σ²_θ = σ²_θ, σ²_ϵ = σ²_ϵ)
        bpar     = BPars(n, par)

        Π_test   = Transition_Matrix2(par, bpar)
        Π_direct = Transition_Probability(par, bpar)

        @test sum((Π_direct - Π_test).^2) < 1e-4
    end
end


text:
"""
Volume Image

This composite type contains volume image information

#### Fields

* `data`: contains the recorded data
* `x`, `y`, `z`, `t` Arrays containing spatial and time information
* `method` AbstractString of method used to compute tomography
* `info`: additional information in dictionary

#### `processing` Fields
The following standard names are used when saving data to the info dictionary.
* `Regularisation`: Regularisation used in tomography
* `NormalisationConstant`: Value used to normalise image to maximum of 1
* `FileName`: Name of file

"""
mutable struct VolumeImage
    data::Array{AbstractFloat,4}
    units::AbstractString
    x::Vector{typeof(1.0u"m")}
    y::Vector{typeof(1.0u"m")}
    z::Vector{typeof(1.0u"m")}
    t::Vector{typeof(1.0u"s")}
    method::AbstractString
    info::Dict
    coord_system::AbstractString

    function VolumeImage(
        data::Array{F,4},
        units::S,
        x::Vector{F},
        y::Vector{F},
        z::Vector{F},
        t::Vector{F},
        method::S,
        info::Dict,
        coord_system::S,
    ) where {F<:AbstractFloat,S<:AbstractString}

        @assert size(data, 1) == length(x)
        @assert size(data, 2) == length(y)
        @assert size(data, 3) == length(z)
        @assert size(data, 4) == length(t)

        new(data, units, x * u"m", y * u"m", z * u"m", t * u"s", method, info, coord_system)
    end

    function VolumeImage(
        data::Array{F,4},
        units::S,
        x::Vector{Met},
        y::Vector{Met},
        z::Vector{Met},
        t::Vector{Sec},
        method::S,
        info::Dict,
        coord_system::S,
    ) where {F<:AbstractFloat,S<:AbstractString,Met<:typeof(1.0u"m"),Sec<:typeof(1.0u"s")}

        @assert size(data, 1) == length(x)
        @assert size(data, 2) == length(y)
        @assert size(data, 3) == length(z)
        @assert size(data, 4) == length(t)

        new(data, units, x, y, z, t, method, info, coord_system)
    end

    function VolumeImage(
        data::Vector{F},
        units::S,
        x::Vector{F},
        y::Vector{F},
        z::Vector{F},
        t::Vector{F},
        method::S,
        info::Dict,
        coord_system::S,
    ) where {F<:AbstractFloat,S<:AbstractString}

        @assert length(x) == length(data)
        @assert length(y) == length(data)
        @assert length(z) == length(data)
        @assert length(t) == length(data)

        newX = sort(unique(x))
        newY = sort(unique(y))
        newZ = sort(unique(z))
        newT = sort(unique(t))

        L = zeros(typeof(data[1]), length(newX), length(newY), length(newZ), length(newT))

        for idx = 1:length(data)
            idxX = something(findfirst(isequal(x[idx]), newX), 0)
            idxY = something(findfirst(isequal(y[idx]), newY), 0)
            idxZ = something(findfirst(isequal(z[idx]), newZ), 0)
            idxT = something(findfirst(isequal(t[idx]), newT), 0)
            L[idxX, idxY, idxZ, idxT] = data[idx]
        end

        new(
            L,
            units,
            newX * u"m",
            newY * u"m",
            newZ * u"m",
            newT * u"s",
            method,
            info,
            coord_system,
        )
    end

end

#
# Basic operations
# ----------------
#

import Base: show

function Base.show(io::IO, vi::VolumeImage)

    println(io, "VolumeImage of method $(vi.method) and units $(vi.units)")
    println(io, "  Spanning x: $(vi.x[1]) : $(vi.x[end])")
    println(io, "  Spanning y: $(vi.y[1]) : $(vi.y[end])")
    println(io, "  Spanning z: $(vi.z[1]) : $(vi.z[end])")
    println(io, "  Spanning t: $(vi.t[1]) : $(vi.t[end])")

    if haskey(vi.info, "Regularisation")
        println(io, "  Regularisation: $(vi.info["Regularisation"])")
    end
    if haskey(vi.info, "NormalisationConstant")
        println(io, "  Image has been normalised")
    end
end


text:
const X_MS_VERSION = "2017-04-17"

# Dates.RFC1123Format with " GMT" in the end
const RFC1123_GMT = Dates.DateFormat("e, dd u yyyy HH:MM:SS G\\MT")

const UTC = TimeZones.tz"UTC"

# In some virtualized environments it is not possible to get the local timezone
const LOCALZONE = try TimeZones.localzone() catch; UTC end

max_stars_repo_name: robertdj/AzureBlob

max_issues_repo_name: robertdj/AzureBlob

max_forks_repo_name: robertdj/AzureBlob


text:
# prepare the sets N_{ij} and M_{ij} from the article
# for further details take a look at our article
function sets_m_n(airdist::Array{Float64,2}, hex::Int64)
    N = Array{Bool,3}(undef,hex,hex,hex) .= 0
    M = Array{Bool,3}(undef,hex,hex,hex) .= 0
    maxdist = Array{Float64,2}(undef,hex,hex) .= 0
    for j = 1:hex
        for v = 1:hex
            if adjacent[j,v] == 1
                for i = 1:hex
                    if airdist[i,v] < airdist[i,j]
                        N[i,j,v] = 1
                    end
                    maxdist[i,j] = max(maxdist[i,j], airdist[i,v])
                end
            end
        end
    end
    for j = 1:hex
        for v = 1:hex
            if adjacent[j,v] == 1
                for i = 1:hex
                    if airdist[i,v] < maxdist[i,j]
                        M[i,j,v] = 1
                    end
                end
            end
        end
    end
    card_n = sum(N, dims=3)[:,:,1]
    card_m = sum(M, dims=3)[:,:,1]
    return  N::Array{Bool,3},
            M::Array{Bool,3},
            card_n::Array{Int64,2},
            card_m::Array{Int64,2}
end


text:
## Quick test of envelope detection

pl = Gadfly.plot(y=leftwavdata[0], Geom.line);

pl |> PDF("/tmp/test.pdf",30cm,20cm)

# now abs and correlate with 2^16 size

leftenv = loadPs3WavEnvelopes(joinpath(datadir,"left_wavdata_corrected.jld2"))

pl = vstack(
Gadfly.plot(y=leftenv[0], Geom.line),
Gadfly.plot(y=leftenv[1], Geom.line),
Gadfly.plot(y=leftenv[2], Geom.line),
Gadfly.plot(y=leftenv[3], Geom.line),
);

pl |> PDF("/tmp/test.pdf",20cm,30cm)

##

@async run(`evince /tmp/test.pdf`)


text:
function diffeq!(du, u, p, t)
    v = Dict{Int64,Float64}()

    v[1] = p[C.V1] * p[C.a] * u[V.ppMEKc] * u[V.ERKc] / ( p[C.Km1] * (1 + u[V.pERKc] / p[C.Km2]) + u[V.ERKc] )
    v[2] = p[C.V2] * p[C.a] * u[V.ppMEKc] * u[V.pERKc] / ( p[C.Km2] * (1 + u[V.ERKc] / p[C.Km1]) + u[V.pERKc] )
    v[3] = p[C.V3] * u[V.pERKc] / ( p[C.Km3] * (1 + u[V.ppERKc] / p[C.Km4]) + u[V.pERKc] )
    v[4] = p[C.V4] * u[V.ppERKc] / ( p[C.Km4] * (1 + u[V.pERKc] / p[C.Km3]) + u[V.ppERKc] )
    v[5] = p[C.V5] * u[V.pERKn] / ( p[C.Km5] * (1 + u[V.ppERKn] / p[C.Km6]) + u[V.pERKn] )
    v[6] = p[C.V6] * u[V.ppERKn] / ( p[C.Km6] * (1 + u[V.pERKn] / p[C.Km5]) + u[V.ppERKn] )
    v[7] = p[C.KimERK] * u[V.ERKc] - p[C.KexERK] * (p[C.Vn] / p[C.Vc]) * u[V.ERKn]
    v[8] = p[C.KimpERK] * u[V.pERKc] - p[C.KexpERK] * (p[C.Vn] / p[C.Vc]) * u[V.pERKn]
    v[9] = p[C.KimppERK] * u[V.ppERKc] - p[C.KexppERK] * (p[C.Vn] / p[C.Vc]) * u[V.ppERKn]
    v[10] = p[C.V10] * u[V.ppERKn]^p[C.n10] / ( p[C.Km10]^p[C.n10] + u[V.ppERKn]^p[C.n10] )
    v[11] = p[C.p11] * u[V.PreduspmRNAn]
    v[12] = p[C.p12] * u[V.duspmRNAc]
    v[13] = p[C.p13] * u[V.duspmRNAc]
    v[14] = p[C.V14] * u[V.ppERKc] * u[V.DUSPc] / ( p[C.Km14] + u[V.DUSPc] )
    v[15] = p[C.V15] * u[V.pDUSPc] / ( p[C.Km15] + u[V.pDUSPc] )
    v[16] = p[C.p16] * u[V.DUSPc]
    v[17] = p[C.p17] * u[V.pDUSPc]
    v[18] = p[C.KimDUSP] * u[V.DUSPc] - p[C.KexDUSP] * (p[C.Vn] / p[C.Vc]) * u[V.DUSPn]
    v[19] = p[C.KimpDUSP] * u[V.pDUSPc] - p[C.KexpDUSP] * (p[C.Vn] / p[C.Vc]) * u[V.pDUSPn]
    v[20] = p[C.V20] * u[V.ppERKn] * u[V.DUSPn] / ( p[C.Km20] + u[V.DUSPn] )
    v[21] = p[C.V21] * u[V.pDUSPn] / ( p[C.Km21] + u[V.pDUSPn] )
    v[22] = p[C.p22] * u[V.DUSPn]
    v[23] = p[C.p23] * u[V.pDUSPn]
    v[24] = p[C.V24] * u[V.ppERKc] * u[V.RSKc] / ( p[C.Km24] + u[V.RSKc] )
    v[25] = p[C.V25] * u[V.pRSKc] / ( p[C.Km25] + u[V.pRSKc] )
    v[26] = p[C.KimRSK] * u[V.pRSKc] - p[C.KexRSK] * (p[C.Vn] / p[C.Vc]) * u[V.pRSKn]
    v[27] = p[C.V27] * u[V.pRSKn] * u[V.CREBn] / ( p[C.Km27] + u[V.CREBn] )
    v[28] = p[C.V28] * u[V.pCREBn] / ( p[C.Km28] + u[V.pCREBn] )
    v[29] = p[C.V29] * u[V.ppERKn] * u[V.Elk1n] / ( p[C.Km29] + u[V.Elk1n] )
    v[30] = p[C.V30] * u[V.pElk1n] / ( p[C.Km30] + u[V.pElk1n] )
    v[31] = p[C.V31] * (u[V.pCREBn] * u[V.pElk1n])^p[C.n31] / ( p[C.Km31]^p[C.n31] + (u[V.pCREBn] * u[V.pElk1n])^p[C.n31] + (u[V.Fn] / p[C.KF31])^p[C.nF31] )
    v[32] = p[C.p32] * u[V.PrecfosmRNAn]
    v[33] = p[C.p33] * u[V.cfosmRNAc]
    v[34] = p[C.p34] * u[V.cfosmRNAc]
    v[35] = p[C.V35] * u[V.ppERKc] * u[V.cFOSc] / ( p[C.Km35] + u[V.cFOSc] )
    v[36] = p[C.V36] * u[V.pRSKc] * u[V.cFOSc] / ( p[C.Km36] + u[V.cFOSc] )
    v[37] = p[C.V37] * u[V.pcFOSc] / ( p[C.Km37] + u[V.pcFOSc] )
    v[38] = p[C.p38] * u[V.cFOSc]
    v[39] = p[C.p39] * u[V.pcFOSc]
    v[40] = p[C.KimFOS] * u[V.cFOSc] - p[C.KexFOS] * (p[C.Vn] / p[C.Vc]) * u[V.cFOSn]
    v[41] = p[C.KimpcFOS] * u[V.pcFOSc] - p[C.KexpcFOS] * (p[C.Vn] / p[C.Vc]) * u[V.pcFOSn]
    v[42] = p[C.V42] * u[V.ppERKn] * u[V.cFOSn] / ( p[C.Km42] + u[V.cFOSn] )
    v[43] = p[C.V43] * u[V.pRSKn] * u[V.cFOSn] / ( p[C.Km43] + u[V.cFOSn] )
    v[44] = p[C.V44] * u[V.pcFOSn] / ( p[C.Km44] + u[V.pcFOSn] )
    v[45] = p[C.p45] * u[V.cFOSn]
    v[46] = p[C.p46] * u[V.pcFOSn]
    v[47] = p[C.p47] * u[V.DUSPn] * u[V.ppERKn] - p[C.m47] * u[V.DUSPn_ppERKn]
    v[48] = p[C.p48] * u[V.DUSPn_ppERKn]
    v[49] = p[C.p49] * u[V.DUSPn] * u[V.pERKn] - p[C.m49] * u[V.DUSPn_pERKn]
    v[50] = p[C.p50] * u[V.DUSPn_pERKn]
    v[51] = p[C.p51] * u[V.DUSPn] * u[V.ERKn] - p[C.m51] * u[V.DUSPn_ERKn]
    v[52] = p[C.p52] * u[V.pDUSPn] * u[V.ppERKn] - p[C.m52] * u[V.pDUSPn_ppERKn]
    v[53] = p[C.p53] * u[V.pDUSPn_ppERKn]
    v[54] = p[C.p54] * u[V.pDUSPn] * u[V.pERKn] - p[C.m54] * u[V.pDUSPn_pERKn]
    v[55] = p[C.p55] * u[V.pDUSPn_pERKn]
    v[56] = p[C.p56] * u[V.pDUSPn] * u[V.ERKn] - p[C.m56] * u[V.pDUSPn_ERKn]
    v[57] = p[C.V57] * u[V.pcFOSn]^p[C.n57] / ( p[C.Km57]^p[C.n57] + u[V.pcFOSn]^p[C.n57] )
    v[58] = p[C.p58] * u[V.PreFmRNAn]
    v[59] = p[C.p59] * u[V.FmRNAc]
    v[60] = p[C.p60] * u[V.FmRNAc]
    v[61] = p[C.p61] * u[V.Fc]
    v[62] = p[C.KimF] * u[V.Fc] - p[C.KexF] * (p[C.Vn] / p[C.Vc]) * u[V.Fn]
    v[63] = p[C.p63] * u[V.Fn]

    if p[C.Ligand] == p[C.EGF]  # EGF=10nM
        if t < 300.
            du[V.ppMEKc] = 0.00258
        elseif t < 600.
            du[V.ppMEKc] = -0.00111
        elseif t < 900.
            du[V.ppMEKc] = -0.000625
        elseif t < 1200.
            du[V.ppMEKc] = -0.000135
        elseif t < 1800.
            du[V.ppMEKc] = -0.000135
        elseif t < 2700.
            du[V.ppMEKc] = -0.0000480
        elseif t < 3600.
            du[V.ppMEKc] = -0.00000852
        elseif t <= 5400.
            du[V.ppMEKc] = -0.00000728
        end
    elseif p[C.Ligand] == p[C.HRG]  # HRG=10nM
        if t < 300.
            du[V.ppMEKc] = 0.00288
        elseif t < 600.
            du[V.ppMEKc] = 0.000451
        elseif t < 900.
            du[V.ppMEKc] = -0.000545
        elseif t < 1200.
            du[V.ppMEKc] = 0.0000522
        elseif t < 1800.
            du[V.ppMEKc] = 0.0000522
        elseif t < 2700.
            du[V.ppMEKc] = 0.0000399
        elseif t < 3600.
            du[V.ppMEKc] = -0.0000500
        elseif t <= 5400.
            du[V.ppMEKc] = -0.0000478
        end
    else
        du[V.ppMEKc] = 0.0
    end

    du[V.CREBn] = -v[27] + v[28]
    du[V.pCREBn] = v[27] - v[28]
    du[V.ERKc] = -v[1] + v[3] - v[7]
    du[V.ERKn] = v[5] + v[7] * (p[C.Vc] / p[C.Vn]) + v[50] - v[51] + v[55] - v[56]
    du[V.pERKc] = v[1] - v[2] - v[3] + v[4] - v[8]
    du[V.pERKn] = -v[5] + v[6] + v[8] * (p[C.Vc] / p[C.Vn]) + v[48] - v[49] + v[53] - v[54]
    du[V.ppERKc] =  v[2] - v[4] - v[9]
    du[V.ppERKn] = -v[6] + v[9] * (p[C.Vc] / p[C.Vn]) - v[47] - v[52]
    du[V.Elk1n] = -v[29] + v[30]
    du[V.pElk1n] = v[29] - v[30]
    du[V.cFOSc] = v[34] - v[35] - v[36] + v[37] - v[38] - v[40]
    du[V.cFOSn] = v[40] * (p[C.Vc] / p[C.Vn]) - v[42] - v[43] + v[44] - v[45]
    du[V.pcFOSc] = v[35] + v[36] - v[37] - v[39] - v[41]
    du[V.pcFOSn] = v[41] * (p[C.Vc] / p[C.Vn]) + v[42] + v[43] - v[44] - v[46]
    du[V.DUSPc] = v[13] - v[14] + v[15] - v[16] - v[18]
    du[V.DUSPn] = v[18] * (p[C.Vc] / p[C.Vn]) - v[20] + v[21] - v[22] - v[47] + v[48] - v[49] + v[50] - v[51]
    du[V.pDUSPc] = v[14] - v[15] - v[17] - v[19]
    du[V.pDUSPn] = v[19] * (p[C.Vc] / p[C.Vn]) + v[20] - v[21] - v[23] - v[52] + v[53] - v[54] + v[55] - v[56]
    du[V.DUSPn_ERKn] = v[51]
    du[V.DUSPn_pERKn] = v[49] - v[50]
    du[V.DUSPn_ppERKn] = v[47] - v[48]
    du[V.pDUSPn_ERKn] = v[56]
    du[V.pDUSPn_pERKn] = v[54] - v[55]
    du[V.pDUSPn_ppERKn] = v[52] - v[53]
    du[V.RSKc] = -v[24] + v[25]
    du[V.pRSKc] = v[24] - v[25] - v[26]
    du[V.pRSKn] = v[26] * (p[C.Vc] / p[C.Vn])
    du[V.PrecfosmRNAn] = v[31] - v[32]
    du[V.PreduspmRNAn] = v[10] - v[11]
    du[V.cfosmRNAc] = v[32] * (p[C.Vn] / p[C.Vc]) - v[33]
    du[V.duspmRNAc] = v[11] * (p[C.Vn] / p[C.Vc]) - v[12]
    du[V.Fc] = v[60] - v[61] - v[62]
    du[V.Fn] = v[62] * (p[C.Vc] / p[C.Vn]) - v[63]
    du[V.FmRNAc] = v[58] * (p[C.Vn] / p[C.Vc]) - v[59]
    du[V.PreFmRNAn] = v[57] - v[58]

end

function param_values()::Vector{Float64}
    p::Vector{Float64} = zeros(C.NUM)

    p[C.V1] = 0.34284837
    p[C.Km1] = 307.0415253
    p[C.V2] = 2.20e-01
    p[C.Km2] = 3.50e+02
    p[C.V3] = 7.20e-01
    p[C.Km3] = 1.60e+02
    p[C.V4] = 6.48e-01
    p[C.Km4] = 6.00e+01
    p[C.V5] = 19.49872346
    p[C.Km5] = 29.94073716
    p[C.V6] = p[C.V5]
    p[C.Km6] = p[C.Km5]
    p[C.KimERK] = 1.20e-02
    p[C.KexERK] = 1.80e-02
    p[C.KimpERK] = 1.20e-02
    p[C.KexpERK] = 1.80e-02
    p[C.KimppERK] = 1.10e-02
    p[C.KexppERK] = 1.30e-02
    p[C.V10] = 29.24109258
    p[C.Km10] = 169.0473748
    p[C.n10] = 3.970849295
    p[C.p11] = 0.000126129
    p[C.p12] = 0.007875765
    p[C.p13] = 0.001245747
    p[C.V14] = 5.636949216
    p[C.Km14] = 34180.48
    p[C.V15] = 2.992346912
    p[C.Km15] = 0.001172165
    p[C.p16] = 2.57e-04
    p[C.p17] = 9.63e-05
    p[C.KimDUSP] = 0.024269764
    p[C.KexDUSP] = 0.070467899
    p[C.KimpDUSP] = p[C.KimDUSP]
    p[C.KexpDUSP] = p[C.KexDUSP]
    p[C.V20] = 0.157678678
    p[C.Km20] = 735598.6967
    p[C.V21] = 0.005648117
    p[C.Km21] = 387.8377182
    p[C.p22] = 2.57e-04
    p[C.p23] = 9.63e-05
    p[C.V24] = 0.550346114
    p[C.Km24] = 29516.06587
    p[C.V25] = 10.09063736
    p[C.Km25] = 0.913939859
    p[C.KimRSK] = 0.025925065
    p[C.KexRSK] = 0.129803956
    p[C.V27] = 19.23118154
    p[C.Km27] = 441.5834425
    p[C.V28] = 6.574759504
    p[C.Km28] = 14.99180922
    p[C.V29] = 0.518529841
    p[C.Km29] = 21312.69109
    p[C.V30] = 13.79479021
    p[C.Km30] = 15.04396629
    p[C.V31] = 0.655214248
    p[C.Km31] = 185.9760682
    p[C.n31] = 1.988003164
    p[C.KF31] = 0.013844393
    p[C.nF31] = 2.800340453
    p[C.p32] = 0.003284434
    p[C.p33] = 0.000601234
    p[C.p34] = 7.65E-05
    p[C.V35] = 8.907637012
    p[C.Km35] = 8562.744184
    p[C.V36] = 0.000597315
    p[C.Km36] = 528.552427
    p[C.V37] = 1.745848179
    p[C.Km37] = 0.070379236
    p[C.p38] = 2.57e-04
    p[C.p39] = 9.63e-05
    p[C.KimFOS] = 0.54528521
    p[C.KexFOS] = 0.133249762
    p[C.KimpcFOS] = p[C.KimFOS]
    p[C.KexpcFOS] = p[C.KexFOS]
    p[C.V42] = 0.909968714
    p[C.Km42] = 3992.061328
    p[C.V43] = 0.076717457
    p[C.Km43] = 1157.116021
    p[C.V44] = 0.078344305
    p[C.Km44] = 0.051168202
    p[C.p45] = 2.57e-04
    p[C.p46] = 9.63e-05
    p[C.p47] = 0.001670815
    p[C.m47] = 15.80783969
    p[C.p48] = 0.686020478
    p[C.p49] = 0.314470502
    p[C.m49] = 2.335459127
    p[C.p50] = 26.59483436
    p[C.p51] = 0.01646825
    p[C.m51] = 9.544308421
    p[C.p52] = p[C.p47]
    p[C.m52] = p[C.m47]
    p[C.p53] = p[C.p48]
    p[C.p54] = p[C.p49]
    p[C.m54] = p[C.m49]
    p[C.p55] = p[C.p50]
    p[C.p56] = p[C.p51]
    p[C.m56] = p[C.m51]
    p[C.V57] = 1.026834758
    p[C.Km57] = 0.637490056
    p[C.n57] = 3.584464176
    p[C.p58] = 0.000270488
    p[C.p59] = 0.001443889
    p[C.p60] = 0.002448164
    p[C.p61] = 3.50E-05
    p[C.KimF] = 0.019898797
    p[C.KexF] = 0.396950616
    p[C.p63] = 4.13E-05

    p[C.a] = 218.6276381
    p[C.Vn] = 0.22
    p[C.Vc] = 0.94

    p[C.EGF] = 0.
    p[C.HRG] = 1.
    p[C.no_ligand] = 2.

    return p
end

function initial_values()::Vector{Float64}
    u0::Vector{Float64} = zeros(V.NUM)

    u0[V.ERKc] = 9.60e+02
    u0[V.RSKc] = 3.53e+02
    u0[V.CREBn] = 1.00e+03
    u0[V.Elk1n] = 1.51e+03

    return u0
end


text:
##

##
"""
    function tux()

TuX
(c) pborys@p-soft.silesia.linux.org.pl

# Example
```jldoctest
julia> cowsay("Talk is cheap. Show me the code.", cow=Cowsay.tux)
 __________________________________
< Talk is cheap. Show me the code. >
 ----------------------------------
   \\
    \\
        .--.
       |o_o |
       |:_/ |
      //   \\ \\
     (|     | )
    /'\\_   _/`\\
    \\___)=(___/
```
"""
function tux(;eyes="oo", tongue="  ", thoughts="\\")
the_cow =
"""
   $thoughts
    $thoughts
        .--.
       |o_o |
       |:_/ |
      //   \\ \\
     (|     | )
    /'\\_   _/`\\
    \\___)=(___/

"""
return the_cow
end


text:
# This file is a part of LegendUpROOTIO.jl, licensed under the MIT License (MIT).

@enum EWFEncScheme EWFEncScheme_kUndefined=0 EWFEncScheme_kDouble=1 EWFEncScheme_kShort=2 EWFEncScheme_kDiffVarInt=3 EWFEncScheme_kUShort=4

@enum EWFType EWFType_kNoType=0 EWFType_kCharge=1 EWFType_kCurrent=2 EWFType_kADC=3 EWFType_kVoltage=4

@enum MGEventType MGEventType_kReal=0 MGEventType_kPulser=1 MGEventType_kMC=2 MGEventType_kUndefined=3 MGEventType_kBaseline=4

@enum EWaveformTag EWaveformTag_kNormal=0 EWaveformTag_kUnderflow=1 EWaveformTag_kOverflow=2 EWaveformTag_kUndefined=3 EWaveformTag_kCorrupted=4

struct MGTDataObjectContent{S<:AbstractString}
    tnamed::TNamedContent{S}
end

function Base.read(io::ROOTIOBuffer, ::Type{MGTDataObjectContent})
    ver = read(io, ROOTClassVersion)
    @assert ver.version == 2
    tnamed = read(io, TNamedContent)
    MGTDataObjectContent(tnamed)
end

function read_mgtclonesarray_header(io::ROOTIOBuffer)
    unknown_ver = read(io, ROOTClassVersion)
    if unknown_ver.version == -1
        # Interpretation of unknown_ver.version == -1 unclear
        @assert String(read(io.buffer, length("\xff\xffMGTClonesArray\0"))) == "\xff\xffMGTClonesArray\0"
    elseif (unknown_ver.version & 0xff00) == 0x8000
        # Interpretation of unknown_ver.version == 0x80.. unclear
        unknown_short = read(io, UInt16)
    else
        throw(ErrorException("Unexpected class version $(unknown_ver.version)"))
    end
    MGTClonesArray_wf_ver = read(io, ROOTClassVersion)
    @assert MGTClonesArray_wf_ver.version == 2
end

function _read_wfsamples_noenc(io::ROOTIOBuffer, T::Type{<:Real}, n::Int)
    samples = Vector{T}(undef, n)
    read!(io.buffer, samples)
    samples .= ntoh.(samples)
    samples
end

@inline function _read_vl_int(io::IO, T::Type{<:Unsigned})
    maxPos = 8 * sizeof(T)
    x::T = 0
    pos::Int = 0
    while true
        (pos >= maxPos) && throw(ErrorException("Overflow during decoding of variable-length encoded number."))
        b = read(io, UInt8)
        x = x | (T(b & 0x7f) << pos)
        if ((b & 0x80) == 0)
            return x
        else
            pos += 7
        end
    end
end

@inline _read_vl_int(io::IO, T::Type{<:Signed}) = _zigzagdec(_read_vl_int(io, unsigned(T)))

@inline _zigzagdec(x::Unsigned) = signed(xor((x >>> 1), (-(x & 1))))

function _read_wfsamples_diffvarint(io::ROOTIOBuffer, T::Type{<:Real}, n::Int)
    T = Int32
    buf = io.buffer
    x::T = 0
    samples = Vector{T}(undef, n)
    for i in eachindex(samples)
        x = _read_vl_int(buf, T) + x
        samples[i] = x
    end
    samples
end

function _read_wfsamples(io::ROOTIOBuffer, enc::EWFEncScheme, n::Integer)
    if enc == EWFEncScheme_kUndefined || enc == EWFEncScheme_kDouble
        _read_wfsamples_noenc(io, Float64, n)
    elseif (enc == EWFEncScheme_kShort)
        _read_wfsamples_noenc(io, Int16, n)
    elseif (enc == EWFEncScheme_kUShort)
        _read_wfsamples_noenc(io, UInt16, n)
    elseif (enc == EWFEncScheme_kDiffVarInt)
        _read_wfsamples_diffvarint(io, Int32, n)
    else
        throw(ErrorException("Sample encoding scheme $enc not supported"))
    end
end

function read_MGTWaveform(io::ROOTIOBuffer)
    MGTWaveform_ver = read(io, ROOTClassVersion)
    @assert MGTWaveform_ver.version == 5
    mgtdataobject = read(io, MGTDataObjectContent)

    fSampFreq = read(io, Float64)u"GHz"
    fTOffset = read(io, Float64)u"ns"
    fWFType = EWFType(read(io, Int32))
    fData_size = Int(read(io, UInt64))
    fWFEncScheme = EWFEncScheme(read(io, Int32))
    fData = _read_wfsamples(io, fWFEncScheme, fData_size)
    fID = read(io, Int32)

    unknown_byte = read(io, UInt8)
    @assert unknown_byte == 0x01

    time_units = u"ns"

    dt = Int32(uconvert(time_units, inv(fSampFreq)) / time_units) * time_units
    t0 = Int32(uconvert(time_units, fTOffset) / time_units) * time_units
    t_idxs = Int32(0):Int32(size(fData,1) - 1)
    t = t0 .+ t_idxs .* dt

    (
        ch = fID,
        wf = RDWaveform(t, fData),
    )
end

function read_tca_MGTWaveform(io::ROOTIOBuffer)
    read_mgtclonesarray_header(io)
    tcahdr = read(io, TClonesArrayHeader)
    @assert tcahdr.fName == "MGTWaveforms"
    @assert tcahdr.fClass == "MGTWaveform;5"

    wfdata_1 = read_MGTWaveform(io)
    ch = [wfdata_1.ch]
    wf_t = StructArray([wfdata_1.wf.time])
    wf_v = VectorOfVectors([wfdata_1.wf.value])

    for i in 2:tcahdr.nobjects
        wfdata = read_MGTWaveform(io)
        push!(ch, wfdata.ch)
        push!(wf_t, wfdata.wf.time)
        push!(wf_v, wfdata.wf.value)
    end

    TypedTables.Table(
        ch = ch,
        wf = StructArray{RDWaveform}((wf_t, wf_v))
    )
end

#=
# Possible digitizer data fields, depending on class:

bool fADCNMinusOneTriggerFlag;
bool fADCNPlusOneTriggerFlag;
bool fCFDCrossingOccurred;
bool fExternalTriggerFlag;
bool fIsInverted;
bool fIsMuVetoed;
bool fIsWrapModeEnabled;
bool fLEDCrossingSignIsNegative;
bool fPileupFlag;
bool fRetriggerFlag;
bool fTriggerFlag;
bool fTTCLTimeoutFlag;
double fClockFrequency;
double fEnergy;
int32_t fCFDPoint1;
int32_t fCFDPoint2;
int fEventNumber;
MGWaveformTag::EWaveformTag fWaveformTag;
size_t fPretrigger;
size_t fTriggerNumber;
uint32_t fNMaxChannels;
uint32_t fBitResolution;
uint32_t fBoardID;
uint32_t fEnergyInitial;
uint32_t fEnergyMax;
uint32_t fFastTriggerCounter;
uint32_t fID;
uint32_t fNChannels;
uint32_t fNofWrapSamples;
uint32_t fWrapStartIndex;
uint64_t fCFDTimeStamp;
uint64_t fIndex;
uint64_t fTimeStamp;
unsigned int fDecimalTimeStamp;
unsigned int fMuVetoSample;
unsigned short fBoardSerialNumber;
=#

# Use fnv1a32 hash of names for enum int values:

function read_MGTVDigitizerData(io::ROOTIOBuffer)
    MGTVDigitizerData_ver = read(io, ROOTClassVersion)
    @assert MGTVDigitizerData_ver.version == 1

    mgtdataobject = read(io, MGTDataObjectContent)

    fEnergy = read(io, Float64)
    fTimeStamp = read(io, UInt64)
    fID = read(io, Int32)
    fIndex = read(io, UInt64)

    (
        fEnergy = fEnergy,
        fTimeStamp = fTimeStamp,
        fID = fID,
        fIndex = fIndex,
    )
end

function read_GETVGERDADigitizerData(io::ROOTIOBuffer)
    GETVGERDADigitizerData_ver = read(io, ROOTClassVersion)
    @assert GETVGERDADigitizerData_ver.version == 2

    parent_content = read_MGTVDigitizerData(io)

    fNMaxChannels = UInt32(4) # meaningful?

    parent_content
end

function read_GETGERDADigitizerData(io::ROOTIOBuffer)
    GETGERDADigitizerData_ver = read(io, ROOTClassVersion)
    @assert GETGERDADigitizerData_ver.version == 1

    parent_content = read_GETVGERDADigitizerData(io)

    fClockFrequency = read(io, Float64)
    fEventNumber = read(io, Int32)
    fPretrigger = read(io, UInt64)
    fTriggerNumber = read(io, UInt64)
    fDecimalTimeStamp = read(io, UInt32)
    fIsMuVetoed = read(io, Bool)
    fMuVetoSample = read(io, UInt32)
    fIsInverted = read(io, Bool)
    fWaveformTag = EWaveformTag(read(io, Int32))

    unknown_byte = read(io, UInt8)
    @assert unknown_byte == 0x01

    content = (
        type = Val(:GETGERDADigitizerData),

        fClockFrequency = fClockFrequency,
        fEventNumber = fEventNumber,
        fPretrigger = fPretrigger,
        fTriggerNumber = fTriggerNumber,
        fDecimalTimeStamp = fDecimalTimeStamp,
        fIsMuVetoed = fIsMuVetoed,
        fMuVetoSample = fMuVetoSample,
        fIsInverted = fIsInverted,
        fWaveformTag = fWaveformTag,
    )
    merge(parent_content, content)
end

function read_tca_MGTVDigitizerData(io::ROOTIOBuffer)
    read_mgtclonesarray_header(io)
    tcahdr = read(io, TClonesArrayHeader)

    # ToDo: Support more digitizer classes

    @assert tcahdr.fName == "GETGERDADigitizerDatas"
    @assert tcahdr.fClass == "GETGERDADigitizerData;1"

    TypedTables.Table([read_GETGERDADigitizerData(io) for i in 1:tcahdr.nobjects])
end

function read_mgtevent(io::ROOTIOBuffer)
    MGTEvent_header = read(io, ROOTClassHeader)
    @assert MGTEvent_header.classname == "MGTEvent"

    version_MGTEvent = read(io, ROOTClassVersion)
    @assert version_MGTEvent.version == 9

    mgtdataobject = read(io, MGTDataObjectContent)

    fMGEventType = read(io, Int32)
    fETotal = read(io, Float64)
    fTime = read(io, Float64)

    fWaveforms = read_tca_MGTWaveform(io)

    fDigitizerData = read_tca_MGTVDigitizerData(io)

    fActiveID = let
        unknown_ver = read(io, ROOTClassVersion)
        @assert unknown_ver.version == -1
        @assert String(read(io.buffer, length("\xff\xffvector<bool>\0"))) == "\xff\xffvector<bool>\0"
        read(io, Vector{Bool})
    end

    fUseAuxWaveformArray = read(io, Bool)

    fAuxWaveforms = read_tca_MGTWaveform(io)

    fEventNumber = read(io, Int32)

    (
        fMGEventType = MGEventType(fMGEventType),
        fETotal = fETotal,
        fTime = fTime,
        fWaveforms = fWaveforms,
        fDigitizerData = fDigitizerData,
        fActiveID = fActiveID,
        fUseAuxWaveformArray = fUseAuxWaveformArray,
        fAuxWaveforms = fAuxWaveforms,
        fEventNumber = fEventNumber,
    )
end

function raw2mgtevent(evt_raw::UpROOT.OpaqueObject{:MGTEvent})
    io = ROOTIOBuffer(evt_raw.data)
    event = read_mgtevent(io)
    @assert eof(io)
    event
end

const _dict_EWaveformTag2daqflags = IdDict(
    EWaveformTag_kNormal => daq_noflags,
    EWaveformTag_kUnderflow => daq_underflow,
    EWaveformTag_kOverflow => daq_overflow,
    # EWaveformTag_kUndefined => ??,
    EWaveformTag_kCorrupted => daq_corrupt,
)

mgdo2legend(x::EWaveformTag) = _dict_EWaveformTag2daqflags[x]

const _dict_MGEventType2EventType = IdDict(
    MGEventType_kReal => evt_real,
    MGEventType_kPulser => evt_pulser,
    MGEventType_kMC => evt_mc,
    MGEventType_kUndefined => evt_undef,
    MGEventType_kBaseline => evt_baseline,
)

mgdo2legend(x::MGEventType) = _dict_MGEventType2EventType[x]

function mgdo2legend(::Val{:GETGERDADigitizerData}, x::TypedTables.Table)
    TypedTables.Table(
        ch = x.fID,

        daqclk = x.fTimeStamp,
        index = x.fIndex, # Meaning?

        psa_energy = x.fEnergy,

        # ignore = x.fClockFrequency, # better kept in metadata/config files
        daqevtno = x.fEventNumber, # Meaning?
        # ignore = x.fPretrigger, # use in waveform
        trigno = x.fTriggerNumber, # Meaning?
        muveto = BitArray(x.fIsMuVetoed),
        muveto_sample = x.fMuVetoSample, # Meaning?
        decimal_timestamp = x.fDecimalTimeStamp, # Meaning?
        inverted = BitArray(x.fIsInverted), #
        #flags = mgdo2legend.(x.fWaveformTag),
    )
end

function mgdo2legend(evt::NamedTuple)
    @argcheck evt.fWaveforms.ch == evt.fDigitizerData.fID
    @argcheck isempty(evt.fAuxWaveforms) || evt.fWaveforms.ch == evt.fAuxWaveforms.ch
    @argcheck evt.fETotal ≈ sum(evt.fDigitizerData.fEnergy)
    # EWFType_kADC

    digidata = mgdo2legend(eltype(evt.fDigitizerData.type)(), evt.fDigitizerData)

    n = length(eachindex(digidata))

    evtdata = (
        evttype = fill(mgdo2legend(evt.fMGEventType), n),
        unixtime = fill(evt.fTime, n),
        evtno = fill(evt.fEventNumber, n), # LEGEND daq files will typically not have this
        # evt.fActiveID, # better kept in metadata/config files
        # evt.fETotal # Don't need to store this
        # evt.fUseAuxWaveformArray # Don't need to store this
    )

    wfdata = (
        waveform_lf = evt.fWaveforms.wf,
        waveform_hf = evt.fAuxWaveforms.wf,
    )

    TypedTables.Table(merge(evtdata, Tables.columns(digidata), wfdata))
end

function mgdo2legend(events::AbstractVector{<:NamedTuple})
    tbl = mgdo2legend(deepcopy(events[firstindex(events)]))
    for i in (firstindex(events) + 1):lastindex(events)
        append!(tbl, mgdo2legend(events[i]))
    end
    tbl
end


text:
int_rules_4_3_9 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 4.3.9*trig^m*(a+b*tan^n+c*tan^(2*n))^p =#
    @apply_utils Antiderivative((~a + ~(b') * tan(~(d') + ~(e') * ~x) ^ ~(n') + ~(c') * tan(~(d') + ~(e') * ~x) ^ ~(n2')) ^ ~(p'), ~x) => (1 / (4 ^ ~p * (~c) ^ ~p)) * Antiderivative((~b + 2 * ~c * tan(~d + ~e * ~x) ^ ~n) ^ (2 * ~p), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n], ~x) && (EqQ(~n2, 2 * ~n) && (EqQ((~b) ^ 2 - 4 * ~a * ~c, 0) && IntegerQ(~p)))
    @apply_utils Antiderivative((~a + ~(b') * cot(~(d') + ~(e') * ~x) ^ ~(n') + ~(c') * cot(~(d') + ~(e') * ~x) ^ ~(n2')) ^ ~(p'), ~x) => (1 / (4 ^ ~p * (~c) ^ ~p)) * Antiderivative((~b + 2 * ~c * cot(~d + ~e * ~x) ^ ~n) ^ (2 * ~p), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n], ~x) && (EqQ(~n2, 2 * ~n) && (EqQ((~b) ^ 2 - 4 * ~a * ~c, 0) && IntegerQ(~p)))
    @apply_utils Antiderivative((~a + ~(b') * tan(~(d') + ~(e') * ~x) ^ ~(n') + ~(c') * tan(~(d') + ~(e') * ~x) ^ ~(n2')) ^ ~p, ~x) => ((~a + ~b * tan(~d + ~e * ~x) ^ ~n + ~c * tan(~d + ~e * ~x) ^ (2 * ~n)) ^ ~p / (~b + 2 * ~c * tan(~d + ~e * ~x) ^ ~n) ^ (2 * ~p)) * Antiderivative((~b + 2 * ~c * tan(~d + ~e * ~x) ^ ~n) ^ (2 * ~p), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n], ~x) && (EqQ(~n2, 2 * ~n) && (EqQ((~b) ^ 2 - 4 * ~a * ~c, 0) && Not(IntegerQ(~p))))
    @apply_utils Antiderivative((~a + ~(b') * cot(~(d') + ~(e') * ~x) ^ ~(n') + ~(c') * cot(~(d') + ~(e') * ~x) ^ ~(n2')) ^ ~p, ~x) => ((~a + ~b * cot(~d + ~e * ~x) ^ ~n + ~c * cot(~d + ~e * ~x) ^ (2 * ~n)) ^ ~p / (~b + 2 * ~c * cot(~d + ~e * ~x) ^ ~n) ^ (2 * ~p)) * Antiderivative((~b + 2 * ~c * cot(~d + ~e * ~x) ^ ~n) ^ (2 * ~p), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n], ~x) && (EqQ(~n2, 2 * ~n) && (EqQ((~b) ^ 2 - 4 * ~a * ~c, 0) && Not(IntegerQ(~p))))
    @apply_utils Antiderivative(1 / (~(a') + ~(b') * tan(~(d') + ~(e') * ~x) ^ ~(n') + ~(c') * tan(~(d') + ~(e') * ~x) ^ ~(n2')), ~x) => Module([q = Rt((~b) ^ 2 - 4 * ~a * ~c, 2)], ((2 * ~c) / q) * Antiderivative(1 / ((~b - q) + 2 * ~c * tan(~d + ~e * ~x) ^ ~n), ~x) - ((2 * ~c) / q) * Antiderivative(1 / (~b + q + 2 * ~c * tan(~d + ~e * ~x) ^ ~n), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n], ~x) && (EqQ(~n2, 2 * ~n) && NeQ((~b) ^ 2 - 4 * ~a * ~c, 0))
    @apply_utils Antiderivative(1 / (~(a') + ~(b') * cot(~(d') + ~(e') * ~x) ^ ~(n') + ~(c') * cot(~(d') + ~(e') * ~x) ^ ~(n2')), ~x) => Module([q = Rt((~b) ^ 2 - 4 * ~a * ~c, 2)], ((2 * ~c) / q) * Antiderivative(1 / ((~b - q) + 2 * ~c * cot(~d + ~e * ~x) ^ ~n), ~x) - ((2 * ~c) / q) * Antiderivative(1 / (~b + q + 2 * ~c * cot(~d + ~e * ~x) ^ ~n), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n], ~x) && (EqQ(~n2, 2 * ~n) && NeQ((~b) ^ 2 - 4 * ~a * ~c, 0))
    @apply_utils Antiderivative(sin(~(d') + ~(e') * ~x) ^ ~m * (~(a') + ~(b') * (~(f') * tan(~(d') + ~(e') * ~x)) ^ ~(n') + ~(c') * (~(f') * tan(~(d') + ~(e') * ~x)) ^ ~(n2')) ^ ~p, ~x) => (~f / ~e) * Subst(Antiderivative(((~x) ^ ~m * (~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ ~p) / ((~f) ^ 2 + (~x) ^ 2) ^ (~m / 2 + 1), ~x), ~x, ~f * tan(~d + ~e * ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~n, ~p], ~x) && (EqQ(~n2, 2 * ~n) && IntegerQ(~m / 2))
    @apply_utils Antiderivative(cos(~(d') + ~(e') * ~x) ^ ~m * (~(a') + ~(b') * (~(f') * cot(~(d') + ~(e') * ~x)) ^ ~(n') + ~(c') * (~(f') * cot(~(d') + ~(e') * ~x)) ^ ~(n2')) ^ ~p, ~x) => (-(~f) / ~e) * Subst(Antiderivative(((~x) ^ ~m * (~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ ~p) / ((~f) ^ 2 + (~x) ^ 2) ^ (~m / 2 + 1), ~x), ~x, ~f * cot(~d + ~e * ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~n, ~p], ~x) && (EqQ(~n2, 2 * ~n) && IntegerQ(~m / 2))
    @apply_utils Antiderivative(sin(~(d') + ~(e') * ~x) ^ ~(m') * (~(a') + ~(b') * tan(~(d') + ~(e') * ~x) ^ ~(n') + ~(c') * tan(~(d') + ~(e') * ~x) ^ ~(n2')) ^ ~(p'), ~x) => Module([g = FreeFactors(cos(~d + ~e * ~x), ~x)], (-g / ~e) * Subst(Antiderivative(((1 - g ^ 2 * (~x) ^ 2) ^ ((~m - 1) / 2) * ExpandToSum(~a * (g * ~x) ^ (2 * ~n) + ~b * (g * ~x) ^ ~n * (1 - g ^ 2 * (~x) ^ 2) ^ (~n / 2) + ~c * (1 - g ^ 2 * (~x) ^ 2) ^ ~n, ~x) ^ ~p) / (g * ~x) ^ (2 * ~n * ~p), ~x), ~x, cos(~d + ~e * ~x) / g)) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (EqQ(~n2, 2 * ~n) && (IntegerQ((~m - 1) / 2) && (IntegerQ(~n / 2) && IntegerQ(~p))))
    @apply_utils Antiderivative(cos(~(d') + ~(e') * ~x) ^ ~(m') * (~(a') + ~(b') * cot(~(d') + ~(e') * ~x) ^ ~(n') + ~(c') * tan(~(d') + ~(e') * ~x) ^ ~(n2')) ^ ~(p'), ~x) => Module([g = FreeFactors(sin(~d + ~e * ~x), ~x)], (g / ~e) * Subst(Antiderivative(((1 - g ^ 2 * (~x) ^ 2) ^ ((~m - 1) / 2) * ExpandToSum(~a * (g * ~x) ^ (2 * ~n) + ~b * (g * ~x) ^ ~n * (1 - g ^ 2 * (~x) ^ 2) ^ (~n / 2) + ~c * (1 - g ^ 2 * (~x) ^ 2) ^ ~n, ~x) ^ ~p) / (g * ~x) ^ (2 * ~n * ~p), ~x), ~x, sin(~d + ~e * ~x) / g)) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (EqQ(~n2, 2 * ~n) && (IntegerQ((~m - 1) / 2) && (IntegerQ(~n / 2) && IntegerQ(~p))))
    @apply_utils Antiderivative(cos(~(d') + ~(e') * ~x) ^ ~m * (~(a') + ~(b') * (~(f') * tan(~(d') + ~(e') * ~x)) ^ ~(n') + ~(c') * (~(f') * tan(~(d') + ~(e') * ~x)) ^ ~(n2')) ^ ~(p'), ~x) => ((~f) ^ (~m + 1) / ~e) * Subst(Antiderivative((~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ ~p / ((~f) ^ 2 + (~x) ^ 2) ^ (~m / 2 + 1), ~x), ~x, ~f * tan(~d + ~e * ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~n, ~p], ~x) && (EqQ(~n2, 2 * ~n) && IntegerQ(~m / 2))
    @apply_utils Antiderivative(sin(~(d') + ~(e') * ~x) ^ ~m * (~(a') + ~(b') * (~(f') * cot(~(d') + ~(e') * ~x)) ^ ~(n') + ~(c') * (~(f') * cot(~(d') + ~(e') * ~x)) ^ ~(n2')) ^ ~(p'), ~x) => (-((~f) ^ (~m + 1)) / ~e) * Subst(Antiderivative((~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ ~p / ((~f) ^ 2 + (~x) ^ 2) ^ (~m / 2 + 1), ~x), ~x, ~f * cot(~d + ~e * ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~n, ~p], ~x) && (EqQ(~n2, 2 * ~n) && IntegerQ(~m / 2))
    @apply_utils Antiderivative(cos(~(d') + ~(e') * ~x) ^ ~m * (~(a') + ~(b') * tan(~(d') + ~(e') * ~x) ^ ~(n') + ~(c') * tan(~(d') + ~(e') * ~x) ^ ~(n2')) ^ ~(p'), ~x) => Module([g = FreeFactors(sin(~d + ~e * ~x), ~x)], (g / ~e) * Subst(Antiderivative((1 - g ^ 2 * (~x) ^ 2) ^ (((~m - 2 * ~n * ~p) - 1) / 2) * ExpandToSum(~c * (~x) ^ (2 * ~n) + ~b * (~x) ^ ~n * (1 - (~x) ^ 2) ^ (~n / 2) + ~a * (1 - (~x) ^ 2) ^ ~n, ~x) ^ ~p, ~x), ~x, sin(~d + ~e * ~x) / g)) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (EqQ(~n2, 2 * ~n) && (IntegerQ((~m - 1) / 2) && (IntegerQ(~n / 2) && IntegerQ(~p))))
    @apply_utils Antiderivative(sin(~(d') + ~(e') * ~x) ^ ~m * (~(a') + ~(b') * cot(~(d') + ~(e') * ~x) ^ ~(n') + ~(c') * cot(~(d') + ~(e') * ~x) ^ ~(n2')) ^ ~(p'), ~x) => Module([g = FreeFactors(cos(~d + ~e * ~x), ~x)], (-g / ~e) * Subst(Antiderivative((1 - g ^ 2 * (~x) ^ 2) ^ (((~m - 2 * ~n * ~p) - 1) / 2) * ExpandToSum(~c * (~x) ^ (2 * ~n) + ~b * (~x) ^ ~n * (1 - (~x) ^ 2) ^ (~n / 2) + ~a * (1 - (~x) ^ 2) ^ ~n, ~x) ^ ~p, ~x), ~x, cos(~d + ~e * ~x) / g)) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (EqQ(~n2, 2 * ~n) && (IntegerQ((~m - 1) / 2) && (IntegerQ(~n / 2) && IntegerQ(~p))))
    @apply_utils Antiderivative(tan(~(d') + ~(e') * ~x) ^ ~(m') * (~(a') + ~(b') * tan(~(d') + ~(e') * ~x) ^ ~(n') + ~(c') * tan(~(d') + ~(e') * ~x) ^ ~(n2')) ^ ~p, ~x) => (1 / (4 ^ ~p * (~c) ^ ~p)) * Antiderivative(tan(~d + ~e * ~x) ^ ~m * (~b + 2 * ~c * tan(~d + ~e * ~x) ^ ~n) ^ (2 * ~p), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~n], ~x) && (EqQ(~n2, 2 * ~n) && (EqQ((~b) ^ 2 - 4 * ~a * ~c, 0) && IntegerQ(~p)))
    @apply_utils Antiderivative(cot(~(d') + ~(e') * ~x) ^ ~(m') * (~(a') + ~(b') * cot(~(d') + ~(e') * ~x) ^ ~(n') + ~(c') * cot(~(d') + ~(e') * ~x) ^ ~(n2')) ^ ~p, ~x) => (1 / (4 ^ ~p * (~c) ^ ~p)) * Antiderivative(cot(~d + ~e * ~x) ^ ~m * (~b + 2 * ~c * cot(~d + ~e * ~x) ^ ~n) ^ (2 * ~p), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~n], ~x) && (EqQ(~n2, 2 * ~n) && (EqQ((~b) ^ 2 - 4 * ~a * ~c, 0) && IntegerQ(~p)))
    @apply_utils Antiderivative(tan(~(d') + ~(e') * ~x) ^ ~(m') * (~(a') + ~(b') * tan(~(d') + ~(e') * ~x) ^ ~(n') + ~(c') * tan(~(d') + ~(e') * ~x) ^ ~(n2')) ^ ~p, ~x) => ((~a + ~b * tan(~d + ~e * ~x) ^ ~n + ~c * tan(~d + ~e * ~x) ^ (2 * ~n)) ^ ~p / (~b + 2 * ~c * tan(~d + ~e * ~x) ^ ~n) ^ (2 * ~p)) * Antiderivative(tan(~d + ~e * ~x) ^ ~m * (~b + 2 * ~c * tan(~d + ~e * ~x) ^ ~n) ^ (2 * ~p), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~n, ~p], ~x) && (EqQ(~n2, 2 * ~n) && (EqQ((~b) ^ 2 - 4 * ~a * ~c, 0) && Not(IntegerQ(~p))))
    @apply_utils Antiderivative(cot(~(d') + ~(e') * ~x) ^ ~(m') * (~(a') + ~(b') * cot(~(d') + ~(e') * ~x) ^ ~(n') + ~(c') * cot(~(d') + ~(e') * ~x) ^ ~(n2')) ^ ~p, ~x) => ((~a + ~b * cot(~d + ~e * ~x) ^ ~n + ~c * cot(~d + ~e * ~x) ^ (2 * ~n)) ^ ~p / (~b + 2 * ~c * cot(~d + ~e * ~x) ^ ~n) ^ (2 * ~p)) * Antiderivative(cot(~d + ~e * ~x) ^ ~m * (~b + 2 * ~c * cot(~d + ~e * ~x) ^ ~n) ^ (2 * ~p), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~n, ~p], ~x) && (EqQ(~n2, 2 * ~n) && (EqQ((~b) ^ 2 - 4 * ~a * ~c, 0) && Not(IntegerQ(~p))))
    @apply_utils Antiderivative(tan(~(d') + ~(e') * ~x) ^ ~(m') * (~(a') + ~(b') * (~(f') * tan(~(d') + ~(e') * ~x)) ^ ~(n') + ~(c') * (~(f') * tan(~(d') + ~(e') * ~x)) ^ ~(n2')) ^ ~p, ~x) => (~f / ~e) * Subst(Antiderivative(((~x / ~f) ^ ~m * (~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ ~p) / ((~f) ^ 2 + (~x) ^ 2), ~x), ~x, ~f * tan(~d + ~e * ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m, ~n, ~p], ~x) && (EqQ(~n2, 2 * ~n) && NeQ((~b) ^ 2 - 4 * ~a * ~c, 0))
    @apply_utils Antiderivative(cot(~(d') + ~(e') * ~x) ^ ~(m') * (~(a') + ~(b') * (~(f') * cot(~(d') + ~(e') * ~x)) ^ ~(n') + ~(c') * (~(f') * cot(~(d') + ~(e') * ~x)) ^ ~(n2')) ^ ~p, ~x) => (-(~f) / ~e) * Subst(Antiderivative(((~x / ~f) ^ ~m * (~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ ~p) / ((~f) ^ 2 + (~x) ^ 2), ~x), ~x, ~f * cot(~d + ~e * ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m, ~n, ~p], ~x) && (EqQ(~n2, 2 * ~n) && NeQ((~b) ^ 2 - 4 * ~a * ~c, 0))
    @apply_utils Antiderivative(cot(~(d') + ~(e') * ~x) ^ ~(m') * (~(a') + ~(b') * tan(~(d') + ~(e') * ~x) ^ ~(n') + ~(c') * tan(~(d') + ~(e') * ~x) ^ ~(n2')) ^ ~(p'), ~x) => (1 / (4 ^ ~p * (~c) ^ ~p)) * Antiderivative(cot(~d + ~e * ~x) ^ ~m * (~b + 2 * ~c * tan(~d + ~e * ~x) ^ ~n) ^ (2 * ~p), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~n], ~x) && (EqQ(~n2, 2 * ~n) && (EqQ((~b) ^ 2 - 4 * ~a * ~c, 0) && IntegerQ(~p)))
    @apply_utils Antiderivative(tan(~(d') + ~(e') * ~x) ^ ~(m') * (~(a') + ~(b') * cot(~(d') + ~(e') * ~x) ^ ~(n') + ~(c') * cot(~(d') + ~(e') * ~x) ^ ~(n2')) ^ ~(p'), ~x) => (1 / (4 ^ ~p * (~c) ^ ~p)) * Antiderivative(tan(~d + ~e * ~x) ^ ~m * (~b + 2 * ~c * cot(~d + ~e * ~x) ^ ~n) ^ (2 * ~p), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~n], ~x) && (EqQ(~n2, 2 * ~n) && (EqQ((~b) ^ 2 - 4 * ~a * ~c, 0) && IntegerQ(~p)))
    @apply_utils Antiderivative(cot(~(d') + ~(e') * ~x) ^ ~(m') * (~(a') + ~(b') * tan(~(d') + ~(e') * ~x) ^ ~(n') + ~(c') * tan(~(d') + ~(e') * ~x) ^ ~(n2')) ^ ~p, ~x) => ((~a + ~b * tan(~d + ~e * ~x) ^ ~n + ~c * tan(~d + ~e * ~x) ^ (2 * ~n)) ^ ~p / (~b + 2 * ~c * tan(~d + ~e * ~x) ^ ~n) ^ (2 * ~p)) * Antiderivative(cot(~d + ~e * ~x) ^ ~m * (~b + 2 * ~c * tan(~d + ~e * ~x) ^ ~n) ^ (2 * ~p), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~n, ~p], ~x) && (EqQ(~n2, 2 * ~n) && (EqQ((~b) ^ 2 - 4 * ~a * ~c, 0) && Not(IntegerQ(~p))))
    @apply_utils Antiderivative(tan(~(d') + ~(e') * ~x) ^ ~(m') * (~(a') + ~(b') * cot(~(d') + ~(e') * ~x) ^ ~(n') + ~(c') * cot(~(d') + ~(e') * ~x) ^ ~(n2')) ^ ~p, ~x) => ((~a + ~b * cot(~d + ~e * ~x) ^ ~n + ~c * cot(~d + ~e * ~x) ^ (2 * ~n)) ^ ~p / (~b + 2 * ~c * cot(~d + ~e * ~x) ^ ~n) ^ (2 * ~p)) * Antiderivative(tan(~d + ~e * ~x) ^ ~m * (~b + 2 * ~c * cot(~d + ~e * ~x) ^ ~n) ^ (2 * ~p), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~n, ~p], ~x) && (EqQ(~n2, 2 * ~n) && (EqQ((~b) ^ 2 - 4 * ~a * ~c, 0) && Not(IntegerQ(~p))))
    @apply_utils Antiderivative(cot(~(d') + ~(e') * ~x) ^ ~(m') * (~(a') + ~(b') * tan(~(d') + ~(e') * ~x) ^ ~n + ~(c') * tan(~(d') + ~(e') * ~x) ^ ~n2) ^ ~(p'), ~x) => Module([g = FreeFactors(cot(~d + ~e * ~x), ~x)], (g / ~e) * Subst(Antiderivative(((g * ~x) ^ (~m - 2 * ~n * ~p) * (~c + ~b * (g * ~x) ^ ~n + ~a * (g * ~x) ^ (2 * ~n)) ^ ~p) / (1 + g ^ 2 * (~x) ^ 2), ~x), ~x, cot(~d + ~e * ~x) / g)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~p], ~x) && (EqQ(~n2, 2 * ~n) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && IntegerQ(~n / 2)))
    @apply_utils Antiderivative(tan(~(d') + ~(e') * ~x) ^ ~(m') * (~(a') + ~(b') * cot(~(d') + ~(e') * ~x) ^ ~n + ~(c') * cot(~(d') + ~(e') * ~x) ^ ~n2) ^ ~(p'), ~x) => Module([g = FreeFactors(tan(~d + ~e * ~x), ~x)], (-g / ~e) * Subst(Antiderivative(((g * ~x) ^ (~m - 2 * ~n * ~p) * (~c + ~b * (g * ~x) ^ ~n + ~a * (g * ~x) ^ (2 * ~n)) ^ ~p) / (1 + g ^ 2 * (~x) ^ 2), ~x), ~x, tan(~d + ~e * ~x) / g)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~p], ~x) && (EqQ(~n2, 2 * ~n) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && IntegerQ(~n / 2)))
    @apply_utils Antiderivative((~A + ~(B') * tan(~(d') + ~(e') * ~x)) * (~a + ~(b') * tan(~(d') + ~(e') * ~x) + ~(c') * tan(~(d') + ~(e') * ~x) ^ 2) ^ ~n, ~x) => (1 / (4 ^ ~n * (~c) ^ ~n)) * Antiderivative((~A + ~B * tan(~d + ~e * ~x)) * (~b + 2 * ~c * tan(~d + ~e * ~x)) ^ (2 * ~n), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~A, ~B], ~x) && (EqQ((~b) ^ 2 - 4 * ~a * ~c, 0) && IntegerQ(~n))
    @apply_utils Antiderivative((~A + ~(B') * cot(~(d') + ~(e') * ~x)) * (~a + ~(b') * cot(~(d') + ~(e') * ~x) + ~(c') * cot(~(d') + ~(e') * ~x) ^ 2) ^ ~n, ~x) => (1 / (4 ^ ~n * (~c) ^ ~n)) * Antiderivative((~A + ~B * cot(~d + ~e * ~x)) * (~b + 2 * ~c * cot(~d + ~e * ~x)) ^ (2 * ~n), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~A, ~B], ~x) && (EqQ((~b) ^ 2 - 4 * ~a * ~c, 0) && IntegerQ(~n))
    @apply_utils Antiderivative((~A + ~(B') * tan(~(d') + ~(e') * ~x)) * (~a + ~(b') * tan(~(d') + ~(e') * ~x) + ~(c') * tan(~(d') + ~(e') * ~x) ^ 2) ^ ~n, ~x) => ((~a + ~b * tan(~d + ~e * ~x) + ~c * tan(~d + ~e * ~x) ^ 2) ^ ~n / (~b + 2 * ~c * tan(~d + ~e * ~x)) ^ (2 * ~n)) * Antiderivative((~A + ~B * tan(~d + ~e * ~x)) * (~b + 2 * ~c * tan(~d + ~e * ~x)) ^ (2 * ~n), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~A, ~B], ~x) && (EqQ((~b) ^ 2 - 4 * ~a * ~c, 0) && Not(IntegerQ(~n)))
    @apply_utils Antiderivative((~A + ~(B') * cot(~(d') + ~(e') * ~x)) * (~a + ~(b') * cot(~(d') + ~(e') * ~x) + ~(c') * cot(~(d') + ~(e') * ~x) ^ 2) ^ ~n, ~x) => ((~a + ~b * cot(~d + ~e * ~x) + ~c * cot(~d + ~e * ~x) ^ 2) ^ ~n / (~b + 2 * ~c * cot(~d + ~e * ~x)) ^ (2 * ~n)) * Antiderivative((~A + ~B * cot(~d + ~e * ~x)) * (~b + 2 * ~c * cot(~d + ~e * ~x)) ^ (2 * ~n), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~A, ~B], ~x) && (EqQ((~b) ^ 2 - 4 * ~a * ~c, 0) && Not(IntegerQ(~n)))
    @apply_utils Antiderivative((~A + ~(B') * tan(~(d') + ~(e') * ~x)) / (~(a') + ~(b') * tan(~(d') + ~(e') * ~x) + ~(c') * tan(~(d') + ~(e') * ~x) ^ 2), ~x) => Module([q = Rt((~b) ^ 2 - 4 * ~a * ~c, 2)], (~B + (~b * ~B - 2 * ~A * ~c) / q) * Antiderivative(1 / Simp(~b + q + 2 * ~c * tan(~d + ~e * ~x), ~x), ~x) + (~B - (~b * ~B - 2 * ~A * ~c) / q) * Antiderivative(1 / Simp((~b - q) + 2 * ~c * tan(~d + ~e * ~x), ~x), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~A, ~B], ~x) && NeQ((~b) ^ 2 - 4 * ~a * ~c, 0)
    @apply_utils Antiderivative((~A + ~(B') * cot(~(d') + ~(e') * ~x)) / (~(a') + ~(b') * cot(~(d') + ~(e') * ~x) + ~(c') * cot(~(d') + ~(e') * ~x) ^ 2), ~x) => Module([q = Rt((~b) ^ 2 - 4 * ~a * ~c, 2)], (~B + (~b * ~B - 2 * ~A * ~c) / q) * Antiderivative(1 / Simp(~b + q + 2 * ~c * cot(~d + ~e * ~x), ~x), ~x) + (~B - (~b * ~B - 2 * ~A * ~c) / q) * Antiderivative(1 / Simp((~b - q) + 2 * ~c * cot(~d + ~e * ~x), ~x), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~A, ~B], ~x) && NeQ((~b) ^ 2 - 4 * ~a * ~c, 0)
    @apply_utils Antiderivative((~A + ~(B') * tan(~(d') + ~(e') * ~x)) * (~(a') + ~(b') * tan(~(d') + ~(e') * ~x) + ~(c') * tan(~(d') + ~(e') * ~x) ^ 2) ^ ~n, ~x) => Antiderivative(ExpandTrig((~A + ~B * tan(~d + ~e * ~x)) * (~a + ~b * tan(~d + ~e * ~x) + ~c * tan(~d + ~e * ~x) ^ 2) ^ ~n, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~A, ~B], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && IntegerQ(~n))
    @apply_utils Antiderivative((~A + ~(B') * cot(~(d') + ~(e') * ~x)) * (~(a') + ~(b') * cot(~(d') + ~(e') * ~x) + ~(c') * cot(~(d') + ~(e') * ~x) ^ 2) ^ ~n, ~x) => Antiderivative(ExpandTrig((~A + ~B * cot(~d + ~e * ~x)) * (~a + ~b * cot(~d + ~e * ~x) + ~c * cot(~d + ~e * ~x) ^ 2) ^ ~n, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~A, ~B], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && IntegerQ(~n))
end

max_stars_repo_path: src/rules/4 Trig functions/4.3 Tangent/4.3.9 trig^m (a+b tan^n+c tan^(2 n))^p.jl


text:
function set_blob_data{T}(data::Array{T}, blob::CuTensorBlob{T}, blob_idx::Int)
  n_fea = get_fea_size(blob)
  ptr = Base.unsafe_convert(Ptr{Void}, blob.ptr.p) + sizeof(T) * n_fea * (blob_idx-1) # note 0-based indexing in CUDA Vector
  CuBLAS.set_vector(length(data), sizeof(T), convert(Ptr{Void},pointer(data)), 1, ptr, 1)
end


text:
const stata_epoch_time = DateTime(1960, 1, 1)
const stata_epoch_date = Date(1960, 1, 1)
const sas_epoch_time = DateTime(1960, 1, 1)
const sas_epoch_date = Date(1960, 1, 1)
const spss_epoch_time = DateTime(1582, 10, 14)

# Reference: Stata documentation
const stata_dt_formats = Dict{String, Tuple{Union{DateTime,Date}, Period}}(
    "%tc" => (stata_epoch_time, Millisecond(1)),
    "%td" => (stata_epoch_date, Day(1)),
    "%tw" => (stata_epoch_date, Week(1)),
    "%tm" => (stata_epoch_date, Month(1)),
    "%tq" => (stata_epoch_date, Month(3)),
    "%th" => (stata_epoch_date, Month(6)),
    "%ty" => (Date(0), Year(1))
)

# Reference: https://github.com/Roche/pyreadstat/blob/master/pyreadstat/_readstat_parser.pyx
const sas_dt_formats = Dict{String, Tuple{Union{DateTime,Date}, Period}}(
    "TIME" => (sas_epoch_time, Second(1)),
    "HHMM" => (sas_epoch_time, Second(1)),
    "TIME20.3" => (sas_epoch_time, Second(1)),
    "DATETIME" => (sas_epoch_time, Second(1)),
    "DATETIME20" => (sas_epoch_time, Second(1)),
    "WEEKDATE" => (sas_epoch_date, Day(1)),
    "MMDDYY" => (sas_epoch_date, Day(1)),
    "DDMMYY" => (sas_epoch_date, Day(1)),
    "YYMMDD" => (sas_epoch_date, Day(1)),
    "DATE" => (sas_epoch_date, Day(1)),
    "DATE9" => (sas_epoch_date, Day(1)),
    "YYMMDD10" => (sas_epoch_date, Day(1))
)

# Reference: https://github.com/Roche/pyreadstat/blob/master/pyreadstat/_readstat_parser.pyx
const spss_dt_formats = Dict{String, Tuple{Union{DateTime,Date}, Period}}(
    "TIME" => (spss_epoch_time, Second(1)),
    "DTIME" => (spss_epoch_time, Second(1)),
    "TIME8" => (spss_epoch_time, Second(1)),
    "TIME5" => (spss_epoch_time, Second(1)),
    "TIME11.2" => (spss_epoch_time, Second(1)),
    "DATETIME" => (spss_epoch_time, Second(1)),
    "DATETIME20" => (spss_epoch_time, Second(1)),
    "DATETIME23.2" => (spss_epoch_time, Second(1)),
    "DATETIME8" => (spss_epoch_time, Second(1)),
    "YMDHMS20" => (spss_epoch_time, Second(1)),
    "DATE" => (spss_epoch_time, Second(1)),
    "ADATE" => (spss_epoch_time, Second(1)),
    "EDATE" => (spss_epoch_time, Second(1)),
    "JDATE" => (spss_epoch_time, Second(1)),
    "SDATE" => (spss_epoch_time, Second(1)),
    "EDATE10" => (spss_epoch_time, Second(1)),
    "DATE8" => (spss_epoch_time, Second(1)),
    "EDATE8" => (spss_epoch_time, Second(1)),
    "DATE11" => (spss_epoch_time, Second(1))
)

const dt_formats = Dict{Val, Dict}(
    Val(:dta) => stata_dt_formats,
    Val(:sav) => spss_dt_formats,
    Val(:por) => spss_dt_formats,
    Val(:sas7bdat) => sas_dt_formats,
    Val(:xport) => sas_dt_formats
)

"""
    parse_datetime(col::Vector, epoch::Union{DateTime,Date}, delta::Period)
    parse_datetime(col::Vector, epoch::Union{DateTime,Date}, delta::Period, missingvalue)

Construct a vector of time values of type `DateTime` or `Date`
by interpreting the elements in `col` as the number of periods passed
since `epoch` with the length of each period being `delta`.
If `missingvalue` is specified,
indices where the elements in `col` are equal to `missingvalue` based on `isequal`
are set to be `missing` no matter what is specified with `missingvalue`.
"""
function parse_datetime(col::Vector, epoch::Union{DateTime,Date}, delta::Period)
    out = Vector{typeof(epoch)}(undef, length(col))
    @inbounds for i in eachindex(col)
        v = col[i]
        out[i] = epoch + v * delta
    end
    return out
end

# Missing value is set to be missing no matter what is specified with missingvalue
function parse_datetime(col::Vector, epoch::Union{DateTime,Date}, delta::Period, missingvalue)
    out = Vector{Union{typeof(epoch), Missing}}(undef, length(col))
    @inbounds for i in eachindex(col)
        v = col[i]
        out[i] = isequal(v, missingvalue) ? missing : epoch + v * delta
    end
    return out
end


text:
include("abstract_map.jl")
include("map_translate.jl")
include("map_affine.jl")
include("map_quadratic.jl")

include("map_from_vectfunc.jl")
include("transform_src_insert_dest_map.jl")


text:
using Tokens
using Test

function testshared()
end

# tests of locate in IOShared
function testlocate()

end


text:
# Various heap implementation

###########################################################
#
#   Heap interface specification
#
#   Each heap is associated with a handle type (H), and
#   a value type v.
#
#   Here, the value type must be comparable, and a handle
#   is an object through which one can refer to a specific
#   node of the heap and thus update its value.
#
#   Each heap type must implement all of the following
#   functions. Here, let h be a heap, i be a handle, and
#   v be a value.
#
#   - length(h)         returns the number of elements
#
#   - isempty(h)        returns whether the heap is
#                       empty
#
#   - push!(h, v)       add a value to the heap
#
#   - sizehint(h)       set size hint to a heap
#
#   - top(h)            return the top value of a heap
#
#   - pop!(h)           removes the top value, and
#                       returns it
#
#  For mutable heaps, it should also support
#
#   - push!(h, v)       adds a value to the heap and
#                       returns a handle to v
#
#   - update!(h, i, v)  updates the value of an element
#                       (referred to by the handle i)
#
###########################################################

# HT: handle type
# VT: value type

abstract AbstractHeap{VT}

abstract AbstractMutableHeap{VT,HT} <: AbstractHeap{VT}

# comparer

immutable LessThan
end

immutable GreaterThan
end

compare(c::LessThan, x, y) = x < y
compare(c::GreaterThan, x, y) = x > y

# heap implementations

include("heaps/binary_heap.jl")
include("heaps/mutable_binary_heap.jl")

# generic functions

function extract_all!{VT}(h::AbstractHeap{VT})
    n = length(h)
    r = Array(VT, n)
    for i = 1 : n
        r[i] = pop!(h)
    end
    r
end