text:
# Collections of ScalarGridData

macro collectionfield(wrapper,nctypes)

  ctype = Symbol[]
  ctypetype = Expr[]
  crtype = Expr[]
  for i in 1:eval(nctypes)
    c = Symbol("C",string(i))
    push!(ctype,c)
    push!(ctypetype,:(::Type{$c}))
    push!(crtype,:($c<:CellType))
  end

  return esc(quote

    # This is the basic function that uses reshape(view()) of an input vector
    # to provide the components of this type
    function (::Type{$wrapper{$(ctype...),NX,NY,T,DT}})(data::AbstractVector{R}) where
                              {$(ctype...),NX,NY,T<:Number,DT,R}
      n0 = 0
      flist = []
      gtype = $wrapper{$(ctype...),NX,NY,R,typeof(data)}
      for ft in _fieldtypes(gtype)
        if ft <: GridData
          ftype = griddatatype(ft)
          dims = size(ft)
          dn = prod(dims)
          fdata = reshape(view(data,n0+1:n0+dn),dims)
          n0 += dn
          push!(flist,ftype{celltype(ft),NX,NY,R,typeof(fdata)}(fdata))
        end
      end
      gtype(data,flist...)
    end

    # This provides the most straightforward constructor, calling upon the constructor
    # above to deal with the split between components
    function $wrapper($(ctypetype...), dualnodedims::Tuple{Int, Int};dtype=Float64) where {$(crtype...)}
        len = 0
        for ft in _fieldtypes($wrapper{$(ctype...),dualnodedims...,dtype})
          if ft <: GridData
            len += prod(size(ft))
          end
        end
        data = zeros(dtype,len)
        $wrapper{$(ctype...),dualnodedims...,dtype,typeof(data)}(data)
    end

    # Create several other basic functions for this type
    @griddata($wrapper,$nctypes)

  end)
end

"""
    Edges

`Edges` is a wrapper for vector-valued data that lie at the faces of either dual cells or
primary cells. `Edges` type data have fields `u` and `v` for the components of the
vector field. These are the normal components of the vector field on the vertical
and horizontal faces of the corresponding cell.

# Constructors
- `Edges(C,dims)` creates a vector field of zeros in cells of type `C` (where `C` is
  either `Dual` or `Primal`), on a grid of dimensions `dims`. Note that `dims`
  represent the number of dual cells on the grid.
- `Edges(C,w)` performs the same construction, but uses existing field data `w`
  of `GridData` type to determine the size of the grid.
-  Adding the `dtype=` keyword allows the data type of the field data to be
  changed. The default is `Float64`, but can be changed to, e.g., `ComplexF64`
"""
struct Edges{C <: CellType, NX, NY, T <: Number, DT} <: VectorGridData{NX,NY,T}
    data::DT
    u::XEdges{C,NX,NY,T}
    v::YEdges{C,NX,NY,T}
end

@collectionfield Edges 1

# Based on number of dual nodes, return the number of edges
edge_inds(::Type{C},   dualnodedims) where {C <: CellType} =
            xedge_inds(C,dualnodedims), yedge_inds(C,dualnodedims)

function Base.show(io::IO, edges::Edges{C, NX, NY, T, DT}) where {C, NX, NY, T, DT}
    nodedims = "(nx = $NX, ny = $NY)"
    udims = "(nx = $(size(edges.u,1)), ny = $(size(edges.u,2)))"
    vdims = "(nx = $(size(edges.v,1)), ny = $(size(edges.v,2)))"
    println(io, "$C edges for a $nodedims cell grid of type $T data")
    println(io, "  Internal u-faces: $udims")
    print(io, "  Internal v-faces: $vdims")
end

function Base.show(io::IO, m::MIME"text/plain", edges::Edges)
    println(io,"$(typeof(edges)) data")
    println(io,"u (in grid orientation)")
    show(io,m,reverse(transpose(edges.u),dims=1))
    println(io)
    println(io,"v (in grid orientation)")
    show(io,m,reverse(transpose(edges.v),dims=1))
end

### Edge gradient ###

"""
    EdgeGradient

`EdgeGradient` is a wrapper for tensor-valued data that lie partly at the nodes of dual cells and
primary cells. `EdgeGradient` type data have fields `dudx`, `dudy`, `dvdx`, `dvdy` for the components of the
tensor field. The diagonal components lie at one set of nodes (e.g. Primal),
and the offdiagonal at the other set (e.g. Dual).

# Constructors
- `EdgeGradient(C,dims)` creates a tensor field of zeros in cells of type `C` (where `C` is
  either `Dual` or `Primal`), on a grid of dimensions `dims`. Note that `dims`
  represent the number of dual cells on the grid.
- `EdgeGradient(C,w)` performs the same construction, but uses existing field data `w`
  of `GridData` type to determine the size of the grid.
-  Adding the `dtype=` keyword allows the data type of the field data to be
  changed. The default is `Float64`, but can be changed to, e.g., `ComplexF64`
"""
struct EdgeGradient{C <: CellType,D <: CellType, NX,NY, T<: Number, DT} <: GridData{NX,NY,T}
  data :: DT
  dudx :: Nodes{C,NX,NY,T}
  dudy :: Nodes{D,NX,NY,T}
  dvdx :: Nodes{D,NX,NY,T}
  dvdy :: Nodes{C,NX,NY,T}
end

@collectionfield EdgeGradient 2

# Only provide one of the cell types -- the other cell type is set automatically
EdgeGradient(::Type{C}, dualnodedims::Tuple{Int, Int};dtype=Float64) where {C <: CellType} =
    EdgeGradient(C,othertype(C),dualnodedims,dtype=dtype)

function Base.show(io::IO, nodes::EdgeGradient{R, S, NX, NY, T, DT}) where {R, S, NX, NY, T, DT}
    nodedims = "(nx = $NX, ny = $NY)"
    udims = "(nx = $(size(nodes.dudx,1)), ny = $(size(nodes.dudx,2)))"
    vdims = "(nx = $(size(nodes.dudy,1)), ny = $(size(nodes.dudy,2)))"
    println(io, "Edge gradient data in a $nodedims cell grid of type $T data")
    println(io, "  du/dx and dv/dy as $R nodes : $udims")
    println(io, "  du/dy and dv/dx as $S nodes : $vdims")
end

function Base.show(io::IO, m::MIME"text/plain", nodes::EdgeGradient)
    println(io,"$(typeof(nodes)) data")
    println(io,"du/dx (in grid orientation)")
    show(io,m,reverse(transpose(nodes.dudx),dims=1))
    println(io)
    println(io,"du/dy (in grid orientation)")
    show(io,m,reverse(transpose(nodes.dudy),dims=1))
    println(io)
    println(io,"dv/dx (in grid orientation)")
    show(io,m,reverse(transpose(nodes.dvdx),dims=1))
    println(io)
    println(io,"dv/dy (in grid orientation)")
    show(io,m,reverse(transpose(nodes.dvdy),dims=1))
end

### Node pair ###

"""
    NodePair

`NodePair` is a wrapper for vector-valued data that lie at the nodes of dual cells and
primal cells. `NodePair` type data have fields `u` and `v` for the components of the
vector field. These are the normal components of a vector field on nodes that
form the faces of a virtual cell centered at one of the faces of the primal cell.

# Constructors
- `NodePair(C,dims)` creates a vector field of zeros in cells of type `C` (where `C` is
  either `Dual` or `Primal`), on a grid of dimensions `dims`. Note that `dims`
  represent the number of dual cells on the grid.
- `NodePair(C,w)` performs the same construction, but uses existing field data `w`
  of `GridData` type to determine the size of the grid.
-  Adding the `dtype=` keyword allows the data type of the field data to be
  changed. The default is `Float64`, but can be changed to, e.g., `ComplexF64`
"""
struct NodePair{C <: CellType,D <: CellType, NX,NY,T, DT} <: GridData{NX,NY,T}
  data :: DT
  u :: Nodes{C,NX,NY,T}
  v :: Nodes{D,NX,NY,T}
end

@collectionfield NodePair 2

# Only provide one of the cell types -- the other cell type is set automatically
NodePair(::Type{C}, dualnodedims::Tuple{Int, Int};dtype=Float64) where {C <: CellType} =
    NodePair(C,othertype(C),dualnodedims,dtype=dtype)

function Base.show(io::IO, nodes::NodePair{R, S, NX, NY, T}) where {R, S, NX, NY, T}
    nodedims = "(nx = $NX, ny = $NY)"
    udims = "(nx = $(size(nodes.u,1)), ny = $(size(nodes.u,2)))"
    vdims = "(nx = $(size(nodes.v,1)), ny = $(size(nodes.v,2)))"
    println(io, "($R,$S) node data pair in a $nodedims cell grid of type $T data")
    println(io, "  Number of $R nodes: $udims")
    println(io, "  Number of $S nodes: $vdims")
end

function Base.show(io::IO, m::MIME"text/plain", nodes::NodePair)
    println(io,"$(typeof(nodes)) data")
    println(io,"u (in grid orientation)")
    show(io,m,reverse(transpose(nodes.u),dims=1))
    println(io)
    println(io,"v (in grid orientation)")
    show(io,m,reverse(transpose(nodes.v),dims=1))
end


text:
using CompileBot

println("tests infer benchmark")

snoop_bench(BotConfig("TestPackage3", tmin =0.0))

ext: jl


text:
# Automatically generated using Clang.jl

const CUTENSOR_MAJOR = 1
const CUTENSOR_MINOR = 0
const CUTENSOR_PATCH = 0

@cenum cutensorAlgo_t::Int32 begin
    CUTENSOR_ALGO_GETT = -4
    CUTENSOR_ALGO_TGETT = -3
    CUTENSOR_ALGO_TTGT = -2
    CUTENSOR_ALGO_DEFAULT = -1
end

@cenum cutensorWorksizePreference_t::UInt32 begin
    CUTENSOR_WORKSPACE_MIN = 1
    CUTENSOR_WORKSPACE_RECOMMENDED = 2
    CUTENSOR_WORKSPACE_MAX = 3
end

@cenum cutensorOperator_t::UInt32 begin
    CUTENSOR_OP_IDENTITY = 1
    CUTENSOR_OP_SQRT = 2
    CUTENSOR_OP_RELU = 8
    CUTENSOR_OP_CONJ = 9
    CUTENSOR_OP_RCP = 10
    CUTENSOR_OP_SIGMOID = 11
    CUTENSOR_OP_TANH = 12
    CUTENSOR_OP_EXP = 22
    CUTENSOR_OP_LOG = 23
    CUTENSOR_OP_ABS = 24
    CUTENSOR_OP_NEG = 25
    CUTENSOR_OP_SIN = 26
    CUTENSOR_OP_COS = 27
    CUTENSOR_OP_TAN = 28
    CUTENSOR_OP_SINH = 29
    CUTENSOR_OP_COSH = 30
    CUTENSOR_OP_ASIN = 31
    CUTENSOR_OP_ACOS = 32
    CUTENSOR_OP_ATAN = 33
    CUTENSOR_OP_ASINH = 34
    CUTENSOR_OP_ACOSH = 35
    CUTENSOR_OP_ATANH = 36
    CUTENSOR_OP_CEIL = 37
    CUTENSOR_OP_FLOOR = 38
    CUTENSOR_OP_ADD = 3
    CUTENSOR_OP_MUL = 5
    CUTENSOR_OP_MAX = 6
    CUTENSOR_OP_MIN = 7
    CUTENSOR_OP_UNKNOWN = 126
end

@cenum cutensorStatus_t::UInt32 begin
    CUTENSOR_STATUS_SUCCESS = 0
    CUTENSOR_STATUS_NOT_INITIALIZED = 1
    CUTENSOR_STATUS_ALLOC_FAILED = 3
    CUTENSOR_STATUS_INVALID_VALUE = 7
    CUTENSOR_STATUS_ARCH_MISMATCH = 8
    CUTENSOR_STATUS_MAPPING_ERROR = 11
    CUTENSOR_STATUS_EXECUTION_FAILED = 13
    CUTENSOR_STATUS_INTERNAL_ERROR = 14
    CUTENSOR_STATUS_NOT_SUPPORTED = 15
    CUTENSOR_STATUS_LICENSE_ERROR = 16
    CUTENSOR_STATUS_CUBLAS_ERROR = 17
    CUTENSOR_STATUS_CUDA_ERROR = 18
    CUTENSOR_STATUS_INSUFFICIENT_WORKSPACE = 19
    CUTENSOR_STATUS_INSUFFICIENT_DRIVER = 20
end

@cenum cutensorComputeType_t::UInt32 begin
    CUTENSOR_R_MIN_16F = 1
    CUTENSOR_C_MIN_16F = 2
    CUTENSOR_R_MIN_32F = 4
    CUTENSOR_C_MIN_32F = 8
    CUTENSOR_R_MIN_64F = 16
    CUTENSOR_C_MIN_64F = 32
    CUTENSOR_R_MIN_8U = 64
    CUTENSOR_R_MIN_32U = 128
    CUTENSOR_R_MIN_8I = 256
    CUTENSOR_R_MIN_32I = 512
end

function cutensorComputeType(T::DataType)
    if T == Float32
        return CUTENSOR_R_MIN_32F
    elseif T == ComplexF32
        return CUTENSOR_C_MIN_32F
    elseif T == Float16
        return CUTENSOR_R_MIN_16F
    elseif T == ComplexF16
        return CUTENSOR_C_MIN_16F
    elseif T == Float64
        return CUTENSOR_R_MIN_64F
    elseif T == ComplexF64
        return CUTENSOR_C_MIN_64F
    elseif T == Int8
        return CUTENSOR_R_MIN_8I
    elseif T == Int32
        return CUTENSOR_R_MIN_32I
    elseif T == UInt8
        return CUTENSOR_R_MIN_8U
    elseif T == UInt32
        return CUTENSOR_R_MIN_32U
    else
        throw(ArgumentError("cutensorComputeType equivalent for input type $T does not exist!"))
    end
end

struct cutensorHandle_t
    fields::NTuple{512, Int64}
end

struct cutensorTensorDescriptor_t
    fields::NTuple{64, Int64}
end

struct cutensorContractionDescriptor_t
    fields::NTuple{256, Int64}
end

struct cutensorContractionPlan_t
    fields::NTuple{640, Int64}
end

struct cutensorContractionFind_t
    fields::NTuple{64, Int64}
end


text:
using Dates
"""
`line = LinearTransform(k, b)`
`y = line(x)`
`x = inverse(line, y)`

Linear function `kx + b`
"""
struct LinearTransform #??? where {T <: Real} #
    b
    k
    inv_k # inverse of k for inverse function speedup
    LinearTransform(b, k) = new(b, k, 1/k)
end

line = LinearTransform(3.0, 10.0)

function (tr::LinearTransform)(x)
    x * tr.k + tr.b
end
function inverse(tr::LinearTransform, x)
    (x - tr.b) * tr.inv_k
end

"""
TimeGrid to represents time axes of regularly sampled signal,
sratring from `timestart` with optional start `offset`
and sampling rate `fs`
"""
struct TimeGrid
    timestart::DateTime
    fs::Float64
    offset::Int
    tr::LinearTransform # {Float64}
    TimeGrid(timestart::DateTime, fs::Union{Int64,Float64}, offset::Int = 0) =
    new(timestart, fs, offset, LinearTransform(1 - offset, fs / 1000))
end

"""
`toindex2 = index_transform(tg1, tg2)`
get function that transforms indexes of one `TimeGrid` to another.
It then can be used to transform indexes:
`i2 = toindex2(i1)`
"""
function index_transform(tg1::TimeGrid, tg2::TimeGrid)
    index_transform(tg1.tr, tg2.tr)
end
function index_transform(tr1::LinearTransform, tr2::LinearTransform)
    kk = tr1.k / tr2.k
    bb = (tr1.b - tr2.b) / tr2.k
    tr = LinearTransform(bb, kk)
    return x->floor(Int, tr(x)) # transform, then floor to integer index
end

# internal reinterpret functions
@inline time2ms(time::Period) = Dates.value(time)
@inline ms2time(msec) = Millisecond(msec)
# ms2time(msec) = DatePeriod(Dates.UTM(msec)) - transforms to absolute DateTime

"""
`ind = tg[time]`
transform relative time (period from `tg.timestart`) to index
"""
function Base.getindex(tg::TimeGrid, time::Period)
    # ind = Int(time2ms(time) * tg.fs ÷ 1000 - tg.offset + 1)
    ind = floor(Int, tg.tr(time2ms(time)))
end

"""
`ind = tg[time]`
transform absolute time to index
"""
function Base.getindex(tg::TimeGrid, time::DateTime)
    # ind = Int(time2ms(time - tg.timestart) * tg.fs ÷ 1000 - tg.offset + 1)
    ind = floor(Int, tg.tr(time2ms(time - tg.timestart)))
end

"""
`time = tg[ind]`
transform index to absolute time
"""
function Base.getindex(tg::TimeGrid, ind::Real)
    # time = ms2time((ind - 1 + tg.offset) * 1000 / tg.fs) + tg.timestart
    time = ms2time(floor(Int, inverse(tg.tr, ind))) + tg.timestart
end

"""
`time = tg[ind, Period]`
transform index to time, relative to the `tg.timestart`
"""
function Base.getindex(tg::TimeGrid, ind::Real, ::Type{Period})
    # time = ms2time((ind - 1 + tg.offset) * 1000 / tg.fs)
    time = ms2time(floor(Int, inverse(tg.tr, ind)))
end


text:
#=
type and methods to track the evolution of a strategy's securities portfolio composition
=#

mutable struct Portfolio
    quantity::Matrix{Float64}
    weight::Matrix{Float64}
    entry_price::Matrix{Float64}
    close_price::Matrix{Float64}
    pnl::Matrix{Float64}
    idx::Vector{<:TimeType}
    function Portfolio(universe::Universe)
        idx = get_overall_index(universe)
        quantity = weight = entry_price = close_price = pnl =
            zeros(Float64, length(idx), length(universe.assets))
        return new(quantity, weight, entry_price, close_price, pnl, idx)
    end
end

#function update_portfolio!(portfolio::Portfolio, order::Order, universe::Universe)
#    i = findfirst(portfolio.idx .> order.time)
#    quantity[i]
#end


text:
export Transform, IdentityTransform, ScaleTransform, ARDTransform, LowRankTransform, FunctionTransform, ChainTransform
export transform

"""
```julia
    transform(t::Transform, X::AbstractMatrix)
    transform(k::Kernel, X::AbstractMatrix)
```
Apply the transfomration `t` or `k.transform` on the input `X`
"""
transform

include("scaletransform.jl")
include("ardtransform.jl")
include("lowranktransform.jl")
include("functiontransform.jl")
include("selecttransform.jl")
include("chaintransform.jl")

"""
IdentityTransform
Return exactly the input
"""
struct IdentityTransform <: Transform end

params(t::IdentityTransform) = nothing
duplicate(t::IdentityTransform,θ) = t

transform(t::IdentityTransform,x::AbstractArray,obsdim::Int=defaultobs) = x #TODO add test

### TODO Maybe defining adjoints could help but so far it's not working

# @adjoint function ScaleTransform(s::T) where {T<:Real}
#     @check_args(ScaleTransform, s, s > zero(T), "s > 0")
#     ScaleTransform{T}(s),Δ->ScaleTransform{T}(Δ)
# end
#
# @adjoint function ScaleTransform(s::A) where {A<:AbstractVector{<:Real}}
#     @check_args(ScaleTransform, s, all(s.>zero(eltype(A))), "s > 0")
#     ScaleTransform{A}(s),Δ->begin; @show Δ,size(Δ); ScaleTransform{A}(Δ); end
# end

# @adjoint transform(t::ScaleTransform{<:AbstractVector{<:Real}},x::AbstractVector{<:Real}) = transform(t,x),Δ->(ScaleTransform(nothing),t.s.*Δ)
#
#     @adjoint transform(t::ScaleTransform{<:AbstractVector{<:Real}},X::AbstractMatrix{<:Real},obsdim::Int) = transform(t,X,obsdim),Δ->begin
#     @show Δ,size(Δ);
#     return (obsdim == 1 ? ScaleTransform()Δ'.*X : ScaleTransform()Δ.*X,transform(t,Δ,obsdim),nothing)
#     end
#
# @adjoint transform(t::ScaleTransform{T},x::AbstractVecOrMat,obsdim::Int) where {T<:Real} = transform(t,x), Δ->(ScaleTransform(one(T)),t.s.*Δ,nothing)


text:
module SmallInts

export Int1, Int2, Int4
export UInt1, UInt2, UInt4

immutable Int1 <: Signed; val::Int8; end
immutable Int2 <: Signed; val::Int8; end
immutable Int4 <: Signed; val::Int8; end
immutable UInt1 <: Unsigned; val::UInt8; end
immutable UInt2 <: Unsigned; val::UInt8; end
immutable UInt4 <: Unsigned; val::UInt8; end

nbits(::Type{Int1}) = 1
nbits(::Type{Int2}) = 2
nbits(::Type{Int4}) = 4
nbits(::Type{UInt1}) = 1
nbits(::Type{UInt2}) = 2
nbits(::Type{UInt4}) = 4

import Base: print, show
import Base: typemin, typemax, rem, convert, promote_rule
import Base: leading_zeros, leading_ones, trailing_zeros, trailing_ones
import Base: ~, &, |, $
import Base: <<, >>, >>>
import Base: <, <=
import Base: +, -, abs, *, div, rem, fld, mod, cld

typealias BaseSigned Union{Int8, Int16, Int32, Int64, Int128}
typealias BaseUnsigned Union{UInt8, UInt16, UInt32, UInt64, UInt128}
typealias BaseInteger Union{BaseSigned, BaseUnsigned}

for I in (:Int1, :Int2, :Int4)
    @eval begin
        print(io::IO, x::$I) = print(io, x.val)
        show(io::IO, x::$I) = show(io, x.val)

        sext(x::$I) = $I(x.val << (8 - nbits($I)) >> (8 - nbits($I)))

        typemin(::Type{$I}) = $I(Int8(-1) << (nbits($I)-1))
        typemax(::Type{$I}) = $I(Int8(1) << (nbits($I)-1) - Int8(1))
        rem(x::$I, ::Type{$I}) = x
        rem(x::$I, ::Type{Bool}) = rem(x.val, Bool)
        rem{T<:BaseInteger}(x::$I, ::Type{T}) = rem(x.val, T)
        # rem{T<:Integer}(x::$I, ::Type{T}) = rem(x.val, T)
        rem(x::Bool, ::Type{$I}) = $I(rem(x, Int8))
        rem{T<:BaseInteger}(x::T, ::Type{$I}) = sext($I(rem(x, Int8)))
        # rem{T<:Integer}(x::T, ::Type{$I}) = sext($I(rem(x, Int8)))
        convert(::Type{$I}, x::$I) = x
        convert(::Type{$I}, x::Bool) = $I(convert(Int8, x))
        function convert{T<:BaseInteger}(::Type{$I}, x::T)
            (x<typemin($I).val) | (x>typemax($I).val) && throw(InexactError())
            rem(x, $I)
        end
        # function convert{T<:Integer}(::Type{$I}, x::T)
        #     (x<typemin($I).val) | (x>typemax($I).val) && throw(InexactError())
        #     rem(x, $I)
        # end
        convert(::Type{Bool}, x::$I) = convert(Bool, x.val)
        convert{T<:BaseInteger}(::Type{T}, x::$I) = convert(T, x.val)
        # convert{T<:Integer}(::Type{T}, x::$I) = convert(T, x.val)
        promote_rule{T<:BaseSigned}(::Type{$I}, ::Type{T}) = T

        leading_zeros(x::$I) =
            leading_zeros(x.val << (8 - nbits($I)) | Int8(-1) >>> nbits($I))
        leading_ones(x::$I) = leading_ones(x.val << (8 - nbits($I)))
        trailing_zeros(x::$I) = trailing_zeros(x.val | Int8(-1) << nbits($I))
        trailing_ones(x::$I) =
            trailing_ones(x.val & ~(Int8(-1) << nbits($I)))

        ~(x::$I) = $I(~x.val)
        (&)(x::$I, y::$I) = $I(x.val & y.val)
        (|)(x::$I, y::$I) = $I(x.val | y.val)
        ($)(x::$I, y::$I) = $I(x.val $ y.val)

        <<(x::$I, y::Int) = rem(x.val << y, $I)
        >>(x::$I, y::Int) = $I(x.val >> y)
        >>>(x::$I, y::Int) = rem((x.val & ~(Int(-1) << nbits($I))) >>> y, $I)

        <=(x::$I, y::$I) = x.val <= y.val
        <(x::$I, y::$I) = x.val < y.val

        +(x::$I) = x
        -(x::$I) = rem(-x.val, $I)
        abs(x::$I) = rem(abs(x.val), $I)
        +(x::$I, y::$I) = rem(x.val + y.val, $I)
        -(x::$I, y::$I) = rem(x.val - y.val, $I)
        *(x::$I, y::$I) = rem(x.val * y.val, $I)
        function div(x::$I, y::$I)
            (x == typemin($I)) & (y == $I(-1)) && throw(DivideError())
            rem(div(x.val, y.val), $I)
        end
        rem(x::$I, y::$I) = rem(rem(x.val, y.val), $I)
        function fld(x::$I, y::$I)
            (x == typemin($I)) & (y == $I(-1)) && throw(DivideError())
            rem(fld(x.val, y.val), $I)
        end
        mod(x::$I, y::$I) = rem(mod(x.val, y.val), $I)
        function cld(x::$I, y::$I)
            (x == typemin($I)) & (y == $I(-1)) && throw(DivideError())
            rem(cld(x.val, y.val), $I)
        end
    end
end

for U in (:UInt1, :UInt2, :UInt4)
    @eval begin
        print(io::IO, x::$U) = print(io, x.val)
        show(io::IO, x::$U) = show(io, x.val)

        mask(::Type{$U}) = (0x01 << nbits($U)) - 0x01
        # zext(x::$U) = x << (8 - nbits($IU) >> (8 - nbits($U))

        typemin(::Type{$U}) = $U(0x00)
        typemax(::Type{$U}) = $U(mask($U))
        rem(x::$U, ::Type{$U}) = x
        rem(x::$U, ::Type{Bool}) = rem(x.val, Bool)
        rem{T<:BaseInteger}(x::$U, ::Type{T}) = rem(x.val, T)
        rem(x::Bool, ::Type{$U}) = $U(rem(x, UInt8))
        rem{T<:BaseInteger}(x::T, ::Type{$U}) = $U(rem(x, UInt8) & mask($U))
        convert(::Type{$U}, x::$U) = x
        convert(::Type{$U}, x::Bool) = $U(convert(UInt8, x))
        function convert{T<:BaseInteger}(::Type{$U}, x::T)
            (x<typemin($U).val) | (x>typemax($U).val) && throw(InexactError())
            rem(x, $U)
        end
        convert(::Type{Bool}, x::$U) = convert(Bool, x.val)
        convert{T<:BaseInteger}(::Type{T}, x::$U) = convert(T, x.val)
        promote_rule{T<:BaseUnsigned}(::Type{$U}, ::Type{T}) = T

        leading_zeros(x::$U) = leading_zeros(x.val) - (8 - nbits($U))
        leading_ones(x::$U) = leading_ones(x.val | ~mask($U)) - (8 - nbits($U))
        trailing_zeros(x::$U) = trailing_zeros(x.val | ~mask($U))
        trailing_ones(x::$U) = trailing_ones(x.val)

        ~(x::$U) = rem(~x.val, $U)
        (&)(x::$U, y::$U) = $U(x.val & y.val)
        (|)(x::$U, y::$U) = $U(x.val | y.val)
        ($)(x::$U, y::$U) = $U(x.val $ y.val)

        <<(x::$U, y::Int) = rem(x.val << y, $U)
        >>(x::$U, y::Int) = $U(x.val >> y)
        >>>(x::$U, y::Int) = >>(x, y)

        <=(x::$U, y::$U) = x.val <= y.val
        <(x::$U, y::$U) = x.val < y.val

        +(x::$U) = x
        -(x::$U) = rem(-x.val, $U)
        abs(x::$U) = x
        +(x::$U, y::$U) = rem(x.val + y.val, $U)
        -(x::$U, y::$U) = rem(x.val - y.val, $U)
        *(x::$U, y::$U) = rem(x.val * y.val, $U)
        div(x::$U, y::$U) = rem(div(x.val, y.val), $U)
        rem(x::$U, y::$U) = rem(rem(x.val, y.val), $U)
        fld(x::$U, y::$U) = rem(fld(x.val, y.val), $U)
        mod(x::$U, y::$U) = rem(mod(x.val, y.val), $U)
        cld(x::$U, y::$U) = rem(cld(x.val, y.val), $U)
    end
end

end


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    @assert precompile(Tuple{CUDA.HostKernel{identity, Tuple{Nothing}},Nothing})
    @assert precompile(Tuple{Type{CuModule},String,Dict{CUDA.CUjit_option_enum, Any}})
    @assert precompile(Tuple{typeof(CUDA.prepare_cuda_call)})
    @assert precompile(Tuple{typeof(GPUCompiler.load_runtime),GPUCompiler.CompilerJob{GPUCompiler.PTXCompilerTarget, CUDA.CUDACompilerParams},LLVM.Context})
    @assert precompile(Tuple{typeof(cufunction),typeof(identity),Type{Tuple{Nothing}}})
    @assert precompile(Tuple{typeof(which(CUDA.pack_arguments,(Function,Vararg{Any, N} where N,)).generator.gen),Any,Any,Any})
end


text:
bar(a,b,x::Vararg{Any,2}) = (a,b,x)
println(bar(1,2,3))
println(bar(1,2,3,4))
println(bar(1,2,3,4,5))


text:
using MakieGallery

function record_selection(
    database, dir;
    selection = MakieGallery.get_unrecorded_examples(MakieGallery.database, dir),
    kwargs...
    )

    db = filter(x -> x.unique_name ∈ selection, MakieGallery.database)

    old_db = copy(database)

    empty!(MakieGallery.database)
    append!(MakieGallery.database, db)

    examples = MakieGallery.record_examples(repo; kwargs...)

    empty!(MakieGallery.database)
    append!(MakieGallery.database, old_db)

    return examples
end

empty!(MakieGallery.plotting_backends)
append!(MakieGallery.plotting_backends, ["Makie"])

# we need to clone the repo - this means git must be installed on the user system.

# Here, we reorder the database, to make it easier to see.
database = MakieGallery.load_database([
    "tutorials.jl",
    "layouting.jl",
    "legends.jl",
    "statsmakie.jl",
    "geomakie.jl",
    "examples2d.jl",
    "examples3d.jl",
    "interactive.jl",
    "documentation.jl",
    "diffeq.jl",
    "implicits.jl",
    "recipes.jl",
    "bigdata.jl",
    "short_tests.jl",
    "attributes.jl",
    "intermediate.jl",
]);

# where is the refimage repo?
repo = get(ENV, "MAKIEGALLERY_REFIMG_PATH", joinpath(homedir(), ".julia", "dev", "MakieReferenceImages"))

# clone if not present
isdir(repo) || run(`git clone --depth=1 https://github.com/JuliaPlots/MakieReferenceImages $repo`)

gallery     = joinpath(repo, "gallery")
recordings  = joinpath(@__DIR__, "test_recordings")
differences = joinpath(@__DIR__, "tested_different")

# MakieGallery.generate_preview(recordings, joinpath(recordings, "index.html"))

# Here, we redefine the order of the refimages, to our preferred order.

preferred_order = abspath.(joinpath.(
    dirname(pathof(MakieGallery)),
    "..",
    "examples",
    [
    "tutorials.jl",
    "statsmakie.jl",
    "examples2d.jl",
    "examples3d.jl",
    "interactive.jl",
    "layouting.jl",
    "legends.jl",
    "diffeq.jl",
    "implicits.jl",
    "geomakie.jl",
    "documentation.jl",
    "recipes.jl",
    "bigdata.jl",
    "attributes.jl",
    "intermediate.jl",
    "short_tests.jl",
    ]
))

sort!(database, by = x -> findfirst(==(x.file), preferred_order));

diff_uids  = Symbol.(readdir(differences))
unrec_uids = MakieGallery.get_unrecorded_examples(MakieGallery.database, gallery)
recorded_examples = Symbol.(readdir(recordings))
gallery_examples  = Symbol.(readdir(gallery))

to_record = setdiff(unrec_uids, gallery_examples)

record_selection(database, recordings; selection = to_record, generate_thumbnail = true)

for uid in union(diff_uids, unrec_uids)
    if uid ∈ recorded_examples
        mkpath(joinpath(gallery, string(uid)))
        cp(
            joinpath(recordings, string(uid)),
            joinpath(gallery, string(uid));
            force = true
        )
    else
        @warn "Example with unique name $uid was not recorded and not listed."
    end
end

# generate `thumb.jpg` for every directory in `recordings`
MakieGallery.generate_thumbnails(gallery)

# generate HTML pages for the Gallery
MakieGallery.gallery_from_recordings(
    repo * "/gallery",
    joinpath(gallery, "index.html");
    hltheme = MakieGallery.Highlights.Themes.DefaultTheme
)

cd(repo)
run(`git stage -A`)
run(`git commit -am "Reference image update"`)
run(`git push`)


text:
@testset "ISInDy" begin

    # Create a test problem
    function simple(u, p, t)
        return [(2.0u[2]^2 - 3.0)/(1.0 + u[1]^2); -u[1]^2/(2.0 + u[2]^2); (1-u[2])/(1+u[3]^2)]
    end

    u0 = [2.37; 1.58; -3.10]
    tspan = (0.0, 10.0)
    prob = ODEProblem(simple, u0, tspan)
    sol = solve(prob, Tsit5(), saveat = 0.1)

    # Create the differential data
    X = sol[:,:]
    DX = similar(X)
    for (i, xi) in enumerate(eachcol(X))
        DX[:, i] = simple(xi, [], 0.0)
    end

    # Create a basis
    @variables u[1:3]
    polys = ModelingToolkit.Operation[]
    # Lots of basis functions
    for i ∈ 0:5
        if i == 0
            push!(polys, u[1]^0)
        end
        for ui in u
            if i > 0
                push!(polys, ui^i)
            end
        end
    end

    basis= Basis(polys, u)

    opt = ADM(1e-3)
    Ψ = ISInDy(X, DX, basis, opt = opt, maxiter = 100, rtol = 0.1)

    # Simulate
    estimator = ODEProblem(dynamics(Ψ), u0, tspan)
    sol_ = solve(estimator, Tsit5(), saveat = 0.1)
    @test sol[:,:] ≈ sol_[:,:]
end


text:
module test_julia_op

using Test

for arg in [(+, +, +) (+, -, -) (-, +, -) (-, -, +)]
    (a, b, c) = nameof.(arg)
    @eval Base.$a(::typeof($b)) = $c
end

@test ==(+-+, -)
@test ==(+-+-+, +)

@test (π < π) == false

@test √4π == 2pi
@test 1e-3 == 0.001
@test √4π ≈ 2pi atol=1e-3

struct G
end

import Base: +
x::G + y::G = 2

@test G() + G() == 2

# Julia PR #29819
if VERSION >= v"1.1.0-DEV.764"
for x in ('a', "a", ['a'], 1, pi)
    @test identity.(x) == identity(x) == x
end
end # if VERSION >= v"1.1.0-DEV.764"

end # module test_julia_op

module test_julia_op2

using Test

f(op::Union{typeof.((+, -, *, /))...}) = 1

@test f(+) == f(-) == f(*) == f(/) == 1
@test_throws MethodError f(\)

end # module test_julia_op2

module test_julia_plus

using Test

module A
    + = 3
end

module B
end

module C
    + = 3
module D
end
end

@test A.:+ == 3
@test B.:+ == Base.:+
@test A.:+ != B.:+
@test C.:+ != C.D.:+

end # module test_julia_plus